source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Panel.ipf @ 849

Last change on this file since 849 was 849, checked in by srkline, 11 years ago

Changes to the real-space modeling to sped up the drawing of cylinders, and to provide a few more examples of scritped, unattended calculations, and saving the 3D matrix.
File size: 22.6 KB
Line 
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
43Proc 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_estTime = 0
54               
55                Variable/G root:FFT_SolventSLD = 0
56                Variable/G root:FFT_delRho = 1e-7                       //multiplier for SLD (other value is 1e-7)
57               
58                FFT_Qmax :=2*pi/FFT_T
59                FFT_QmaxReal := FFT_Qmax/2
60                FFT_DQ := pi/(FFT_N*FFT_T)
61                //empirical fit (cubic) of time vs N=50 to N=256
62                FFT_estTime := 0.56 - 0.0156*FFT_N + 0.000116*FFT_N^2 + 8e-7*FFT_N^3
63//              FFT_estTime := FFT_N/128
64               
65                FFT_Panel()
66        endif
67End
68
69
70Proc FFT_Panel()
71        PauseUpdate; Silent 1           // building window...
72        NewPanel /W=(1452,44,1768,531)/K=1 as "FFT_Panel"
73        DoWindow/C FFT_Panel
74        SetDrawLayer UserBack
75        DrawLine 5,68,311,68
76        DrawLine 5,142,311,142
77        DrawLine 5,250,311,250
78        SetVariable FFTSetVar_0,pos={7,7},size={150,15},title="Cells per edge (N)"
79        SetVariable FFTSetVar_0,limits={50,512,2},value= FFT_N,live= 1
80        SetVariable FFTSetVar_1,pos={7,33},size={150,15},title="Length per Cell (T)"
81        SetVariable FFTSetVar_1,limits={1,500,0.2},value= FFT_T,live= 1
82        SetVariable FFTSetVar_2,pos={183,26},size={120,15},title="Real Qmax"
83        SetVariable FFTSetVar_2,limits={0,0,0},value= FFT_QmaxReal,noedit= 1,live= 1
84        SetVariable FFTSetVar_3,pos={183,48},size={120,15},title="delta Q (A)"
85        SetVariable FFTSetVar_3,limits={0,0,0},value= FFT_DQ,noedit= 1,live= 1
86        Button FFTButton_0,pos={15,79},size={90,20},proc=FFT_MakeMatrixButtonProc,title="Make Matrix"
87        Button FFTButton_1,pos={14,157},size={90,20},proc=FFTMakeGizmoButtonProc,title="Make Gizmo"
88        Button FFTButton_2,pos={14,187},size={100,20},proc=FFTDrawSphereButtonProc,title="Draw Sphere"
89        Button FFTButton_3,pos={14,265},size={70,20},proc=DoTheFFT_ButtonProc,title="Do FFT"
90        Button FFTButton_4,pos={180,264},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot FFT Results"
91        Button FFTButton_5,pos={13,218},size={120,20},proc=FFTDrawZCylinderButtonProc,title="Draw Cylinder"
92//      Button FFTButton_6,pos={134,79},size={90,20},proc=FFTEraseMatrixButtonProc,title="Erase Matrix"
93        Button FFTButton_6a,pos={160,79},size={60,20},proc=FFTSaveMatrixButtonProc,title="Save"
94        Button FFTButton_6b,pos={240,79},size={60,20},proc=FFTLoadMatrixButtonProc,title="Load"
95        Button FFTButton_7,pos={13,329},size={130,20},proc=FFT_BinnedSpheresButtonProc,title="Do Binned Debye"
96        Button FFTButton_7a,pos={180,329},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot Binned Results"
97
98        Button FFTButton_8,pos={13,297},size={130,20},proc=FFT_AltiSpheresButtonProc,title="Do Debye Spheres"
99        Button FFTButton_8a,pos={180,297},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot Debye Results"
100
101        Button FFTButton_14,pos={13,360},size={130,20},proc=FFT_BinnedSLDButtonProc,title="Do Binned SLD"
102        Button FFTButton_14a,pos={180,360},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot SLD Results"
103
104        SetVariable FFTSetVar_4,pos={201,4},size={100,15},title="FFT time(s)"
105        SetVariable FFTSetVar_4,limits={0,0,0},value= FFT_estTime,noedit= 1,live= 1,format="%d"
106        Button FFTButton_9,pos={200,400},size={100,20},proc=FFT_Get2DSlice,title="Get 2D Slice"
107        Button FFTButton_10,pos={169,156},size={130,20},proc=FFT_TransposeMat,title="Transpose Matrix"
108        Button FFTButton_11,pos={169,189},size={130,20},proc=FFT_RotateMat,title="Rotate Matrix"
109        Button FFTButton_12,pos={168,219},size={130,20},proc=FFT_AddRotatedObject,title="Add Rotated Obj"
110        Button FFTButton_12,disable=2           // hide this button
111        Button FFTButton_13,pos={14,109},size={120,20},proc=FFTFillSolventMatrixProc,title="Solvent Matrix"
112        SetVariable FFTSetVar_5,pos={155,111},size={150,15},title="Solvent SLD (10^-7)"
113        SetVariable FFTSetVar_5,limits={-99,99,1},value= FFT_SolventSLD,live= 1
114        Button FFTButton_15,pos={209,430},size={90,20},proc=Interp2DSliceButton,title="Interp 2D"
115        Button FFTButton_16,pos={14,460},size={70,20},proc=FFTHelpButton,title="Help"
116EndMacro
117
118// Save a matrix wave, plus the N, T, and solvent values in the wave note for reloading
119Function FFTSaveMatrixButtonProc(ba) : ButtonControl
120        STRUCT WMButtonAction &ba
121       
122        String win = ba.win
123
124        switch (ba.eventCode)
125                case 2:
126                        // click code here
127                        String fileStr=""
128                        SaveMyMatrix(fileStr)
129                       
130                        break
131        endswitch
132
133        return 0
134End
135
136// this will wave as Igor Binary, so be sure to use the ".ibw extension.
137// - this could possibly be enforced, but that's maybe not necessary at this stage.
138//
139Function SaveMyMatrix(fileStr)
140        String fileStr
141       
142        WAVE mat=root:mat
143        NVAR FFT_T = root:FFT_T
144        NVAR FFT_N = root:FFT_N
145        NVAR FFT_SolventSLD = root:FFT_SolventSLD
146        String str=""
147        sprintf str,"FFT_T=%g;FFT_N=%d;FFT_SolventSLD=%d;",FFT_T,FFT_N,FFT_SolventSLD
148        Note mat,str
149        Save/C/P=home mat as fileStr    //will ask for a file name, save as Igor Binary
150        Note/K mat                      //kill wave note on exiting since I don't properly update this anywhere else
151                       
152        return(0)
153end
154
155
156// load in a previously saved matrix, and reset FFT_N, FFT_T and solvent
157// from the wave note when saved
158Function FFTLoadMatrixButtonProc(ba) : ButtonControl
159        STRUCT WMButtonAction &ba
160       
161        String win = ba.win
162
163        switch (ba.eventCode)
164                case 2:
165                        // click code here
166                        String fileStr=""
167                        ReloadMatrix(fileStr)
168                       
169                        break
170        endswitch
171
172        return 0
173End
174
175
176Function ReloadMatrix(fileStr)
177        String fileStr
178       
179                LoadWave/M/O/W/P=home           fileStr         //will ask for a file, Igor Binary format is assumed here
180                if(V_flag == 0)
181                        return(0)               //user cancel
182                endif
183               
184                String str
185                str=note(mat)
186                NVAR FFT_T = root:FFT_T
187                NVAR FFT_N = root:FFT_N
188                NVAR FFT_SolventSLD = root:FFT_SolventSLD
189               
190                FFT_T = NumberByKey("FFT_T", str, "=" ,";")
191                FFT_N = NumberByKey("FFT_N", str, "=" ,";")
192                FFT_SolventSLD = NumberByKey("FFT_SolventSLD", str, "=" ,";")
193
194// if I got bad values, put in default values                   
195                if(numtype(FFT_T) != 0 )
196                        FFT_T = 5
197                endif
198                if(numtype(FFT_N) != 0 )
199                        FFT_N = DimSize(mat,0)
200                endif
201                if(numtype(FFT_SolventSLD) != 0 )
202                        FFT_SolventSLD = 0
203                endif                   
204               
205                Print "Loaded matrix parameters = ",str
206                Execute "NumberOfPoints()"
207                       
208        return(0)
209end
210
211// utility function that reads the wave note from a binary file
212//
213Function/S LoadNoteFunc(PName,FName[,FileRef])
214        String PName, FName
215        Variable FileRef
216       
217        Variable noteStart, noteLength, version, dependLength
218        String noteStr, typeStr = ".ibw"
219        if (ParamIsDefault(FileRef))
220                Open/R/P=$PName/T=typeStr fileRef, as FName     //open the file
221        endif
222        FSetPos fileRef, 0
223        FBinRead/F=2 fileRef, version
224       
225       
226        if (version == 5)
227       
228                FSetPos fileRef, 4
229                Make/N=(3)/I/Free SizeWave
230                FBinRead FileRef,SizeWave
231                noteStart = SizeWave[0]
232                DependLength = SizeWave[1]
233                NoteLength = SizeWave[2]
234                noteStart += dependLength+64
235               
236        elseif (version == 2)
237       
238                FBinRead/F=3 fileRef, noteStart
239//              FBinRead/F=4 fileRef, dependLength
240                FBinRead/F=3 fileRef, noteLength
241                noteStart += 16
242               
243        else
244       
245                if (ParamIsDefault(FileRef))
246                        Close(FileRef)          //close the file
247                endif
248                return ""
249       
250        endif
251        if (!NoteLength)
252                if (ParamIsDefault(FileRef))
253                        Close(FileRef)          //close the file
254                endif
255                return("")
256        endif
257        FSetPos fileRef, noteStart
258        NoteStr = PadString("",NoteLength,0)
259        FBinRead FileRef,NoteStr
260       
261        if (ParamIsDefault(FileRef))
262                Close(FileRef)          //close the file
263        endif
264        return noteStr
265       
266End //LoadNoteFunc
267
268
269Function FFTHelpButton(ba) : ButtonControl
270        STRUCT WMButtonAction &ba
271       
272        String win = ba.win
273
274        switch (ba.eventCode)
275                case 2:
276                        // click code here
277                        DisplayHelpTopic/Z/K=1 "Real-Space Modeling of SANS Data"
278                        if(V_flag !=0)
279                                DoAlert 0,"The Real-Space Modeling Help file could not be found"
280                        endif
281                        break
282        endswitch
283
284        return 0
285End
286
287
288Function Interp2DSliceButton(ba) : ButtonControl
289        STRUCT WMButtonAction &ba
290
291        switch( ba.eventCode )
292                case 2: // mouse up
293                        // click code here
294                        String folderStr=""
295                        Prompt folderStr, "Pick a 2D data folder"               // Set prompt for x param
296                        DoPrompt "Enter data folder", folderStr
297                        if (V_Flag || strlen(folderStr) == 0)
298                                return -1                                                               // User canceled, null string entered
299                        endif                   
300       
301                        Interpolate2DSliceToData(folderStr)
302                       
303                        //display the 2D plane
304                        DoWindow FFT_Interp2D_log
305                        if(V_flag == 0)
306                                Execute "Display2DInterpSlice_log()"
307                        endif
308                       
309                        break
310        endswitch
311
312        return 0
313End
314
315
316
317Function FFT_Get2DSlice(ctrlName) : ButtonControl
318        String ctrlName
319               
320        WAVE/Z dmat = root:dmat
321        if(WaveExists(dmat)==1)
322                get2DSlice(dmat)
323        endif
324       
325        //display the 2D plane
326        DoWindow FFT_Slice2D
327        if(V_flag == 0)
328                Execute "Display2DSlice()"
329        endif
330       
331        //display the 2D plane
332        DoWindow FFT_Slice2D_log
333        if(V_flag == 0)
334                Execute "Display2DSlice_log()"
335        endif
336       
337        //display the 1D binning of the 2D plane
338        DoWindow FFT_Slice1D
339        if(V_flag == 0)
340                Execute "Slice2_1D()"
341        endif
342       
343        return(0)       
344End
345
346Proc Display2DSlice()
347        PauseUpdate; Silent 1           // building window...
348        Display /W=(1038,44,1404,403)
349        DoWindow/C FFT_Slice2D
350        AppendImage/T detPlane
351        ModifyImage detPlane ctab= {*,*,YellowHot,0}
352        ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14
353        ModifyGraph mirror=2
354        ModifyGraph nticks=4
355        ModifyGraph minor=1
356        ModifyGraph fSize=9
357        ModifyGraph standoff=0
358        ModifyGraph tkLblRot(left)=90
359        ModifyGraph btLen=3
360        ModifyGraph tlOffset=-2
361        SetAxis/A/R left
362End
363
364Proc Display2DInterpSlice_log()
365        PauseUpdate; Silent 1           // building window...
366        Display /W=(1038,44,1404,403)
367        DoWindow/C FFT_Interp2D_log
368        AppendImage/T interp2DSlice_log
369        ModifyImage interp2DSlice_log ctab= {*,*,YellowHot,0}
370        ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14
371        ModifyGraph mirror=2
372        ModifyGraph nticks=4
373        ModifyGraph minor=1
374        ModifyGraph fSize=9
375        ModifyGraph standoff=0
376        ModifyGraph tkLblRot(left)=90
377        ModifyGraph btLen=3
378        ModifyGraph tlOffset=-2
379        SetAxis/A/R left
380End
381
382Proc Display2DSlice_log()
383        PauseUpdate; Silent 1           // building window...
384        Display /W=(1038,44,1404,403)
385        DoWindow/C FFT_Slice2D_log
386        AppendImage/T logP
387        ModifyImage logP ctab= {*,*,YellowHot,0}
388        ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14
389        ModifyGraph mirror=2
390        ModifyGraph nticks=4
391        ModifyGraph minor=1
392        ModifyGraph fSize=9
393        ModifyGraph standoff=0
394        ModifyGraph tkLblRot(left)=90
395        ModifyGraph btLen=3
396        ModifyGraph tlOffset=-2
397        SetAxis/A/R left
398End
399Proc Slice2_1D()
400        PauseUpdate; Silent 1           // building window...
401        Display /W=(1034,425,1406,763) iBin_2d vs qBin_2d
402        DoWindow/C FFT_Slice1D
403        ModifyGraph mode=4
404        ModifyGraph marker=19
405        ModifyGraph msize=2
406        ModifyGraph grid=1
407        ModifyGraph log=1
408        ModifyGraph mirror=2
409        Legend
410EndMacro
411
412
413Function FFT_TransposeMat(ctrlName) : ButtonControl
414        String ctrlName
415
416        Variable mode=1
417        Prompt mode,"Transform XYZ to:",popup,"XZY;ZXY;ZYX;YZX;YXZ;"
418        DoPrompt "Transform mode",mode
419        if (V_Flag)
420                return 0                                                                        // user canceled
421        endif
422       
423        fFFT_TransposeMat(mode)
424       
425        return (0)
426End
427
428// starting in XYZ
429// mode 1;2;3;4;5;
430//correspond to
431// "XZY;ZXY;ZYX;YZX;YXZ;"
432//
433Function fFFT_TransposeMat(mode)
434        Variable mode
435
436        WAVE/Z mat=mat
437        ImageTransform /G=(mode) transposeVol mat
438        WAVE M_VolumeTranspose=M_VolumeTranspose
439        Duplicate/O M_VolumeTranspose mat
440       
441        return(0)       
442End
443
444Function FFT_RotateMat(ctrlName) : ButtonControl
445        String ctrlName
446       
447        Variable degree=45,sense=1
448        Prompt degree, "Degrees of rotation around Z-axis:"
449//      Prompt sense, "Direction of rotation:",popup,"CW;CCW;"
450        DoPrompt "Enter paramters for rotation", degree
451       
452        if (V_Flag)
453                return 0                                                                        // user canceled
454        endif
455       
456        fFFT_RotateMat(degree)
457        return(0)
458End
459// note that the rotation is not perfect. if the rotation produces an
460// odd number of pixels, then the object will "walk" one pixel. Also, some
461// small artifacts are introduced, simply due to the voxelization of the object
462// as it is rotated. rotating a cylinder 10 then 350 shows a few extra "bumps" on
463// the surface, but the calculated scattering is virtually identical.
464//
465// these issues, may be correctable, if needed
466//
467Function fFFT_RotateMat(degree)
468        Variable degree
469       
470        Variable fill=0
471       
472        WAVE mat = mat
473        Variable dx,dy,dz,nx,ny,nz
474        dx = DimSize(mat,0)
475        dy = DimSize(mat,1)
476        dz = DimSize(mat,2)
477               
478        //? the /W and /C flags seem to be special cases for 90 deg rotations
479        ImageRotate /A=(degree)/E=(fill) mat            //mat is not overwritten, unless /O
480       
481        nx = DimSize(M_RotatedImage,0)
482        ny = DimSize(M_RotatedImage,1)
483        nz = DimSize(M_RotatedImage,2)
484//      Print "rotated dims = ",dx,dy,dz,nx,ny,nz
485        Variable delx,dely,odd=0
486        delx = nx-dx
487        dely = ny-dy
488       
489        if(mod(delx,2) != 0)            //sometimes the new x,y dims are odd!
490                odd = 1
491        endif
492//      delx = (delx + odd)/2
493//      dely = (dely + odd)/2
494        delx = trunc(delx/2)
495        dely = trunc(dely/2)
496
497        //+odd removes an extra point if there is one
498        Duplicate/O/R=[delx+odd,nx-delx-1][dely+odd,ny-dely-1][0,nz-1] M_RotatedImage mat
499
500// - not sure why the duplicate operation changes the start and delta of mat, but I
501// need to reset the start and delta to be sure that the display is correct, and that the scaling
502// is read correctly later
503        SetScale/P x 0,1,"", mat
504        SetScale/P y 0,1,"", mat
505       
506        nx = DimSize(mat,0)
507        ny = DimSize(mat,1)
508        nz = DimSize(mat,2)
509//      Print "redim = ",dx,dy,dz,nx,ny,nz
510       
511        KillWaves/Z M_RotatedImage
512       
513End
514
515// also look for ImageTransform operations that will allow translation of the objects.
516// then simple objects could be oriented @0,0,0, and translated to the correct position (and added)
517//
518Function FFT_AddRotatedObject(ctrlName) : ButtonControl
519        String ctrlName
520
521        Print "Not yet implemented"
522End
523
524
525Function FFT_MakeMatrixButtonProc(ctrlName) : ButtonControl
526        String ctrlName
527
528        NVAR nn=root:FFT_N
529        if(mod(nn, 2 ) ==1)             //force an even number for FFT
530                nn +=1
531        endif
532        MakeMatrix("mat",nn,nn,nn)
533End
534
535Function FFTMakeGizmoButtonProc(ctrlName) : ButtonControl
536        String ctrlName
537       
538        DoWindow/F Gizmo_VoxelMat
539        if(V_flag==0)
540                Execute "Gizmo_VoxelMat()"
541        endif
542End
543
544Function FFTDrawSphereButtonProc(ctrlName)  : ButtonControl
545        String ctrlName
546       
547        Execute "FFTDrawSphereProc()"
548End
549
550Proc FFTDrawSphereProc(matStr,rad,xc,yc,zc,fill)
551        String matStr="mat"
552        Variable rad=25,xc=50,yc=50,zc=50,fill=10
553        Prompt matStr,"the wave"                //,popup,WaveList("*",";","")
554        Prompt rad,"enter real radius (A)"
555        Prompt xc,"enter the X-center"
556        Prompt yc,"enter the Y-center"
557        Prompt zc,"enter the Z-center"
558        Prompt fill,"fill SLD value"
559       
560        Variable grid=root:FFT_T
561       
562        FillSphereRadius($matStr,grid,rad,xc,yc,zc,fill)
563       
564End
565
566Function FFTDrawZCylinderButtonProc(ctrlName)  : ButtonControl
567        String ctrlName
568       
569        Execute "FFTDrawCylinder()"
570End
571
572Proc FFTDrawCylinder(direction,matStr,rad,len,xc,yc,zc,fill)
573        String direction
574        String matStr="mat"
575        Variable rad=25,len=300,xc=50,yc=50,zc=50,fill=10
576        Prompt direction, "Direction", popup "X;Y;Z;"
577        Prompt matStr,"the wave"                //,popup,WaveList("*",";","")
578        Prompt rad,"enter real radius (A)"
579        Prompt len,"enter length (A)"
580        Prompt xc,"enter the X-center"
581        Prompt yc,"enter the Y-center"
582        Prompt zc,"enter the Z-center"
583        Prompt fill,"fill SLD value"
584       
585       
586        Variable grid=root:FFT_T
587       
588        if(cmpstr(direction,"X")==0)
589                FillXCylinder($matStr,grid,rad,xc,yc,zc,len,fill)
590        endif
591        if(cmpstr(direction,"Y")==0)
592                FillYCylinder($matStr,grid,rad,xc,yc,zc,len,fill)
593        endif
594        if(cmpstr(direction,"Z")==0)
595                FillZCylinder($matStr,grid,rad,xc,yc,zc,len,fill)
596        endif
597       
598End
599
600
601Function DoTheFFT_ButtonProc(ctrlName) : ButtonControl
602        String ctrlName
603
604        Execute "DoFFT()"
605End
606
607Function FFT_PlotResultsButtonProc(ctrlName) : ButtonControl
608        String ctrlName
609
610        Variable first=0
611        DoWindow/F FFT_IQ
612        if(V_flag==0)
613                first = 1
614                Display /W=(295,44,627,302)
615                DoWindow/C FFT_IQ
616        Endif
617       
618        // append the desired data, if it's not already there
619        // FFTButton_4 = FFT            = iBin
620        // FFTButton_7a = binned = _XOP
621        // FFTButton_8a = Debye = _full
622        // FFTButton_14a = SLD = _SLD
623        strswitch(ctrlName)     
624                case "FFTButton_4":
625                        if(!isTraceOnGraph("iBin","FFT_IQ") && exists("iBin")==1)               //only append if it's not already there
626                                AppendToGraph /W=FFT_IQ iBin vs qBin
627                                ModifyGraph mode=4,marker=19,msize=2,rgb(iBin)=(65535,0,0)
628                        endif
629                        break
630                case "FFTButton_7a":
631                        if(!isTraceOnGraph("ival_XOP","FFT_IQ") && exists("ival_XOP")==1)               //only append if it's not already there
632                                AppendToGraph /W=FFT_IQ ival_XOP vs qval_XOP
633                                ModifyGraph mode=4,marker=19,msize=2,rgb(ival_XOP)=(1,12815,52428)
634                        endif           
635                        break
636                case "FFTButton_8a":
637                        if(!isTraceOnGraph("ival_full","FFT_IQ") && exists("ival_full")==1)             //only append if it's not already there
638                                AppendToGraph /W=FFT_IQ ival_full vs qval_full
639                                ModifyGraph mode=4,marker=19,msize=2,rgb(ival_full)=(0,0,0)
640                        endif           
641                        break
642                case "FFTButton_14a":
643                        if(!isTraceOnGraph("ival_SLD","FFT_IQ") && exists("ival_SLD")==1)               //only append if it's not already there
644                                AppendToGraph /W=FFT_IQ ival_SLD vs qval_SLD
645                                ModifyGraph mode=4,marker=19,msize=2,rgb(ival_SLD)=(2,39321,1)
646                        endif           
647                        break
648        endswitch
649       
650        if(first)
651                ModifyGraph mode=4
652                ModifyGraph marker=19
653                ModifyGraph msize=2
654                ModifyGraph gaps=0
655                ModifyGraph grid=1
656                ModifyGraph log=1
657                ModifyGraph mirror=2
658                Legend
659        endif
660       
661        return(0)
662End
663
664Function isTraceOnGraph(traceStr,winStr)
665        String traceStr,winStr
666       
667        Variable isOn=0
668        String str
669        str = TraceNameList(winStr,";",1)               //only normal traces
670        isOn = strsearch(str,traceStr,0)                //is the trace there?
671        isOn = isOn == -1 ? 0 : 1                       // return 0 if not there, 1 if there
672
673        return(isOn)
674End
675
676Function FFTEraseMatrixButtonProc(ctrlName) : ButtonControl
677        String ctrlName
678       
679        Wave mat=root:mat
680        FastOp mat=0
681        return(0)
682End
683
684Function FFTFillSolventMatrixProc(ctrlName) : ButtonControl
685        String ctrlName
686       
687        Wave mat=root:mat
688        NVAR val=root:FFT_SolventSLD
689        mat=val
690        return(0)
691End
692
693Function FFT_AltiSpheresButtonProc(ctrlName) : ButtonControl
694        String ctrlName
695
696        Execute "DoSpheresCalcFFTPanel()"
697End
698
699Function FFT_BinnedSpheresButtonProc(ctrlName) : ButtonControl
700        String ctrlName
701
702        Execute "DoBinnedSpheresCalcFFTPanel()"
703End
704
705Function FFT_BinnedSLDButtonProc(ctrlName) : ButtonControl
706        String ctrlName
707
708        Execute "DoBinnedSLDCalcFFTPanel()"
709End
710
711
712
713
714
715/////UTILITIES
716
717// inverts the values in the matrix 0<->1
718Function InvertMatrixFill(mat)
719        Wave mat
720       
721        mat = (mat==1) ? 0 : 1
722End
723
724//overwrites any existing matrix
725// matrix is byte, to save space
726//forces the name to be "mat"
727//
728// switched to signed byte
729//
730Function MakeMatrix(nam,xd,yd,zd)
731        String nam
732        Variable xd,yd,zd
733       
734        nam="mat"
735        Print "Matrix has been created and named \"mat\""
736//      Make/O/U/B/N=(xd,yd,zd) $nam
737        Make/O/B/N=(xd,yd,zd) $nam
738End
739
740
741// calculate the average center-to-center distance between points
742// assumes evenly spaced on a cubic grid
743//
744Proc Center_to_Center(np)
745        Variable np = 100
746       
747        Variable Nedge = root:FFT_N
748        Variable Tscale = root:FFT_T
749       
750        Variable davg
751       
752        davg = (Nedge*Tscale)^3 / np
753        davg = davg^(1/3)
754        Print "Average separation (A) = ",davg
755       
756End
757
758// calculate the number of points required on a grid
759// to yield a given average center-to-center distance between points
760// assumes evenly spaced on a cubic grid
761//
762// davg and Tscale are the same units, Angstroms
763//
764Proc Davg_to_Np(davg)
765        Variable davg = 400
766       
767        Variable Nedge = root:FFT_N
768        Variable Tscale = root:FFT_T
769       
770        Variable np
771       
772        np = (Nedge*Tscale)^3 / davg^3
773        Print "Number of points required = ",np
774       
775End
776
777
778
779
780
781
782// The matrix is not necessarily 0|1, this reports the number of filled voxels
783// - needed to estimate the time required for the AltiVec_Spheres calculation
784Proc NumberOfPoints()
785       
786        Print "Number of points = ",NumberOfPoints_Occ(root:mat)
787        Print "Fraction occupied = ",VolumeFraction_Occ(root:mat)
788        Print "Overall Cube Edge [A] = ",root:FFT_T * root:FFT_N
789        Print "Found values in matrix = ",ListOfValues(root:mat)
790       
791End
792
793Function NumberOfPoints_Occ(m)
794        Wave m
795       
796        Variable num = NonZeroValues(m)
797
798        return(num)
799End
800
801
802Function VolumeFraction_Occ(m)
803        Wave m
804       
805        Variable num = NonZeroValues(m)
806        Variable dim = DimSize(m,0)
807
808        return(num/dim^3)
809End
810
811//it's a byte wave, so I can't use the trick of setting the zero values to NaN
812// since NaN can't be expressed as byte. So make it binary 0|1
813//
814// - the assumption here is that the defined solvent value is not part of the
815// particle volume.
816//
817Function NonZeroValues(m)
818        Wave m
819       
820        Variable num
821        NVAR val = root:FFT_SolventSLD
822       
823//      Variable t1=ticks,dim
824//      dim = DimSize(m, 0 )            //assume NxNxN
825        Duplicate/O m,mz
826       
827        MultiThread mz = (m[p][q] != val) ? 1 : 0
828       
829        WaveStats/Q/M=1 mz              // NaN and Inf are not reported in V_npnts
830       
831        num = V_npnts*V_avg
832       
833        KillWaves/Z mz
834//      Print "Number of points = ",num
835//      Print "Fraction occupied = ",num/V_npnts
836       
837//      Print "time = ",(ticks-t1)/60.15
838        return(num)
839End
840
841//
842// return a list of the different values of the voxels in the matrix
843//
844Function/S ListOfValues(m)
845        Wave m
846       
847        String list=""
848        Variable done
849
850        Duplicate/O m,mz
851       
852        done=0
853        do
854                WaveStats/Q/M=1 mz              // NaN and Inf are not reported in V_npnts
855                if(V_max == V_min)
856                        list += num2str(V_min) + ";"
857                        done = 1
858                else
859                        list += num2str(V_max) + ";"
860                        MultiThread mz = mz[p][q] == V_max ? V_min : mz[p][q]           // replace the max with min                     
861                endif
862        while(!done)   
863       
864//      Print "Found values in matrix = ",list
865        KillWaves/Z mz
866
867        return(list)
868End
869
870
871
872// returns estimate in seconds
873//
874//              updated for quad-core iMac, 2010
875//
876// type 3 = binned distances and SLD
877// type 2 = binned distances
878// type 1 is rather meaningless
879// type 0 = is the Debye, double sum (multithreaded)
880//
881//
882// - types 2 and 3, the binned methods are inherently AAO, so there is no significant
883// dependence on the number of q-values (unless it's  >> 1000). So values reported are
884// for 100 q-values.
885//
886// There is a function Timing_Method(type) in FFT_ConcentratedSpheres.ipf to automate some of this
887//
888Function EstimatedTime(nx,nq,type)
889        Variable nx,nq,type
890       
891        Variable est=0
892       
893        if(type==0)             // "full" XOP
894                est = 0.4 + 1.622e-5*nx+4.86e-7*nx^2
895                est /= 100
896                est *= nq
897        endif
898       
899        if(type==1)             //Igor function (much slower, 9x)
900                est = (nx^2)*0.0000517          //empirical values (seconds) for igor code, 100 q-values
901                est /= 100                                      // per q now...
902                est *= nq
903        endif
904       
905        if(type==2)             //binned distances, slow parts XOPed
906                est = 0.0680 + 2.94e-6*nx+4.51e-9*nx^2                  //with threading
907               
908//              est /= 100                                      // per q now...
909//              est *= nq
910        endif
911       
912        if(type==3)             //binned distances AND sld, slow parts XOPed
913                est = 0.576 + 4.22e-7*nx+1.76e-8*nx^2           //with threading
914               
915//              est /= 100                                      // per q now...
916//              est *= nq
917        endif
918       
919        return(est)
920End
Note: See TracBrowser for help on using the repository browser.