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

Last change on this file since 838 was 838, checked in by srkline, 10 years ago

adding two new help files for real space modeling and description of the 2D resolution function

cleaning up the FFT routines for addition as a beta operation

File size: 18.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3// utility functions and procedures for displaying information
4// setting the matrix, and doing the calculations
5//
6
7
8
9///////////
10//
11// put the SLD multiplier on the panel somewhere - only 2 values are allowed - so use
12// radio buttons, or just display the value of the global, make the change at the command
13// line -- or keep the multiplier at 10^-7 and always use 2 digits (10^-7 is probably
14// better to use, just to allow more flexibility from the start).
15
16// TO DO:
17// -- I need to change a lot of routines (most notably Gizmo) to take "10" as the default SLD
18//               rather than "1"
19//
20////////////
21
22//
23
24//***********
25//For a 3D slicer view of dmat:
26//
27//Duplicate/O dmat dmatView
28//
29//Open a New 3D voxelgram. Don't set any of the levels
30//Open the 3D Slicer under the Gizmo menu
31//Make the X, Y, and Z slices
32//(Yellow Hot seems to be the best)
33//
34//// for a better view
35//dmatView = log(dmat)
36//
37//// if it all goes blank after the log transform get rid of the INF
38//dmatView = numtype(dmatView)== 0 ? dmatView : 0
39//*************
40
41
42///PANEL
43///// panel procedures
44
45Proc Init_FFT()
46        DoWindow/F FFT_Panel
47        if(V_flag==0)
48                SetDataFolder root:
49                //init the globals
50                Variable/G root:FFT_N=128
51                Variable/G root:FFT_T=5
52                Variable/G root:FFT_Qmax = 0
53                Variable/G root:FFT_QmaxReal = 0
54                Variable/G root:FFT_DQ=0
55                Variable/G root:FFT_estTime = 0
56               
57                Variable/G root:FFT_SolventSLD = 0
58                Variable/G root:FFT_delRho = 1e-7                       //multiplier for SLD (other value is 1e-7)
59               
60                FFT_Qmax :=2*pi/FFT_T
61                FFT_QmaxReal := FFT_Qmax/2
62                FFT_DQ := 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
69End
70
71Proc FFT_Panel()
72        PauseUpdate; Silent 1           // building window...
73        NewPanel /W=(1452,44,1768,531)/K=1 as "FFT_Panel"
74        DoWindow/C FFT_Panel
75        SetDrawLayer UserBack
76        DrawLine 5,68,311,68
77        DrawLine 5,142,311,142
78        DrawLine 5,250,311,250
79        SetVariable FFTSetVar_0,pos={7,7},size={150,15},title="Cells per edge (N)"
80        SetVariable FFTSetVar_0,limits={50,512,2},value= FFT_N,live= 1
81        SetVariable FFTSetVar_1,pos={7,33},size={150,15},title="Length per Cell (T)"
82        SetVariable FFTSetVar_1,limits={1,500,0.2},value= FFT_T,live= 1
83        SetVariable FFTSetVar_2,pos={183,26},size={120,15},title="Real Qmax"
84        SetVariable FFTSetVar_2,limits={0,0,0},value= FFT_QmaxReal,noedit= 1,live= 1
85        SetVariable FFTSetVar_3,pos={183,48},size={120,15},title="delta Q (A)"
86        SetVariable FFTSetVar_3,limits={0,0,0},value= FFT_DQ,noedit= 1,live= 1
87        Button FFTButton_0,pos={15,79},size={90,20},proc=FFT_MakeMatrixButtonProc,title="Make Matrix"
88        Button FFTButton_1,pos={14,157},size={90,20},proc=FFTMakeGizmoButtonProc,title="Make Gizmo"
89        Button FFTButton_2,pos={14,187},size={100,20},proc=FFTDrawSphereButtonProc,title="Draw Sphere"
90        Button FFTButton_3,pos={14,265},size={70,20},proc=DoTheFFT_ButtonProc,title="Do FFT"
91        Button FFTButton_4,pos={180,264},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot FFT Results"
92        Button FFTButton_5,pos={13,218},size={120,20},proc=FFTDrawZCylinderButtonProc,title="Draw Cylinder"
93        Button FFTButton_6,pos={134,79},size={90,20},proc=FFTEraseMatrixButtonProc,title="Erase Matrix"
94        Button FFTButton_7,pos={13,329},size={130,20},proc=FFT_BinnedSpheresButtonProc,title="Do Binned Debye"
95        Button FFTButton_7a,pos={180,329},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot Binned Results"
96
97        Button FFTButton_8,pos={13,297},size={130,20},proc=FFT_AltiSpheresButtonProc,title="Do Debye Spheres"
98        Button FFTButton_8a,pos={180,297},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot Debye Results"
99
100        Button FFTButton_14,pos={13,360},size={130,20},proc=FFT_BinnedSLDButtonProc,title="Do Binned SLD"
101        Button FFTButton_14a,pos={180,360},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot SLD Results"
102
103        SetVariable FFTSetVar_4,pos={201,4},size={100,15},title="FFT time(s)"
104        SetVariable FFTSetVar_4,limits={0,0,0},value= FFT_estTime,noedit= 1,live= 1,format="%d"
105        Button FFTButton_9,pos={200,400},size={100,20},proc=FFT_Get2DSlice,title="Get 2D Slice"
106        Button FFTButton_10,pos={169,156},size={130,20},proc=FFT_TransposeMat,title="Transpose Matrix"
107        Button FFTButton_11,pos={169,189},size={130,20},proc=FFT_RotateMat,title="Rotate Matrix"
108        Button FFTButton_12,pos={168,219},size={130,20},proc=FFT_AddRotatedObject,title="Add Rotated Obj"
109        Button FFTButton_12,disable=2           // hide this button
110        Button FFTButton_13,pos={14,109},size={120,20},proc=FFTFillSolventMatrixProc,title="Solvent Matrix"
111        SetVariable FFTSetVar_5,pos={155,111},size={150,15},title="Solvent SLD (10^-7)"
112        SetVariable FFTSetVar_5,limits={-99,99,1},value= FFT_SolventSLD,live= 1
113        Button FFTButton_15,pos={209,430},size={90,20},proc=Interp2DSliceButton,title="Interp 2D"
114        Button FFTButton_16,pos={14,460},size={70,20},proc=FFTHelpButton,title="Help"
115EndMacro
116
117Function FFTHelpButton(ba) : ButtonControl
118        STRUCT WMButtonAction &ba
119       
120        String win = ba.win
121
122        switch (ba.eventCode)
123                case 2:
124                        // click code here
125                        DisplayHelpTopic/Z/K=1 "Real-Space Modeling of SANS Data"
126                        if(V_flag !=0)
127                                DoAlert 0,"The Real-Space Modeling Help file could not be found"
128                        endif
129                        break
130        endswitch
131
132        return 0
133End
134
135
136Function Interp2DSliceButton(ba) : ButtonControl
137        STRUCT WMButtonAction &ba
138
139        switch( ba.eventCode )
140                case 2: // mouse up
141                        // click code here
142                        String folderStr=""
143                        Prompt folderStr, "Pick a 2D data folder"               // Set prompt for x param
144                        DoPrompt "Enter data folder", folderStr
145                        if (V_Flag || strlen(folderStr) == 0)
146                                return -1                                                               // User canceled, null string entered
147                        endif                   
148       
149                        Interpolate2DSliceToData(folderStr)
150                       
151                        //display the 2D plane
152                        DoWindow FFT_Interp2D_log
153                        if(V_flag == 0)
154                                Execute "Display2DInterpSlice_log()"
155                        endif
156                       
157                        break
158        endswitch
159
160        return 0
161End
162
163
164
165Function FFT_Get2DSlice(ctrlName) : ButtonControl
166        String ctrlName
167               
168        WAVE/Z dmat = root:dmat
169        if(WaveExists(dmat)==1)
170                get2DSlice(dmat)
171        endif
172       
173        //display the 2D plane
174        DoWindow FFT_Slice2D
175        if(V_flag == 0)
176                Execute "Display2DSlice()"
177        endif
178       
179        //display the 2D plane
180        DoWindow FFT_Slice2D_log
181        if(V_flag == 0)
182                Execute "Display2DSlice_log()"
183        endif
184       
185        //display the 1D binning of the 2D plane
186        DoWindow FFT_Slice1D
187        if(V_flag == 0)
188                Execute "Slice2_1D()"
189        endif
190       
191        return(0)       
192End
193
194Proc Display2DSlice()
195        PauseUpdate; Silent 1           // building window...
196        Display /W=(1038,44,1404,403)
197        DoWindow/C FFT_Slice2D
198        AppendImage/T detPlane
199        ModifyImage detPlane ctab= {*,*,YellowHot,0}
200        ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14
201        ModifyGraph mirror=2
202        ModifyGraph nticks=4
203        ModifyGraph minor=1
204        ModifyGraph fSize=9
205        ModifyGraph standoff=0
206        ModifyGraph tkLblRot(left)=90
207        ModifyGraph btLen=3
208        ModifyGraph tlOffset=-2
209        SetAxis/A/R left
210End
211
212Proc Display2DInterpSlice_log()
213        PauseUpdate; Silent 1           // building window...
214        Display /W=(1038,44,1404,403)
215        DoWindow/C FFT_Interp2D_log
216        AppendImage/T interp2DSlice_log
217        ModifyImage interp2DSlice_log ctab= {*,*,YellowHot,0}
218        ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14
219        ModifyGraph mirror=2
220        ModifyGraph nticks=4
221        ModifyGraph minor=1
222        ModifyGraph fSize=9
223        ModifyGraph standoff=0
224        ModifyGraph tkLblRot(left)=90
225        ModifyGraph btLen=3
226        ModifyGraph tlOffset=-2
227        SetAxis/A/R left
228End
229
230Proc Display2DSlice_log()
231        PauseUpdate; Silent 1           // building window...
232        Display /W=(1038,44,1404,403)
233        DoWindow/C FFT_Slice2D_log
234        AppendImage/T logP
235        ModifyImage logP ctab= {*,*,YellowHot,0}
236        ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14
237        ModifyGraph mirror=2
238        ModifyGraph nticks=4
239        ModifyGraph minor=1
240        ModifyGraph fSize=9
241        ModifyGraph standoff=0
242        ModifyGraph tkLblRot(left)=90
243        ModifyGraph btLen=3
244        ModifyGraph tlOffset=-2
245        SetAxis/A/R left
246End
247Proc Slice2_1D()
248        PauseUpdate; Silent 1           // building window...
249        Display /W=(1034,425,1406,763) iBin_2d vs qBin_2d
250        DoWindow/C FFT_Slice1D
251        ModifyGraph mode=4
252        ModifyGraph marker=19
253        ModifyGraph msize=2
254        ModifyGraph grid=1
255        ModifyGraph log=1
256        ModifyGraph mirror=2
257        Legend
258EndMacro
259
260
261Function FFT_TransposeMat(ctrlName) : ButtonControl
262        String ctrlName
263
264        Variable mode=1
265        Prompt mode,"Transform XYZ to:",popup,"XZY;ZXY;ZYX;YZX;YXZ;"
266        DoPrompt "Transform mode",mode
267        if (V_Flag)
268                return 0                                                                        // user canceled
269        endif
270       
271        fFFT_TransposeMat(mode)
272       
273        return (0)
274End
275
276// starting in XYZ
277// mode 1;2;3;4;5;
278//correspond to
279// "XZY;ZXY;ZYX;YZX;YXZ;"
280//
281Function fFFT_TransposeMat(mode)
282        Variable mode
283
284        WAVE/Z mat=mat
285        ImageTransform /G=(mode) transposeVol mat
286        WAVE M_VolumeTranspose=M_VolumeTranspose
287        Duplicate/O M_VolumeTranspose mat
288       
289        return(0)       
290End
291
292Function FFT_RotateMat(ctrlName) : ButtonControl
293        String ctrlName
294       
295        Variable degree=45,sense=1
296        Prompt degree, "Degrees of rotation around Z-axis:"
297//      Prompt sense, "Direction of rotation:",popup,"CW;CCW;"
298        DoPrompt "Enter paramters for rotation", degree
299       
300        if (V_Flag)
301                return 0                                                                        // user canceled
302        endif
303       
304        fFFT_RotateMat(degree)
305        return(0)
306End
307// note that the rotation is not perfect. if the rotation produces an
308// odd number of pixels, then the object will "walk" one pixel. Also, some
309// small artifacts are introduced, simply due to the voxelization of the object
310// as it is rotated. rotating a cylinder 10 then 350 shows a few extra "bumps" on
311// the surface, but the calculated scattering is virtually identical.
312//
313// these issues, may be correctable, if needed
314//
315Function fFFT_RotateMat(degree)
316        Variable degree
317       
318        Variable fill=0
319       
320        WAVE mat = mat
321        Variable dx,dy,dz,nx,ny,nz
322        dx = DimSize(mat,0)
323        dy = DimSize(mat,1)
324        dz = DimSize(mat,2)
325               
326        //? the /W and /C flags seem to be special cases for 90 deg rotations
327        ImageRotate /A=(degree)/E=(fill) mat            //mat is not overwritten, unless /O
328       
329        nx = DimSize(M_RotatedImage,0)
330        ny = DimSize(M_RotatedImage,1)
331        nz = DimSize(M_RotatedImage,2)
332//      Print "rotated dims = ",dx,dy,dz,nx,ny,nz
333        Variable delx,dely,odd=0
334        delx = nx-dx
335        dely = ny-dy
336       
337        if(mod(delx,2) != 0)            //sometimes the new x,y dims are odd!
338                odd = 1
339        endif
340//      delx = (delx + odd)/2
341//      dely = (dely + odd)/2
342        delx = trunc(delx/2)
343        dely = trunc(dely/2)
344
345        //+odd removes an extra point if there is one
346        Duplicate/O/R=[delx+odd,nx-delx-1][dely+odd,ny-dely-1][0,nz-1] M_RotatedImage mat
347
348// - not sure why the duplicate operation changes the start and delta of mat, but I
349// need to reset the start and delta to be sure that the display is correct, and that the scaling
350// is read correctly later
351        SetScale/P x 0,1,"", mat
352        SetScale/P y 0,1,"", mat
353       
354        nx = DimSize(mat,0)
355        ny = DimSize(mat,1)
356        nz = DimSize(mat,2)
357//      Print "redim = ",dx,dy,dz,nx,ny,nz
358       
359        KillWaves/Z M_RotatedImage
360       
361End
362
363// also look for ImageTransform operations that will allow translation of the objects.
364// then simple objects could be oriented @0,0,0, and translated to the correct position (and added)
365//
366Function FFT_AddRotatedObject(ctrlName) : ButtonControl
367        String ctrlName
368
369        Print "Not yet implemented"
370End
371
372
373Function FFT_MakeMatrixButtonProc(ctrlName) : ButtonControl
374        String ctrlName
375
376        NVAR nn=root:FFT_N
377        if(mod(nn, 2 ) ==1)             //force an even number for FFT
378                nn +=1
379        endif
380        MakeMatrix("mat",nn,nn,nn)
381End
382
383Function FFTMakeGizmoButtonProc(ctrlName) : ButtonControl
384        String ctrlName
385       
386        DoWindow/F Gizmo_VoxelMat
387        if(V_flag==0)
388                Execute "Gizmo_VoxelMat()"
389        endif
390End
391
392Function FFTDrawSphereButtonProc(ctrlName)  : ButtonControl
393        String ctrlName
394       
395        Execute "FFTDrawSphereProc()"
396End
397
398Proc FFTDrawSphereProc(matStr,rad,xc,yc,zc,fill)
399        String matStr="mat"
400        Variable rad=25,xc=50,yc=50,zc=50,fill=1
401        Prompt matStr,"the wave"                //,popup,WaveList("*",";","")
402        Prompt rad,"enter real radius (A)"
403        Prompt xc,"enter the X-center"
404        Prompt yc,"enter the Y-center"
405        Prompt zc,"enter the Z-center"
406        Prompt fill,"1= fill, 0 = erase"
407       
408        Variable grid=root:FFT_T
409       
410        FillSphereRadius($matStr,grid,rad,xc,yc,zc,fill)
411       
412End
413
414Function FFTDrawZCylinderButtonProc(ctrlName)  : ButtonControl
415        String ctrlName
416       
417        Execute "FFTDrawCylinder()"
418End
419
420Proc FFTDrawCylinder(direction,matStr,rad,len,xc,yc,zc,fill)
421        String direction
422        String matStr="mat"
423        Variable rad=25,len=300,xc=50,yc=50,zc=50,fill=1
424        Prompt direction, "Direction", popup "X;Y;Z;"
425        Prompt matStr,"the wave"                //,popup,WaveList("*",";","")
426        Prompt rad,"enter real radius (A)"
427        Prompt len,"enter length (A)"
428        Prompt xc,"enter the X-center"
429        Prompt yc,"enter the Y-center"
430        Prompt zc,"enter the Z-center"
431        Prompt fill,"1= fill, 0 = erase"
432       
433       
434        Variable grid=root:FFT_T
435       
436        if(cmpstr(direction,"X")==0)
437                FillXCylinder($matStr,grid,rad,xc,yc,zc,len,fill)
438        endif
439        if(cmpstr(direction,"Y")==0)
440                FillYCylinder($matStr,grid,rad,xc,yc,zc,len,fill)
441        endif
442        if(cmpstr(direction,"Z")==0)
443                FillZCylinder($matStr,grid,rad,xc,yc,zc,len,fill)
444        endif
445       
446End
447
448
449Function DoTheFFT_ButtonProc(ctrlName) : ButtonControl
450        String ctrlName
451
452        Execute "DoFFT()"
453End
454
455Function FFT_PlotResultsButtonProc(ctrlName) : ButtonControl
456        String ctrlName
457
458        Variable first=0
459        DoWindow/F FFT_IQ
460        if(V_flag==0)
461                first = 1
462                Display /W=(295,44,627,302)
463                DoWindow/C FFT_IQ
464        Endif
465       
466        // append the desired data, if it's not already there
467        // FFTButton_4 = FFT            = iBin
468        // FFTButton_7a = binned = _XOP
469        // FFTButton_8a = Debye = _full
470        // FFTButton_14a = SLD = _SLD
471        strswitch(ctrlName)     
472                case "FFTButton_4":
473                        if(!isTraceOnGraph("iBin","FFT_IQ") && exists("iBin")==1)               //only append if it's not already there
474                                AppendToGraph /W=FFT_IQ iBin vs qBin
475                                ModifyGraph mode=4,marker=19,msize=2,rgb(iBin)=(65535,0,0)
476                        endif
477                        break
478                case "FFTButton_7a":
479                        if(!isTraceOnGraph("ival_XOP","FFT_IQ") && exists("ival_XOP")==1)               //only append if it's not already there
480                                AppendToGraph /W=FFT_IQ ival_XOP vs qval_XOP
481                                ModifyGraph mode=4,marker=19,msize=2,rgb(ival_XOP)=(1,12815,52428)
482                        endif           
483                        break
484                case "FFTButton_8a":
485                        if(!isTraceOnGraph("ival_full","FFT_IQ") && exists("ival_full")==1)             //only append if it's not already there
486                                AppendToGraph /W=FFT_IQ ival_full vs qval_full
487                                ModifyGraph mode=4,marker=19,msize=2,rgb(ival_full)=(0,0,0)
488                        endif           
489                        break
490                case "FFTButton_14a":
491                        if(!isTraceOnGraph("ival_SLD","FFT_IQ") && exists("ival_SLD")==1)               //only append if it's not already there
492                                AppendToGraph /W=FFT_IQ ival_SLD vs qval_SLD
493                                ModifyGraph mode=4,marker=19,msize=2,rgb(ival_SLD)=(2,39321,1)
494                        endif           
495                        break
496        endswitch
497       
498        if(first)
499                ModifyGraph mode=4
500                ModifyGraph marker=19
501                ModifyGraph msize=2
502                ModifyGraph gaps=0
503                ModifyGraph grid=1
504                ModifyGraph log=1
505                ModifyGraph mirror=2
506                Legend
507        endif
508       
509        return(0)
510End
511
512Function isTraceOnGraph(traceStr,winStr)
513        String traceStr,winStr
514       
515        Variable isOn=0
516        String str
517        str = TraceNameList("winStr",";",1)             //only normal traces
518        isOn = stringmatch(str,traceStr)                //must be an exact match
519
520        return(isOn)
521End
522
523Function FFTEraseMatrixButtonProc(ctrlName) : ButtonControl
524        String ctrlName
525       
526        Wave mat=root:mat
527        FastOp mat=0
528        return(0)
529End
530
531Function FFTFillSolventMatrixProc(ctrlName) : ButtonControl
532        String ctrlName
533       
534        Wave mat=root:mat
535        NVAR val=root:FFT_SolventSLD
536        mat=val
537        return(0)
538End
539
540Function FFT_AltiSpheresButtonProc(ctrlName) : ButtonControl
541        String ctrlName
542
543        Execute "DoSpheresCalcFFTPanel()"
544End
545
546Function FFT_BinnedSpheresButtonProc(ctrlName) : ButtonControl
547        String ctrlName
548
549        Execute "DoBinnedSpheresCalcFFTPanel()"
550End
551
552Function FFT_BinnedSLDButtonProc(ctrlName) : ButtonControl
553        String ctrlName
554
555        Execute "DoBinnedSLDCalcFFTPanel()"
556End
557
558
559
560
561
562/////UTILITIES
563
564// inverts the values in the matrix 0<->1
565Function InvertMatrixFill(mat)
566        Wave mat
567       
568        mat = (mat==1) ? 0 : 1
569End
570
571//overwrites any existing matrix
572// matrix is byte, to save space
573//forces the name to be "mat"
574//
575// switched to signed byte
576//
577Function MakeMatrix(nam,xd,yd,zd)
578        String nam
579        Variable xd,yd,zd
580       
581        nam="mat"
582        Print "Matrix has been created and named \"mat\""
583//      Make/O/U/B/N=(xd,yd,zd) $nam
584        Make/O/B/N=(xd,yd,zd) $nam
585End
586
587
588// calculate the average center-to-center distance between points
589// assumes evenly spaced on a cubic grid
590//
591Proc Center_to_Center(np)
592        Variable np = 100
593       
594        Variable Nedge = root:FFT_N
595        Variable Tscale = root:FFT_T
596       
597        Variable davg
598       
599        davg = (Nedge*Tscale)^3 / np
600        davg = davg^(1/3)
601        Print "Average separation (A) = ",davg
602       
603End
604
605// calculate the number of points required on a grid
606// to yield a given average center-to-center distance between points
607// assumes evenly spaced on a cubic grid
608//
609// davg and Tscale are the same units, Angstroms
610//
611Proc Davg_to_Np(davg)
612        Variable davg = 400
613       
614        Variable Nedge = root:FFT_N
615        Variable Tscale = root:FFT_T
616       
617        Variable np
618       
619        np = (Nedge*Tscale)^3 / davg^3
620        Print "Number of points required = ",np
621       
622End
623
624// The matrix is not necessarily 0|1, this reports the number of filled voxels
625// - needed to estimate the time required for the AltiVec_Spheres calculation
626Proc NumberOfPoints()
627       
628        Print "Number of points = ",NumberOfPoints_Occ(root:mat)
629        Print "Fraction occupied = ",VolumeFraction_Occ(root:mat)
630        Print "Overall Cube Edge [A] = ",root:FFT_T * root:FFT_N
631       
632End
633
634Function NumberOfPoints_Occ(m)
635        Wave m
636       
637        Variable num = NonZeroValues(m)
638
639        return(num)
640End
641
642
643Function VolumeFraction_Occ(m)
644        Wave m
645       
646        Variable num = NonZeroValues(m)
647        Variable dim = DimSize(m,0)
648
649        return(num/dim^3)
650End
651
652//it's a byte wave, so I can't use the trick of setting the zero values to NaN
653// since NaN can't be expressed as byte. So make it binary 0|1
654//
655// - the assumption here is that the defined solvent value is not part of the
656// particle volume.
657//
658Function NonZeroValues(m)
659        Wave m
660       
661        Variable num
662        NVAR val = root:FFT_SolventSLD
663       
664//      Variable t1=ticks,dim
665//      dim = DimSize(m, 0 )            //assume NxNxN
666        Duplicate/O m,mz
667       
668        MultiThread mz = (m[p][q] != val) ? 1 : 0
669       
670        WaveStats/Q/M=1 mz              // NaN and Inf are not reported in V_npnts
671       
672        num = V_npnts*V_avg
673       
674        KillWaves/Z mz
675//      Print "Number of points = ",num
676//      Print "Fraction occupied = ",num/V_npnts
677       
678//      Print "time = ",(ticks-t1)/60.15
679        return(num)
680End
681
682// returns estimate in seconds
683//
684//              updated for quad-core iMac, 2010
685//
686// type 3 = binned distances and SLD
687// type 2 = binned distances
688// type 1 is rather meaningless
689// type 0 = is the Debye, double sum (multithreaded)
690//
691//
692// - types 2 and 3, the binned methods are inherently AAO, so there is no significant
693// dependence on the number of q-values (unless it's  >> 1000). So values reported are
694// for 100 q-values.
695//
696// There is a function Timing_Method(type) in FFT_ConcentratedSpheres.ipf to automate some of this
697//
698Function EstimatedTime(nx,nq,type)
699        Variable nx,nq,type
700       
701        Variable est=0
702       
703        if(type==0)             // "full" XOP
704                est = 0.4 + 1.622e-5*nx+4.86e-7*nx^2
705                est /= 100
706                est *= nq
707        endif
708       
709        if(type==1)             //Igor function (much slower, 9x)
710                est = (nx^2)*0.0000517          //empirical values (seconds) for igor code, 100 q-values
711                est /= 100                                      // per q now...
712                est *= nq
713        endif
714       
715        if(type==2)             //binned distances, slow parts XOPed
716                est = 0.0680 + 2.94e-6*nx+4.51e-9*nx^2                  //with threading
717               
718//              est /= 100                                      // per q now...
719//              est *= nq
720        endif
721       
722        if(type==3)             //binned distances AND sld, slow parts XOPed
723                est = 0.576 + 4.22e-7*nx+1.76e-8*nx^2           //with threading
724               
725//              est /= 100                                      // per q now...
726//              est *= nq
727        endif
728       
729        return(est)
730End
Note: See TracBrowser for help on using the repository browser.