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

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

fixing the shape filling functions to explicitly take fill as a parameter, and make the default fill value = 10 == 10 e-7 A-2

More additions to the RealSpeca? help file

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=10
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,"fill SLD value"
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=10
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,"fill SLD value"
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.