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

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

changes to FFT routines to clean things up for a beta release at the next startup. No functionality changes, just cleaning up the operation, interface, and menu items.

Some important changes to the polarization routines. New equations for calculating the coefficient matrix. Now appears to be correct. Also, proper proportions are used when adding multiple files together to the matrix.

Change to the Model Docs is an updated reference.

WorkFileUtils? was broken when I changed to use linear_data exclusively. This has now been fixed and it operates as expected.

PackageLoader? now has menu items (under Macros) for the Polarization routines. These may later be moved to the SANS menu.

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