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

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

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