-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path3DSimulationsv5.ipf
6805 lines (6176 loc) · 292 KB
/
3DSimulationsv5.ipf
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=1 // Use modern global access method.
#include "opticalconstantsDB"
#include <ImageSlider>
#include <3DWaveDisplay>
structure ThreeDSystem
// parameter inputs
variable rot // boolean saying if we are rotating the system or not (90 degrees is always included)
variable num // resolution of simulation (in each of the long axes of the film)
variable thickness // voxels of thickness along the beam direction
variable voxelsize // the size of each pixel in the simulation (in nanometers)
variable size // the dominant size of features (ie the radius of the sphere, if we're simulating spheres)
variable materialnum // number of materials
string paramstring //specific parameters for creating density matricies according to the model
string modelname //store the model name somewhere
string path // the path to put important things for this calculation
string name // name of simulation
variable movie // are we doing a movie
// morphology outputs, alignment inputs
//technically not needed by simulations, but nice to have around for visualizations
wave density1 // the density of material 1 output from the morphology module, this will be changed into
wave density2 // optional density of material 2 (implies 3 materials)
wave density3 // implies 4 materials
wave density4 // implies 5 materials
wave density5 // cannot be created by morphology, but only as the remnant of the other materials (to make total density everywhere 1)
// Inputs for alignment calculations
string materialalignmentstring // list of materials alignment (for creating alignment)
// Morphology or alignment Outputs, Scattering Inputs
wave m1 //4 dimensional density of material 1 in (x direction, y direction, z direction, unaligned) sum = relative concentration sum(between -1 - 1 for x, y and z)^2 + unaligned = total material density
wave m2 //4 dimensional density of material 2 in (x direction, y direction, z direction, unaligned) sum = relative concentration sum(between -1 - 1 for x, y and z)^2 + unaligned = total material density
wave m3 //4 dimensional density of material 3 in (x direction, y direction, z direction, unaligned) sum = relative concentration sum(between -1 - 1 for x, y and z)^2 + unaligned = total material density
wave m4 //4 dimensional density of material 4 in (x direction, y direction, z direction, unaligned) sum = relative concentration sum(between -1 - 1 for x, y and z)^2 + unaligned = total material density
wave m5 //4 dimensional density of material 5 in (x direction, y direction, z direction, unaligned) sum = relative concentration sum(between -1 - 1 for x, y and z)^2 + unaligned = total material density
// parameter inputs for scattering
string materials // list of materials
string efield // two components which add to 1 which indicate the unit direction of the efield in the (ydirection, z direction)
// x-ray propogation is always in the x direction
// Internal Variables
wave /t logbook // keep a log of progress in the simulation
variable progressbar
variable timer
variable timesofar
// Scattering internal calculations
wave /d/c n //a single complex number for each material which will be set at each energy step, and used in calculations
wave /d /c px // induced polarization in the x direction for real space, x,y,z
wave /d /c pxFFT // fourier transform of above, in qx, qy, qz space
wave /d /c py // same in y direction
wave /d /c pyFFT// fourier transform of above, in qx, qy, qz space
wave /d /c pz // same in z direction
wave /d /c pzFFT// fourier transform of above, in qx, qy, qz space
wave /d /c Escatx, Escaty, Escatz //
wave /d EscatSqr // the sum of the magsquared values as described above, before interpolation
// scattering outputs for each step
wave/d scatter3D // 2d scattering simulation from 3D data
wave/d int3D // radially summed intensity
wave/d int3Dp0 // radially averaged intensity
wave/d int3Dpara // radially summed intensity in direction of electric field
wave/d int3Dp0para // radially averaged intensity in direction of electric field
wave/d int3Dperp // radially summed intensity in direction normal to the electric field
wave/d int3Dp0perp // radially averaged intensity in direction normal to the electric field
// scattering outputs saved (if chosen) and output for all energies
wave/d Scatter3DSave // Full 2D Scattering Scattering Pattern Saved vs Energy
//wave/d/c PZ3DSave // Full 3D PZ field Saved vs Energy
wave /d int3DvsEn // same but for the 3D calculated data
wave /d ratio3DvsEn // same but for the 3D calculated data
wave /d Para3DvsEn // same but for the 3D calculated data
wave /d Perp3DvsEn // same but for the 3D calculated data
// internal scattering variables
variable step // current energy step (integer for storing data)
variable en // current energy = 1239.84197 / wavelength = ( 1239.84197/(2*pi)) * wavevector
variable wavelength // current wavelength = 1239.84197 / energy = 2*pi / wavevector (in nanometers)
variable k // current wavevector magnitude (2*pi / wavelength) = (2 * pi / 1239.84197 ) * Energy (in inverse nanometers)
// current code has the X-rays propotaging in the x direction
endstructure
function addtologbook(s3d,lognote)
struct ThreeDSystem &s3d
string lognote
s3d.timesofar +=stopmstimer(s3d.timer)/10^6
s3d.timer = startmstimer
variable loc = dimsize(s3d.logbook,0)
insertpoints /M=0 loc,1, s3d.logbook
s3d.logbook[loc][0] = time2str(s3d.timesofar)
s3d.logbook[loc][1] = lognote
make /t /free timew = {time2str(s3d.timesofar)}
make /t /free notew = {lognote}
ListBox progresslistbox win=SIM_STATUS,row=max(0,loc-8)
updateprogress(s3d)
pathinfo save3dscattering
if(v_flag)
Save/A=2/J/M="\n"/O/F/p=Save3DScattering timew,notew as s3d.name+".log"
endif
end
function updateprogress(s3d)
struct ThreeDSystem &s3d
nvar pbar = $(s3d.path+"progressbar")
pbar = s3d.progressbar
doupdate
end
function updatephase(s3d,phase)
struct ThreeDSystem &s3d
string phase
SetDrawLayer/w=SIM_Status /k UserBack
SetDrawEnv/w=SIM_Status fsize= 16,textxjust= 1,textyjust= 1
DrawText/w=SIM_Status 425,15,"Simulation: "+s3d.name+" using "+ s3d.modelname + " morphology. Progress displayed : " + phase
end
function model3D(modelname,voxelsize,sizescale,resolution,thickness,paramstring,materiallist,materialalignmentstring,efield,energymin,energymax,energynum[movie, save3d,moviepath,enwave,rotatesys,skipsim,outputcy])
// electric field is always in the y z plane
// modelname - a string with the name of the model accepted models are specificed below
// morphology parameters
// sizsescale - basic parameter of all models, incidating general size scale
// resolution - basic parameter of all models, indicating the grid size of the system (the pixel is the basic size)
// paramstring - a model specific string containing all parameters needed to create the system
//alignment parameters
// materiallist - a list of all the materials in the system. the number of materials needs to match the model requirements seperated by ","
// materialalignment - a list the length of material list (material1alignment;material2alignment ...)
//each materialalignment item is a comma seperated list of alignment parameters
//( alignment (0 faceon, 1 edgeon, other = none) if none, then whatever is produced by morphology is passed through (this may include alignment) ,
// volfrac (negative means as calculated by morphology) ,
// ratio of entropy cost of alignment vs interface alignment)
//scattering parameters
// efield - a string containing a vector (y component, z component)
// energymin - the minimum x-ray energy to calculate
// energymax - the maximum x-ray energy to calculate
// energynum - the number of energysteps
//simulation options (optional)
//movie = make a movie of the process?
// save3d - save the full scattering pattern at each energy
// moviepath - avoid a popup if you supply a valid path for the movie
// enwave - overwrites energy min, max and num and just calculated the energy values in this wave, if it exists
// boolean to rotate the system - greatly increases the quality o he calculation, can be done in he fuure by roating the efield, rather than the elaborate way it is done now
string modelname,paramstring,materiallist,materialalignmentstring,efield, moviepath
variable sizescale, resolution, energymin,energymax,energynum, movie, save3d, thickness,voxelsize, rotatesys, skipsim, outputcy
wave/z enwave // this will override the energy start, step etc, with arbitrary energy steps
struct ThreeDSystem s3D
variable result, timer1
for(timer1=0;timer1<11;timer1+=1)
s3d.timesofar = stopmstimer(timer1)
endfor
s3d.voxelsize = voxelsize
s3d.timesofar = 0
s3d.timer = startmstimer
s3d.thickness = thickness
s3d.num = resolution
s3d.size = sizescale
s3d.materials = materiallist
s3d.efield = efield
s3d.rot = paramisdefault(rotatesys)? 0 : rotatesys
s3d.modelname = modelname
s3d.paramstring = paramstring
s3d.movie = movie
s3d.materialalignmentstring = materialalignmentstring
s3d.materialnum = itemsinlist(materiallist,",")
s3d.progressbar = 0
s3d.name = getdatafolder(0)
string dfolder = getdatafolder(1)
make /n=(0,2)/o /t Progresslist
wave /t s3d.logbook = Progresslist
variable /g progressbar
s3d.path = getdatafolder(1)
killwindow/z SIM_Status
NewPanel /k=1/n=SIM_Status/W=(800,767,1650,1014) as getdatafolder(0) + " Progress"
updatephase(s3d,"Starting Up")
ValDisplay valdispProgress,win=SIM_Status,limits={0,100,0},barmisc={10,1},highColor= (16385,16388,65535),pos={37.00,39.00},size={747.00,41.00},bodyWidth=747
ValDisplay valdispProgress,win=SIM_Status,frame=4,value=#(dfolder+"progressbar")
ListBox progresslistbox,win=SIM_STATUS,pos={18.00,98.00},size={796.00,134.00},listWave=$(dfolder+"Progresslist"),widths={57,688}
if(paramisdefault(moviepath) || strlen(moviepath)<3)
open /D/F=".mov" fileref as getdatafolder(0)+".mov"
//close fileref
moviepath = S_filename
endif
newpath /o/z/q Save3DScattering, parseFilePath(1,moviepath,":",1,0)
if(!v_flag)
Save/J/M="\n"/O/F/p=Save3DScattering s3d.logbook as s3d.name+".log"
endif
addtologbook(s3d,"Starting Up - Using " + s3d.modelname + " Morphology")
addtologbook(s3d,"parameters: " + s3d.paramstring)
addtologbook(s3d,"Voxelsize: " + num2str(s3d.voxelsize))
addtologbook(s3d,"thickness: " + num2str(thickness))
addtologbook(s3d,"resolution: " + num2str(resolution))
addtologbook(s3d,"sizescale: " + num2str(sizescale))
addtologbook(s3d,"materials: " + materiallist)
addtologbook(s3d,"efield: along Z")
addtologbook(s3d,"modelname: " + modelname)
addtologbook(s3d,"movie: " + num2str(movie))
addtologbook(s3d,"Output for CyRSoXS input: " + num2str(outputcy))
addtologbook(s3d,"Skip the Simulation: " + num2str(skipsim))
addtologbook(s3d,"Material Slignment String: " + materialalignmentstring)
string cmdstr = "Model3D("
cmdstr += "\"" + modelname + "\","
cmdstr += num2str(voxelsize) + ","
cmdstr += num2str(sizescale) + ","
cmdstr += num2str(resolution) + ","
cmdstr += num2str(thickness) + ","
cmdstr += "\"" + paramstring + "\","
cmdstr += "\"" + materiallist + "\","
cmdstr += "\"" + materialalignmentstring + "\","
cmdstr += "\"" + efield + "\","
cmdstr += num2str(energymin) + ","
cmdstr += num2str(energymax) + ","
cmdstr += num2str(energynum) + ","
if(!paramisdefault(movie))
cmdstr += "movie=" + num2str(movie) + ","
endif
if(!paramisdefault(save3d))
cmdstr += "save3d=" + num2str(save3d) + ","
endif
if(!paramisdefault(moviepath))
cmdstr += "\r\tmoviepath=\"" + moviepath + "\",\r\t "
endif
if(!paramisdefault(enwave))
cmdstr += "enwave=" + nameofwave(enwave) + ","
endif
if(!paramisdefault(rotatesys))
cmdstr += "rotatesys=" + num2str(rotatesys) + ","
endif
if(!paramisdefault(outputcy))
cmdstr += "outputcy=" + num2str(outputcy) + ","
endif
if(!paramisdefault(skipsim))
cmdstr += "skipsim=" + num2str(skipsim) + ","
endif
cmdstr = removeending(cmdstr,",")
cmdstr += ")"
addtologbook(s3d,"Command to run again: " + replacestring("\r\t",cmdstr,""))
print/len=1000 cmdstr
if(paramisdefault(enwave))
make /n=(energynum)/o enwave=energymin + p * (energymax-energymin)/(energynum-1)
else
if(waveexists(enwave))
energynum = numpnts(enwave)
energymin = wavemin(enwave)
energymax = wavemax(enwave)
else
make /n=(energynum)/o enwave=energymin + p * (energymax-energymin)/(energynum-1)
endif
endif
duplicate/o enwave, enwavedisp
insertpoints 0,1,enwavedisp
enwavedisp[0]=enwave[0]-.1
// set up the display and movie output
Variable DebugEnab
Variable err
dowindow /k simulation_layout
execute("Simulation_Layout()")
if(s3d.movie)
DebuggerOptions
DebugEnab = V_debugOnError //check for debug on error
if (DebugEnab) //if it is on,
DebuggerOptions debugOnError=0 //turn it off
endif
try
closemovie;AbortOnRTE
catch
err = GetRTError(1)
addtologbook(s3d," - Creating new movie")
endtry
if(DebugEnab)
DebuggerOptions debugOnError=1
endif
SavePict/O/WIN=Simulation_Layout /E=-5/P=_PictGallery_/w=(0,0,800,800) as "SimPict"
variable fileref
//print moviepath
newmovie /z /PICT=SimPict /O as moviepath
if(v_flag)
s3d.movie=0
addtologbook(s3d,"movie could not be started")
else
addtologbook(s3d,"movie started")
endif
endif
// try to call the chosen morphology creator
if(exists("model3d_"+modelname)==6)
updatephase(s3d,"Building Morphology")
funcref model3d_existing creatfunc=$("model3d_"+modelname)
creatfunc(s3d)
else
print "no recognized model could be created"
return -1
endif
variable alignmentincluded
if(exists("special_"+modelname)==6) // this is a unique function which if exists, will give some extra information about the model, in this case wether alignment is created
// in the model, or it needs to be calculated - with modern models it is almos always generated within the model creation process
funcref special_existing specfunc=$("special_"+modelname)
alignmentincluded = stringmatch(specfunc(),"*IncludesAlignment*")
endif
//calculate alignment of system if needed - some morphologies already do this (code above checkes for this)
if( stringmatch(materialalignmentstring,"none") )
//Print "Time : "+time2str(s3d.timesofar) +" - Loading Existing material Alignment"
addtologbook(s3d,"Loading Existing Material Alignment")
updatephase(s3d,"Loading Alignment")
loadexistingmaterialalignment(s3d)
elseif( !alignmentincluded )
//Print "Time : "+time2str(s3d.timesofar) +" - Calculating alignment in system"
updatephase(s3d,"Calculating Alignment")
addtologbook(s3d,"Calculating alignment in system")
Align3Dsystem(s3d)
endif
addtologbook(s3d,"Alignment Loaded - Checking System Density for consistancy")
// check the m1-m5 matrices to make sure they make sense
updatephase(s3d,"Checking Morphology and Alignment")
if(sum3dsystem(s3d)<0)
addtologbook(s3d,"warning : System Creation failed to conserve density throughout system")
//Print "Time : "+time2str(s3d.timesofar) +" - warning : System Creation failed to conserve density throughout system"
elseif(sum3dsystem(s3d)>0)
addtologbook(s3d,"Error : System Creation failed no density matrices are found")
//print "Time : "+time2str(s3d.timesofar) +" - Error : System Creation failed no density matrices are found"
return -1
else
addtologbook(s3d,"Check of system density checks out")
//Print "Time : "+time2str(s3d.timesofar) +" - Check of system density checks out"
endif
variable enstep = (energymax-energymin)/(energynum-1)
if(outputcy)
// At this point we have all the values that we need. we could save the model now to disk, or load a model from disk at this point
// data to save to model file
// m1-m5 (if they exist)
// first three dimensions are the aligned component - should go into mat_(1-5)_alignment vectors field
// fourth dimension is the unaligned component should go into mat_(1-5)_unaligned scalar field
writeconfig(s3d,replacestring(".mov",replacestring(".mp4",moviepath,":"),":"),energymin, energymax,enstep)
variable h5file, groupid
HDF5CreateFile /p=cyrsoxspath/O h5file as parsefilepath(0,replacestring("mov",replacestring("mp4",moviepath,"hd5"),"hd5"),":",1,0)
/// this next chunk is the depricated vector morphology representation. This is now replaced with angular representation
HDF5CreateGroup h5file , "vector_morphology" , groupID
newdatafolder /o/s hd5output
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2),3) /d Mat_1_alignment = s3d.m1[p][q][r][2-t] // switching dimension order for CyRSoXS beta 2
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_1_alignment, groupiD
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d Mat_1_unaligned = s3d.m1[p][q][r][3]
hdf5saveData /GZIP = {9,1} /LAYO={2,64,64,64} /MAXD={256,2048,2048} Mat_1_unaligned, groupiD
if(waveexists(s3d.m2))
make /o /n=(dimsize(s3d.m2,0),dimsize(s3d.m2,1),dimsize(s3d.m2,2),3) /d Mat_2_alignment = s3d.m2[p][q][r][2-t]
make /o /n=(dimsize(s3d.m2,0),dimsize(s3d.m2,1),dimsize(s3d.m2,2)) /d Mat_2_unaligned = s3d.m2[p][q][r][3]
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_2_alignment, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,64,64,64} /MAXD={256,2048,2048} Mat_2_unaligned, groupiD
if(waveexists(s3d.m3))
make /o /n=(dimsize(s3d.m3,0),dimsize(s3d.m3,1),dimsize(s3d.m3,2),3) /d Mat_3_alignment = s3d.m3[p][q][r][2-t]
make /o /n=(dimsize(s3d.m3,0),dimsize(s3d.m3,1),dimsize(s3d.m3,2)) /d Mat_3_unaligned = s3d.m3[p][q][r][3]
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_3_alignment, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,64,64,64} /MAXD={256,2048,2048} Mat_3_unaligned, groupiD
if(waveexists(s3d.m4))
make /o /n=(dimsize(s3d.m4,0),dimsize(s3d.m4,1),dimsize(s3d.m4,2),3) /d Mat_4_alignment = s3d.m4[p][q][r][2-t]
make /o /n=(dimsize(s3d.m4,0),dimsize(s3d.m4,1),dimsize(s3d.m4,2)) /d Mat_4_unaligned = s3d.m4[p][q][r][3]
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_4_alignment, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,64,64,64} /MAXD={256,2048,2048} Mat_4_unaligned, groupiD
if(waveexists(s3d.m5))
make /o /n=(dimsize(s3d.m5,0),dimsize(s3d.m5,1),dimsize(s3d.m5,2),3) /d Mat_5_alignment = s3d.m5[p][q][r][2-t]
make /o /n=(dimsize(s3d.m5,0),dimsize(s3d.m5,1),dimsize(s3d.m5,2)) /d Mat_5_unaligned = s3d.m5[p][q][r][3]
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_5_alignment, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64} /MAXD={256,2048,2048} Mat_5_unaligned, groupiD
endif
endif
endif
endif
setdatafolder ::
// this is the new unified morphology representation (for now, both are produced for backwards compatability)
HDF5CreateGroup h5file , "morphology" , groupID
newdatafolder /o/s hd5output
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d Mat_1_S = sqrt(s3d.m1[p][q][r][0]^2 + s3d.m1[p][q][r][1]^2 + s3d.m1[p][q][r][2]^2)
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d Mat_1_Phi = atan(s3d.m1[p][q][r][0]/s3d.m1[p][q][r][1])
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d Mat_1_Theta = acos(s3d.m1[p][q][r][2]/Mat_1_S[p][q][r])
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d Mat_1_vfrac = s3d.m1[p][q][r][3] + Mat_1_S[p][q][r]
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_1_S, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_1_Phi, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_1_Theta, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_1_vfrac, groupiD
if(waveexists(s3d.m2))
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d Mat_2_S = sqrt(s3d.m2[p][q][r][0]^2 + s3d.m2[p][q][r][1]^2 + s3d.m2[p][q][r][2]^2)
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d Mat_2_Phi = atan(s3d.m2[p][q][r][0]/s3d.m2[p][q][r][1])
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d Mat_2_Theta = acos(s3d.m2[p][q][r][2]/Mat_2_S[p][q][r])
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d Mat_2_vfrac = s3d.m2[p][q][r][3] + Mat_2_S[p][q][r]
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_2_S, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_2_Phi, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_2_Theta, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} Mat_2_vfrac, groupiD
if(waveexists(s3d.m3))
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d mat_3_S = sqrt(s3d.m3[p][q][r][0]^2 + s3d.m3[p][q][r][1]^2 + s3d.m3[p][q][r][2]^2)
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d mat_3_Phi = atan(s3d.m3[p][q][r][0]/s3d.m3[p][q][r][1])
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d mat_3_Theta = acos(s3d.m3[p][q][r][2]/mat_3_S[p][q][r])
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d mat_3_vfrac = s3d.m3[p][q][r][3] + mat_3_S[p][q][r]
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} mat_3_S, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} mat_3_Phi, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} mat_3_Theta, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} mat_3_vfrac, groupiD
if(waveexists(s3d.m4))
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d mat_4_S = sqrt(s3d.m4[p][q][r][0]^2 + s3d.m4[p][q][r][1]^2 + s3d.m4[p][q][r][2]^2)
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d mat_4_Phi = atan(s3d.m4[p][q][r][0]/s3d.m4[p][q][r][1])
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d mat_4_Theta = acos(s3d.m4[p][q][r][2]/mat_4_S[p][q][r])
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d mat_4_vfrac = s3d.m4[p][q][r][3] + mat_4_S[p][q][r]
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} mat_4_S, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} mat_4_Phi, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} mat_4_Theta, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} mat_4_vfrac, groupiD
if(waveexists(s3d.m5))
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d mat_5_S = sqrt(s3d.m5[p][q][r][0]^2 + s3d.m2[p][q][r][1]^2 + s3d.m2[p][q][r][2]^2)
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d mat_5_Phi = atan(s3d.m2[p][q][r][0]/s3d.m2[p][q][r][1])
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d mat_5_Theta = acos(s3d.m2[p][q][r][2]/mat_5_S[p][q][r])
make /o /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1),dimsize(s3d.m1,2)) /d mat_5_vfrac = s3d.m2[p][q][r][3] + mat_5_S[p][q][r]
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} mat_5_S, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} mat_5_Phi, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} mat_5_Theta, groupiD
hdf5saveData /GZIP = {9,1} /LAYO={2,1,64,64,1} /MAXD={256,2048,2048,3} mat_5_vfrac, groupiD
endif
endif
endif
endif
newdatafolder /o/s igor_paramaters
// parameters to save to parameter file (*unique to Igor)
variable /g igorrotation = s3d.rot // boolean saying if we are rotating the system or not (90 degrees is always included) * - this should be universal
variable /g igornum = s3d.num // resolution of simulation (in each of the long axes of the film) * - the data itself should have this
variable /g igorthickness = s3d.thickness // voxels of thickness along the beam direction * - the data itself should have this
variable /g igorvoxelsize = s3d.voxelsize // the size of each pixel in the simulation (in nanometers) * - can be changed trivially, it is just a scaling parameter
variable /g igormaterialnum = s3d.materialnum // number of materials * - can be read from materials list
string /g igorparamstring = s3d.paramstring //specific parameters for creating density matricies according to the model * - used for morphology creation
string /g igormodelname = s3d.modelname //store the model name somewhere * - used for morphology creation
string /g igorpath = s3d.path // the path to put important things for this calculation
string /g igorname = s3d.name // name of simulation
variable /g igormovie = s3d.movie // are we doing a movie * - for visualization and some feedback
string /g igormaterials = s3d.materials // list of materials
string /g igorefield = s3d.efield // direction of efield - in the future, if combined with rotation, this doesn't need to be defined, maybe - all orientations will be calculated
setdatafolder ::
newdatafolder /o/s morphology_variables
make /o /n=3 film_normal = {1,0,0}
variable /g voxel_size_nm = s3d.voxelsize
string /g morphology_creator = "EG_"+s3d.modelname
string /g version = secs2date(datetime,-2)
string /g creation_date = secs2date(datetime,-2) +" " + secs2time(datetime,1)
string /g name = s3d.name
setdatafolder ::
HDF5CreateGroup h5file , "igor_parameters" , groupID
HDF5SaveGroup /L=7 /O /R igor_paramaters , h5file , "igor_parameters"
HDF5CreateGroup h5file , "morphology_variables" , groupID
HDF5SaveGroup /L=7 /O /R morphology_variables , h5file , "morphology_variables"
HDF5CloseFile h5file
setdatafolder ::
killdatafolder /z hd5output
endif
if(skipsim) // we will not be doing the simulation, we are done
addtologbook(s3d,"No Simulation will be done")
if(s3d.movie)
addtologbook(s3d,"Finishing Movie, and closing file")
//Print "Time : "+time2str(s3d.timesofar) +" - Finishing Movie, and closing file"
DebuggerOptions
DebugEnab = V_debugOnError //check for debug on error
if (DebugEnab) //if it is on,
DebuggerOptions debugOnError=0 //turn it off
endif
try
closemovie;AbortOnRTE
catch
err = GetRTError(1)
addtologbook(s3d,"No Frames added, no movie created!")
endtry
if(DebugEnab)
DebuggerOptions debugOnError=1
endif
endif
updatephase(s3d,"Complete")
s3d.progressbar=100
addtologbook(s3d,"Complete.")
return 1
endif
// make waves to store integrations and ratios
make/d/o/n=(floor(s3d.num/sqrt(2)),energynum) int3DvsEn=0,ratio3DvsEn=0, para3dvsen, perp3dvsen
make/d/o/n=(s3d.num,s3d.num,energynum) scatter3DSave
//make/c/d/o/n=(s3d.thickness,s3d.num,s3d.num,energynum) PZ3DSave
wave s3d.scatter3DSave=scatter3DSave, s3d.int3DvsEn=int3DvsEn, s3d.ratio3DvsEn=ratio3DvsEn, s3d.para3dvsen=para3dvsen, s3d.perp3dvsen=perp3dvsen
//wave /c s3d.pz3dsave = pz3dsave
setscale /p x,pi/s3d.num,pi/s3d.num, s3d.int3DvsEn,s3d.ratio3DvsEn, s3d.para3dvsen, s3d.perp3dvsen // in physics units of q
setscale /i x, -pi/s3d.voxelsize,pi/s3d.voxelsize, s3d.scatter3DSave
setscale /i y, -pi/s3d.voxelsize,pi/s3d.voxelsize, s3d.scatter3DSave
setscale /i y, energymin, energymax, s3d.int3DvsEn,s3d.ratio3DvsEn, s3d.para3dvsen, s3d.perp3dvsen // only valid if not using a energy wave // best to not used it, and loot at the wave
setscale /i z, energymin, energymax, s3d.scatter3DSave // only valid if not using a wave, best to just use the wave
//setscale /p x,dimoffset(s3d.density1,0), dimdelta(s3d.density1,0), PZ3DSave
//setscale /p y,dimoffset(s3d.density1,1), dimdelta(s3d.density1,1), PZ3DSave
//setscale /p z,dimoffset(s3d.density1,2), dimdelta(s3d.density1,2), PZ3DSave
variable en, j
variable angle=0, anglestep=2 //(1 degree steps from 0 to 180) - 180 steps
if(!s3d.rot)
anglestep = 90
endif
variable offsetx, offsety // offset for rotation to get rid of speckle
s3d.step=0
updatephase(s3d,"Calculating Scattering")
// expand the M matrices (this will allow rotation, although it assumes periodic type boundary conditions, otherwise this will add scattering artifacts.)
make /n=(dimsize(s3d.m1,0),dimsize(s3d.m1,1)*1.8,dimsize(s3d.m1,2)*1.8,dimsize(s3d.m1,3))/o m1ext = s3d.m1[p][mod(q+.9*dimsize(s3d.m1,1),dimsize(s3d.m1,1))][mod(r+.9*dimsize(s3d.m1,2),dimsize(s3d.m1,2))][s]
duplicate/o s3d.m1, m1save
if(s3d.materialnum>1)
make /n=(dimsize(s3d.m2,0),dimsize(s3d.m2,1)*1.8,dimsize(s3d.m2,2)*1.8,dimsize(s3d.m2,3))/o m2ext = s3d.m2[p][mod(q+.9*dimsize(s3d.m2,1),dimsize(s3d.m2,1))][mod(r+.9*dimsize(s3d.m2,2),dimsize(s3d.m2,2))][s]
duplicate/o s3d.m2, m2save
endif
if(s3d.materialnum>2)
make /n=(dimsize(s3d.m3,0),dimsize(s3d.m3,1)*1.8,dimsize(s3d.m3,2)*1.8,dimsize(s3d.m3,3))/o m3ext = s3d.m3[p][mod(q+.9*dimsize(s3d.m3,1),dimsize(s3d.m3,1))][mod(r+.9*dimsize(s3d.m3,2),dimsize(s3d.m3,2))][s]
duplicate/o s3d.m3, m3save
endif
if(s3d.materialnum>3)
make /n=(dimsize(s3d.m4,0),dimsize(s3d.m4,1)*1.8,dimsize(s3d.m4,2)*1.8,dimsize(s3d.m4,3))/o m4ext = s3d.m4[p][mod(q+.9*dimsize(s3d.m4,1),dimsize(s3d.m4,1))][mod(r+.9*dimsize(s3d.m4,2),dimsize(s3d.m4,2))][s]
duplicate/o s3d.m4, m4save
endif
if(s3d.materialnum>4)
make /n=(dimsize(s3d.m5,0),dimsize(s3d.m5,1)*1.8,dimsize(s3d.m5,2)*1.8,dimsize(s3d.m5,3))/o m5ext = s3d.m5[p][mod(q+.9*dimsize(s3d.m5,1),dimsize(s3d.m5,1))][mod(r+.9*dimsize(s3d.m5,2),dimsize(s3d.m5,2))][s]
duplicate/o s3d.m5, m5save
endif
createsimlayout() // open the live display for the scattering simulation
for( j=0; j <dimsize(enwave,0) ; j+=1) // main loop through energies
en=enwave[j]
//s3d.progressbar = 100*(en-energymin) / (energymax-energymin)
s3d.progressbar = 100*(j) / (dimsize(enwave,0))
s3d.en = en
s3d.step = j
s3d.wavelength = 1239.84197/en
s3d.k = 2*pi/s3d.wavelength
addtologbook(s3d,"calculating polarization for " + num2str(s3d.en) + " eV")
//Print "Time : "+time2str(s3d.timesofar) +" - Creating Indexwave for energy : " + num2str(s3d.en)
// loop through polarization and scattering calculations as the system is rotated
// this assumes that the enlarged systems we have created were from periodic boundary conditions. Unfortunately once rotated, the boundaries cannot maintain periodicity, and will introduce some artifacts\
// this is far outweighed by the benefit particulary when looking at anisotropy
// For Future development**** instead of rotating the system, rotate the efield, and then rotate the scattering pattern back. it is very important that sums are only made with efield pointing in the same direction
// otherwise all anisotropy will be removed (simulating circular polarization)
for(angle=0;angle<180;angle+=anglestep)
offsetx = 0//floor(enoise(dimsize(s3d.m1,1)/10)) // this was rotating and offsetting the center randomly, I found that the offsetting didn't matter
offsety = 0//floor(enoise(dimsize(s3d.m1,2)/10))
sysrotate(m1ext,s3d.m1,angle,offsetx,offsety)
if(s3d.materialnum>1)
sysrotate(m2ext,s3d.m2,angle,offsetx,offsety)
endif
if(s3d.materialnum>2)
sysrotate(m3ext,s3d.m3,angle,offsetx,offsety)
endif
if(s3d.materialnum>3)
sysrotate(m4ext,s3d.m4,angle,offsetx,offsety)
endif
if(s3d.materialnum>4)
sysrotate(m5ext,s3d.m5,angle,offsetx,offsety)
endif
result = CalculatePolarizationWave(s3d) // This is the main calculation, taking efield and m matrices with npara and nperp and producing the induced polarization within each voxel
if(result<0)
addtologbook(s3d,"quitting because of critical errors in calculating polarization")
//print "quitting because of critical errors in calculating polarization"
return -1
elseif(result>0 && s3d.en == energymin && angle==0)
addtologbook(s3d,"warning - no aligned optical constants for "+stringfromlist(result-1, s3d.materials,",")+" were found, using normal Optical Constants")
//print "warning - no aligned optical constants for "+stringfromlist(result-1, s3d.materials,",")+" were found, using normal Optical Constants"
endif
// addtologbook(s3d,"Angle="+num2str(angle))
//Print "Time : "+time2str(s3d.timesofar) +" - Creating 3D Scattering Pattern"
Calculate3DScatter(s3d) // this projects the induced polarization to the far field, simulating the scattering which emits from this polarization field
if(angle==0)
duplicate/free s3d.scatter3d, scatter3dsum
else
scatter3dsum+= s3d.scatter3D
endif
if(angle/10 == round(angle/10) && angle + anglestep < 180)
if(angle>0)
s3d.scatter3D = scatter3Dsum/(angle/anglestep) // add up the scattering from different rotations
else
RadialIntegratesystem(s3d)
killwindow /z simulation_layout
createsimlayout()
TextBox/w=Simulation_Layout/C/N=text5/F=0/A=RT/X=51.85/Y=37.15/G=(16386,65535,16385) "\\Z24"+num2str(s3d.en)+" eV"
SetDrawLayer /w=RatioIntensity /k userfront
SetDrawEnv /w=RatioIntensity ycoord= left,linefgc= (65280,0,0),dash= 1,linethick= 3.00;DelayUpdate
DrawLine/w=RatioIntensity 0,s3d.en,1,s3d.en
SetDrawLayer /w=ScatteringIntensity /k userfront
SetDrawEnv /w=ScatteringIntensity ycoord= left,linefgc= (65280,0,0),dash= 1,linethick= 3.00;DelayUpdate
DrawLine/w=ScatteringIntensity 0,s3d.en,1,s3d.en
endif
RadialIntegratesystem(s3d)
ModifyGraph/w=ParaPerpInt hbFill(int3Dpara)=5
doupdate
dowindow /f simulation_layout
if(s3d.movie)
SavePict/O/WIN=Simulation_Layout /E=-5/P=_PictGallery_/w=(0,0,800,800) as "SimPict"
AddMovieFrame/PICT=SimPict
endif
endif
endfor
s3d.scatter3D = scatter3Dsum/(180/anglestep)
//imagerotate /o/H scatter3Dsum
//s3d.scatter3D += scatter3Dsum/(180/anglestep)
//imagerotate /o/V scatter3Dsum
//s3d.scatter3D += scatter3Dsum/(180/anglestep)
//imagerotate /o/H scatter3Dsum
//s3d.scatter3D += scatter3Dsum/(180/anglestep)
RadialIntegratesystem(s3d)
StoreIntegrations(s3d)
// visualizescatter(s3d) // updates displays of scattering profiles
//(re)create layout of scattering for an energy
//CreateSimLayout()
// update the layout (projections, 2D scattering pattern, label, and cuts)
//doupdate
//doupdate
if(s3d.movie)
killwindow /z simulation_layout
// make the display
createsimlayout()
TextBox/w=Simulation_Layout/C/N=text5/F=0/A=RT/X=51.85/Y=37.15/G=(16386,65535,16385) "\\Z24"+num2str(s3d.en)+" eV"
SetDrawLayer /w=RatioIntensity /k userfront
SetDrawEnv /w=RatioIntensity ycoord= left,linefgc= (65280,0,0),dash= 1,linethick= 3.00;DelayUpdate
DrawLine/w=RatioIntensity 0,s3d.en,1,s3d.en
SetDrawLayer /w=ScatteringIntensity /k userfront
SetDrawEnv /w=ScatteringIntensity ycoord= left,linefgc= (65280,0,0),dash= 1,linethick= 3.00;DelayUpdate
DrawLine/w=ScatteringIntensity 0,s3d.en,1,s3d.en
SavePict/O/WIN=Simulation_Layout /E=-5/P=_PictGallery_/w=(0,0,800,800) as "SimPict"
dowindow /f simulationlayout
ModifyGraph/w=ParaPerpInt hbFill(int3Dpara)=2
doupdate
AddMovieFrame/PICT=SimPict
AddMovieFrame/PICT=SimPict
AddMovieFrame/PICT=SimPict
AddMovieFrame/PICT=SimPict
AddMovieFrame/PICT=SimPict
AddMovieFrame/PICT=SimPict
AddMovieFrame/PICT=SimPict
AddMovieFrame/PICT=SimPict
endif
s3d.step +=1
endfor
if(s3d.movie)
addtologbook(s3d,"Finishing Movie, and closing file")
//Print "Time : "+time2str(s3d.timesofar) +" - Finishing Movie, and closing file"
DebuggerOptions
DebugEnab = V_debugOnError //check for debug on error
if (DebugEnab) //if it is on,
DebuggerOptions debugOnError=0 //turn it off
endif
try
closemovie;AbortOnRTE
catch
err = GetRTError(1)
addtologbook(s3d,"No Frames added, no movie created!")
endtry
if(DebugEnab)
DebuggerOptions debugOnError=1
endif
endif
updatephase(s3d,"Complete")
s3d.progressbar=100
addtologbook(s3d,"Complete.")
//killwindow /z SIM_Status
//Print "Time : "+time2str(s3d.timesofar) +" - Complete."
//display 2D rato and integrations vs q and energy
end
function createsimlayout()
//dowindow /k Simulation_Layout
debuggeroptions debugonerror=0
ScatteringIntensitydisp()
RatioIntensitydisp()
ScatteringPatterndisp()
ParaPerpIntdisp()
DeltaProjectiondisp()
BetaProjectiondisp()
dowindow Simulation_Layout
if(!v_flag)
execute("Simulation_Layout()")
endif
dowindow /f Simulation_Layout
doupdate
end
function align3dsystem(s3d) // given a morphology of only scalar density, produce an alignment field
struct ThreeDSystem &s3d
string matstr = s3d.materialalignmentstring
string mat,alignment,alignmentwidth, alignmentstrength
variable nummaterials = itemsinlist(matstr) , i
variable align, width, strength
variable volfrac
if(nummaterials==2)
mat = stringfromlist(0,matstr)
if(stringmatch(mat,"*Face-on*") || stringmatch(mat, "1*"))
align = 1
elseif(stringmatch(mat, "*Edge-on*") || stringmatch(mat, "0*"))
align=0
else
//don't do the next calculation
align=-1
endif
if(align==-1)
//don't do the calculation for m1
make /o/n=(dimsize(s3d.density1,0),dimsize(s3d.density1,1),dimsize(s3d.density1,2),4) s3d.m1=0
s3d.m1[][][][3] = s3d.density1[p][q][r]
else
splitstring/e="^[^,]*,([^,]*),([^,]*)" mat, alignmentwidth, alignmentstrength
strength = str2num(alignmentstrength)
width = str2num(alignmentwidth)
s3d.timesofar += stopmstimer(s3d.timer)/1e6
s3d.timer = startmstimer
wave alignmentw = createalignmentdensity(s3d.density1,width,strength, align, s3d)
//createbinarysystem(s3d.density1,width,strength, 1, 0, 1- volfrac, 0, 1,s=s)
duplicate/o alignmentw, s3d.m1
//s3d.m1 *=s3d.density1[p][q][r]
s3d.m1 *=sqrt(s3d.density1[p][q][r])
s3d.m1[][][][3] *=sqrt(s3d.density1[p][q][r])
endif
mat = stringfromlist(1,matstr)
if(stringmatch( mat,"*Face-on*") || stringmatch(mat, "1*") )
align = 1
elseif(stringmatch( mat,"*Edge-on*") || stringmatch(mat, "0*") )
align=0
else
duplicate/o s3d.m1, s3d.m2
s3d.m2[][][][0,2]=0
s3d.m2[][][][3] = 1 - s3d.m1[p][q][r][0]^2 - s3d.m1[p][q][r][1]^2 - s3d.m1[p][q][r][2]^2 -s3d.m1[p][q][r][3] // all density that isn't accounted for in m1 is unaligned m2
align=-1
endif
if(align==-1)
// we've already created both m1 and m2, so we're done!
else
splitstring/e="^[^,]*,([^,]*),([^,]*)" mat, alignmentwidth, alignmentstrength
strength = str2num(alignmentstrength)
width = str2num(alignmentwidth)
volfrac = mean(s3d.density1)
s3d.timesofar += stopmstimer(s3d.timer)/1e6
s3d.timer = startmstimer
wave alignmentw = createalignmentdensity(s3d.density1,width,strength, align, s3d)
// wave alignmentw = root:packages:ScatterSim3D:test:sn
//createbinarysystem(s3d.density1,width,strength, 1, 0, 1- volfrac, 0, 1,s=s)
duplicate/o alignmentw, s3d.m2
s3d.m2 *=sqrt(1-s3d.density1[p][q][r])
s3d.m2[][][][3] *=sqrt(1-s3d.density1[p][q][r])
endif
else
for(i=0;i<nummaterials;i+=1)
mat = stringfromlist(i,matstr)
// make density and m matricies for all dimensions
if(i==0)
make /o/n=(dimsize(s3d.density1,0),dimsize(s3d.density1,1),dimsize(s3d.density1,2),4) s3d.m1=0
s3d.m1[][][][3] = s3d.density1[p][q][r]
elseif(i==1)
duplicate /o s3d.m1,s3d.m2
s3d.m2 = 0
if(i< nummaterials - 1)
s3d.m2[][][][3] = s3d.density2[p][q][r]
else
//this is the last material, so density wave doesn't exist, need to calculate it
s3d.m2[][][][3] = 1 - s3d.m1[p][q][r][0]^2 - s3d.m1[p][q][r][1]^2 - s3d.m1[p][q][r][2]^2 -s3d.m1[p][q][r][3]
duplicate /o s3d.density1,s3d.density2
s3d.density2[][][] = s3d.m2[p][q][r][3]
endif
elseif(i==2)
duplicate /o s3d.m1,s3d.m3
s3d.m3 = 0
if(i< nummaterials - 1)
s3d.m3[][][][3] = s3d.density3[p][q][r]
else
//this is the last material, so density wave doesn't exist, need to calculate it
s3d.m3[][][][3] = 1
s3d.m3[][][][3] -= s3d.m1[p][q][r][0]^2 + s3d.m1[p][q][r][1]^2 + s3d.m1[p][q][r][2]^2 + s3d.m1[p][q][r][3]
s3d.m3[][][][3] -= s3d.m2[p][q][r][0]^2 + s3d.m2[p][q][r][1]^2 + s3d.m2[p][q][r][2]^2 + s3d.m2[p][q][r][3]
duplicate /o s3d.density1,s3d.density3
s3d.density3[][][] = s3d.m3[p][q][r][3]
endif
elseif(i==3)
duplicate /o s3d.m1,s3d.m4
s3d.m4 = 0
if(i< nummaterials - 1)
s3d.m4[][][][3] = s3d.density4[p][q][r]
else
//this is the last material, so density wave doesn't exist, need to calculate it
s3d.m4[][][][3] = 1
s3d.m4[][][][3] -= s3d.m1[p][q][r][0]^2 + s3d.m1[p][q][r][1]^2 + s3d.m1[p][q][r][2]^2 + s3d.m1[p][q][r][3]
s3d.m4[][][][3] -= s3d.m2[p][q][r][0]^2 + s3d.m2[p][q][r][1]^2 + s3d.m2[p][q][r][2]^2 + s3d.m2[p][q][r][3]
s3d.m4[][][][3] -= s3d.m3[p][q][r][0]^2 + s3d.m3[p][q][r][1]^2 + s3d.m3[p][q][r][2]^2 + s3d.m3[p][q][r][3]
duplicate /o s3d.density1,s3d.density4
s3d.density4[][][] = s3d.m4[p][q][r][3]
endif
elseif(i==4)
duplicate /o s3d.m1,s3d.m4
s3d.m4 = 0
if(i< nummaterials - 1)
s3d.m5[][][][3] = s3d.density5[p][q][r] // this shouldn't exist, we are limited to 5 materials so 4 densities in the current code
else
//this is the last material, so density wave doesn't exist, need to calculate it
s3d.m5[][][][3] = 1
s3d.m5[][][][3] -= s3d.m1[p][q][r][0]^2 + s3d.m1[p][q][r][1]^2 + s3d.m1[p][q][r][2]^2 + s3d.m1[p][q][r][3]
s3d.m5[][][][3] -= s3d.m2[p][q][r][0]^2 + s3d.m2[p][q][r][1]^2 + s3d.m2[p][q][r][2]^2 + s3d.m2[p][q][r][3]
s3d.m5[][][][3] -= s3d.m3[p][q][r][0]^2 + s3d.m3[p][q][r][1]^2 + s3d.m3[p][q][r][2]^2 + s3d.m3[p][q][r][3]
s3d.m5[][][][3] -= s3d.m4[p][q][r][0]^2 + s3d.m4[p][q][r][1]^2 + s3d.m4[p][q][r][2]^2 + s3d.m4[p][q][r][3]
duplicate /o s3d.density1,s3d.density5
s3d.density5[][][] = s3d.m5[p][q][r][3]
endif
endif
if(stringmatch("*Face-on*", mat) || stringmatch(mat, "1*"))
align = 1
elseif(stringmatch("*Edge-on*", mat) || stringmatch(mat, "0*"))
align=0
else
//don't do the calculation for this material, the previously made density and m matrix is sufficient
continue // move to next material
endif
splitstring/e="^[^,]*,([^,]*),([^,]*)" mat, alignmentwidth, alignmentstrength
strength = str2num(alignmentstrength)
width = str2num(alignmentwidth)
volfrac = mean(s3d.density1)
s3d.timesofar += stopmstimer(s3d.timer)/1e6
s3d.timer = startmstimer
if(i==0)
wave alignmentw = createalignmentdensity(s3d.density1,width,strength, align, s3d)
duplicate/o alignmentw, s3d.m1
s3d.m1 *=sqrt(s3d.density1[p][q][r])
s3d.m1[][][][3] *=sqrt(s3d.density1[p][q][r]) // last dimension is density not sqrt density, so it needs an extra multiplication
elseif(i==1)
wave alignmentw = createalignmentdensity(s3d.density2,width,strength, align, s3d)
duplicate/o alignmentw, s3d.m2
s3d.m2 *=sqrt(s3d.density2[p][q][r])
s3d.m2[][][][3] *=sqrt(s3d.density2[p][q][r])
elseif(i==2)
wave alignmentw = createalignmentdensity(s3d.density3,width,strength, align, s3d)
duplicate/o alignmentw, s3d.m3
s3d.m3 *=sqrt(s3d.density3[p][q][r])
s3d.m3[][][][3] *=sqrt(s3d.density3[p][q][r])
elseif(i==3)
wave alignmentw = createalignmentdensity(s3d.density4,width,strength, align, s3d)
duplicate/o alignmentw, s3d.m4
s3d.m4 *=sqrt(s3d.density4[p][q][r])
s3d.m4[][][][3] *=sqrt(s3d.density4[p][q][r])
elseif(i==4)
wave alignmentw = createalignmentdensity(s3d.density5,width,strength, align, s3d)
duplicate/o alignmentw, s3d.m5
s3d.m5 *=sqrt(s3d.density5[p][q][r])
s3d.m5[][][][3] *=sqrt(s3d.density5[p][q][r])
endif
endfor
endif
end
function storeintegrations(s3d)
struct ThreeDSystem &s3d
s3d.int3DvsEn[][s3d.step] = s3d.int3d[p] * x^2
s3d.ratio3DvsEn[][s3d.step] = (s3d.int3dpara(x) - s3d.int3dperp(x)) / (s3d.int3dpara(x) + s3d.int3dperp(x))
s3d.scatter3dSave[][][s3d.step] = s3d.scatter3d(x)(y)
s3d.perp3Dvsen[][s3d.step] = s3d.int3dperp(x)
s3d.para3Dvsen[][s3d.step] = s3d.int3dpara(x)
//s3d.pz3dSave[][][][s3d.step] = s3d.pz[p][q][r]
end
function radialintegratesystem(s3d)
//radially integrate the 2D and 3D transforms and populate the corresponding waves in s (see structure declaration for details)
struct ThreeDSystem &s3d
variable ey =0// str2num(stringfromlist(0,s3d.efield,","))
variable ez =1// str2num(stringfromlist(1,s3d.efield,","))
variable pe = atan(ey/ez)*180/pi
variable pa = atan(-ez/ey)*180/pi
wave s3d.int3dp0 = radialintegratehdf(s3d.scatter3d,0,90,0,pi/s3d.size,floor(s3d.num/sqrt(2)),"int3dp0")
wave s3d.int3dp0para = radialintegratehdf(s3d.scatter3d,pa-20,pa+20,0,pi/s3d.size,floor(s3d.num/sqrt(2)),"int3dp0para")
wave s3d.int3dp0perp = radialintegratehdf(s3d.scatter3d,pe-20,pe+20,0,pi/s3d.size,floor(s3d.num/sqrt(2)),"int3dp0perp")
duplicate/o s3d.int3dp0, s3d.int3d
duplicate/o s3d.int3dp0para, s3d.int3dpara
duplicate/o s3d.int3dp0perp, s3d.int3dperp
// s3d.int3d *= x^2
// s3d.int3dpara *= x^2
// s3d.int3dperp *= x^2
end
function Calculate3DScatter(s3d) // turn the polarization fields into scattering fields (real space to fourier space)
struct ThreeDSystem &s3d
variable thicknum = dimsize(s3d.px,0) , num = s3d.num
if(num/2 != round(num/2) )
num +=1
endif
if(thicknum/2 != round(thicknum/2) )
thicknum +=1
endif
wave px = s3d.px, py=s3d.py, pz=s3d.pz
// each element is fourier transformed seperatly, the different components only add incoherently
FFT /WINF=Hanning /pad={1*thicknum,1*num,1*num} /Dest=pxfft px
FFT /WINF=Hanning /pad={1*thicknum,1*num,1*num} /Dest=pyfft py
FFT /WINF=Hanning /pad={1*thicknum,1*num,1*num} /Dest=pzfft pz
wave/c s3d.pxfft = pxfft, s3d.pyfft=pyfft, s3d.pzfft=pzfft
//make /n=(dimsize(s3d.pxfft,0),dimsize(s3d.pxfft,1),dimsize(s3d.pxfft,2),9) /d/o s3d.qtensor
make /n=(dimsize(s3d.pxfft,0),dimsize(s3d.pxfft,1),dimsize(s3d.pxfft,2)) /d/o s3d.EscatSqr =0
wave EscatSqr=s3d.EscatSqr
setscale /i x, -pi/s3d.voxelsize, pi/s3d.voxelsize, pxfft, pyfft, pzfft, EscatSqr // changing from math units (1/voxelsize) to physics units (2pi/voxelsize)
setscale /i y, -pi/s3d.voxelsize, pi/s3d.voxelsize, pxfft, pyfft, pzfft, EscatSqr
setscale /i z, -pi/s3d.voxelsize, pi/s3d.voxelsize, pxfft, pyfft, pzfft, EscatSqr // all three dimensions have the same range. The x dimension may have a lower number of pixels, but they cover the same range
variable k = s3d.k
multithread EscatSqr = magsqr(x * (2 * k + x) * pxfft[p][q][r] + y * (k + x) * pyfft[p][q][r] + z * (k + x) * pzfft[p][q][r]) // Equation 8 from paper (qx -> -qx ends up equivalent somehow - checked with mathematica)
multithread EscatSqr += magsqr(y * (k + x)* pxfft[p][q][r]+ (y^2 - k^2) * pyfft[p][q][r] + y * z * pzfft[p][q][r])
multithread EscatSqr += magsqr(z * (k + x)* pxfft[p][q][r]+ y * z* pyfft[p][q][r]+ (z^2 - k^2 ) * pzfft[p][q][r])
make /d/n=(s3d.num,s3d.num)/o s3d.scatter3D // this will hold the final scattering result
wave scatter3d = s3d.scatter3d
setscale /i x, -2*pi/s3d.voxelsize, 2*pi/s3d.voxelsize, s3d.scatter3D // changing from math units (1/voxelsize) to physics units (2pi/voxelsize)
setscale /i y, -2*pi/s3d.voxelsize, 2*pi/s3d.voxelsize, s3d.scatter3D
multithread scatter3D = k^2 > x^2 + y^2 ? interp3d(EscatSqr,-k + sqrt(k^2 - x^2 -y^2),x,y) : nan // projection to the EWalds sphere
duplicate /o s3d.EscatSqr, pxreal
//pxreal = magsqr(s3d.pxfft)
imagetransform/meth=2 xprojection, pxreal
duplicate/o m_xprojection, pxfftproj
return 0
end
function CalculatePolarizationWave(s3d)
struct ThreeDSystem &s3d
variable result
//create polarization waves
make /o/c/d/n=(dimsize(s3d.m1,0),s3d.num,s3d.num) s3d.px = 0, s3d.py = 0, s3d.pz = 0
setscale /p x, dimoffset(s3d.m1,0), dimdelta(s3d.m1,0), s3d.px, s3d.py, s3d.pz
setscale /p y, dimoffset(s3d.m1,1), dimdelta(s3d.m1,1), s3d.px, s3d.py, s3d.pz