1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | |
---|
3 | // |
---|
4 | // utility functions and procedures for displaying information |
---|
5 | // setting the matrix, and doing the calculations |
---|
6 | // |
---|
7 | |
---|
8 | |
---|
9 | // |
---|
10 | // TO DO: |
---|
11 | // x- I need to change a lot of routines (most notably Gizmo) to take "10" as the default SLD |
---|
12 | // rather than "1" |
---|
13 | // |
---|
14 | // |
---|
15 | // -- incorporate utility function that reads the wave note from a binary file |
---|
16 | // Function/S LoadNoteFunc(PName,FName[,FileRef]) - use this as an "inspector?" |
---|
17 | //////////// |
---|
18 | |
---|
19 | // |
---|
20 | |
---|
21 | //*********** |
---|
22 | //For a 3D slicer view of dmat: |
---|
23 | // |
---|
24 | //Duplicate/O dmat dmatView |
---|
25 | // |
---|
26 | //Open a New 3D voxelgram. Don't set any of the levels |
---|
27 | //Open the 3D Slicer under the Gizmo menu |
---|
28 | //Make the X, Y, and Z slices |
---|
29 | //(Yellow Hot seems to be the best) |
---|
30 | // |
---|
31 | //// for a better view |
---|
32 | //dmatView = log(dmat) |
---|
33 | // |
---|
34 | //// if it all goes blank after the log transform get rid of the INF |
---|
35 | //dmatView = numtype(dmatView)== 0 ? dmatView : 0 |
---|
36 | // |
---|
37 | //************* |
---|
38 | |
---|
39 | |
---|
40 | ///PANEL |
---|
41 | ///// panel procedures |
---|
42 | |
---|
43 | Proc Init_FFT() |
---|
44 | DoWindow/F FFT_Panel |
---|
45 | if(V_flag==0) |
---|
46 | SetDataFolder root: |
---|
47 | //init the globals |
---|
48 | Variable/G root:FFT_N=128 |
---|
49 | Variable/G root:FFT_T=5 |
---|
50 | Variable/G root:FFT_Qmax = 0 |
---|
51 | Variable/G root:FFT_QmaxReal = 0 |
---|
52 | Variable/G root:FFT_DQ=0 |
---|
53 | Variable/G root:FFT_Qmin=0 |
---|
54 | Variable/G root:FFT_estTime = 0 |
---|
55 | |
---|
56 | Variable/G root:FFT_SolventSLD = 0 |
---|
57 | Variable/G root:FFT_delRho = 1e-7 //multiplier for SLD (other value is 1e-7) |
---|
58 | |
---|
59 | FFT_Qmax :=2*pi/FFT_T |
---|
60 | FFT_QmaxReal := FFT_Qmax/2 |
---|
61 | FFT_DQ := pi/(FFT_N*FFT_T) |
---|
62 | FFT_Qmin := 2*pi/(FFT_N*FFT_T) |
---|
63 | //empirical fit (cubic) of time vs N=50 to N=256 |
---|
64 | FFT_estTime := 0.56 - 0.0156*FFT_N + 0.000116*FFT_N^2 + 8e-7*FFT_N^3 |
---|
65 | // FFT_estTime := FFT_N/128 |
---|
66 | |
---|
67 | FFT_Panel() |
---|
68 | endif |
---|
69 | End |
---|
70 | |
---|
71 | |
---|
72 | Proc FFT_Panel() |
---|
73 | PauseUpdate; Silent 1 // building window... |
---|
74 | NewPanel /W=(1452,44,1768,531)/K=1 as "FFT_Panel" |
---|
75 | DoWindow/C FFT_Panel |
---|
76 | SetDrawLayer UserBack |
---|
77 | DrawLine 5,68,311,68 |
---|
78 | DrawLine 5,142,311,142 |
---|
79 | DrawLine 5,250,311,250 |
---|
80 | SetVariable FFTSetVar_0,pos={7,7},size={150,15},title="Cells per edge (N)" |
---|
81 | SetVariable FFTSetVar_0,limits={50,512,2},value= FFT_N,live= 1 |
---|
82 | SetVariable FFTSetVar_1,pos={7,27},size={150,15},title="Length per Cell (T)" |
---|
83 | SetVariable FFTSetVar_1,limits={1,5000,0.2},value= FFT_T,live= 1 |
---|
84 | SetVariable FFTSetVar_2,pos={183,7},size={120,15},title="Real Qmax" |
---|
85 | SetVariable FFTSetVar_2,limits={0,0,0},value= FFT_QmaxReal,noedit= 1,live= 1 |
---|
86 | SetVariable FFTSetVar_3,pos={183,47},size={120,15},title="delta Q (A)" |
---|
87 | SetVariable FFTSetVar_3,limits={0,0,0},value= FFT_DQ,noedit= 1,live= 1 |
---|
88 | SetVariable FFTSetVar_6,pos={183,27},size={120,15},title="Real Qmin (A)" |
---|
89 | SetVariable FFTSetVar_6,limits={0,0,0},value= FFT_Qmin,noedit= 1,live= 1 |
---|
90 | Button FFTButton_0,pos={15,79},size={90,20},proc=FFT_MakeMatrixButtonProc,title="Make Matrix" |
---|
91 | |
---|
92 | Button FFTButton_3,pos={14,265},size={70,20},proc=DoTheFFT_ButtonProc,title="Do FFT" |
---|
93 | Button FFTButton_4,pos={180,264},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot FFT Results" |
---|
94 | |
---|
95 | Button FFTButton_1,pos={14,150},size={90,20},proc=FFTMakeGizmoButtonProc,title="Make Gizmo" |
---|
96 | Button FFTButton_2,pos={14,175},size={100,20},proc=FFTDrawSphereButtonProc,title="Draw Sphere" |
---|
97 | Button FFTButton_5,pos={14,200},size={130,20},proc=FFTDrawZCylinderButtonProc,title="Draw XYZ Cylinder" |
---|
98 | Button FFTButton_20,pos={14,225},size={130,20},proc=FFTDrawRotCylinderButton,title="Draw Rot Cylinder" |
---|
99 | |
---|
100 | // Button FFTButton_6,pos={134,79},size={90,20},proc=FFTEraseMatrixButtonProc,title="Erase Matrix" |
---|
101 | Button FFTButton_6a,pos={140,79},size={50,20},proc=FFTSaveMatrixButtonProc,title="Save" |
---|
102 | Button FFTButton_6b,pos={200,79},size={50,20},proc=FFTLoadMatrixButtonProc,title="Load" |
---|
103 | Button FFTButton_6c,pos={260,79},size={50,20},proc=FFT_AddMatrixButtonProc,title="Add" |
---|
104 | Button FFTButton_7,pos={13,329},size={130,20},proc=FFT_BinnedSpheresButtonProc,title="Do Binned Debye" |
---|
105 | Button FFTButton_7a,pos={180,329},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot Binned Results" |
---|
106 | |
---|
107 | Button FFTButton_8,pos={13,297},size={130,20},proc=FFT_AltiSpheresButtonProc,title="Do Debye Spheres" |
---|
108 | Button FFTButton_8a,pos={180,297},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot Debye Results" |
---|
109 | |
---|
110 | Button FFTButton_14,pos={13,360},size={130,20},proc=FFT_BinnedSLDButtonProc,title="Do Binned SLD" |
---|
111 | Button FFTButton_14a,pos={180,360},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot SLD Results" |
---|
112 | |
---|
113 | SetVariable FFTSetVar_4,pos={7,47},size={100,15},title="FFT time(s)" |
---|
114 | SetVariable FFTSetVar_4,limits={0,0,0},value= FFT_estTime,noedit= 1,live= 1,format="%g" |
---|
115 | Button FFTButton_9,pos={200,400},size={100,20},proc=FFT_Get2DSlice,title="Get 2D Slice" |
---|
116 | |
---|
117 | Button FFTButton_19,pos={168,150},size={130,20},proc=FFT_ChangeMatrixValuesButton,title="Replace Voxels" |
---|
118 | Button FFTButton_12,pos={168,175},size={130,20},proc=FFT_ReplaceSolventButton,title="Replace Solvent" |
---|
119 | Button FFTButton_11,pos={169,200},size={130,20},proc=FFT_RotateMat,title="Rotate Matrix" |
---|
120 | Button FFTButton_10,pos={169,225},size={130,20},proc=FFT_TransposeMat,title="Transpose Matrix" |
---|
121 | |
---|
122 | Button FFTButton_13,pos={14,109},size={120,20},proc=FFTFillSolventMatrixProc,title="Solvent Matrix" |
---|
123 | SetVariable FFTSetVar_5,pos={155,111},size={150,15},title="Solvent SLD (10^-7)" |
---|
124 | SetVariable FFTSetVar_5,limits={-99,99,1},value= FFT_SolventSLD,live= 1 |
---|
125 | Button FFTButton_15,pos={209,430},size={90,20},proc=Interp2DSliceButton,title="Interp 2D" |
---|
126 | Button FFTButton_16,pos={14,460},size={70,20},proc=FFTHelpButton,title="Help" |
---|
127 | |
---|
128 | Button FFTButton_17,pos={13,400},size={120,20},proc=FFT_Iso2USANS,title="Iso to USANS" |
---|
129 | Button FFTButton_18,pos={13,430},size={120,20},proc=FFT_Aniso2USANS,title="Aniso to USANS" |
---|
130 | |
---|
131 | EndMacro |
---|
132 | |
---|
133 | // Save a matrix wave, plus the N, T, and solvent values in the wave note for reloading |
---|
134 | Function FFTDrawRotCylinderButton(ba) : ButtonControl |
---|
135 | STRUCT WMButtonAction &ba |
---|
136 | |
---|
137 | String win = ba.win |
---|
138 | |
---|
139 | switch (ba.eventCode) |
---|
140 | case 2: |
---|
141 | // click code here |
---|
142 | Execute "FFTDrawRotCylinder()" |
---|
143 | |
---|
144 | break |
---|
145 | endswitch |
---|
146 | |
---|
147 | return 0 |
---|
148 | End |
---|
149 | |
---|
150 | // Save a matrix wave, plus the N, T, and solvent values in the wave note for reloading |
---|
151 | Function FFTSaveMatrixButtonProc(ba) : ButtonControl |
---|
152 | STRUCT WMButtonAction &ba |
---|
153 | |
---|
154 | String win = ba.win |
---|
155 | |
---|
156 | switch (ba.eventCode) |
---|
157 | case 2: |
---|
158 | // click code here |
---|
159 | String fileStr="" |
---|
160 | SaveMyMatrix(fileStr) |
---|
161 | |
---|
162 | break |
---|
163 | endswitch |
---|
164 | |
---|
165 | return 0 |
---|
166 | End |
---|
167 | |
---|
168 | // this will wave as Igor Binary, so be sure to use the ".ibw extension. |
---|
169 | // - this could possibly be enforced, but that's maybe not necessary at this stage. |
---|
170 | // |
---|
171 | Function SaveMyMatrix(fileStr) |
---|
172 | String fileStr |
---|
173 | |
---|
174 | WAVE mat=root:mat |
---|
175 | NVAR FFT_T = root:FFT_T |
---|
176 | NVAR FFT_N = root:FFT_N |
---|
177 | NVAR FFT_SolventSLD = root:FFT_SolventSLD |
---|
178 | String str="" |
---|
179 | sprintf str,"FFT_T=%g;FFT_N=%d;FFT_SolventSLD=%d;",FFT_T,FFT_N,FFT_SolventSLD |
---|
180 | Note mat,str |
---|
181 | Save/C/P=home mat as fileStr //will ask for a file name if fileStr="" save as Igor Binary |
---|
182 | Note/K mat //kill wave note on exiting since I don't properly update this anywhere else |
---|
183 | |
---|
184 | return(0) |
---|
185 | end |
---|
186 | |
---|
187 | |
---|
188 | // load in a previously saved matrix, and reset FFT_N, FFT_T and solvent |
---|
189 | // from the wave note when saved |
---|
190 | Function FFTLoadMatrixButtonProc(ba) : ButtonControl |
---|
191 | STRUCT WMButtonAction &ba |
---|
192 | |
---|
193 | String win = ba.win |
---|
194 | |
---|
195 | switch (ba.eventCode) |
---|
196 | case 2: |
---|
197 | // click code here |
---|
198 | String fileStr="" |
---|
199 | ReloadMatrix(fileStr) |
---|
200 | |
---|
201 | break |
---|
202 | endswitch |
---|
203 | |
---|
204 | return 0 |
---|
205 | End |
---|
206 | |
---|
207 | // /H flag on the LoadWave command severs the connection with the binary file |
---|
208 | // -- this seems to be important - otherwise I get odd results and the wave (on disk) can change! |
---|
209 | // |
---|
210 | Function ReloadMatrix(fileStr) |
---|
211 | String fileStr |
---|
212 | |
---|
213 | LoadWave/H/M/O/W/P=home fileStr //will ask for a file, Igor Binary format is assumed here |
---|
214 | if(V_flag == 0) |
---|
215 | return(0) //user cancel |
---|
216 | endif |
---|
217 | |
---|
218 | String str |
---|
219 | str=note(mat) |
---|
220 | NVAR FFT_T = root:FFT_T |
---|
221 | NVAR FFT_N = root:FFT_N |
---|
222 | NVAR FFT_SolventSLD = root:FFT_SolventSLD |
---|
223 | |
---|
224 | FFT_T = NumberByKey("FFT_T", str, "=" ,";") |
---|
225 | FFT_N = NumberByKey("FFT_N", str, "=" ,";") |
---|
226 | FFT_SolventSLD = NumberByKey("FFT_SolventSLD", str, "=" ,";") |
---|
227 | |
---|
228 | // if I got bad values, put in default values |
---|
229 | if(numtype(FFT_T) != 0 ) |
---|
230 | FFT_T = 5 |
---|
231 | endif |
---|
232 | if(numtype(FFT_N) != 0 ) |
---|
233 | FFT_N = DimSize(mat,0) |
---|
234 | endif |
---|
235 | if(numtype(FFT_SolventSLD) != 0 ) |
---|
236 | FFT_SolventSLD = 0 |
---|
237 | endif |
---|
238 | |
---|
239 | ColorizeGizmo() |
---|
240 | |
---|
241 | Print "Loaded matrix parameters = ",str |
---|
242 | Execute "NumberOfPoints()" |
---|
243 | |
---|
244 | return(0) |
---|
245 | end |
---|
246 | |
---|
247 | // load in a previously saved matrix, and reset FFT_N, FFT_T and solvent |
---|
248 | // from the wave note when saved |
---|
249 | Function FFT_AddMatrixButtonProc(ba) : ButtonControl |
---|
250 | STRUCT WMButtonAction &ba |
---|
251 | |
---|
252 | String win = ba.win |
---|
253 | |
---|
254 | switch (ba.eventCode) |
---|
255 | case 2: |
---|
256 | // click code here |
---|
257 | String fileStr="" |
---|
258 | LoadAndAddMatrix(fileStr) |
---|
259 | |
---|
260 | break |
---|
261 | endswitch |
---|
262 | |
---|
263 | return 0 |
---|
264 | End |
---|
265 | |
---|
266 | // This function will load and add the matrix to whatever is currently present |
---|
267 | // - it is checked that the N, T, and SolventSLD are the same |
---|
268 | // |
---|
269 | // - it currently SUMS the voxels - so if there is overlap, you get a new value |
---|
270 | // |
---|
271 | Function LoadAndAddMatrix(fileStr) |
---|
272 | String fileStr |
---|
273 | |
---|
274 | String toLoadStr="" |
---|
275 | |
---|
276 | // if the passed in file name is null, pick a file |
---|
277 | if(strlen(fileStr)==0) |
---|
278 | toLoadStr = DoOpenFileDialog("Pick the binary file to load") |
---|
279 | if(strlen(toLoadStr)==0) |
---|
280 | return(0) |
---|
281 | endif |
---|
282 | endif |
---|
283 | fileStr=toLoadStr |
---|
284 | |
---|
285 | // make sure that N and T are correct, just read in the note |
---|
286 | String noteStr=LoadNoteFunc("home",fileStr) |
---|
287 | String abortStr |
---|
288 | |
---|
289 | // current values |
---|
290 | NVAR FFT_T = root:FFT_T |
---|
291 | NVAR FFT_N = root:FFT_N |
---|
292 | NVAR FFT_SolventSLD = root:FFT_SolventSLD |
---|
293 | |
---|
294 | // possible set to load |
---|
295 | Variable test_T,test_N,test_SolventSLD |
---|
296 | test_T = NumberByKey("FFT_T", noteStr, "=" ,";") |
---|
297 | test_N = NumberByKey("FFT_N", noteStr, "=" ,";") |
---|
298 | test_SolventSLD = NumberByKey("FFT_SolventSLD", noteStr, "=" ,";") |
---|
299 | |
---|
300 | // if I got bad values, warn and abort |
---|
301 | if(FFT_T != test_T) |
---|
302 | abortStr = "Current T = "+num2str(FFT_T)+" and Selected T = "+num2str(test_T)+", aborting load" |
---|
303 | Abort abortStr |
---|
304 | endif |
---|
305 | if(FFT_N != test_N) |
---|
306 | abortStr = "Current N = "+num2str(FFT_N)+" and Selected N = "+num2str(test_N)+", aborting load" |
---|
307 | Abort abortStr |
---|
308 | endif |
---|
309 | if(FFT_SolventSLD != test_SolventSLD) |
---|
310 | abortStr = "Current SolventSLD = "+num2str(FFT_SolventSLD)+" and Selected SolventSLD = "+num2str(test_SolventSLD)+", aborting load" |
---|
311 | Abort abortStr |
---|
312 | endif |
---|
313 | |
---|
314 | // OK, then load and add the matrix |
---|
315 | // the loaded matrix is "mat" as usual, so make a copy of the current first |
---|
316 | Duplicate/O mat tmpMat |
---|
317 | LoadWave/H/M/O/W/P=home fileStr //will ask for a file, Igor Binary format is assumed here |
---|
318 | if(V_flag == 0) |
---|
319 | return(0) //user cancel |
---|
320 | endif |
---|
321 | |
---|
322 | Wave mat=root:mat |
---|
323 | FastOp mat = mat + tmpMat |
---|
324 | |
---|
325 | KillWaves/Z tmpMat |
---|
326 | |
---|
327 | Print "Loaded matrix parameters = ",noteStr |
---|
328 | Execute "NumberOfPoints()" |
---|
329 | |
---|
330 | ColorizeGizmo() |
---|
331 | |
---|
332 | return(0) |
---|
333 | end |
---|
334 | |
---|
335 | // utility function that reads the wave note from a binary file |
---|
336 | // where did I get this from? Igor Exchange? I didn't write this myself... |
---|
337 | // |
---|
338 | Function/S LoadNoteFunc(PName,FName[,FileRef]) |
---|
339 | String PName, FName |
---|
340 | Variable FileRef |
---|
341 | |
---|
342 | Variable noteStart, noteLength, version, dependLength |
---|
343 | String noteStr, typeStr = ".ibw" |
---|
344 | if (ParamIsDefault(FileRef)) |
---|
345 | Open/R/P=$PName/T=typeStr fileRef, as FName //open the file |
---|
346 | endif |
---|
347 | FSetPos fileRef, 0 |
---|
348 | FBinRead/F=2 fileRef, version |
---|
349 | |
---|
350 | |
---|
351 | if (version == 5) |
---|
352 | |
---|
353 | FSetPos fileRef, 4 |
---|
354 | Make/N=(3)/I/Free SizeWave |
---|
355 | FBinRead FileRef,SizeWave |
---|
356 | noteStart = SizeWave[0] |
---|
357 | DependLength = SizeWave[1] |
---|
358 | NoteLength = SizeWave[2] |
---|
359 | noteStart += dependLength+64 |
---|
360 | |
---|
361 | elseif (version == 2) |
---|
362 | |
---|
363 | FBinRead/F=3 fileRef, noteStart |
---|
364 | // FBinRead/F=4 fileRef, dependLength |
---|
365 | FBinRead/F=3 fileRef, noteLength |
---|
366 | noteStart += 16 |
---|
367 | |
---|
368 | else |
---|
369 | |
---|
370 | if (ParamIsDefault(FileRef)) |
---|
371 | Close(FileRef) //close the file |
---|
372 | endif |
---|
373 | return "" |
---|
374 | |
---|
375 | endif |
---|
376 | if (!NoteLength) |
---|
377 | if (ParamIsDefault(FileRef)) |
---|
378 | Close(FileRef) //close the file |
---|
379 | endif |
---|
380 | return("") |
---|
381 | endif |
---|
382 | FSetPos fileRef, noteStart |
---|
383 | NoteStr = PadString("",NoteLength,0) |
---|
384 | FBinRead FileRef,NoteStr |
---|
385 | |
---|
386 | if (ParamIsDefault(FileRef)) |
---|
387 | Close(FileRef) //close the file |
---|
388 | endif |
---|
389 | return noteStr |
---|
390 | |
---|
391 | End //LoadNoteFunc |
---|
392 | |
---|
393 | |
---|
394 | Function FFTHelpButton(ba) : ButtonControl |
---|
395 | STRUCT WMButtonAction &ba |
---|
396 | |
---|
397 | String win = ba.win |
---|
398 | |
---|
399 | switch (ba.eventCode) |
---|
400 | case 2: |
---|
401 | // click code here |
---|
402 | DisplayHelpTopic/Z/K=1 "Real-Space Modeling of SANS Data" |
---|
403 | if(V_flag !=0) |
---|
404 | DoAlert 0,"The Real-Space Modeling Help file could not be found" |
---|
405 | endif |
---|
406 | break |
---|
407 | endswitch |
---|
408 | |
---|
409 | return 0 |
---|
410 | End |
---|
411 | |
---|
412 | |
---|
413 | Function Interp2DSliceButton(ba) : ButtonControl |
---|
414 | STRUCT WMButtonAction &ba |
---|
415 | |
---|
416 | switch( ba.eventCode ) |
---|
417 | case 2: // mouse up |
---|
418 | // click code here |
---|
419 | String folderStr="" |
---|
420 | Prompt folderStr, "Pick a 2D data folder",popup,W_DataSetPopupList() // Set prompt for x param |
---|
421 | DoPrompt "Enter data folder", folderStr |
---|
422 | if (V_Flag || strlen(folderStr) == 0) |
---|
423 | return -1 // User canceled, null string entered |
---|
424 | endif |
---|
425 | |
---|
426 | Interpolate2DSliceToData(folderStr) |
---|
427 | |
---|
428 | //display the 2D plane |
---|
429 | DoWindow FFT_Interp2D_log |
---|
430 | if(V_flag == 0) |
---|
431 | Execute "Display2DInterpSlice_log()" |
---|
432 | endif |
---|
433 | |
---|
434 | break |
---|
435 | endswitch |
---|
436 | |
---|
437 | return 0 |
---|
438 | End |
---|
439 | |
---|
440 | |
---|
441 | |
---|
442 | Function FFT_Get2DSlice(ctrlName) : ButtonControl |
---|
443 | String ctrlName |
---|
444 | |
---|
445 | WAVE/Z dmat = root:dmat |
---|
446 | if(WaveExists(dmat)==1) |
---|
447 | get2DSlice(dmat) |
---|
448 | endif |
---|
449 | |
---|
450 | //display the 2D plane |
---|
451 | DoWindow FFT_Slice2D |
---|
452 | if(V_flag == 0) |
---|
453 | Execute "Display2DSlice()" |
---|
454 | endif |
---|
455 | |
---|
456 | //display the 2D plane |
---|
457 | DoWindow FFT_Slice2D_log |
---|
458 | if(V_flag == 0) |
---|
459 | Execute "Display2DSlice_log()" |
---|
460 | endif |
---|
461 | |
---|
462 | //display the 1D binning of the 2D plane |
---|
463 | DoWindow FFT_Slice1D |
---|
464 | if(V_flag == 0) |
---|
465 | Execute "Slice2_1D()" |
---|
466 | endif |
---|
467 | |
---|
468 | return(0) |
---|
469 | End |
---|
470 | |
---|
471 | Proc Display2DSlice() |
---|
472 | PauseUpdate; Silent 1 // building window... |
---|
473 | Display /W=(1038,44,1404,403) |
---|
474 | DoWindow/C FFT_Slice2D |
---|
475 | AppendImage/T detPlane |
---|
476 | ModifyImage detPlane ctab= {*,*,YellowHot,0} |
---|
477 | ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14 |
---|
478 | ModifyGraph mirror=2 |
---|
479 | ModifyGraph nticks=4 |
---|
480 | ModifyGraph minor=1 |
---|
481 | ModifyGraph fSize=9 |
---|
482 | ModifyGraph standoff=0 |
---|
483 | ModifyGraph tkLblRot(left)=90 |
---|
484 | ModifyGraph btLen=3 |
---|
485 | ModifyGraph tlOffset=-2 |
---|
486 | // SetAxis/A/R left |
---|
487 | End |
---|
488 | |
---|
489 | Proc Display2DInterpSlice_log() |
---|
490 | PauseUpdate; Silent 1 // building window... |
---|
491 | Display /W=(1038,44,1404,403) |
---|
492 | DoWindow/C FFT_Interp2D_log |
---|
493 | AppendImage/T interp2DSlice_log |
---|
494 | ModifyImage interp2DSlice_log ctab= {*,*,YellowHot,0} |
---|
495 | ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14 |
---|
496 | ModifyGraph mirror=2 |
---|
497 | ModifyGraph nticks=4 |
---|
498 | ModifyGraph minor=1 |
---|
499 | ModifyGraph fSize=9 |
---|
500 | ModifyGraph standoff=0 |
---|
501 | ModifyGraph tkLblRot(left)=90 |
---|
502 | ModifyGraph btLen=3 |
---|
503 | ModifyGraph tlOffset=-2 |
---|
504 | // SetAxis/A/R left |
---|
505 | End |
---|
506 | |
---|
507 | Proc Display2DSlice_log() |
---|
508 | PauseUpdate; Silent 1 // building window... |
---|
509 | Display /W=(1038,44,1404,403) |
---|
510 | DoWindow/C FFT_Slice2D_log |
---|
511 | AppendImage/T logP |
---|
512 | ModifyImage logP ctab= {*,*,YellowHot,0} |
---|
513 | ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14 |
---|
514 | ModifyGraph mirror=2 |
---|
515 | ModifyGraph nticks=4 |
---|
516 | ModifyGraph minor=1 |
---|
517 | ModifyGraph fSize=9 |
---|
518 | ModifyGraph standoff=0 |
---|
519 | ModifyGraph tkLblRot(left)=90 |
---|
520 | ModifyGraph btLen=3 |
---|
521 | ModifyGraph tlOffset=-2 |
---|
522 | // SetAxis/A/R left |
---|
523 | End |
---|
524 | Proc Slice2_1D() |
---|
525 | PauseUpdate; Silent 1 // building window... |
---|
526 | Display /W=(1034,425,1406,763) iBin_2d vs qBin_2d |
---|
527 | DoWindow/C FFT_Slice1D |
---|
528 | ModifyGraph mode=4 |
---|
529 | ModifyGraph marker=19 |
---|
530 | ModifyGraph msize=2 |
---|
531 | ModifyGraph grid=1 |
---|
532 | ModifyGraph log=1 |
---|
533 | ModifyGraph mirror=2 |
---|
534 | Legend |
---|
535 | EndMacro |
---|
536 | |
---|
537 | |
---|
538 | Function FFT_TransposeMat(ctrlName) : ButtonControl |
---|
539 | String ctrlName |
---|
540 | |
---|
541 | Variable mode=1 |
---|
542 | Prompt mode,"Transform XYZ to:",popup,"XZY;ZXY;ZYX;YZX;YXZ;" |
---|
543 | DoPrompt "Transform mode",mode |
---|
544 | if (V_Flag) |
---|
545 | return 0 // user canceled |
---|
546 | endif |
---|
547 | |
---|
548 | fFFT_TransposeMat(mode) |
---|
549 | |
---|
550 | return (0) |
---|
551 | End |
---|
552 | |
---|
553 | // starting in XYZ |
---|
554 | // mode 1;2;3;4;5; |
---|
555 | //correspond to |
---|
556 | // "XZY;ZXY;ZYX;YZX;YXZ;" |
---|
557 | // |
---|
558 | Function fFFT_TransposeMat(mode) |
---|
559 | Variable mode |
---|
560 | |
---|
561 | WAVE/Z mat=mat |
---|
562 | ImageTransform /G=(mode) transposeVol mat |
---|
563 | WAVE M_VolumeTranspose=M_VolumeTranspose |
---|
564 | Duplicate/O M_VolumeTranspose mat |
---|
565 | |
---|
566 | return(0) |
---|
567 | End |
---|
568 | |
---|
569 | Function FFT_RotateMat(ctrlName) : ButtonControl |
---|
570 | String ctrlName |
---|
571 | |
---|
572 | Variable angleX=45,angleY=0,angleZ=0 |
---|
573 | Prompt angleX, "Degrees of rotation around X-axis:" |
---|
574 | Prompt angleY, "Degrees of rotation around Y-axis:" |
---|
575 | Prompt angleZ, "Degrees of rotation around Z-axis:" |
---|
576 | DoPrompt "Enter angles for rotation", angleX,angleY,angleZ |
---|
577 | |
---|
578 | if (V_Flag) |
---|
579 | return 0 // user canceled |
---|
580 | endif |
---|
581 | |
---|
582 | XYZRotate(angleX,angleY,angleZ) |
---|
583 | |
---|
584 | |
---|
585 | //////////// old way, only one angle |
---|
586 | // Variable degree=45,sense=1 |
---|
587 | // Prompt degree, "Degrees of rotation around Z-axis:" |
---|
588 | //// Prompt sense, "Direction of rotation:",popup,"CW;CCW;" |
---|
589 | // DoPrompt "Enter parameters for rotation", degree |
---|
590 | // |
---|
591 | // if (V_Flag) |
---|
592 | // return 0 // user canceled |
---|
593 | // endif |
---|
594 | |
---|
595 | |
---|
596 | // old way using ImageRotate that interpolates, and is only around Z-axis |
---|
597 | // fFFT_RotateMat(degree) |
---|
598 | return(0) |
---|
599 | End |
---|
600 | |
---|
601 | // note that the rotation is not perfect. if the rotation produces an |
---|
602 | // odd number of pixels, then the object will "walk" one pixel. Also, some |
---|
603 | // small artifacts are introduced, simply due to the voxelization of the object |
---|
604 | // as it is rotated. rotating a cylinder 10 then 350 shows a few extra "bumps" on |
---|
605 | // the surface, but the calculated scattering is virtually identical. |
---|
606 | // |
---|
607 | // these issues, may be correctable, if needed |
---|
608 | // |
---|
609 | Function fFFT_RotateMat(degree) |
---|
610 | Variable degree |
---|
611 | |
---|
612 | Variable fill=0 |
---|
613 | NVAR solventSLD = root:FFT_SolventSLD |
---|
614 | fill = solventSLD |
---|
615 | |
---|
616 | WAVE mat = root:mat |
---|
617 | Variable dx,dy,dz,nx,ny,nz |
---|
618 | dx = DimSize(mat,0) |
---|
619 | dy = DimSize(mat,1) |
---|
620 | dz = DimSize(mat,2) |
---|
621 | |
---|
622 | //? the /W and /C flags seem to be special cases for 90 deg rotations |
---|
623 | ImageRotate /A=(degree)/E=(fill) mat //mat is not overwritten, unless /O |
---|
624 | |
---|
625 | nx = DimSize(M_RotatedImage,0) |
---|
626 | ny = DimSize(M_RotatedImage,1) |
---|
627 | nz = DimSize(M_RotatedImage,2) |
---|
628 | // Print "rotated dims = ",dx,dy,dz,nx,ny,nz |
---|
629 | Variable delx,dely,odd=0 |
---|
630 | delx = nx-dx |
---|
631 | dely = ny-dy |
---|
632 | |
---|
633 | if(mod(delx,2) != 0) //sometimes the new x,y dims are odd! |
---|
634 | odd = 1 |
---|
635 | endif |
---|
636 | // delx = (delx + odd)/2 |
---|
637 | // dely = (dely + odd)/2 |
---|
638 | delx = trunc(delx/2) |
---|
639 | dely = trunc(dely/2) |
---|
640 | |
---|
641 | //+odd removes an extra point if there is one |
---|
642 | Duplicate/O/R=[delx+odd,nx-delx-1][dely+odd,ny-dely-1][0,nz-1] M_RotatedImage mat |
---|
643 | |
---|
644 | // - not sure why the duplicate operation changes the start and delta of mat, but I |
---|
645 | // need to reset the start and delta to be sure that the display is correct, and that the scaling |
---|
646 | // is read correctly later |
---|
647 | SetScale/P x 0,1,"", mat |
---|
648 | SetScale/P y 0,1,"", mat |
---|
649 | |
---|
650 | nx = DimSize(mat,0) |
---|
651 | ny = DimSize(mat,1) |
---|
652 | nz = DimSize(mat,2) |
---|
653 | // Print "redim = ",dx,dy,dz,nx,ny,nz |
---|
654 | |
---|
655 | KillWaves/Z M_RotatedImage |
---|
656 | |
---|
657 | End |
---|
658 | |
---|
659 | // also look for ImageTransform operations that will allow translation of the objects. |
---|
660 | // then simple objects could be oriented @0,0,0, and translated to the correct position (and added) |
---|
661 | // |
---|
662 | Function FFT_AddRotatedObject(ctrlName) : ButtonControl |
---|
663 | String ctrlName |
---|
664 | |
---|
665 | Print "Not yet implemented" |
---|
666 | End |
---|
667 | |
---|
668 | Function FFT_ChangeMatrixValuesButton(ctrlName) |
---|
669 | String ctrlName |
---|
670 | |
---|
671 | Execute "ChangeMatrixValues()" |
---|
672 | |
---|
673 | end |
---|
674 | |
---|
675 | Function FFT_ReplaceSolventButton(ctrlname) |
---|
676 | String ctrlName |
---|
677 | |
---|
678 | Execute "ReplaceSolvent()" |
---|
679 | end |
---|
680 | |
---|
681 | Function FFT_MakeMatrixButtonProc(ctrlName) : ButtonControl |
---|
682 | String ctrlName |
---|
683 | |
---|
684 | NVAR nn=root:FFT_N |
---|
685 | if(mod(nn, 2 ) ==1) //force an even number for FFT |
---|
686 | nn +=1 |
---|
687 | endif |
---|
688 | MakeMatrix("mat",nn,nn,nn) |
---|
689 | End |
---|
690 | |
---|
691 | Function FFTMakeGizmoButtonProc(ctrlName) : ButtonControl |
---|
692 | String ctrlName |
---|
693 | |
---|
694 | DoWindow/F Gizmo_VoxelMat |
---|
695 | if(V_flag==0) |
---|
696 | Execute "Gizmo_VoxelMat()" |
---|
697 | endif |
---|
698 | ColorizeGizmo() |
---|
699 | End |
---|
700 | |
---|
701 | Function FFTDrawSphereButtonProc(ctrlName) : ButtonControl |
---|
702 | String ctrlName |
---|
703 | |
---|
704 | Execute "FFTDrawSphereProc()" |
---|
705 | End |
---|
706 | |
---|
707 | Proc FFTDrawSphereProc(matStr,rad,xc,yc,zc,fill,periodic) |
---|
708 | String matStr="mat" |
---|
709 | Variable rad=25,xc=50,yc=50,zc=50,fill=10,periodic=1 |
---|
710 | Prompt matStr,"the wave" //,popup,WaveList("*",";","") |
---|
711 | Prompt rad,"enter real radius (A)" |
---|
712 | Prompt xc,"enter the X-center" |
---|
713 | Prompt yc,"enter the Y-center" |
---|
714 | Prompt zc,"enter the Z-center" |
---|
715 | Prompt fill,"fill SLD value" |
---|
716 | Prompt periodic,"enter 1 for periodic, 0 for non-periodic fill" |
---|
717 | |
---|
718 | Variable grid=root:FFT_T |
---|
719 | |
---|
720 | if(periodic) |
---|
721 | FillSphereRadiusPeriodic($matStr,grid,rad,xc,yc,zc,fill) |
---|
722 | else |
---|
723 | FillSphereRadius($matStr,grid,rad,xc,yc,zc,fill) |
---|
724 | endif |
---|
725 | End |
---|
726 | |
---|
727 | Function FFTDrawZCylinderButtonProc(ctrlName) : ButtonControl |
---|
728 | String ctrlName |
---|
729 | |
---|
730 | Execute "FFTDrawCylinder()" |
---|
731 | End |
---|
732 | |
---|
733 | Proc FFTDrawCylinder(direction,matStr,rad,len,xc,yc,zc,fill) |
---|
734 | String direction |
---|
735 | String matStr="mat" |
---|
736 | Variable rad=25,len=300,xc=50,yc=50,zc=50,fill=10 |
---|
737 | Prompt direction, "Direction", popup "X;Y;Z;" |
---|
738 | Prompt matStr,"the wave" //,popup,WaveList("*",";","") |
---|
739 | Prompt rad,"enter real radius (A)" |
---|
740 | Prompt len,"enter length (A)" |
---|
741 | Prompt xc,"enter the X-center" |
---|
742 | Prompt yc,"enter the Y-center" |
---|
743 | Prompt zc,"enter the Z-center" |
---|
744 | Prompt fill,"fill SLD value" |
---|
745 | |
---|
746 | |
---|
747 | Variable grid=root:FFT_T |
---|
748 | |
---|
749 | if(cmpstr(direction,"X")==0) |
---|
750 | FillXCylinder($matStr,grid,rad,xc,yc,zc,len,fill) |
---|
751 | endif |
---|
752 | if(cmpstr(direction,"Y")==0) |
---|
753 | FillYCylinder($matStr,grid,rad,xc,yc,zc,len,fill) |
---|
754 | endif |
---|
755 | if(cmpstr(direction,"Z")==0) |
---|
756 | FillZCylinder($matStr,grid,rad,xc,yc,zc,len,fill) |
---|
757 | endif |
---|
758 | |
---|
759 | End |
---|
760 | |
---|
761 | |
---|
762 | Function DoTheFFT_ButtonProc(ctrlName) : ButtonControl |
---|
763 | String ctrlName |
---|
764 | |
---|
765 | Calc_IQ_FFT() |
---|
766 | // Execute "DoFFT()" |
---|
767 | End |
---|
768 | |
---|
769 | Function FFT_PlotResultsButtonProc(ctrlName) : ButtonControl |
---|
770 | String ctrlName |
---|
771 | |
---|
772 | Variable first=0 |
---|
773 | DoWindow/F FFT_IQ |
---|
774 | if(V_flag==0) |
---|
775 | first = 1 |
---|
776 | Display /W=(295,44,627,302) |
---|
777 | DoWindow/C FFT_IQ |
---|
778 | Endif |
---|
779 | |
---|
780 | // append the desired data, if it's not already there |
---|
781 | // FFTButton_4 = FFT = iBin |
---|
782 | // FFTButton_7a = binned = _XOP |
---|
783 | // FFTButton_8a = Debye = _full |
---|
784 | // FFTButton_14a = SLD = _SLD |
---|
785 | // 17 = iso USANS |
---|
786 | // 18 = Anisotropic USANS |
---|
787 | // |
---|
788 | strswitch(ctrlName) |
---|
789 | case "FFTButton_4": |
---|
790 | if(!isTraceOnGraph("iBin","FFT_IQ") && exists("iBin")==1) //only append if it's not already there |
---|
791 | AppendToGraph /W=FFT_IQ iBin vs qBin |
---|
792 | ModifyGraph mode=4,marker=19,msize=2,rgb(iBin)=(65535,0,0) |
---|
793 | endif |
---|
794 | break |
---|
795 | case "FFTButton_7a": |
---|
796 | if(!isTraceOnGraph("ival_XOP","FFT_IQ") && exists("ival_XOP")==1) //only append if it's not already there |
---|
797 | AppendToGraph /W=FFT_IQ ival_XOP vs qval_XOP |
---|
798 | ModifyGraph mode=4,marker=19,msize=2,rgb(ival_XOP)=(1,12815,52428) |
---|
799 | endif |
---|
800 | break |
---|
801 | case "FFTButton_8a": |
---|
802 | if(!isTraceOnGraph("ival_full","FFT_IQ") && exists("ival_full")==1) //only append if it's not already there |
---|
803 | AppendToGraph /W=FFT_IQ ival_full vs qval_full |
---|
804 | ModifyGraph mode=4,marker=19,msize=2,rgb(ival_full)=(0,0,0) |
---|
805 | endif |
---|
806 | break |
---|
807 | case "FFTButton_14a": |
---|
808 | if(!isTraceOnGraph("ival_SLD","FFT_IQ") && exists("ival_SLD")==1) //only append if it's not already there |
---|
809 | AppendToGraph /W=FFT_IQ ival_SLD vs qval_SLD |
---|
810 | ModifyGraph mode=4,marker=19,msize=2,rgb(ival_SLD)=(2,39321,1) |
---|
811 | endif |
---|
812 | break |
---|
813 | case "FFTButton_14a": |
---|
814 | if(!isTraceOnGraph("ival_SLD","FFT_IQ") && exists("ival_SLD")==1) //only append if it's not already there |
---|
815 | AppendToGraph /W=FFT_IQ ival_SLD vs qval_SLD |
---|
816 | ModifyGraph mode=4,marker=19,msize=2,rgb(ival_SLD)=(2,39321,1) |
---|
817 | endif |
---|
818 | break |
---|
819 | case "FFTButton_17": |
---|
820 | if(!isTraceOnGraph("FFT_iUSANS_i","FFT_IQ") && exists("FFT_iUSANS_i")==1) //only append if it's not already there |
---|
821 | AppendToGraph /W=FFT_IQ FFT_iUSANS_i vs FFT_iUSANS_q |
---|
822 | ModifyGraph mode=4,marker=19,msize=2,rgb(FFT_iUSANS_i)=(39321,1,31457) |
---|
823 | endif |
---|
824 | break |
---|
825 | case "FFTButton_18": |
---|
826 | if(!isTraceOnGraph("FFT_aUSANS_i","FFT_IQ") && exists("FFT_aUSANS_i")==1) //only append if it's not already there |
---|
827 | AppendToGraph /W=FFT_IQ FFT_aUSANS_i vs FFT_aUSANS_q |
---|
828 | ModifyGraph mode=4,marker=19,msize=2,rgb(FFT_aUSANS_i)=(52428,34958,1) |
---|
829 | endif |
---|
830 | break |
---|
831 | |
---|
832 | |
---|
833 | endswitch |
---|
834 | |
---|
835 | if(first) |
---|
836 | ModifyGraph mode=4 |
---|
837 | ModifyGraph marker=19 |
---|
838 | ModifyGraph msize=2 |
---|
839 | ModifyGraph gaps=0 |
---|
840 | ModifyGraph grid=1 |
---|
841 | ModifyGraph log=1 |
---|
842 | ModifyGraph mirror=2 |
---|
843 | Legend |
---|
844 | endif |
---|
845 | |
---|
846 | return(0) |
---|
847 | End |
---|
848 | |
---|
849 | Function isTraceOnGraph(traceStr,winStr) |
---|
850 | String traceStr,winStr |
---|
851 | |
---|
852 | Variable isOn=0 |
---|
853 | String str |
---|
854 | str = TraceNameList(winStr,";",1) //only normal traces |
---|
855 | isOn = strsearch(str,traceStr,0) //is the trace there? |
---|
856 | isOn = isOn == -1 ? 0 : 1 // return 0 if not there, 1 if there |
---|
857 | |
---|
858 | return(isOn) |
---|
859 | End |
---|
860 | |
---|
861 | Function FFTEraseMatrixButtonProc(ctrlName) : ButtonControl |
---|
862 | String ctrlName |
---|
863 | |
---|
864 | Wave mat=root:mat |
---|
865 | FastOp mat=0 |
---|
866 | return(0) |
---|
867 | End |
---|
868 | |
---|
869 | Function FFTFillSolventMatrixProc(ctrlName) : ButtonControl |
---|
870 | String ctrlName |
---|
871 | |
---|
872 | Wave mat=root:mat |
---|
873 | NVAR val=root:FFT_SolventSLD |
---|
874 | FastOp mat=(val) |
---|
875 | return(0) |
---|
876 | End |
---|
877 | |
---|
878 | Function FFT_AltiSpheresButtonProc(ctrlName) : ButtonControl |
---|
879 | String ctrlName |
---|
880 | |
---|
881 | Execute "DoSpheresCalcFFTPanel()" |
---|
882 | End |
---|
883 | |
---|
884 | Function FFT_BinnedSpheresButtonProc(ctrlName) : ButtonControl |
---|
885 | String ctrlName |
---|
886 | |
---|
887 | Execute "DoBinnedSpheresCalcFFTPanel()" |
---|
888 | End |
---|
889 | |
---|
890 | Function FFT_BinnedSLDButtonProc(ctrlName) : ButtonControl |
---|
891 | String ctrlName |
---|
892 | |
---|
893 | Execute "DoBinnedSLDCalcFFTPanel()" |
---|
894 | End |
---|
895 | |
---|
896 | Function FFT_Iso2USANS(ctrlName) : ButtonControl |
---|
897 | String ctrlName |
---|
898 | |
---|
899 | Execute "Isotropic_FFT_to_USANS()" |
---|
900 | FFT_PlotResultsButtonProc(ctrlName) |
---|
901 | End |
---|
902 | |
---|
903 | Function FFT_Aniso2USANS(ctrlName) : ButtonControl |
---|
904 | String ctrlName |
---|
905 | |
---|
906 | Execute "Anisotropic_FFT_to_USANS()" |
---|
907 | FFT_PlotResultsButtonProc(ctrlName) |
---|
908 | End |
---|
909 | |
---|
910 | |
---|
911 | |
---|
912 | |
---|
913 | /////UTILITIES |
---|
914 | |
---|
915 | // inverts the values in the matrix 0<->1 |
---|
916 | Function InvertMatrixFill(mat) |
---|
917 | Wave mat |
---|
918 | |
---|
919 | mat = (mat==1) ? 0 : 1 |
---|
920 | End |
---|
921 | |
---|
922 | // replaces specified values |
---|
923 | Proc ChangeMatrixValues(old,new) |
---|
924 | Variable old,new |
---|
925 | |
---|
926 | mat = (mat==old) ? new : mat |
---|
927 | // sequence of steps to get the gizmo to update the display correctly |
---|
928 | RemoveFromGizmo/N=Gizmo_VoxelMat object=Voxelgram0 |
---|
929 | RemoveFromGizmo/N=Gizmo_VoxelMat displayItem=axes0 |
---|
930 | AppendToGizmo/N=Gizmo_VoxelMat voxelgram=root:mat,name=voxelgram0 |
---|
931 | ModifyGizmo/N=Gizmo_VoxelMat setDisplayList=-1, object=voxelgram0 |
---|
932 | ModifyGizmo/N=Gizmo_VoxelMat setDisplayList=-1, object=axes0 //so that the axes are drawn last |
---|
933 | ModifyGizmo ModifyObject=voxelgram0 property={ pointSize,3} |
---|
934 | |
---|
935 | |
---|
936 | ColorizeGizmo() |
---|
937 | End |
---|
938 | |
---|
939 | |
---|
940 | |
---|
941 | // replaces the solvent value and updates the global |
---|
942 | Proc ReplaceSolvent(newSolv) |
---|
943 | Variable newSolv |
---|
944 | |
---|
945 | Variable solv = root:FFT_SolventSLD |
---|
946 | |
---|
947 | mat = (mat==solv) ? newSolv : mat |
---|
948 | |
---|
949 | root:FFT_SolventSLD = newSolv |
---|
950 | |
---|
951 | // sequence of steps to get the gizmo to update the display correctly |
---|
952 | RemoveFromGizmo/N=Gizmo_VoxelMat object=Voxelgram0 |
---|
953 | RemoveFromGizmo/N=Gizmo_VoxelMat displayItem=axes0 |
---|
954 | AppendToGizmo/N=Gizmo_VoxelMat voxelgram=root:mat,name=voxelgram0 |
---|
955 | ModifyGizmo/N=Gizmo_VoxelMat setDisplayList=-1, object=voxelgram0 |
---|
956 | ModifyGizmo/N=Gizmo_VoxelMat setDisplayList=-1, object=axes0 //so that the axes are drawn last |
---|
957 | ModifyGizmo ModifyObject=voxelgram0 property={ pointSize,3} |
---|
958 | |
---|
959 | |
---|
960 | ColorizeGizmo() |
---|
961 | End |
---|
962 | |
---|
963 | |
---|
964 | //overwrites any existing matrix |
---|
965 | // matrix is byte, to save space |
---|
966 | //forces the name to be "mat" |
---|
967 | // |
---|
968 | // switched to signed byte |
---|
969 | // |
---|
970 | Function MakeMatrix(nam,xd,yd,zd) |
---|
971 | String nam |
---|
972 | Variable xd,yd,zd |
---|
973 | |
---|
974 | nam="mat" |
---|
975 | Print "Matrix has been created and named \"mat\"" |
---|
976 | // Make/O/U/B/N=(xd,yd,zd) $nam |
---|
977 | Make/O/B/N=(xd,yd,zd) $nam |
---|
978 | End |
---|
979 | |
---|
980 | |
---|
981 | // calculate the average center-to-center distance between points |
---|
982 | // assumes evenly spaced on a cubic grid |
---|
983 | // |
---|
984 | Proc Center_to_Center(np) |
---|
985 | Variable np = 100 |
---|
986 | |
---|
987 | Variable Nedge = root:FFT_N |
---|
988 | Variable Tscale = root:FFT_T |
---|
989 | |
---|
990 | Variable davg |
---|
991 | |
---|
992 | davg = (Nedge*Tscale)^3 / np |
---|
993 | davg = davg^(1/3) |
---|
994 | Print "Average separation (A) = ",davg |
---|
995 | |
---|
996 | End |
---|
997 | |
---|
998 | // calculate the number of points required on a grid |
---|
999 | // to yield a given average center-to-center distance between points |
---|
1000 | // assumes evenly spaced on a cubic grid |
---|
1001 | // |
---|
1002 | // davg and Tscale are the same units, Angstroms |
---|
1003 | // |
---|
1004 | Proc Davg_to_Np(davg) |
---|
1005 | Variable davg = 400 |
---|
1006 | |
---|
1007 | Variable Nedge = root:FFT_N |
---|
1008 | Variable Tscale = root:FFT_T |
---|
1009 | |
---|
1010 | Variable np |
---|
1011 | |
---|
1012 | np = (Nedge*Tscale)^3 / davg^3 |
---|
1013 | Print "Number of points required = ",np |
---|
1014 | |
---|
1015 | End |
---|
1016 | |
---|
1017 | |
---|
1018 | |
---|
1019 | // The matrix is not necessarily 0|1, this reports the number of filled voxels |
---|
1020 | // - needed to estimate the time required for the AltiVec_Spheres calculation |
---|
1021 | Proc NumberOfPoints() |
---|
1022 | |
---|
1023 | Print "Number of points = ",NumberOfPoints_Occ(root:mat) |
---|
1024 | Print "Fraction occupied = ",VolumeFraction_Occ(root:mat) |
---|
1025 | Print "Overall Cube Edge [A] = ",root:FFT_T * root:FFT_N |
---|
1026 | Print "Found values in matrix = ",ListOfValues(root:mat) |
---|
1027 | |
---|
1028 | End |
---|
1029 | |
---|
1030 | Function NumberOfPoints_Occ(m) |
---|
1031 | Wave m |
---|
1032 | |
---|
1033 | Variable num = NonZeroValues(m) |
---|
1034 | |
---|
1035 | return(num) |
---|
1036 | End |
---|
1037 | |
---|
1038 | |
---|
1039 | Function VolumeFraction_Occ(m) |
---|
1040 | Wave m |
---|
1041 | |
---|
1042 | Variable num = NonZeroValues(m) |
---|
1043 | Variable dim = DimSize(m,0) |
---|
1044 | |
---|
1045 | return(num/dim^3) |
---|
1046 | End |
---|
1047 | |
---|
1048 | //it's a byte wave, so I can't use the trick of setting the zero values to NaN |
---|
1049 | // since NaN can't be expressed as byte. So make it binary 0|1 |
---|
1050 | // |
---|
1051 | // - the assumption here is that the defined solvent value is not part of the |
---|
1052 | // particle volume. |
---|
1053 | // |
---|
1054 | Function NonZeroValues(m) |
---|
1055 | Wave m |
---|
1056 | |
---|
1057 | Variable num |
---|
1058 | NVAR val = root:FFT_SolventSLD |
---|
1059 | |
---|
1060 | // Variable t1=ticks,dim |
---|
1061 | // dim = DimSize(m, 0 ) //assume NxNxN |
---|
1062 | Duplicate/O m,mz |
---|
1063 | |
---|
1064 | MultiThread mz = (m[p][q][r] != val) ? 1 : 0 |
---|
1065 | |
---|
1066 | WaveStats/Q/M=1 mz // NaN and Inf are not reported in V_npnts |
---|
1067 | |
---|
1068 | num = V_npnts*V_avg |
---|
1069 | |
---|
1070 | KillWaves/Z mz |
---|
1071 | // Print "Number of points = ",num |
---|
1072 | // Print "Fraction occupied = ",num/V_npnts |
---|
1073 | |
---|
1074 | // Print "time = ",(ticks-t1)/60.15 |
---|
1075 | return(num) |
---|
1076 | End |
---|
1077 | |
---|
1078 | // |
---|
1079 | // return a list of the different values of the voxels in the matrix |
---|
1080 | // |
---|
1081 | Function/S ListOfValues(m) |
---|
1082 | Wave m |
---|
1083 | |
---|
1084 | String list="" |
---|
1085 | Variable done |
---|
1086 | |
---|
1087 | Duplicate/O m,mz |
---|
1088 | |
---|
1089 | done=0 |
---|
1090 | do |
---|
1091 | WaveStats/Q/M=1 mz // NaN and Inf are not reported in V_npnts |
---|
1092 | if(V_max == V_min) |
---|
1093 | list += num2str(V_min) + ";" |
---|
1094 | done = 1 |
---|
1095 | else |
---|
1096 | list += num2str(V_max) + ";" |
---|
1097 | MultiThread mz = mz[p][q] == V_max ? V_min : mz[p][q] // replace the max with min |
---|
1098 | endif |
---|
1099 | while(!done) |
---|
1100 | |
---|
1101 | // Print "Found values in matrix = ",list |
---|
1102 | KillWaves/Z mz |
---|
1103 | |
---|
1104 | return(list) |
---|
1105 | End |
---|
1106 | |
---|
1107 | |
---|
1108 | |
---|
1109 | // returns estimate in seconds |
---|
1110 | // |
---|
1111 | // updated for quad-core iMac, 2010 |
---|
1112 | // |
---|
1113 | // type 3 = binned distances and SLD |
---|
1114 | // type 2 = binned distances |
---|
1115 | // type 1 is rather meaningless |
---|
1116 | // type 0 = is the Debye, double sum (multithreaded) |
---|
1117 | // |
---|
1118 | // |
---|
1119 | // - types 2 and 3, the binned methods are inherently AAO, so there is no significant |
---|
1120 | // dependence on the number of q-values (unless it's >> 1000). So values reported are |
---|
1121 | // for 100 q-values. |
---|
1122 | // |
---|
1123 | // There is a function Timing_Method(type) in FFT_ConcentratedSpheres.ipf to automate some of this |
---|
1124 | // |
---|
1125 | Function EstimatedTime(nx,nq,type) |
---|
1126 | Variable nx,nq,type |
---|
1127 | |
---|
1128 | Variable est=0 |
---|
1129 | |
---|
1130 | if(type==0) // "full" XOP |
---|
1131 | est = 0.4 + 1.622e-5*nx+4.86e-7*nx^2 |
---|
1132 | est /= 100 |
---|
1133 | est *= nq |
---|
1134 | endif |
---|
1135 | |
---|
1136 | if(type==1) //Igor function (much slower, 9x) |
---|
1137 | est = (nx^2)*0.0000517 //empirical values (seconds) for igor code, 100 q-values |
---|
1138 | est /= 100 // per q now... |
---|
1139 | est *= nq |
---|
1140 | endif |
---|
1141 | |
---|
1142 | if(type==2) //binned distances, slow parts XOPed |
---|
1143 | est = 0.0680 + 2.94e-6*nx+4.51e-9*nx^2 //with threading |
---|
1144 | |
---|
1145 | // est /= 100 // per q now... |
---|
1146 | // est *= nq |
---|
1147 | endif |
---|
1148 | |
---|
1149 | if(type==3) //binned distances AND sld, slow parts XOPed |
---|
1150 | est = 0.576 + 4.22e-7*nx+1.76e-8*nx^2 //with threading |
---|
1151 | |
---|
1152 | // est /= 100 // per q now... |
---|
1153 | // est *= nq |
---|
1154 | endif |
---|
1155 | |
---|
1156 | return(est) |
---|
1157 | End |
---|
1158 | |
---|
1159 | ////////////// my functions to rotate the matrix in a XYZ coordinate system |
---|
1160 | // |
---|
1161 | // |
---|
1162 | // definitely not generic, is expecting NxNxN volume |
---|
1163 | // |
---|
1164 | // not very friendly in that it "clips" anything that rotates out of the volume |
---|
1165 | // |
---|
1166 | // does translate the center of the box to 000, rotates, then translates back |
---|
1167 | // |
---|
1168 | // friendly in the sense that the rotated matrix is the same size as the original. |
---|
1169 | // --this is important for my final application (FFT) |
---|
1170 | // |
---|
1171 | // does no interpolation of values, so be sure to keep a copy of the original |
---|
1172 | // -- multiple rotation steps are going to make a mess of things. |
---|
1173 | // |
---|
1174 | // The multi axis rotation is done as one step, and probably violates every conventional |
---|
1175 | // coordinate system. The rotation is applied as RxRyRz, but this could easily be changed |
---|
1176 | // |
---|
1177 | // I just want it to be correct, so speed was not an issue. |
---|
1178 | // -- it's nested for loops. |
---|
1179 | // -- it's working with the full matrix, even when 99% is empty. |
---|
1180 | // |
---|
1181 | // 20 NOV 2013 SRK |
---|
1182 | // |
---|
1183 | // |
---|
1184 | |
---|
1185 | // |
---|
1186 | // mat is the input volume |
---|
1187 | // rotVol is the output rotated volume |
---|
1188 | Function XYZRotate(angleX,angleY,angleZ) |
---|
1189 | Variable angleX,angleY,angleZ |
---|
1190 | |
---|
1191 | NVAR FFT_N=root:FFT_N |
---|
1192 | WAVE mat=root:mat |
---|
1193 | Variable dist=FFT_N/2 |
---|
1194 | |
---|
1195 | |
---|
1196 | // convert the NxNxN into 3xN xyz locations + wave of "w" values named "values" |
---|
1197 | fVolumeToXYZTriplet(mat,"trip") |
---|
1198 | Wave trip=root:trip |
---|
1199 | |
---|
1200 | // translate to get the center of the xyz values to 0,0,0 |
---|
1201 | fTranslateCoordinate(trip,dist) //subtracts dist |
---|
1202 | |
---|
1203 | // do the rotation as a matrix multiplication |
---|
1204 | // putting zero is no rotation around that axis |
---|
1205 | // the triplet wave "trip" is overwritten with the output |
---|
1206 | DoRotation(trip,angleX,angleY,angleZ) |
---|
1207 | // Wave rotated=root:rotated |
---|
1208 | |
---|
1209 | // translate back to a 0->N based coordinate |
---|
1210 | // fTranslateCoordinate(rotated,-dist) |
---|
1211 | fTranslateCoordinate(trip,-dist) |
---|
1212 | Wave values=root:values |
---|
1213 | |
---|
1214 | // convert the triplet back to a volume |
---|
1215 | // this CLIPS anything that has rotated out of the NxNxN volume |
---|
1216 | // fXYZTripletToVolume(rotated,values,"rotVol",FFT_N) |
---|
1217 | fXYZTripletToVolume(trip,values,"rotVol",FFT_N) |
---|
1218 | |
---|
1219 | // clean up by killng the extra waves that were generated |
---|
1220 | // |
---|
1221 | KillWaves/Z trip,rotated,values |
---|
1222 | |
---|
1223 | Wave rotVol=root:rotVol |
---|
1224 | mat=rotVol |
---|
1225 | |
---|
1226 | return(0) |
---|
1227 | End |
---|
1228 | |
---|
1229 | |
---|
1230 | |
---|
1231 | Function fVolumeToXYZTriplet(matrixWave, outputName) |
---|
1232 | Wave matrixWave |
---|
1233 | String outputName |
---|
1234 | |
---|
1235 | Variable dimx=DimSize(matrixWave,0) |
---|
1236 | Variable dimy=DimSize(matrixWave,1) |
---|
1237 | Variable dimz=DimSize(matrixWave,2) |
---|
1238 | Variable rows=dimx*dimy*dimz |
---|
1239 | Make/O/N=(3,rows) $outputName |
---|
1240 | Make/O/N=(rows) values |
---|
1241 | WAVE TripletWave= $outputName |
---|
1242 | Wave values=values |
---|
1243 | |
---|
1244 | |
---|
1245 | Variable ii,jj,kk,count=0 |
---|
1246 | Variable xVal,yVal,zval |
---|
1247 | for(kk=0;kk<dimz;kk+=1) // kk is z (layer) |
---|
1248 | zval=kk |
---|
1249 | for(jj=0;jj<dimy;jj+=1) // jj is y (column) |
---|
1250 | yVal=jj |
---|
1251 | for(ii=0;ii<dimx;ii+=1) // ii is x (row) |
---|
1252 | xVal=ii |
---|
1253 | TripletWave[0][count]=xVal |
---|
1254 | TripletWave[1][count]=yVal |
---|
1255 | TripletWave[2][count]=zval |
---|
1256 | values[count]=matrixWave[ii][jj][kk] // value at [row][col][lay] |
---|
1257 | count+=1 |
---|
1258 | endfor |
---|
1259 | endfor |
---|
1260 | endfor |
---|
1261 | |
---|
1262 | return(0) |
---|
1263 | End |
---|
1264 | |
---|
1265 | |
---|
1266 | Function fXYZTripletToVolume(triplet, values, outputName, outputDim) |
---|
1267 | Wave triplet,values |
---|
1268 | String outputName |
---|
1269 | Variable outputDim |
---|
1270 | |
---|
1271 | Variable numPt=DimSize(triplet,1) |
---|
1272 | |
---|
1273 | Variable num = outputDim |
---|
1274 | |
---|
1275 | Make/O/B/N=(num,num,num) $outputName |
---|
1276 | WAVE newVol= $outputName |
---|
1277 | |
---|
1278 | FastOp newVol = 0 |
---|
1279 | |
---|
1280 | Variable ii,jj,kk,count=0 |
---|
1281 | Variable xVal,yVal,zval |
---|
1282 | Variable xOK, yOK, zOK |
---|
1283 | |
---|
1284 | |
---|
1285 | for(ii=0;ii<numPt;ii+=1) |
---|
1286 | xval = round(triplet[0][ii]) |
---|
1287 | yval = round(triplet[1][ii]) |
---|
1288 | zval = round(triplet[2][ii]) |
---|
1289 | |
---|
1290 | // round and keep in bounds (returns truth) |
---|
1291 | xOK = inRange(xval,0,num-1) |
---|
1292 | yOK = inRange(yval,0,num-1) |
---|
1293 | zOK = inRange(zval,0,num-1) |
---|
1294 | |
---|
1295 | if(xOK && yOK && zOK) |
---|
1296 | newVol[xval][yval][zval] = values[ii] |
---|
1297 | endif |
---|
1298 | |
---|
1299 | endfor |
---|
1300 | |
---|
1301 | |
---|
1302 | return(0) |
---|
1303 | End |
---|
1304 | |
---|
1305 | // if val < lo or > hi, bad val |
---|
1306 | // fi both of these pass, pt is OK |
---|
1307 | ThreadSafe Static Function inRange(val, lo, hi) |
---|
1308 | Variable val, lo, hi |
---|
1309 | |
---|
1310 | if(val < lo) |
---|
1311 | return (0) |
---|
1312 | endif |
---|
1313 | if(val > hi) |
---|
1314 | return (0) |
---|
1315 | endif |
---|
1316 | |
---|
1317 | return(1) |
---|
1318 | |
---|
1319 | End |
---|
1320 | |
---|
1321 | // Rotation is applied in the order Rx Ry Rz |
---|
1322 | Function DoRotation(triplet,angleX,angleY,angleZ) |
---|
1323 | Wave triplet |
---|
1324 | Variable angleX,angleY,angleZ |
---|
1325 | |
---|
1326 | Variable thetaX,thetaY,thetaZ |
---|
1327 | thetaX = angleX*pi/180 // convert degrees to radians |
---|
1328 | thetaY = angleY*pi/180 // convert degrees to radians |
---|
1329 | thetaZ = angleZ*pi/180 // convert degrees to radians |
---|
1330 | |
---|
1331 | Make/O/D/N=(3,3) Rx,Ry,Rz |
---|
1332 | Rx=0 |
---|
1333 | Ry=0 |
---|
1334 | Rz=0 |
---|
1335 | |
---|
1336 | Rx[0][0] = 1 |
---|
1337 | Rx[1][1] = cos(thetaX) |
---|
1338 | Rx[1][2] = -sin(thetaX) |
---|
1339 | Rx[2][1] = sin(thetaX) |
---|
1340 | Rx[2][2] = cos(thetaX) |
---|
1341 | |
---|
1342 | Ry[0][0] = cos(thetaY) |
---|
1343 | Ry[0][2] = sin(thetaY) |
---|
1344 | Ry[1][1] = 1 |
---|
1345 | Ry[2][0] = -sin(thetaY) |
---|
1346 | Ry[2][2] = cos(thetaY) |
---|
1347 | |
---|
1348 | Rz[0][0] = cos(thetaZ) |
---|
1349 | Rz[0][1] = -sin(thetaZ) |
---|
1350 | Rz[1][0] = sin(thetaZ) |
---|
1351 | Rz[1][1] = cos(thetaZ) |
---|
1352 | Rz[2][2] = 1 |
---|
1353 | |
---|
1354 | |
---|
1355 | // MatrixOp/O rotated = Rx x Ry x Rz x triplet |
---|
1356 | MatrixOp/O triplet = Rx x Ry x Rz x triplet |
---|
1357 | |
---|
1358 | |
---|
1359 | return(0) |
---|
1360 | end |
---|
1361 | |
---|
1362 | |
---|
1363 | // the rotation matrix, as I copied from wikipedia, is a rotation around (0,0,0), not |
---|
1364 | // the center of the gizmo plot |
---|
1365 | Function fTranslateCoordinate(trip,dist) |
---|
1366 | Wave trip |
---|
1367 | Variable dist |
---|
1368 | |
---|
1369 | // MatrixOp/O trip = trip - dist |
---|
1370 | trip = trip - dist |
---|
1371 | |
---|
1372 | return(0) |
---|
1373 | End |
---|
1374 | |
---|