Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

VASP scripts

bubbleplot3 was not created by me, see comments in code
  • Loading branch information...
commit bf7153bbba68ee6ac46336b16d45786da1ef3521 0 parents
@cheniffer authored
Showing with 3,905 additions and 0 deletions.
  1. +22 −0 .gitattributes
  2. +163 −0 .gitignore
  3. 0  CONTCAR.txt
  4. +172 −0 POSCAR
  5. +158 −0 VASP_plot.asv
  6. +162 −0 VASP_plot.m
  7. +170 −0 boo.txt
  8. +120 −0 bubbleplot3.m
  9. +47 −0 center.m
  10. +25 −0 check_oszicar.asv
  11. +25 −0 check_oszicar.m
  12. +43 −0 check_outcar.m
  13. BIN  chgcar.mat
  14. +94 −0 contcar.m
  15. BIN  copper.mat
  16. BIN  fig code/benzene_pd.fig
  17. +11 −0 fig code/blah.m
  18. BIN  fig code/light_on_off.fig
  19. BIN  fig code/o2.fig
  20. BIN  fig code/origin_plots.mat
  21. BIN  fig code/photo_curr.fig
  22. BIN  fig code/photo_curr.mat
  23. +66 −0 fig code/polyfit0.m
  24. BIN  fig code/smoothed_photocurr_sub_dark.fig
  25. +9 −0 fig code/temp_fit.asv
  26. +12 −0 fig code/temp_fit.m
  27. BIN  figs/bubbles.fig
  28. BIN  figs/iso2at219_1x1x1.fig
  29. BIN  figs/iso2at219_2x2x1.fig
  30. BIN  figs/iso2at219_2x2x2 2.fig
  31. BIN  figs/iso2at219_2x2x2.fig
  32. BIN  figs/test.fig
  33. +23 −0 fscript.m
  34. BIN  iso_chain.mat
  35. BIN  isochain.fig
  36. +21 −0 lateral_offset.m
  37. +81 −0 lateral_script.m
  38. +102 −0 live/daqstopbutton.m
  39. +57 −0 live/demo_tmwidgets.m
  40. +28 −0 live/license.txt
  41. +289 −0 live/specgramscope.m
  42. +61 −0 live/specgramscope_documentation.m
  43. +295 −0 live/spectrumscope.m
  44. +71 −0 live/spectrumscope_documentation.m
  45. +308 −0 live/stripchart.m
  46. +74 −0 live/stripchart_documentation.m
  47. +164 −0 live/thermometer.m
  48. +55 −0 live/thermometer_documentation.m
  49. +37 −0 make_pos.m
  50. BIN  new_model.mat
  51. +235 −0 peakfinder.m
  52. +20 −0 plot_all.m
  53. +72 −0 process.m
  54. +72 −0 process_chgcar.asv
  55. +71 −0 process_chgcar.m
  56. +34 −0 script.m
  57. +55 −0 supercell.asv
  58. +120 −0 supercell.m
  59. +4 −0 temp.asv
  60. +4 −0 temp.m
  61. +10 −0 tile_script.m
  62. +116 −0 view.asv
  63. +127 −0 view_cell.m
22 .gitattributes
@@ -0,0 +1,22 @@
+# Auto detect text files and perform LF normalization
+* text=auto
+
+# Custom for Visual Studio
+*.cs diff=csharp
+*.sln merge=union
+*.csproj merge=union
+*.vbproj merge=union
+*.fsproj merge=union
+*.dbproj merge=union
+
+# Standard to msysgit
+*.doc diff=astextplain
+*.DOC diff=astextplain
+*.docx diff=astextplain
+*.DOCX diff=astextplain
+*.dot diff=astextplain
+*.DOT diff=astextplain
+*.pdf diff=astextplain
+*.PDF diff=astextplain
+*.rtf diff=astextplain
+*.RTF diff=astextplain
163 .gitignore
@@ -0,0 +1,163 @@
+#################
+## Eclipse
+#################
+
+*.pydevproject
+.project
+.metadata
+bin/
+tmp/
+*.tmp
+*.bak
+*.swp
+*~.nib
+local.properties
+.classpath
+.settings/
+.loadpath
+
+# External tool builders
+.externalToolBuilders/
+
+# Locally stored "Eclipse launch configurations"
+*.launch
+
+# CDT-specific
+.cproject
+
+# PDT-specific
+.buildpath
+
+
+#################
+## Visual Studio
+#################
+
+## Ignore Visual Studio temporary files, build results, and
+## files generated by popular Visual Studio add-ons.
+
+# User-specific files
+*.suo
+*.user
+*.sln.docstates
+
+# Build results
+[Dd]ebug/
+[Rr]elease/
+*_i.c
+*_p.c
+*.ilk
+*.meta
+*.obj
+*.pch
+*.pdb
+*.pgc
+*.pgd
+*.rsp
+*.sbr
+*.tlb
+*.tli
+*.tlh
+*.tmp
+*.vspscc
+.builds
+*.dotCover
+
+## TODO: If you have NuGet Package Restore enabled, uncomment this
+#packages/
+
+# Visual C++ cache files
+ipch/
+*.aps
+*.ncb
+*.opensdf
+*.sdf
+
+# Visual Studio profiler
+*.psess
+*.vsp
+
+# ReSharper is a .NET coding add-in
+_ReSharper*
+
+# Installshield output folder
+[Ee]xpress
+
+# DocProject is a documentation generator add-in
+DocProject/buildhelp/
+DocProject/Help/*.HxT
+DocProject/Help/*.HxC
+DocProject/Help/*.hhc
+DocProject/Help/*.hhk
+DocProject/Help/*.hhp
+DocProject/Help/Html2
+DocProject/Help/html
+
+# Click-Once directory
+publish
+
+# Others
+[Bb]in
+[Oo]bj
+sql
+TestResults
+*.Cache
+ClientBin
+stylecop.*
+~$*
+*.dbmdl
+Generated_Code #added for RIA/Silverlight projects
+
+# Backup & report files from converting an old project file to a newer
+# Visual Studio version. Backup files are not needed, because we have git ;-)
+_UpgradeReport_Files/
+Backup*/
+UpgradeLog*.XML
+
+
+
+############
+## Windows
+############
+
+# Windows image file caches
+Thumbs.db
+
+# Folder config file
+Desktop.ini
+
+
+#############
+## Python
+#############
+
+*.py[co]
+
+# Packages
+*.egg
+*.egg-info
+dist
+build
+eggs
+parts
+bin
+var
+sdist
+develop-eggs
+.installed.cfg
+
+# Installer logs
+pip-log.txt
+
+# Unit test / coverage reports
+.coverage
+.tox
+
+#Translations
+*.mo
+
+#Mr Developer
+.mr.developer.cfg
+
+# Mac crap
+.DS_Store
0  CONTCAR.txt
No changes.
172 POSCAR
@@ -0,0 +1,172 @@
+temp
+1.0
+0 15.3 0
+-19.8753 1.275 0
+0 0 15.3
+Cu C H N
+111 32 16 4
+Selective Dynamics
+Cartesian
+-0.95161 1.2569 15.2995 T T T
+-0.95383 3.8026 15.2863 T T T
+-0.97685 6.4074 15.1533 T T T
+-0.95828 9.0358 15.2928 T T T
+-0.95228 11.5851 15.2989 T T T
+-3.184 5.1207 15.2854 T T T
+-5.4063 1.2744 15.293 T T T
+-5.4055 3.8479 15.2899 T T T
+-5.4075 6.4237 15.2826 T T T
+-5.4017 8.9935 15.2956 T T T
+-5.40546 11.568 15.2928 T T T
+-17.9047 15.3329 13.1507 T T T
+-17.9176 2.57573 13.1751 T T T
+-17.9177 5.12397 13.1904 T T T
+-17.9241 7.68501 13.2223 T T T
+-17.921 10.2384 13.2027 T T T
+-17.9087 12.7885 13.1825 T T T
+-0.238742 0.0270389 13.1411 T T T
+-0.242283 2.56242 13.1743 T T T
+-0.250246 5.11123 13.1275 T T T
+-0.248035 7.69806 13.1156 T T T
+-0.241627 10.2558 13.1678 T T T
+-0.236595 12.7902 13.1287 T T T
+-2.45034 1.30344 13.1738 T T T
+-2.46286 3.84805 13.1457 T T T
+-2.47584 6.40951 13.1168 T T T
+-2.4578 8.97239 13.1895 T T T
+-2.44718 11.5181 13.1693 T T T
+-2.44235 14.0567 13.1377 T T T
+-4.65625 15.3304 13.1729 T T T
+-4.66636 2.57769 13.1794 T T T
+-4.6784 5.13264 13.183 T T T
+-4.66931 7.68952 13.1981 T T T
+-4.66158 10.2356 13.1701 T T T
+-4.6484 12.7789 13.1466 T T T
+-6.87981 1.30296 13.1747 T T T
+-6.88781 3.85472 13.1771 T T T
+-6.88918 6.40439 13.1742 T T T
+-6.88132 8.95455 13.1715 T T T
+-6.8695 11.4977 13.1541 T T T
+-6.86478 14.045 13.1744 T T T
+-9.07771 15.3102 13.1629 T T T
+-9.0889 2.55866 13.1425 T T T
+-9.09311 5.10941 13.1624 T T T
+-9.09529 7.66451 13.1563 T T T
+-9.08924 10.2113 13.1435 T T T
+-9.07683 12.7616 13.1659 T T T
+-11.2846 1.27898 13.1628 T T T
+-11.2884 3.8325 13.1691 T T T
+-11.2951 6.3832 13.137 T T T
+-11.292 8.93376 13.1658 T T T
+-11.2853 11.4879 13.1589 T T T
+-11.2833 14.0341 13.1614 T T T
+-13.4868 15.3064 13.162 T T T
+-13.4921 2.55586 13.1683 T T T
+-13.4868 5.10555 13.1719 T T T
+-13.4897 7.66345 13.1683 T T T
+-13.4925 10.2125 13.1661 T T T
+-13.4868 12.7579 13.1627 T T T
+-15.6939 1.2921 13.1763 T T T
+-15.7004 3.84227 13.1789 T T T
+-15.7079 6.39133 13.1922 T T T
+-15.7009 8.94033 13.1878 T T T
+-15.697 11.4912 13.1718 T T T
+-15.6931 14.0408 13.1712 T T T
+-10.7287 11.5898 17.5376 T T T
+-0.95161 1.2569 12.0992 T T T
+-0.95383 3.8026 12.0926 T T T
+-18.6003 1.275 15.3 T T T
+-18.5962 3.838 15.3242 T T T
+-18.5864 6.4025 15.3163 T T T
+-18.6174 8.9905 15.3595 T T T
+-18.5962 11.5571 15.3149 T T T
+-18.6006 14.119 15.3 T T T
+-0.93336 14.025 15.3 T T T
+-3.1881 2.5496 15.3022 T T T
+-3.1701 7.7172 15.312 T T T
+-3.18777 10.2915 15.3026 T T T
+-3.17779 12.8458 15.3077 T T T
+-3.14173 15.3 15.3 T T T
+-5.35009 14.025 15.3 T T T
+-7.6242 2.5878 15.3036 T T T
+-7.6245 5.1553 15.3039 T T T
+-7.6237 7.7213 15.3023 T T T
+-7.62166 10.2896 15.3 T T T
+-7.55846 12.75 15.3 T T T
+-7.55846 15.3 15.3 T T T
+-9.7668 1.275 15.3 T T T
+-9.7668 3.825 15.3 T T T
+-9.7668 6.375 15.3 T T T
+-9.7668 8.925 15.3 T T T
+-9.7668 11.475 15.3 T T T
+-9.7668 14.025 15.3 T T T
+-11.9752 2.55 15.3 T T T
+-11.9752 5.1 15.3 T T T
+-11.9752 7.65 15.3 T T T
+-11.9752 10.2 15.3 T T T
+-11.9752 12.75 15.3 T T T
+-11.9752 15.3 15.3 T T T
+-14.1836 1.275 15.3 T T T
+-14.1836 3.825 15.3 T T T
+-14.1836 6.375 15.3 T T T
+-14.1836 8.925 15.3 T T T
+-14.1836 11.475 15.3 T T T
+-14.1836 14.025 15.3 T T T
+-16.3919 2.55 15.3 T T T
+-16.3919 5.1 15.3 T T T
+-16.3919 7.65 15.3 T T T
+-16.3919 10.2 15.3 T T T
+-16.3919 12.75 15.3 T T T
+-16.3919 15.3 15.3 T T T
+-13.7339 6.82259 18.1583 T T T
+-16.0849 7.36001 18.958 T T T
+-16.9847 6.2821 19.223 T T T
+-16.5399 4.9361 19.2325 T T T
+-17.4025 3.90181 19.5435 T T T
+-18.8629 7.91422 19.2902 T T T
+-18.3751 6.56218 19.471 T T T
+-19.2139 5.47862 19.8301 T T T
+-18.7407 4.17747 19.8761 T T T
+-16.5016 8.72148 19.073 T T T
+-15.5891 9.80465 18.9716 T T T
+-16.0203 11.1105 19.1094 T T T
+-17.8926 8.99899 19.2662 T T T
+-18.2954 10.3478 19.4178 T T T
+-17.3815 11.3839 19.3529 T T T
+-1.32951 6.11638 18.9679 T T T
+-10.2633 6.35546 18.2007 T T T
+-7.85554 6.17481 19.0006 T T T
+-7.07576 7.35068 19.2018 T T T
+-7.64529 8.65017 19.1377 T T T
+-6.86734 9.77559 19.3151 T T T
+-5.06243 5.92001 19.4586 T T T
+-5.65946 7.22407 19.4377 T T T
+-4.89502 8.40527 19.6185 T T T
+-5.48188 9.65361 19.5621 T T T
+-7.29114 4.87086 19.1013 T T T
+-8.0622 3.68927 18.9335 T T T
+-7.47677 2.44278 19.0228 T T T
+-5.88385 4.74809 19.3486 T T T
+-5.31743 3.4494 19.438 T T T
+-6.09321 2.32043 19.2835 T T T
+-2.61054 6.11597 19.2663 T T T
+-15.4921 4.72438 19.0156 T T T
+-17.035 2.87434 19.5518 T T T
+-0.374807 4.40602 20.0919 T T T
+-19.4109 3.36803 20.1685 T T T
+-14.532 9.59425 18.8046 T T T
+-15.3017 11.9295 19.0389 T T T
+-19.3489 10.5554 19.6076 T T T
+-17.7159 12.4138 19.4908 T T T
+-8.71639 8.75104 18.9561 T T T
+-7.32709 10.7645 19.2655 T T T
+-3.82544 8.3167 19.813 T T T
+-4.87498 10.5491 19.706 T T T
+-9.13177 3.77451 18.7382 T T T
+-8.0874 1.54746 18.8936 T T T
+-4.24884 3.36902 19.6402 T T T
+-5.63751 1.33167 19.3591 T T T
+-14.8189 7.08153 18.5727 T T T
+-0.281853 6.9998 19.0829 T T T
+-9.15907 6.29424 18.6336 T T T
+-3.7346 5.70679 19.5975 T T T
158 VASP_plot.asv
@@ -0,0 +1,158 @@
+function coord_list = VASP_plot(supercell)
+
+e_table = ['C '; 'N '; 'H '; 'Cu'];
+r_table = [.73; .71; .31; 1.28];
+c_table = [[0,1,0]; [1,0,0]; [0,0,1]; [0,1,1]; [1,0,1]; [1,1,0]; [1,1,1]];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% GET COORDINATES FROM CONTCAR %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+fid=fopen('CONTCAR.txt');
+unit_vectors = [];
+coord_list = [];
+restricted = [];
+counter = 0;
+
+fgetl(fid); % skip some lines
+fgetl(fid);
+
+%% build the unit cell vectors
+for i=1:3
+ temp = fgetl(fid);
+ [v1, temp] = strtok(temp); %x
+ [v2, temp] = strtok(temp); %y
+ [v3, temp] = strtok(temp); %z
+ unit_vectors = [unit_vectors; str2num(v1) str2num(v2) str2num(v3)];
+end
+
+%% get the coordinates and scale them by the unit vectors
+while(~feof(fid))
+ line = fgetl(fid);
+ if strmatch('Direct',line)==1
+ coord = fgetl(fid);
+ [length,width] = size(coord);
+ while (width>10)
+ formatted_coord = strrep(coord, '-', ' -'); %make spacing correct for delim
+ formatted_coord = strrep(formatted_coord, ' ', ' '); %make spacing correct for delim
+ %formatted_coord = strrep(formatted_coord, 'T', ' 1'); %change to int
+ %formatted_coord = strrep(formatted_coord, 'F', ' 0'); %change to int
+ [x,y,z,r1,r2,r3] = strread(formatted_coord, '%f%f%f%c%c%c', 'delimiter', ' ');
+ coord_list = [coord_list; unit_vectors(1,:)*x + unit_vectors(2,:)*y + unit_vectors(3,:)*z]; %multiply by unit vectors to scale correctly
+ restricted = [restricted; ' ' r1 ' ' r2 ' ' r3]; %store T/F restrictions
+ coord = fgetl(fid);
+ [length,width] = size(coord);
+ end
+ end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% WRITE NEW POSCAR FILE %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+[length,width] =size(coord_list);
+poscar = [num2str(coord_list) restricted];
+
+fileID = fopen('POSCAR2.txt','w');
+fid3=fopen('POSCAR.txt'); %open poscar
+l = fgetl(fid3);
+fprintf(fileID,'%s 2\n', l);
+for i=2:6 % copy these lines
+ l = fgetl(fid3);
+ fprintf(fileID,'%s\n', l);
+end
+fprintf(fileID,'Selective Dynamics \n');
+fprintf(fileID,'Cartesian \n');
+for i=1:length
+ fprintf(fileID, '%s\n', poscar(i,:));
+end
+fclose(fileID);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% EXTENSION TO SUPERCELL %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+counter1 = 0;
+counter2 = 0;
+counter3 = 0;
+if max(supercell)>1
+ for i=1:supercell(1)-1
+ temp_coord_list = [coord_list(1:length,1)-i*min(coord_list(:,1)) coord_list(1:length, 2:3)];
+ coord_list = [coord_list; temp_coord_list];
+ counter1 = counter1+1;
+ end
+ [length,width] =size(coord_list);
+ for j=1:supercell(2)-1
+ temp_coord_list = [coord_list(1:length,1) coord_list(1:length,2)+j*max(coord_list(:,2)) coord_list(1:end,3)];
+ coord_list = [coord_list; temp_coord_list];
+ counter2 = counter2+1;
+ end
+ [length,width] =size(coord_list);
+ for k=1:supercell(3)-1
+ temp_coord_list = [coord_list(1:length,1:2) coord_list(1:length,3)+k*max(coord_list(:,3))];
+ coord_list = [coord_list; temp_coord_list];
+ counter3 = counter3+1;
+ end
+end
+counter = (counter1+1)*(counter2+1)*(counter3+1)-1;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% GET RADII FROM POTCAR AND POSCAR %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%grab elements from potcar
+fid2=fopen('POTCAR.txt');
+elements = [];
+while(~feof(fid2))
+ line = fgetl(fid2);
+ [len,wid] = size(line);
+ if wid>6
+ if strmatch('PAW_PBE',line(2:8))
+ [e, junk] = strtok(line(9:end));
+ [junk,check_size] = size(e);
+ if check_size<2
+ e = [e ' '];
+ end
+ elements = [elements; e];
+ end
+ end
+end
+
+% elements
+
+%%grab radii from poscar
+[n_e, junk] = size(elements);
+radii = [];
+colors = [];
+for i=1:n_e
+ [temp,l] = strtok(l);
+ r = r_table(find(e_table==elements(i),1));
+ radii = [radii r.*ones(1,str2num(temp))];
+ rgb = r*[c_table(mod(i,7),1)*ones(str2num(temp),1) c_table(mod(i,7),2)*ones(str2num(temp),1) c_table(mod(i,7),3)*ones(str2num(temp),1)];
+ colors = [colors; rgb];
+end
+
+temp_radii = radii;
+temp_colors = colors;
+
+for i=1:counter
+ radii = [radii temp_radii];
+ colors = [colors; temp_colors];
+end
+
+fclose('all');
+
+%%%%%%%%%%%%%%%%%%%%%
+%% PLOT EVERYTHING %%
+%%%%%%%%%%%%%%%%%%%%%
+
+% output = [coord_list(:,1), coord_list(:,2), coord_list(:,3), transpose(radii), colors]
+% size(output)
+
+figure()
+bubbleplot3(coord_list(:,1), coord_list(:,2), coord_list(:,3), radii, colors)
+camlight right; lighting phong; view(60,30);
+% figure()
+% scatter3(coord_list(:,1), coord_list(:,2), coord_list(:,3))
+
+end
162 VASP_plot.m
@@ -0,0 +1,162 @@
+function coord_list = VASP_plot(supercell)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% GET COORDINATES FROM CONTCAR %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+e_table = ['C '; 'N '; 'H '; 'Cu'];
+r_table = [.73; .71; .31; 1.28];
+c_table = [[0,1,0]; [1,0,0]; [0,0,1]; [0,1,1]; [1,0,1]; [1,1,0]; [1,1,1]];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% GET COORDINATES FROM CONTCAR %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+fid=fopen('CONTCAR.txt');
+unit_vectors = [];
+coord_list = [];
+restricted = [];
+counter = 0;
+
+fgetl(fid); % skip some lines
+fgetl(fid);
+
+%% build the unit cell vectors
+for i=1:3
+ temp = fgetl(fid);
+ [v1, temp] = strtok(temp); %x
+ [v2, temp] = strtok(temp); %y
+ [v3, temp] = strtok(temp); %z
+ unit_vectors = [unit_vectors; str2num(v1) str2num(v2) str2num(v3)];
+end
+
+%% get the coordinates and scale them by the unit vectors
+while(~feof(fid))
+ line = fgetl(fid);
+ if strmatch('Direct',line)==1
+ coord = fgetl(fid);
+ [length,width] = size(coord);
+ while (width>10)
+ formatted_coord = strrep(coord, '-', ' -'); %make spacing correct for delim
+ formatted_coord = strrep(formatted_coord, ' ', ' '); %make spacing correct for delim
+ %formatted_coord = strrep(formatted_coord, 'T', ' 1'); %change to int
+ %formatted_coord = strrep(formatted_coord, 'F', ' 0'); %change to int
+ [x,y,z,r1,r2,r3] = strread(formatted_coord, '%f%f%f%c%c%c', 'delimiter', ' ');
+ coord_list = [coord_list; unit_vectors(1,:)*x + unit_vectors(2,:)*y + unit_vectors(3,:)*z]; %multiply by unit vectors to scale correctly
+ restricted = [restricted; ' ' r1 ' ' r2 ' ' r3]; %store T/F restrictions
+ coord = fgetl(fid);
+ [length,width] = size(coord);
+ end
+ end
+end
+
+[length,width] =size(coord_list);
+poscar = [num2str(coord_list) restricted];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% WRITE NEW POSCAR FILE %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+fileID = fopen('POSCAR2.txt','w');
+fid3=fopen('POSCAR.txt'); %open poscar
+l = fgetl(fid3);
+fprintf(fileID,'%s 2\n', l);
+for i=2:6 % copy these lines
+ l = fgetl(fid3);
+ fprintf(fileID,'%s\n', l);
+end
+fprintf(fileID,'Selective Dynamics \n');
+fprintf(fileID,'Cartesian \n');
+for i=1:length
+ fprintf(fileID, '%s\n', poscar(i,:));
+end
+fclose(fileID);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% EXTENSION TO SUPERCELL %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+counter1 = 0;
+counter2 = 0;
+counter3 = 0;
+if max(supercell)>1
+ for i=1:supercell(1)-1
+ temp_coord_list = [coord_list(1:length,1)-i*min(coord_list(:,1)) coord_list(1:length, 2:3)];
+ coord_list = [coord_list; temp_coord_list];
+ counter1 = counter1+1;
+ end
+ [length,width] =size(coord_list);
+ for j=1:supercell(2)-1
+ temp_coord_list = [coord_list(1:length,1) coord_list(1:length,2)+j*max(coord_list(:,2)) coord_list(1:end,3)];
+ coord_list = [coord_list; temp_coord_list];
+ counter2 = counter2+1;
+ end
+ [length,width] =size(coord_list);
+ for k=1:supercell(3)-1
+ temp_coord_list = [coord_list(1:length,1:2) coord_list(1:length,3)+k*max(coord_list(:,3))];
+ coord_list = [coord_list; temp_coord_list];
+ counter3 = counter3+1;
+ end
+end
+counter = (counter1+1)*(counter2+1)*(counter3+1)-1;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% GET RADII FROM POTCAR AND POSCAR %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%grab elements from potcar
+fid2=fopen('POTCAR.txt');
+elements = [];
+while(~feof(fid2))
+ line = fgetl(fid2);
+ [len,wid] = size(line);
+ if wid>6
+ if strmatch('PAW_PBE',line(2:8))
+ [e, junk] = strtok(line(9:end));
+ [junk,check_size] = size(e);
+ if check_size<2
+ e = [e ' '];
+ end
+ elements = [elements; e];
+ end
+ end
+end
+
+% elements
+
+%%grab radii from poscar
+[n_e, junk] = size(elements);
+radii = [];
+colors = [];
+for i=1:n_e
+ [temp,l] = strtok(l);
+ r = r_table(find(e_table==elements(i),1));
+ radii = [radii r.*ones(1,str2num(temp))];
+ rgb = r*[c_table(mod(i,7),1)*ones(str2num(temp),1) c_table(mod(i,7),2)*ones(str2num(temp),1) c_table(mod(i,7),3)*ones(str2num(temp),1)];
+ colors = [colors; rgb];
+end
+
+temp_radii = radii;
+temp_colors = colors;
+
+for i=1:counter
+ radii = [radii temp_radii];
+ colors = [colors; temp_colors];
+end
+
+fclose('all');
+
+%%%%%%%%%%%%%%%%%%%%%
+%% PLOT EVERYTHING %%
+%%%%%%%%%%%%%%%%%%%%%
+
+% output = [coord_list(:,1), coord_list(:,2), coord_list(:,3), transpose(radii), colors]
+% size(output)
+
+figure()
+bubbleplot3(coord_list(:,1), coord_list(:,2), coord_list(:,3), radii, colors)
+camlight right; lighting phong; view(60,30);
+% figure()
+% scatter3(coord_list(:,1), coord_list(:,2), coord_list(:,3))
+
+end
170 boo.txt
@@ -0,0 +1,170 @@
+iso2at2l9 2
+1.
+0. 15.3 0.
+-19.8753 1.275 0.
+0. 0. 15.3
+110 32 16 4
+Selective Dynamics
+Cartesian
+ -1.38304 7.48376 1.80637 T T T
+ -2.24969 5.0981 1.73302 T T T
+ 0.0203784 0.00425982 -0.0941033 T T T
+ 0.0186515 2.5484 -0.0149547 T T T
+-0.0320651 5.12635 0.107165 T T T
+ 0.0633685 7.70543 -0.132111 T T T
+ 0.0194848 10.222 -0.038992 T T T
+ 0.0205167 12.7586 -0.114824 T T T
+ -2.16554 1.26611 -0.13285 T T T
+ -2.13198 3.80087 -0.277456 T T T
+ -2.14959 6.4295 15.0159 T T T
+ -2.16806 8.94978 -0.079781 T T T
+ -2.18355 11.4963 -0.112526 T T T
+ -2.17888 14.0313 -0.141566 T T T
+ -4.36017 2.56214 -0.125364 T T T
+ -4.32979 5.10579 -0.0964857 T T T
+ -4.3605 7.6687 -0.0849354 T T T
+ -4.38139 10.2223 -0.108674 T T T
+ -4.37901 12.7681 -0.135923 T T T
+ -4.37182 15.3111 -0.148121 T T T
+ -6.5728 1.28148 -0.124898 T T T
+ -6.56542 3.83016 -0.0998559 T T T
+ -6.56629 6.3781 -0.102931 T T T
+ -6.57752 8.9327 -0.0915743 T T T
+ -6.59022 11.4834 -0.110949 T T T
+ -6.58124 14.0356 -0.131274 T T T
+ -8.78786 2.55141 -0.104491 T T T
+ -8.78479 5.09978 -0.0992059 T T T
+ -8.78408 7.65398 -0.0982009 T T T
+ -8.78712 10.206 -0.0959065 T T T
+ -8.79404 12.7574 -0.106209 T T T
+ -8.79083 15.3065 -0.109974 T T T
+ -11.0049 1.28873 -0.105455 T T T
+ -11.0015 3.83931 -0.10181 T T T
+ -10.9965 6.38863 -0.0988328 T T T
+ -11.0012 8.9423 -0.0946417 T T T
+ -11.0055 11.4964 -0.0940508 T T T
+ -11.0024 14.042 -0.0936663 T T T
+ -13.2243 2.56362 -0.0835642 T T T
+ -13.2224 5.11344 -0.0909856 T T T
+ -13.2124 7.65552 -0.0962605 T T T
+ -13.2143 10.2101 -0.0882736 T T T
+ -13.2205 12.7655 -0.0898119 T T T
+ -13.223 15.3073 -0.0745089 T T T
+ -15.4377 1.29063 -0.084301 T T T
+ -15.4392 3.83842 -0.0659205 T T T
+ -15.4393 6.38276 -0.074799 T T T
+ -15.4212 8.93044 -0.0787159 T T T
+ -15.422 11.487 -0.0927341 T T T
+ -15.4341 14.0413 -0.0989481 T T T
+ -17.6497 2.57353 -0.0661135 T T T
+ -17.6588 5.1283 -0.0219986 T T T
+ -17.6449 7.66943 -0.0785924 T T T
+ -17.6224 10.2282 -0.0783746 T T T
+ -17.6385 12.7653 -0.0925731 T T T
+ -17.6534 15.3228 -0.105011 T T T
+ -19.1528 15.3138 13.056 T T T
+ -19.1556 2.56601 13.1037 T T T
+ -19.1864 5.12379 13.2029 T T T
+ -19.168 7.65354 13.1329 T T T
+ -19.1333 10.2336 13.0826 T T T
+ -19.146 12.7613 13.0876 T T T
+ -1.46614 15.2996 13.0007 T T T
+ -1.46133 2.52896 12.9695 T T T
+ -1.47821 5.11308 12.9251 T T T
+ -1.46064 7.69439 12.9595 T T T
+ -1.46322 10.2298 13.0752 T T T
+ -1.46775 12.7654 13.0112 T T T
+ -3.668 1.28733 13.0221 T T T
+ -3.67717 3.83289 12.9962 T T T
+ -3.67986 6.40089 13.0157 T T T
+ -3.66547 8.95228 13.0772 T T T
+ -3.66935 11.493 13.0356 T T T
+ -3.67019 14.0384 12.9996 T T T
+ -5.88214 15.3082 13.0278 T T T
+ -5.87841 2.5618 13.0478 T T T
+ -5.87613 5.10737 13.0901 T T T
+ -5.8738 7.66389 13.0962 T T T
+ -5.87649 10.2099 13.063 T T T
+ -5.8828 12.7527 13.042 T T T
+ -8.08967 1.27745 13.0549 T T T
+ -8.08774 3.83056 13.0727 T T T
+ -8.0887 6.37588 13.081 T T T
+ -8.08156 8.93034 13.0867 T T T
+ -8.09596 11.4859 13.0728 T T T
+ -8.10098 14.026 13.0494 T T T
+ -10.3074 15.3171 13.0647 T T T
+ -10.3038 2.5675 13.0756 T T T
+ -10.3036 5.11902 13.0741 T T T
+ -10.3004 7.66494 13.08 T T T
+ -10.3052 10.2138 13.0903 T T T
+ -10.3082 12.7743 13.0756 T T T
+ -12.524 1.29214 13.0915 T T T
+ -12.5199 3.83889 13.086 T T T
+ -12.5164 6.38632 13.0795 T T T
+ -12.5164 8.93281 13.0842 T T T
+ -12.5218 11.4873 13.0916 T T T
+ -12.5267 14.0501 13.0922 T T T
+ -14.7384 15.321 13.0931 T T T
+ -14.7403 2.57027 13.1024 T T T
+ -14.7373 5.10739 13.105 T T T
+ -14.7374 7.6515 13.0988 T T T
+ -14.7283 10.206 13.0926 T T T
+ -14.7272 12.7612 13.0878 T T T
+ -16.9472 1.3059 13.098 T T T
+ -16.9632 3.85727 13.1234 T T T
+ -16.9604 6.38874 13.1123 T T T
+ -16.942 8.94133 13.1068 T T T
+ -16.9296 11.4944 13.0861 T T T
+ -16.9438 14.0424 13.0868 T T T
+ -12.3543 8.4237 4.56015 T T T
+ -14.9022 8.32534 4.33517 T T T
+ -15.5966 7.09128 4.25538 T T T
+ -14.9261 5.83593 4.3118 T T T
+ -15.6389 4.66439 4.19328 T T T
+ -17.6801 8.36631 3.98164 T T T
+ -17.0278 7.11249 4.07781 T T T
+ -17.7253 5.87711 3.96522 T T T
+ -17.0466 4.6825 4.01953 T T T
+ -15.5596 9.57577 4.24891 T T T
+ -14.8548 10.8117 4.30135 T T T
+ -15.5253 12.0059 4.18666 T T T
+ -16.9885 9.60156 4.07495 T T T
+ -17.6477 10.8551 3.96634 T T T
+ -16.9316 12.0282 4.02289 T T T
+ -0.227992 7.33499 3.26535 T T T
+ -11.0541 5.61669 4.54341 T T T
+ -8.50903 5.75167 4.34845 T T T
+ -7.81016 6.98201 4.26525 T T T
+ -8.47535 8.24119 4.31318 T T T
+ -7.76214 9.40827 4.1904 T T T
+ -5.73487 5.69775 3.98876 T T T
+ -6.38047 6.95728 4.08426 T T T
+ -5.67653 8.18905 3.96648 T T T
+ -6.35336 9.38656 4.0161 T T T
+ -7.85599 4.49608 4.27611 T T T
+ -8.56182 3.26146 4.34235 T T T
+ -7.89309 2.06285 4.2321 T T T
+ -6.42939 4.46624 4.09474 T T T
+ -5.77394 3.21081 3.98731 T T T
+ -6.48772 2.03897 4.05358 T T T
+ -3.35102 5.37323 3.20681 T T T
+ -13.8418 5.81279 4.43856 T T T
+ -15.1149 3.70832 4.22616 T T T
+ -18.8083 5.88979 3.82825 T T T
+ -17.5914 3.74485 3.91673 T T T
+ -13.7709 10.7922 4.42302 T T T
+ -14.9685 12.944 4.21783 T T T
+ -18.7301 10.8796 3.83742 T T T
+ -17.4528 12.984 3.93215 T T T
+ -9.55872 8.26757 4.44937 T T T
+ -8.28142 10.3661 4.22373 T T T
+ -4.59447 8.16934 3.8293 T T T
+ -5.80151 10.324 3.9153 T T T
+ -9.64286 3.28175 4.47629 T T T
+ -8.44893 1.12458 4.27643 T T T
+ -4.69224 3.18274 3.84632 T T T
+ -5.9685 1.08288 3.96091 T T T
+ -13.5334 8.32434 4.46932 T T T
+ -19.0163 8.41191 3.7061 T T T
+ -9.87562 5.74687 4.47161 T T T
+ -4.40238 5.64453 3.68446 T T T
120 bubbleplot3.m
@@ -0,0 +1,120 @@
+function handles = bubbleplot3(x,y,z,r,varargin)
+%BUBBLEPLOT3 A simple 3D-bubbleplot.
+% BUBBLEPLOT3() is a three-dimensional bubbleplot.
+%
+% BUBBLEPLOT3(x,y,z,r), where x, y, z and r are four vectors of the
+% same length, plots bubbles of radii r in 3-space with centers at
+% the points whose coordinates are the elements of x, y and z. If r
+% is a matrix of size numel(x)x3, BUBBLEPLOT3 produces ellipsoids with
+% centers x(i),y(i),z(i) and radii r(i,1), r(i,2) and r(i,3).
+%
+% BUBBLEPLOT3(x,y,z,r,c), where c is a rgb-triplet array (in [0,1])
+% with numel(x) rows, plots bubbles with colours specified by c.
+%
+% BUBBLEPLOT3(x,y,z,r,c,Alpha), where Alpha is a scalar with value from
+% 0.0 to 1.0, plots bubbles with FaceAlpha Alpha.
+%
+% BUBBLEPLOT3(x,y,z,r,c,Alpha,n,m), where m and n are scalar values that
+% decides the size of the arrays used to render the bubbles.
+% The largest radius in the set is rendered with (n+1)x(n+1) points.
+% To increase efficiency, the number of rendering points used is
+% decreasing linearly with the radius but is never rendered with
+% fewer points than (m+1)x(m+1).
+%
+% BUBBLEPLOT3(x,y,z,r,c,Alpha,n,m,'PropertyName',PropertyValue,...)
+% Any property-value pair setting valid for a SURFACE object can be
+% passed as optional parameters.
+%
+% BUBBLEPLOT3 returns a column vector of handles to surface objects.
+%
+% Example: Three bubbles using your standard colormap
+%
+% x=[0 1 -1];
+% y=[-1 0 1];
+% z=[1 -1 0];
+% r=[.3 .6 .9]
+% bubbleplot3(x,y,z,r,[],[],[],[],'Tag','MyFirstBubbleplot3Bubbles')
+% shading interp; camlight right; lighting phong; view(60,30);
+%
+% See also scatter3 surface surf plot3 scatter
+
+
+% Author: Peter (PB) Bodin
+% Email: pbodin@kth.se
+% Created: 06-Aug-2005 14:02:27
+
+
+msgstruct = nargchk(0,1,nargout,'struct');
+error(msgstruct);
+if nargin > 8
+ [xargs{1:nargin-8}]=varargin{5:end};
+else
+ xargs={};
+end
+if nargin <8 || isempty(varargin{4}) || varargin{4}<7
+ M=10;
+else
+ M=ceil(varargin{4});
+end
+
+if nargin <7 || isempty(varargin{3})
+ N=20;
+else
+ N=floor(varargin{3});
+end
+
+% Use 3 (N+1)x(N+1) arrays to generate the sphere with the largest radii.
+% There´s no need to render the smaller spheres/ellipsoids as detailed as
+% the largest one, scale them linearly against R-max. Let the smallest N
+% be M (default 10).
+
+N=ceil(r/norm(r,2)*N/(max(r/norm(r,2))+eps));
+N(N<M)=M;
+h=zeros(1,numel(x));
+
+% Make sure that x is a column vector, makes it easier to check if the user
+% supplied bubble or ellipsoid radiis for x,y,z
+x=x(:);
+
+% Set up equal axis, set DoubleBuffer on and add a grid
+axis equal;set(gcf,'DoubleBuffer','on');grid on;
+for k=1:numel(x)
+ if size(r)==[numel(x) 3]
+ [X{k},Y{k},Z{k}]=ellipsoid(x(k),y(k),z(k),r(k,1),r(k,2),r(k,3),N(k));
+ else
+ [X{k},Y{k},Z{k}]=ellipsoid(x(k),y(k),z(k),r(k),r(k),r(k),N(k));
+ end
+ h(k)=surface(X{k},Y{k},Z{k},xargs{:});
+end
+
+% Set LineStyle to 'none'
+[lstyle{1:numel(x),1}]=deal('none');
+set(h,{'LineStyle'},lstyle);
+
+% If a color specification is available ...
+if numel(varargin)>0 && ~isempty(varargin{1})
+ try
+ set(h,{'FaceColor'},mat2cell(varargin{1},ones(size(varargin{1},1),1),3));
+ catch
+ [errmsg,errid]=lasterr;
+ error(errmsg);
+ end
+end
+
+% If Alpha is specified ...
+if numel(varargin)>1 && ~isempty(varargin{2})
+ [Alpha_{1:numel(x),1}]=deal(varargin{2});
+ try
+ set(h,{'FaceAlpha'},Alpha_);
+ catch
+ [errmsg,errid]=lasterr;
+ error(errmsg);
+ end
+end
+
+view(3);
+
+% outputs?
+if nargout>0
+ handles=h';
+end
47 center.m
@@ -0,0 +1,47 @@
+function [pos] = center(Pa,Pb,Pc,Type)
+% CENTER calculates and shows the orthocenter, circumcenter, barycenter and
+% incenter of a triangle, given their vertex's coordinates Pa, Pb and Pc
+%
+% Example: center([0 0.5 0], [1 0 0], [1 1 3], 'orthocenter')
+%
+% Made by: Ing. Gustavo Morales, University of Carabobo, Venezuela.
+% 09/14/09
+%
+Pa = Pa(:); Pb = Pb(:); Pc = Pc(:); % Converting to column vectors (if needed)
+AB = Pb - Pa; AC = Pc - Pa; BC = Pc - Pb; % Side vectors
+switch Type
+ case 'incenter'%
+ uab = AB./norm(AB); uac = AC./norm(AC); ubc = BC./norm(BC); uba = -uab;
+ L1 = uab + uac; L2 = uba + ubc; % directors
+ P21 = Pb - Pa;
+ P1 = Pa;
+ case 'barycenter'
+ L1 = (Pb + Pc)/2 -Pa; L2 = (Pa + Pc)/2 - Pb; % directors
+ P21 = Pb - Pa;
+ P1 = Pa;
+ case 'circumcenter'
+ N = cross(AC,AB);
+ L1 = cross(AB,N); L2 = cross(BC,N); % directors
+ P21 = (Pc - Pa)/2;
+ P1 = (Pa + Pb)/2;
+ case 'orthocenter'
+ N = cross(AC,AB);
+ L1 = cross(N,BC); L2 = cross(AC,N); % directors
+ P21 = Pb - Pa;
+ P1 = Pa;
+ otherwise
+ error('Unknown Center Type');
+end
+ML = [L1 -L2]; % Coefficient Matrix
+lambda = ML\P21; % Solving the linear system
+pos = P1 + lambda(1)*L1; % Line Equation evaluated at lambda(1)
+X = [Pa(1); Pb(1); Pc(1)]; Y = [Pa(2); Pb(2); Pc(2)];
+if numel(Pa)== 3 % Tridimensional Case
+ Z = [Pa(3); Pb(3); Pc(3)];
+ %patch(X,Y,Z,'b','FaceAlpha',0.5); hold on;
+ %plot3(pos(1),pos(2),pos(3),'o','Color','r'); view(3);
+else % Bidimensional case
+ %patch(X,Y,'b','FaceAlpha',0.5); hold on;
+ %plot(pos(1),pos(2),'o','Color','r');
+end
+ %grid on; daspect([1 1 1]);
25 check_oszicar.asv
@@ -0,0 +1,25 @@
+function dE = check_oszicar(path)
+
+file = [path '\OSZICAR'];
+fid=fopen(file);
+flag1 = ' ';
+flag2 = ' F=';
+counter = 1;
+
+while (~feof(fid))
+ flag = [flag1 num2str(counter) flag2];
+ line = fgetl(fid);
+ if strmatch(flag,line(1:7))==1
+ [t, r] = strtok(line);
+ for i = 1:6
+ [t,r] = strtok(r)
+ end
+ dE(counter) = str2num(r(3:end));
+ counter = counter + 1;
+ end
+end
+
+figure()
+plot(dE)
+
+end
25 check_oszicar.m
@@ -0,0 +1,25 @@
+function dE = check_oszicar(path)
+
+file = [path '\OSZICAR'];
+fid=fopen(file);
+flag1 = ' ';
+flag2 = ' F=';
+counter = 1;
+
+while (~feof(fid))
+ flag = [flag1 num2str(counter) flag2];
+ line = fgetl(fid);
+ if strmatch(flag,line(1:7))==1
+ [t, r] = strtok(line);
+ for i = 1:6
+ [t,r] = strtok(r)
+ end
+ dE(counter) = str2num(r(3:end));
+ counter = counter + 1;
+ end
+end
+
+figure()
+semilogy(abs(dE))
+
+end
43 check_outcar.m
@@ -0,0 +1,43 @@
+function f = check_outcar(path)
+
+%% using Cartesian
+[xyz, uv] = process(path, 'p');
+all = view_cell(path, xyz, uv, [1 1 1]);
+xyz_axes = [min(all(:,1))-3,max(all(:,1))+3, min(all(:,2))-3, max(all(:,2))+3, min(all(:,3))-3, max(all(:,3))+3]
+[l,w] = size(all)
+
+xyz_axes = [-30,0,0,16,0,12]
+
+file = [path '\OUTCAR'];
+fid=fopen(file);
+flag = ' POSITION TOTAL-FORCE (eV/Angst)';
+match = 0;
+coord_list = all;
+clf('reset')
+counter = 0;
+
+
+while (~feof(fid))
+ line = fgetl(fid)
+ if strmatch(flag,line)==1
+ counter = counter + 1;
+ fgetl(fid);
+ for i = 1:l
+ coord = fgetl(fid);
+ [x,r] = strtok(coord);
+ [y,r] = strtok(r);
+ [z,r] = strtok(r);
+ coord_list(i,1:3) = [str2num(x) str2num(y) str2num(z)];
+ end
+ clf()
+ plot_all(coord_list);
+ axis(xyz_axes);
+ %camlight right; lighting phong; view(90,30);
+ f(counter) = getframe
+ end
+end
+
+% figure()
+% plot_all(coord_list)
+
+end
BIN  chgcar.mat
Binary file not shown
94 contcar.m
@@ -0,0 +1,94 @@
+function [coord_list, restricted, radii, colors] = contcar()
+
+fid=fopen('CONTCAR.txt');
+unit_vectors = [];
+coord_list = [];
+restricted = [];
+counter = 0;
+
+fgetl(fid); % skip some lines
+fgetl(fid);
+
+%% build the unit cell vectors
+for i=1:3
+ temp = fgetl(fid);
+ [v1, temp] = strtok(temp); %x
+ [v2, temp] = strtok(temp); %y
+ [v3, temp] = strtok(temp); %z
+ unit_vectors = [unit_vectors; str2num(v1) str2num(v2) str2num(v3)];
+end
+
+%% get the coordinates and scale them by the unit vectors
+while(~feof(fid))
+ line = fgetl(fid);
+ if strmatch('Direct',line)==1
+ coord = fgetl(fid);
+ [length,width] = size(coord);
+ while (width>10)
+ formatted_coord = strrep(coord, '-', ' -'); %make spacing correct for delim
+ formatted_coord = strrep(formatted_coord, ' ', ' '); %make spacing correct for delim
+ %formatted_coord = strrep(formatted_coord, 'T', ' 1'); %change to int
+ %formatted_coord = strrep(formatted_coord, 'F', ' 0'); %change to int
+ [x,y,z,r1,r2,r3] = strread(formatted_coord, '%f%f%f%c%c%c', 'delimiter', ' ');
+ coord_list = [coord_list; unit_vectors(1,:)*x + unit_vectors(2,:)*y + unit_vectors(3,:)*z]; %multiply by unit vectors to scale correctly
+ restricted = [restricted; ' ' r1 ' ' r2 ' ' r3]; %store T/F restrictions
+ coord = fgetl(fid);
+ [length,width] = size(coord);
+ end
+ end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% GET RADII FROM POTCAR AND POSCAR %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+[length,width] =size(coord_list);
+poscar = [num2str(coord_list) restricted];
+
+fid3=fopen('POSCAR.txt'); %open poscar
+l = fgetl(fid3);
+for i=2:6 % copy these lines
+ l = fgetl(fid3);
+end
+
+%%grab elements from potcar
+fid2=fopen('POTCAR.txt');
+elements = [];
+while(~feof(fid2))
+ line = fgetl(fid2);
+ [len,wid] = size(line);
+ if wid>6
+ if strmatch('PAW_PBE',line(2:8))
+ [e, junk] = strtok(line(9:end));
+ [junk,check_size] = size(e);
+ if check_size<2
+ e = [e ' '];
+ end
+ elements = [elements; e];
+ end
+ end
+end
+
+% elements
+
+%%grab radii from poscar
+[n_e, junk] = size(elements);
+radii = [];
+colors = [];
+for i=1:n_e
+ [temp,l] = strtok(l);
+ r = r_table(find(e_table==elements(i),1));
+ radii = [radii r.*ones(1,str2num(temp))];
+ rgb = r*[c_table(mod(i,7),1)*ones(str2num(temp),1) c_table(mod(i,7),2)*ones(str2num(temp),1) c_table(mod(i,7),3)*ones(str2num(temp),1)];
+ colors = [colors; rgb];
+end
+
+temp_radii = radii;
+temp_colors = colors;
+
+for i=1:counter
+ radii = [radii temp_radii];
+ colors = [colors; temp_colors];
+end
+
+end
BIN  copper.mat
Binary file not shown
BIN  fig code/benzene_pd.fig
Binary file not shown
11 fig code/blah.m
@@ -0,0 +1,11 @@
+vis_sub = dark_and_vis;
+pot = dark_and_vis(:,1);
+vis_sub(:,1) = dark_and_vis(:,4*2)-dark_and_vis(:,8*2);
+vis_sub(:,2) = dark_and_vis(:,2*2)-dark_and_vis(:,5*2);
+vis_sub(:,3) = dark_and_vis(:,1*2)-dark_and_vis(:,6*2);
+vis_sub(:,4) = dark_and_vis(:,3*2)-dark_and_vis(:,7*2);
+figure()
+hold on
+for i=1:4
+ plot(pot,vis_sub(:,i))
+end
BIN  fig code/light_on_off.fig
Binary file not shown
BIN  fig code/o2.fig
Binary file not shown
BIN  fig code/origin_plots.mat
Binary file not shown
BIN  fig code/photo_curr.fig
Binary file not shown
BIN  fig code/photo_curr.mat
Binary file not shown
66 fig code/polyfit0.m
@@ -0,0 +1,66 @@
+function [p,S] = polyfit0(x,y,n)
+%POLYFIT Fit polynomial to data.
+% POLYFIT0(X,Y,N) finds the coefficients of a polynomial P(X) of
+% degree N that fits the data, P(X(I))~=Y(I), in a least-squares sense.
+% HOWEVER THIS VERSION OF POLYFIT HAS A CURVE THAT PASSES THRU 0
+% POLYFITO HAS NO CONSTANT TERM
+%
+% [P,S] = POLYFIT0(X,Y,N) returns the polynomial coefficients P and a
+% structure S for use with POLYVAL to obtain error estimates on
+% predictions. If the errors in the data, Y, are independent normal
+% with constant variance, POLYVAL will produce error bounds which
+% contain at least 50% of the predictions.
+%
+% The structure S contains the Cholesky factor of the Vandermonde
+% matrix (R), the degrees of freedom (df), and the norm of the
+% residuals (normr) as fields.
+%
+% See also POLY, POLYVAL, ROOTS.
+
+% J.N. Little 4-21-85, 8-23-86; CBM, 12-27-91 BAJ, 5-7-93.
+% Copyright (c) 1984-98 by The MathWorks, Inc.
+% $Revision: 5.9 $ $Date: 1997/11/21 23:40:57 $
+% POLYFIT0 CREATED BY AJUTAN MARCH 6-2000
+
+% The regression problem is formulated in matrix format as:
+%
+% y = V*p or
+%
+% 3 2
+% y = [x x x ] [p3
+% p2
+% p1]
+%
+%
+% where the vector p contains the coefficients to be found. For a
+% 7th order polynomial, matrix V would be:
+%
+% V = [x.^7 x.^6 x.^5 x.^4 x.^3 x.^2 x ];
+
+if ~isequal(size(x),size(y))
+ error('X and Y vectors must be the same size.')
+end
+
+x = x(:);
+y = y(:);
+
+% Construct Vandermonde matrix.
+V(:,n+1) = ones(length(x),1);
+for j = n:-1:1
+ V(:,j) = x.*V(:,j+1);
+end
+% STRIP OFF LAST COL OF ONES TO REMOVE CONSTANT TERM
+V(:,n+1)=[];
+
+% Solve least squares problem, and save the Cholesky factor.
+[Q,R] = qr(V,0);
+p = R\(Q'*y); % Same as p = V\y;
+r = y - V*p;
+p = p.'; % Polynomial coefficients are row vectors by convention.
+
+% S is a structure containing three elements: the Cholesky factor of the
+% Vandermonde matrix, the degrees of freedom and the norm of the residuals.
+
+S.R = R;
+S.df = length(y) - (n+1);
+S.normr = norm(r);
BIN  fig code/smoothed_photocurr_sub_dark.fig
Binary file not shown
9 fig code/temp_fit.asv
@@ -0,0 +1,9 @@
+fit = zeros(100,5);
+x_space = linspace(-1,5);
+counter = 0;
+figure()
+for x = 1:4
+ [coeff,temp] = polyfit0(o2(:,x+counter), o2(:,x+counter+1),1)
+ fit(:,x) = x_space*coeff;
+ counter = counter+1
+end
12 fig code/temp_fit.m
@@ -0,0 +1,12 @@
+o2_fit = zeros(100,5);
+x_space = linspace(-1,5);
+counter = 0;
+figure()
+hold on
+for x = 1:5
+ [coeff,temp] = polyfit0(o2(1:end-1,x+counter), o2(1:end-1,x+counter+1),1)
+ plot(o2(:,x+counter), o2(:,x+counter+1), 'o');
+ o2_fit(:,x) = x_space*coeff;
+ plot(x_space, o2_fit(:,x));
+ counter = counter+1
+end
BIN  figs/bubbles.fig
Binary file not shown
BIN  figs/iso2at219_1x1x1.fig
Binary file not shown
BIN  figs/iso2at219_2x2x1.fig
Binary file not shown
BIN  figs/iso2at219_2x2x2 2.fig
Binary file not shown
BIN  figs/iso2at219_2x2x2.fig
Binary file not shown
BIN  figs/test.fig
Binary file not shown
23 fscript.m
@@ -0,0 +1,23 @@
+function all = fscript(p,sc)
+
+%p = 'C:\Users\Chen\Documents\dca\isochain_loose';
+format long
+[xyz, res, uv] = process(p, 'c');
+figure()
+all = view_cell(p, xyz, uv, sc);
+%%edits%%
+% figure()
+% i = find(all(:,3)>-5&all(:,3)<5);
+% new_all = all(i,:);
+% figure()
+% plot_all(new_all)
+
+i = find(all(:,3)>10&all(:,3)<25);
+new_all = all(i,:);
+figure()
+plot_all(new_all)
+
+%make_pos(all)
+%supercell(p, new_all, uv, supercell)
+
+end
BIN  iso_chain.mat
Binary file not shown
BIN  isochain.fig
Binary file not shown
21 lateral_offset.m
@@ -0,0 +1,21 @@
+function [closest, lateral, dir] = lateral_offset(i, radius, all)
+format short
+everything = all(find(all(:,4)==radius),:);
+[row,col] = find(all(:,4)==radius);
+index = find(all(:,4)==radius);
+everything(:,1) = everything(:,1) - all(i,1);
+everything(:,2) = everything(:,2) - all(i,2);
+everything(:,3) = everything(:,3) - all(i,3);
+v = [(sum(everything(:,1:3).^2,2)).^.5 row];
+[w,j] = sort(v,1,'ascend');
+closest = v(j(2:4),2)
+c = center(all(closest(1),1:3),all(closest(2),1:3),all(closest(3),1:3), 'incenter');
+dir = (transpose(c)-all(i,1:3))*-1;
+offset = sum((transpose(c)-all(i,1:3)).^2)^.5;
+lateral = sum((transpose(c(1:2))-all(i,1:2)).^2)^.5;
+
+all(closest,5:6) = .2;
+figure();
+plot_all(all);
+
+end
81 lateral_script.m
@@ -0,0 +1,81 @@
+ld = zeros(4,4);
+counter = 1;
+
+%%isochain_loose
+sc = [1,1,2];
+p = 'C:\Users\Chen\Documents\dca\isochain_loose';
+format long
+[xyz, res, uv] = process(p, 'c');
+figure()
+all = view_cell(p, xyz, uv, sc);
+i = find(all(:,3)>10&all(:,3)<25);
+new_all = all(i,:);
+figure()
+plot_all(new_all)
+test = new_all;
+test = [new_all(1:65,:); new_all(69:end,:)]
+test = [new_all(66,:); test]
+[c,l,d] = lateral_offset(1,1.28,test);
+ld(counter,:) = [l d];
+counter = counter + 1;
+test(1,:) = new_all(67,:);
+[c,l,d] = lateral_offset(1,1.28,test);
+ld(counter,:) = [l d];
+counter = counter + 1;
+test(1,:) = new_all(68,:);
+[c,l,d] = lateral_offset(1,1.28,test);
+ld(counter,:) = [l d];
+counter = counter + 1;
+
+%%isochain_float
+sc = [1,1,2];
+p = 'C:\Users\Chen\Documents\dca\isochain_float';
+[xyz, res, uv] = process(p, 'c');
+figure()
+all = view_cell(p, xyz, uv, sc);
+i = find(all(:,3)>10&all(:,3)<25);
+new_all = all(i,:);
+figure()
+plot_all(new_all)
+test = new_all;
+test = [new_all(1:65,:); new_all(69:end,:)]
+test = [new_all(66,:); test]
+[c,l,d] = lateral_offset(1,1.28,test);
+ld(counter,:) = [l d];
+counter = counter + 1;
+test(1,:) = new_all(67,:);
+[c,l,d] = lateral_offset(1,1.28,test);
+ld(counter,:) = [l d];
+counter = counter + 1;
+test(1,:) = new_all(68,:);
+[c,l,d] = lateral_offset(1,1.28,test);
+ld(counter,:) = [l d];
+counter = counter + 1;
+
+%% isotrigonal float
+sc = [2,2,1];
+p = 'C:\Users\Chen\Documents\dca\isotrigonal_float';
+[xyz, res, uv] = process(p, 'c');
+figure()
+all = view_cell(p, xyz, uv, sc);
+test = [all(1:526,:); all(528:end,:)];
+test = [all(527,:); test];
+[c,l,d] = lateral_offset(1,1.28,test);
+ld(counter,:) = [l d];
+counter = counter + 1;
+
+
+%% isotrigonal
+sc = [2,2,1];
+p = 'C:\Users\Chen\Documents\dca\isotrigonal';
+[xyz, res, uv] = process(p, 'c');
+figure()
+all = view_cell(p, xyz, uv, sc);
+test = [all(1,:); all(3:end,:)];
+test = [all(2,:); test];
+[c,l,d] = lateral_offset(1,1.28,test);
+ld(counter,:) = [l d];
+counter = counter + 1;
+
+'isochain_loose, isochain_float, isotrigonal, isotrigonal_float'
+ld
102 live/daqstopbutton.m
@@ -0,0 +1,102 @@
+function varargout = daqstopbutton(fig,obj,varargin);
+%DAQSTOPBUTTON Add a stop/start button to a data acquisition application
+%
+% DAQSTOPBUTTON(FIG,OBJ) adds a start/stop button to figure FIG. This
+% button can be used to start and stop data acquisition object OBJ.
+% DAQSTOPBBUTTON will also delete OBJ when FIG is closed (i.e., it sets
+% FIG's CloseRequestFcn to delete the object)
+%
+% DAQSTOPBUTTON(FIG,OBJ,'P1','V1','P2','V2', ...) specifies Property-Value
+% pairs for configuring properties of the start/stop button. Any valid
+% property of a togglebutton can be specified.
+%
+% DAQSTOPBUTTON(HBUTTON,OBJ) turns existing uicontrol HBUTTON into a
+% start/stop button.
+%
+% HBUTTON = DAQSTOPBUTTON(...) returns a handle to the start/stop button
+%
+% Note: This button does not "listen" to determine if your object changes
+% state. For instance, if you have a finite number of samples
+% acquired after starting, the button will not reset automatically
+% when your object stops running. Don't despair, since it is easy
+% enough to do on your own!
+%
+% Example:
+% fh = figure; % Create a figure
+% ai = analoginput('winsound'); % Create an input object
+% addchannel(ai,1); % Add a channel
+% set(ai,'TriggerRepeat',inf); % Configure to run infinitely
+% set(ai,'TimerFcn','plot(peekdata(ai,500))'); % Each timer event will plot recent data
+% hButton = daqstopbutton(fh,ai); % Add the stopbutton
+
+% Scott Hirsch
+% shirsch@mathworks.com
+% Copyright 2003 The MathWorks, Inc.
+
+%Error checking
+msg = nargchk(2,inf,nargin);
+error(msg)
+
+if ~all(isvalid(obj))
+ error('Second input argument must be valid daq object')
+end;
+
+if ~ishandle(fig)
+ error('First input argument must be handle to a figure or a uicontrol');
+end;
+
+% Check if user input handle to figure or handle to button
+switch get(fig,'Type')
+ case 'figure'
+
+ % Create the button
+ hButton = uicontrol(fig,'style','togglebutton',varargin{:});
+ case 'uicontrol'
+ hButton = fig;
+ fig = get(hButton,'Parent');
+ % In R14, it's possible that fig would return a handle to a panel, not a figure
+ if ~strcmp(get(fig,'Type'),'figure')
+ fig = get(fig,'Parent');
+ end;
+end;
+
+%Check current state of the object
+val = strcmp(obj.Running,'On');
+
+if val==1 %Already running
+ set(hButton,'String','Stop');
+ set(hButton,'Value',1)
+else
+ set(hButton,'String','Start');
+ set(hButton,'Value',0)
+end;
+
+set(hButton,'Callback',{@localStartStopObject,obj})
+
+% Configure the CloseRequestFcn of the figure
+setappdata(fig,'DaqStopButtonObject',obj)
+cr = get(fig,'CloseRequestFcn');
+cr = ['obj=getappdata(gcf,''DaqStopButtonObject'');stop(obj);delete(obj);' cr];
+set(fig,'CloseRequestFcn',cr);
+
+% Return a handle to the button
+if nargout
+ varargout{1} = hButton;
+end;
+
+
+function localStartStopObject(hButton,action,obj)
+% Callback for the start/stop button
+val = get(hButton,'Value');
+if val==1 %Pushed in
+ set(hButton,'String','Stop');
+ if all(isvalid(obj))
+ start(obj)
+ end;
+else
+ set(hButton,'String','Start');
+ if all(isvalid(obj))
+ stop(obj)
+ end;
+end;
+
57 live/demo_tmwidgets.m
@@ -0,0 +1,57 @@
+function demo_tmwidgets
+%% DEMO_TMWIDGETS Demonstration of the TM Widgets with Analog Input
+
+% Copyright 2004 The MathWorks, Inc
+
+%% Initialize scopes
+% Initialize the scopes by creating an axes for each scope and calling the
+% initialization routines to configure the scopes.
+Fs = 44100; % Sample rate
+Nfft = 1024; % FFT Length
+ax1 = subplot(411); % Strip chart
+ax2 = subplot(4,2,[3 5 7]); % Spectrumscope
+ax3 = subplot(4,2,[4 6 8]); % Specgramscope
+stripchart(ax1,Fs,1);ylim(ax1,[-.01 .01])
+spectrumscope(ax2,Fs,Nfft);
+specgramscope(ax3,Fs,Nfft);
+
+%% Initialize Hardware
+% Initialize the sound card for a single channel of continuous acquisition.
+% Use a very fast timer to have the display run as fast as possible.
+% Program TimerFcn to send current data to scopes. Program the figure's
+% CloseRequestFcn to stop and destroy the analoginput object. These two
+% functions are included as nested functions inside the main function.
+% This gives them visibility into all data from the main function,
+% including axes handles.
+ai = analoginput('winsound');
+addchannel(ai,1);
+set(ai,'SampleRate',Fs)
+set(ai,'TimerPeriod',.01);
+set(ai,'TriggerRepeat',inf)
+set(ai,'TimerFcn',@aicallback);
+set(gcf,'CloseRequestFcn',@figclosefcn);
+
+%% Start Hardware
+start(ai)
+
+ % AICALLBACK (Nested Function)
+ % Execute upon Timer event. Update scopes with current data.
+ function aicallback(ai,event)
+ if get(ai,'SamplesAvailable') < get(ai,'SamplesPerTrigger'), return, end;
+ d=peekdata(ai,ai.SamplesPerTrigger);
+ stripchart(ax1,d);
+ spectrumscope(ax2,d);
+ specgramscope(ax3,d)
+ end % END AICALLBACK
+
+ % FIGCLOSEFCN (Nested Function)
+ % Execute upon closing figure. Stop and delete AI objects.
+ function figclosefcn(obj,event)
+ try
+ stop(ai);
+ delete(ai);
+ end
+ closereq
+ end % END FIGCLOSEFCN
+
+end % END DEMO_TMWIDGETS (Main function)
28 live/license.txt
@@ -0,0 +1,28 @@
+Copyright (c) 2004, The MathWorks, Inc.
+Copyright (c) 1997, Christophe COUVREUR
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in
+ the documentation and/or other materials provided with the distribution
+ * Neither the name of the The MathWorks, Inc. nor the names
+ of its contributors may be used to endorse or promote products derived
+ from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.
289 live/specgramscope.m
@@ -0,0 +1,289 @@
+function varargout = specgramscope(varargin)
+% SPECGRAMSCOPE Live updating spectrogram display
+%
+% STEP 1: Initialize the scope
+% SPECGRAMSCOPE(FS,NFFT) initializes a spectrogram scope in the current axes.
+% This spectrogram scope will compute and displays the NFFT-point FFT of a vector
+% signal with sample rate FS Hz. It will show the last 10 FFTs
+%
+% STEP 2: Update the scope
+% SPECGRAMSCOPE(S) updates the spectrogram scope in the current axes with the
+% FFT of vector S. The scope should first be initialized as above with
+% sample rate and FFT length. If not, the sample rate will be 1 Hz and the FFT
+% length will be the length of S. Differences between the length of S and
+% the specified FFT length are handled the same as MATLAB's built-in FFT
+% function (i.e., zero-padding or truncation, as appropriate).
+%
+% SPECGRAMSCOPE(FS,NFFT,NTRACES) initializes a spectrogram scope in the
+% current axes with NTRACES traces. This is how many records to keep in
+% time domain. Default = 10
+%
+% SPECGRAMSCOPE(HAX, ...) defines the scope in specified axes HAX instead of GCA. i.e.,
+% SPECGRAMSCOPE(HAX,FS,NFFT) initializes axes HAX as a spectrogram scope, and
+% SPECGRAMSCOPE(HAX,S) updates axes HAX with vector S.
+%
+% HAX = SPECGRAMSCOPE(...) returns a handle to the axes initialized by the
+% spectrogram scope. This is useful if you allow SPECGRAMSCOPE to create an
+% axes for you, and want to be able to easily reference the axes for
+% updates. The surface created by SPECGRAMSCOPE all have the tag
+% 'SpecgramScope'. If you would like to manually modify the properties of
+% these lines, their handles can be found by:
+%
+% HAX = SPECGRAMSCOPE(...);
+% HSurf = findobj(HAX,'Tag','SpecgramScope');
+%
+% Example
+% %% Initialize data
+% Fs = 16384;
+% Nfft = 2048;
+% t = (0:1:Nfft-1)'/Fs;
+% fo = logspace(3.5,3.7); % Range of fundamental frequencies
+% s1 = sin(2*pi*t*fo) + .1*rand(Nfft,length(fo));
+%
+%
+% %% Initialize scope
+% specgramscope(Fs,Nfft,30);
+% view([103 30])
+%
+% %% Update scope
+% for ii = 1:length(fo)
+% specgramscope(s1(:,ii));
+% drawnow;pause(.01);
+% end;
+
+% Scott Hirsch 5-05.
+% shirsch@mathworks.com
+% Copyright 2004-2005 The MathWorks, Inc.
+
+%% Parse input arguments
+% Decision tree:
+% + Initialize or update?
+% o If update -> OK
+% o If initialize -> Axes specified, or use GCA?
+
+error(nargchk(1,4,nargin))
+NTracesDefault = 10;
+
+%% Initialize or update?
+% If first or second input argument is not a scalar, it must be data - i.e. we are
+% updating
+
+if prod(size(varargin{1})) > 1 | prod(size(varargin{2})) > 1 % Update
+ action = 'update';
+
+ if nargin==1 % Use current axes
+ hAxes = gca;
+ data = varargin{1};
+ else
+ hAxes = varargin{1}; % Axes was specified
+ data = varargin{2};
+ end;
+
+ % If the user has not initialized this scope, do it for them
+ parms = getappdata(hAxes,'SpecgramScopeParameters');
+
+ % Ensure that scope has been initialized
+ if isempty(parms)
+ % Use default values
+ Fs = 1;
+ data = rowmajor(data);
+ Nfft = length(data);
+ NTraces = NTracesDefault; % Number of time histories
+ feval(mfilename,hAxes,Fs,Nfft,NTraces); % This recursive call will initialize the scope
+ % Get the new parameter structure
+ parms = getappdata(hAxes,'SpecgramScopeParameters');
+ end;
+
+
+
+else % Initialize
+ action = 'init';
+
+ if ~isaxes(varargin{1}) % Easy mode, no handle passed in
+ % Use current axes
+ hAxes = gca;
+ Fs = varargin{1};
+ Nfft = varargin{2};
+ if nargin==3
+ NTraces = varargin{3};
+ else
+ NTraces = NTracesDefault;
+ end;
+
+ else % Expert mode, passed handle in
+ hAxes = varargin{1};
+ Fs = varargin{2};
+ Nfft = varargin{3};
+ if nargin==4
+ NTraces = varargin{4};
+ else
+ NTraces = NTracesDefault;
+ end;
+ end;
+end;
+
+%% Dole out the work
+%
+switch action
+ case 'init' % Initialize
+
+ % Build structure to internally pass information
+ parms.Fs = Fs; % Sample Rate
+ parms.NTraces = NTraces; % Number of records in time
+ parms.hAxes = hAxes; % Handle to axes
+ parms.Nfft = Nfft; % FFT Block size
+
+ % Store parameter structure
+ setappdata(hAxes,'SpecgramScopeParameters',parms);
+
+ localInitScope(parms) % Initialize scope
+
+ case 'update' % Update
+ parms = getappdata(hAxes,'SpecgramScopeParameters');
+
+ % Error checking
+ % Ensure that scope has been initialized. This shouldn't slip
+ % through to here.
+ if isempty(parms)
+ error(['The spectrogram scope must first be initialized ' ...
+ 'with the sample rate: specgramscope(hAxes,Fs)']);
+ end;
+
+ % Force data to be in columns. Allow for multiple columns. This will
+ % error if data actually has more channels than samples.
+ data = rowmajor(data);
+
+ % Check that the number of columns corresponds to the number of lines
+ nc = size(data,2); % Number of columns
+ if nc ~= 1
+ error(['spectrogram scope requires a single column vector of data']);
+ end;
+
+ localUpdateScope(data,parms) % Update the scope
+end;
+
+% Return appropriate output argument
+if nargout
+ varargout{1} = parms.hAxes;
+end;
+
+
+% ***********************************************************************
+% Initialize the Scope
+function localInitScope(parms)
+
+% Set axes
+f = (0:parms.Nfft/2-1)*parms.Fs/parms.Nfft;
+f = f(:);
+
+% Add surface
+[X,Y] = meshgrid(-parms.NTraces+1:0,f);
+parms.hSurf = surf(parms.hAxes,X,Y,NaN*ones(length(f),parms.NTraces), ...
+ 'Tag','SpecgramScope');
+
+set(parms.hAxes, ...
+ 'XLim',[-parms.NTraces+1 0], ...
+ 'YLim',[0 f(end)]);%, ...
+% 'YDir','reverse'); % Reverse frequency axes direction
+shading(parms.hAxes,'interp');
+
+setappdata(parms.hAxes,'SpecgramScopeParameters',parms);
+
+%% Get handle to the figure
+% Turn doublebuffer on to eliminate flickering
+hFig = get(parms.hAxes,'Parent');
+
+% In R14, it's possible that hFig would return a handle to a panel, not a
+% figure
+if ~strcmp(get(hFig,'Type'),'figure')
+ hFig = get(hFig,'Parent');
+end;
+
+%%
+% Label the plot.
+% There's a bug in R13 when creating xlabel and ylabel with direct
+% parenting - the alignment gets all messed up. Instead, make hAx current
+% axes
+ca = gca;
+set(hFig,'CurrentAxes',parms.hAxes);
+xlabel('History')
+ylabel('Frequency (Hz)');
+zlabel('Magnitude (dB)');
+set(hFig,'CurrentAxes',ca);
+view([103 30])
+
+%%
+% Turn doublebuffer on to eliminate flickering
+set(hFig,'DoubleBuffer','on');
+
+% ***********************************************************************
+% Update the plot.
+function localUpdateScope(data,parms)
+
+[f,mag] = localfft(data,parms);
+
+% Dynamically modify Magnitude axis as we go. Expand, but don't shrink.
+maxM=max(mag(:));
+minM=min(mag(:));
+yax2=get(parms.hAxes,'YLim');
+if minM<yax2(1),
+ yax2(1)=minM;
+end
+if maxM>yax2(2),
+ yax2(2)=maxM;
+end
+set(parms.hAxes,'YLim',yax2)
+
+
+% Update the plot
+hSurf = parms.hSurf;
+zd = get(hSurf,'ZData');
+zd = [mag zd(:,1:end-1)];
+set(hSurf,'ZData',zd,'CData',zd)
+
+
+% set(parms.hLine, 'XData', f(:,1), 'YData', mag(:,1));
+% set(parms.hLine, {'YData'}, Mag');
+
+% Note: It looks like it's faster to update one line at a time
+% in a loop than to update with a cell array
+
+% ***********************************************************************
+% Calculate the fft of the data.
+function [f, mag] = localfft(data,parms)
+
+% Calculate the fft of the data.
+xfft = 2/parms.Nfft*fft(data,parms.Nfft);
+
+% Avoid taking the log of 0.
+xfft(xfft == 0) = 1e-17;
+
+% Compute magnitude, dB
+mag = 20*log10(abs(xfft(1:parms.Nfft/2,:)));
+
+f = (0:length(mag)-1)*parms.Fs/parms.Nfft;
+f = f(:);
+
+% ***********************************************************************
+% Utility - isaxes
+function truefalse = isaxes(h);
+% ISAXES(H) True if H is a handle to a valid axes
+
+truefalse = 0; % Start false
+if ishandle(h)
+ if strcmp('axes',get(h,'Type'))
+ truefalse = 1;
+ end;
+end;
+
+% ***********************************************************************
+% Utility - rowmajor
+function data = rowmajor(data);
+% Force data to be row major. i.e. more rows than columns
+
+[nr,nc] = size(data);
+if nc>nr
+ data = data';
+ [nr,nc] = size(data);
+end;
+
61 live/specgramscope_documentation.m
@@ -0,0 +1,61 @@
+%% SPECGRAMSCOPE Live updating spectrogram display
+% SPECGRAMSCOPE makes it fairly easy to include a spectrogram scope in your
+% real-time data acquisition and analysis application. You feed
+% spectrumscope your data, and it plots the FFT - simple enough! It takes
+% 2 steps to use SPECGRAMSCOPE. First, you initialize the scope with basic
+% information needed for the FFT (sample rate, fft length, and history length).
+% After that, all you need to do is pass your data to the scope.
+%
+% This documentation starts with the simplest syntax for the two steps,
+% then provides a few more advanced options.
+
+%% STEP 1: Initialize the scope
+% SPECGRAMSCOPE(FS,NFFT) initializes a spectrogram scope in the current axes.
+% This spectrogram scope will compute and displays the NFFT-point FFT of a vector
+% signal with sample rate FS Hz. It will show the last 10 FFTs
+
+%% STEP 2: Update the scope
+% SPECGRAMSCOPE(S) updates the spectrogram scope in the current axes with the
+% FFT of vector S. The scope should first be initialized as above with
+% sample rate and FFT length. If not, the sample rate will be 1 Hz and the FFT
+% length will be the length of S. Differences between the length of S and
+% the specified FFT length are handled the same as MATLAB's built-in FFT
+% function (i.e., zero-padding or truncation, as appropriate).
+
+%% Specifying the axes to locate the scope
+% SPECGRAMSCOPE(HAX, ...) defines the scope in specified axes HAX instead of GCA. i.e.,
+% SPECGRAMSCOPE(HAX,FS,NFFT) initializes axes HAX as a spectrogram scope, and
+% SPECGRAMSCOPE(HAX,S) updates axes HAX with vector S.
+%
+
+%% Returning a handle to the scope
+% HAX = SPECGRAMSCOPE(...) returns a handle to the axes initialized by the
+% spectrogram scope. This is useful if you allow SPECGRAMSCOPE to create an
+% axes for you, and want to be able to easily reference the axes for
+% updates. The surface created by SPECGRAMSCOPE all have the tag
+% 'SpecgramScope'. If you would like to manually modify the properties of
+% these lines, their handles can be found by:
+%
+% HAX = SPECGRAMSCOPE(...);
+% HSurf = findobj(HAX,'Tag','SpecgramScope');
+
+%% Example
+%% Initialize data
+Fs = 16384;
+Nfft = 2048;
+t = (0:1:Nfft-1)'/Fs;
+fo = logspace(3.5,3.7); % Range of fundamental frequencies
+s1 = sin(2*pi*t*fo) + .1*rand(Nfft,length(fo));
+
+
+%% Initialize scope
+specgramscope(Fs,Nfft,30);
+view([103 30])
+
+%% Update scope
+for ii = 1:length(fo)
+ specgramscope(s1(:,ii));
+ drawnow;pause(.01);
+end;
+
+% Copyright 2004 The MathWorks, Inc
295 live/spectrumscope.m
@@ -0,0 +1,295 @@
+function varargout = spectrumscope(varargin)
+% SPECTRUMSCOPE Compute and display a real-time FFT
+%
+% SPECTRUMSCOPE makes it fairly easy to include a spectrum scope in your
+% real-time data acquisition and analysis application. You feed
+% spectrumscope your data, and it plots the FFT - simple enough! It takes
+% 2 steps to use SPECTRUMSCOPE. First, you initialize the scope with basic
+% information needed for the FFT (sample rate, fft length, and number of traces).
+% After that, all you need to do is pass your data to the scope.
+%
+% This documentation starts with the simplest syntax for the two steps,
+% then provides a few more advanced options.
+%
+% Step 1: Initialize the scope
+% SPECTRUMSCOPE(FS,NFFT) initializes a spectrum scope in the current axes.
+% This spectrum scope will compute and displays the NFFT-point FFT of a vector
+% signal with sample rate FS Hz.
+%
+% Step 2: Update the scope
+% SPECTRUMSCOPE(S) updates the spectrum scope in the current axes with the
+% FFT of vector S. The scope should first be initialized as above with
+% sample rate and FFT length. If not, the sample rate will be 1 Hz and the FFT
+% length will be the length of S. Differences between the length of S and
+% the specified FFT length are handled the same as MATLAB's built-in FFT
+% function (i.e., zero-padding or truncation, as appropriate).
+%
+% SPECTRUMSCOPE(FS,NFFT,NTRACES) initializes a spectrum scope in the
+% current axes with NTRACES traces. A trace is a single line on the scope;
+% typically one will display one trace per channel of data. The default
+% for NTRACES is 1. To update a spectrum scope with multiple traces,
+% SPECTRUMSCOPE(S) must specify a matrix S with shorter dimension length =
+% NTRACES. SPECTRUMSCOPE computes the FFT along the longer dimension,
+% assuming the shorter dimension corresponds to traces. I know that this
+% differs from FFT (which always defaults to applying to each column), but
+% I find this deviation convenient.
+%
+% SPECTRUMSCOPE(HAX, ...) defines the scope in specified axes HAX instead of GCA. i.e.,
+% SPECTRUMSCOPE(HAX,FS,NFFT) initializes axes HAX as a spectrum scope, and
+% SPECTRUMSCOPE(HAX,S) updates axes HAX with vector S.
+%
+% HAX = SPECTRUMSCOPE(...) returns a handle to the axes initialized by the
+% spectrum scope. This is useful if you allow SPECTRUMSCOPE to create an
+% axes for you, and want to be able to easily reference the axes for
+% updates. The lines created by SPECTRUMSCOPE all have the tag
+% 'SpectrumScope'. If you would like to manually modify the properties of
+% these lines, their handles can be found by:
+%
+% HAX = SPECTRUMSCOPE(...);
+% HLINE = findobj(HAX,'Tag','SpectrumScope');
+%
+% Example
+% %% Create data
+% Fs = 1024;
+% Nfft = 2048;
+% t = (0:1:Nfft-1)'/Fs;
+% fo = 100:5:300; % Range of fundamental frequencies
+% s1 = sin(2*pi*t*fo);
+%
+% %% Initialize scope
+% spectrumscope(Fs,Nfft);
+%
+% %% Update scope
+% for ii = 1:length(fo)
+% spectrumscope(s1(:,ii));
+% drawnow;pause(.01);
+% end;
+
+% Scott Hirsch 2-25-04
+% shirsch@mathworks.com
+% Copyright 2004 The MathWorks, Inc.
+
+%% Parse input arguments
+% Decision tree:
+% + Initialize or update?
+% o If update -> OK
+% o If initialize -> Axes specified, or use GCA?
+
+error(nargchk(1,4,nargin))
+
+%% Initialize or update?
+% If first or second input argument is not a scalar, it must be data - i.e. we are
+% updating
+
+if prod(size(varargin{1})) > 1 | prod(size(varargin{2})) > 1 % Update
+ action = 'update';
+
+ if nargin==1 % Use current axes
+ hAxes = gca;
+ data = varargin{1};
+ else
+ hAxes = varargin{1}; % Axes was specified
+ data = varargin{2};
+ end;
+
+ % If the user has not initialized this scope, do it for them
+ parms = getappdata(hAxes,'SpectrumScopeParameters');
+
+ % Ensure that scope has been initialized
+ if isempty(parms)
+ % Use default values
+ Fs = 1;
+ data = rowmajor(data);
+ [Nfft,NTraces] = size(data);
+ feval(mfilename,hAxes,Fs,Nfft,NTraces); % This recursive call will initialize the scope
+ % Get the new parameter structure
+ parms = getappdata(hAxes,'SpectrumScopeParameters');
+ end;
+
+
+
+else % Initialize
+ action = 'init';
+
+ if ~isaxes(varargin{1}) % Easy mode, no handle passed in
+ % Use current axes
+ hAxes = gca;
+ Fs = varargin{1};
+ Nfft = varargin{2};
+ if nargin==3
+ NTraces = varargin{3};
+ else
+ NTraces = 1;
+ end;
+
+ else % Expert mode, passed handle in
+ hAxes = varargin{1};
+ Fs = varargin{2};
+ Nfft = varargin{3};
+ if nargin==4
+ NTraces = varargin{4};
+ else
+ NTraces = 1;
+ end;
+ end;
+end;
+
+%% Dole out the work
+%
+switch action
+ case 'init' % Initialize
+
+ % Build structure to internally pass information
+ parms.Fs = Fs; % Sample Rate
+ parms.NTraces = NTraces; % Number of lines in plot
+ parms.hAxes = hAxes; % Handle to axes
+ parms.Nfft = Nfft; % FFT Block size
+
+ % Store parameter structure
+ setappdata(hAxes,'SpectrumScopeParameters',parms);
+
+ localInitScope(parms) % Initialize scope
+
+ case 'update' % Update
+ parms = getappdata(hAxes,'SpectrumScopeParameters');
+
+ % Error checking
+ % Ensure that scope has been initialized. This shouldn't slip
+ % through to here.
+ if isempty(parms)
+ error(['The spectrum scope must first be initialized ' ...
+ 'with the sample rate: spectrumscope(hAxes,Fs)']);
+ end;
+
+ % Force data to be in columns. Allow for multiple columns. This will
+ % error if data actually has more channels than samples.
+ data = rowmajor(data);
+
+ % Check that the number of columns corresponds to the number of lines
+ nc = size(data,2); % Number of columns
+ if nc ~= parms.NTraces
+ error(['Size mismatch. You initialized spectrumscope with ' num2str(parms.NTraces) ...
+ ' lines, but just passed in ' num2str(nc) ' channels of data. These' ...
+ ' numbers must be the same.']);
+ end;
+
+ localUpdateScope(data,parms) % Update the scope
+end;
+
+% Return appropriate output argument
+if nargout
+ varargout{1} = parms.hAxes;
+end;
+
+
+% ***********************************************************************
+% Initialize the Scope
+function localInitScope(parms)
+
+% Set axes
+f = (0:parms.Nfft/2-1)*parms.Fs/parms.Nfft;
+f = f(:);
+
+% Add line(s)
+parms.hLine = plot(f,NaN*ones(length(f),parms.NTraces), ...
+ 'Tag','SpectrumScope', ...
+ 'Parent',parms.hAxes);
+set(parms.hAxes,'XLim',[0 parms.Fs/2]);
+setappdata(parms.hAxes,'SpectrumScopeParameters',parms);
+
+%% Get handle to the figure
+% Turn doublebuffer on to eliminate flickering
+hFig = get(parms.hAxes,'Parent');
+
+% In R14, it's possible that hFig would return a handle to a panel, not a
+% figure
+if ~strcmp(get(hFig,'Type'),'figure')
+ hFig = get(hFig,'Parent');
+end;
+
+%%
+% Label the plot.
+% There's a bug in R13 when creating xlabel and ylabel with direct
+% parenting - the alignment gets all messed up. Instead, make hAx current
+% axes
+ca = gca;
+set(hFig,'CurrentAxes',parms.hAxes);
+xlabel('Frequency (Hz)');
+ylabel('Magnitude (dB)');
+set(hFig,'CurrentAxes',ca);
+
+%%
+% Turn doublebuffer on to eliminate flickering
+set(hFig,'DoubleBuffer','on');
+
+% ***********************************************************************
+% Update the plot.
+function localUpdateScope(data,parms)
+
+[f,mag] = localfft(data,parms);
+
+% Dynamically modify Magnitude axis as we go. Expand, but don't shrink.
+maxM=max(mag(:));
+minM=min(mag(:));
+yax2=get(parms.hAxes,'YLim');
+if minM<yax2(1),
+ yax2(1)=minM;
+end
+if maxM>yax2(2),
+ yax2(2)=maxM;
+end
+set(parms.hAxes,'YLim',yax2)
+
+
+hLine = parms.hLine;
+
+% Update the plot
+for ii=1:parms.NTraces
+% Mag{ii} = mag(:,ii);
+ set(hLine(ii), 'YData', mag(:,ii));
+end;
+% set(parms.hLine, 'XData', f(:,1), 'YData', mag(:,1));
+% set(parms.hLine, {'YData'}, Mag');
+
+% Note: It looks like it's faster to update one line at a time
+% in a loop than to update with a cell array
+
+% ***********************************************************************
+% Calculate the fft of the data.
+function [f, mag] = localfft(data,parms)
+
+% Calculate the fft of the data.
+xfft = 2/parms.Nfft*fft(data,parms.Nfft);
+
+% Avoid taking the log of 0.
+xfft(xfft == 0) = 1e-17;
+
+% Compute magnitude, dB
+mag = 20*log10(abs(xfft(1:parms.Nfft/2,:)));
+
+f = (0:length(mag)-1)*parms.Fs/parms.Nfft;
+f = f(:);
+
+% ***********************************************************************
+% Utility - isaxes
+function truefalse = isaxes(h);
+% ISAXES(H) True if H is a handle to a valid axes
+
+truefalse = 0; % Start false
+if ishandle(h)
+ if strcmp('axes',get(h,'Type'))
+ truefalse = 1;
+ end;
+end;
+
+% ***********************************************************************
+% Utility - rowmajor
+function data = rowmajor(data);
+% Force data to be row major. i.e. more rows than columns
+
+[nr,nc] = size(data);
+if nc>nr
+ data = data';
+ [nr,nc] = size(data);
+end;
+
71 live/spectrumscope_documentation.m
@@ -0,0 +1,71 @@
+%% SPECTRUMSCOPE Compute and display a real-time FFT
+% SPECTRUMSCOPE makes it fairly easy to include a spectrum scope in your
+% real-time data acquisition and analysis application. You feed
+% spectrumscope your data, and it plots the FFT - simple enough! It takes
+% 2 steps to use SPECTRUMSCOPE. First, you initialize the scope with basic
+% information needed for the FFT (sample rate, fft length, and number of traces).
+% After that, all you need to do is pass your data to the scope.
+%
+% This documentation starts with the simplest syntax for the two steps,
+% then provides a few more advanced options.
+
+%% Step 1: Initialize the scope
+% SPECTRUMSCOPE(FS,NFFT) initializes a spectrum scope in the current axes.
+% This spectrum scope will compute and displays the NFFT-point FFT of a vector
+% signal with sample rate FS Hz.
+
+%% Step 2: Update the scope
+% SPECTRUMSCOPE(S) updates the spectrum scope in the current axes with the
+% FFT of vector S. The scope should first be initialized as above with
+% sample rate and FFT length. If not, the sample rate will be 1 Hz and the FFT
+% length will be the length of S. Differences between the length of S and
+% the specified FFT length are handled the same as MATLAB's built-in FFT
+% function (i.e., zero-padding or truncation, as appropriate).
+
+%% Including multiple traces
+% SPECTRUMSCOPE(FS,NFFT,NTRACES) initializes a spectrum scope in the
+% current axes with NTRACES traces. A trace is a single line on the scope;
+% typically one will display one trace per channel of data. The default
+% for NTRACES is 1. To update a spectrum scope with multiple traces,
+% SPECTRUMSCOPE(S) must specify a matrix S with shorter dimension length =
+% NTRACES. SPECTRUMSCOPE computes the FFT along the longer dimension,
+% assuming the shorter dimension corresponds to traces. I know that this
+% differs from FFT (which always defaults to applying to each column), but
+% I find this deviation convenient.
+
+%% Specifying the axes to locate the scope
+% SPECTRUMSCOPE(HAX, ...) defines the scope in specified axes HAX instead of GCA. i.e.,
+% SPECTRUMSCOPE(HAX,FS,NFFT) initializes axes HAX as a spectrum scope, and
+% SPECTRUMSCOPE(HAX,S) updates axes HAX with vector S.
+
+%% Returning a handle to the scope
+% HAX = SPECTRUMSCOPE(...) returns a handle to the axes initialized by the
+% spectrum scope. This is useful if you allow SPECTRUMSCOPE to create an
+% axes for you, and want to be able to easily reference the axes for
+% updates. The lines created by SPECTRUMSCOPE all have the tag
+% 'SpectrumScope'. If you would like to manually modify the properties of
+% these lines, their handles can be found by:
+%
+% HAX = SPECTRUMSCOPE(...);
+% HLINE = findobj(HAX,'Tag','SpectrumScope');
+
+%% Example
+% Create data
+Fs = 1024;
+Nfft = 2048;
+t = (0:1:Nfft-1)'/Fs;
+fo = 100:5:300; % Range of fundamental frequencies
+s1 = sin(2*pi*t*fo);
+
+%%
+% Initialize scope
+spectrumscope(Fs,Nfft);
+