source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_ConcentratedSpheres.ipf @ 798

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

I'm not sure this is a great idea, but I'm putting the FFT / Debye sphere work that I have completed into SVN. It's all really rough - the math, I believe is correct, but the interface if really, really rough. But it's not going to develop without help.
File size: 9.2 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3
4
5
6Proc concSphereLoop()
7
8// constant for all steps
9        root:FFT_T = 5
10        root:FFT_N = 128
11        root:FFT_SolventSLD = 0
12       
13// always start fresh
14        FFT_MakeMatrixButtonProc("")
15        FFTEraseMatrixButtonProc("")
16       
17       
18//
19// this block tests the number of passes needed for a "good" average
20//
21//      testConcSpheres(10,20,0,5,"_a5")
22//      testConcSpheres(10,20,0,10,"_a10")
23//      testConcSpheres(10,20,0,20,"_a20")
24//      testConcSpheres(10,20,0,50,"_a50")
25//      testConcSpheres(10,20,0,100,"_a100")
26       
27//
28// this block tests the concentration
29//
30//      testConcSpheres(100,20,0,20,"_b")
31//      testConcSpheres(200,20,0,20,"_c")
32//      testConcSpheres(300,20,0,20,"_d")
33//      testConcSpheres(350,20,0,20,"_e")
34
35//      testConcSpheres(600,20,0,20,"_f")
36//      testConcSpheres(800,20,0,20,"_g")
37//      testConcSpheres(1200,20,0,20,"_h")
38//      testConcSpheres(1600,20,0,20,"_i")
39//
40//      testConcSpheres(2000,20,0,20,"_j")
41//      testConcSpheres(2600,20,0,20,"_k")
42//      testConcSpheres(3200,20,0,20,"_l")
43////
44//// this block tests the concentration and polydispersity
45////
46
47//      testConcSpheres(84,20,0.25,20,"_bp")
48//      testConcSpheres(168,20,0.25,20,"_cp")
49//      testConcSpheres(253,20,0.25,20,"_dp")
50//      testConcSpheres(295,20,0.25,20,"_ep")
51//
52//      testConcSpheres(505,20,0.25,20,"_fp")
53//      testConcSpheres(674,20,0.25,20,"_gp")
54//      testConcSpheres(1011,20,0.25,20,"_hp")
55//      testConcSpheres(1348,20,0.25,20,"_ip")
56
57        testConcSpheres(1685,20,0.25,20,"_jp")
58        testConcSpheres(2190,20,0.25,20,"_kp")
59        testConcSpheres(2696,20,0.25,20,"_lp")
60
61
62       
63//      // this block simply is to save the gizmo as a wave
64//
65//      DoWindow/F Gizmo_VoxelMat
66//     
67//      testConcSpheres(100,20,0.25,1,"_dum")
68//      DoUpdate/W=Gizmo_VoxelMat
69//      ExportGizmo wave as "p_phi_0p012"
70//     
71//      testConcSpheres(200,20,0.25,1,"_dum")
72//      DoUpdate/W=Gizmo_VoxelMat
73//      ExportGizmo wave as "p_phi_0p024"
74//     
75//      testConcSpheres(300,20,0.25,1,"_dum")
76//      DoUpdate/W=Gizmo_VoxelMat
77//      ExportGizmo wave as "p_phi_0p036"
78//     
79//      testConcSpheres(350,20,0.25,1,"_dum")
80//      DoUpdate/W=Gizmo_VoxelMat
81//      ExportGizmo wave as "p_phi_0p041"
82//     
83//     
84//
85//      testConcSpheres(600,20,0.25,1,"_dum")
86//      DoUpdate/W=Gizmo_VoxelMat
87//      ExportGizmo wave as "p_phi_0p070"
88//     
89//      testConcSpheres(800,20,0.25,1,"_dum")
90//      DoUpdate/W=Gizmo_VoxelMat
91//      ExportGizmo wave as "p_phi_0p094"
92//     
93//      testConcSpheres(1200,20,0.25,1,"_dum")
94//      DoUpdate/W=Gizmo_VoxelMat
95//      ExportGizmo wave as "p_phi_0p141"
96//     
97//      testConcSpheres(1600,20,0.25,1,"_dum")
98//      DoUpdate/W=Gizmo_VoxelMat
99//      ExportGizmo wave as "p_phi_0p188"
100//     
101//
102//
103//      testConcSpheres(2000,20,0.25,1,"_dum")
104//      DoUpdate/W=Gizmo_VoxelMat
105//      ExportGizmo wave as "p_phi_0p233"
106//     
107//      testConcSpheres(2600,20,0.25,1,"_dum")
108//      DoUpdate/W=Gizmo_VoxelMat
109//      ExportGizmo wave as "p_phi_0p299"
110//     
111//      testConcSpheres(3200,20,0.25,1,"_dum")
112//      DoUpdate/W=Gizmo_VoxelMat
113//      ExportGizmo wave as "p_phi_0p366"
114////   
115///////////////////////////     
116end
117
118
119//
120// Note that there is still an issue with the algorithm that skews towards
121// smaller radii as the concentration increases (for the polydisperse case)
122// if a sphere fails (and the larger ones are more likely to fail), a new position
123// AND a new sphere radius is selected
124//
125//
126// with the double do loop modification, now there is the possibility that it
127// will try forever to fill a sphere, but there was always that possibility before.
128//
129//
130Function TestConcSpheres(nSph,rad,pd,nPass,tagStr)
131        Variable nSph,rad,pd,nPass
132        String tagStr
133       
134        Variable ii,np,frac,jj,tmprad,t1
135        Variable row,col,lay,xt,yt,zt,err,fails
136        String printStr
137        Wave m=mat
138        NVAR grid=root:FFT_T
139        NVAR Nedge = root:FFT_N
140       
141        row=DimSize(m,0)               
142        col=DimSize(m,1)
143        lay=DimSize(m,2)
144       
145        Make/O/D/N=(nSph) failures,radius,skipped
146        Variable skipCount=0
147        failures = 0
148        radius = 0
149        skipped = 0
150       
151        t1 = ticks
152        for(ii=0;ii<npass;ii+=1)                //number of averaging passes
153                //fill the spheres into the box
154                m=0
155                skipCount = 0
156                failures = 0
157                radius = 0
158                skipped = 0
159               
160                for(jj=0;jj<nSph;jj+=1)         //number of spheres
161                       
162                        // pick a sphere radius
163                        if(pd !=0 )
164                                tmprad = Rad + gnoise(pd*Rad)
165                        else
166                                tmprad = rad
167                        endif
168                       
169                        fails = -1
170                        //now find a place for this sphere
171                        do
172                                fails += 1
173                                //find an unfilled voxel
174                                do
175                                        xt=trunc(abs(enoise(row)))              //distr betw (0,npt)
176                                        yt=trunc(abs(enoise(col)))
177                                        zt=trunc(abs(enoise(lay)))
178                                while(m[xt][yt][zt] == 1)
179                       
180                                //try to put the sphere there, and keep trying forever
181                                err = FillSphereRadiusNoOverlap(m,grid,tmprad,xt,yt,zt,1)
182                               
183//                              if(fails == 10)
184//                                      Print "failed 10x on tmprad = ",tmprad
185//                              endif
186//                              if(fails == 100)
187//                                      Print "failed 100x on tmprad = ",tmprad
188//                              endif
189                                if(fails == 1000)
190                                        Print "failed 1000x on tmprad, skipping = ",tmprad
191                                        skipped[jj] = 1
192                                        skipCount += 1
193                                        err = 0
194                                endif
195                        while(err==1)
196                        failures[jj] = fails
197                        radius[jj] = tmprad
198                                               
199                        if(mod(jj, 100 )        == 0)
200                                Print "sphere jj done, pass = ",jj,ii+1
201                                Print "time = ",(ticks-t1)/60.15
202                        endif
203                        if(jj > 2000 && mod(jj, 10 )    == 0)
204                                Print "sphere jj done, pass = ",jj,ii+1
205                                Print "time = ",(ticks-t1)/60.15
206                        endif
207                endfor
208                // spheres have been placed, do the calculation
209                ParseMatrix3D_rho(m)   
210               
211                Execute "DoFFT()"
212                sprintf printStr,"completed pass %d of %d \r",ii+1,npass
213                Print printStr
214                if(ii==0)
215                        Duplicate/O iBin, $("iBin"+tagStr),$("qBin"+tagStr),$("sBin"+tagStr),isq
216                        wave ib=$("iBin"+tagStr)
217                        wave qb=$("qBin"+tagStr)
218                        wave sb=$("sBin"+tagStr)
219                        Wave iBin=iBin
220                        Wave qbin=qBin
221                        Wave isq=isq
222                        qb=qBin
223                        isq = ib*ib
224                else
225                        ib += iBin
226                        isq += iBin*iBin
227                endif
228        endfor
229
230        ib /= npass
231        isq /= npass
232       
233        sb = sqrt( (isq-ib^2)/(npass - 1) )             // <I^2> - <I>^2
234       
235        Variable vol
236        NVAR delRho = root:FFT_delRho
237        np = NonZeroValues(m)
238        frac = np/Nedge^3
239//      delRho = 1e-6
240        vol = 4*pi/3*(rad)^3            //individual sphere volume based on the voxel equivalent
241       
242        Print "vol = ",vol
243        Print "frac = ",frac
244       
245//      ib *= delrho*delrho
246//      ib *= 1e8
247//      ib /= vol
248//      ib *= frac
249//      ib /= (Nedge)^3
250       
251        nSph -= skipCount
252       
253        String nStr
254        sprintf nstr,"T=%d;N=%d;Npass=%d;NSpheres=%d;Rad=%g;PD=%g;VolFrac=%g;",grid,Nedge,nPass,nSph,rad,pd,frac
255        Note ib, nStr
256
257       
258End
259
260
261Proc FillStatistics()
262
263        Variable num
264       
265        Wavestats/Q skipped
266        Print "Number skipped = ",V_sum
267       
268        Wavestats/Q radius
269        Print "average, sdev = ",V_avg,V_sdev
270        fails_vs_radius()
271       
272
273End
274
275// for testing of the timing
276// type is 0 | 2 | 3
277//
278Proc Timing_Method(type)
279        Variable type           // type of calculation, not used right now
280
281// constant for all steps
282        root:FFT_T = 5
283        root:FFT_N = 128
284        root:FFT_SolventSLD = 0
285        Variable num,qMin,qMax,grid,xc,yc,zc,rad,fill
286        num=100
287        qMin = 0.004
288        qMax = 0.4
289        grid = root:FFT_T
290        xc = 64
291        yc = 64
292        zc = 64
293       
294// always start fresh
295        FFT_MakeMatrixButtonProc("")
296        FFTEraseMatrixButtonProc("")
297       
298//// type 0
299        fill = 1
300        rad = 100
301        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
302        DoSpheresCalcFFTPanel(num,qMin,qMax)
303        FFTEraseMatrixButtonProc("")
304
305        rad = 80
306        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
307        DoSpheresCalcFFTPanel(num,qMin,qMax)
308        FFTEraseMatrixButtonProc("")
309
310        rad = 70
311        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
312        DoSpheresCalcFFTPanel(num,qMin,qMax)
313        FFTEraseMatrixButtonProc("")
314
315        rad = 50
316        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
317        DoSpheresCalcFFTPanel(num,qMin,qMax)
318        FFTEraseMatrixButtonProc("")
319
320        rad = 30
321        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
322        DoSpheresCalcFFTPanel(num,qMin,qMax)
323        FFTEraseMatrixButtonProc("")
324
325       
326//      type 2
327        fill = 1
328        rad = 150
329        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
330        DoBinnedSpheresCalcFFTPanel(num,qMin,qMax)
331        FFTEraseMatrixButtonProc("")
332       
333        rad = 120
334        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
335        DoBinnedSpheresCalcFFTPanel(num,qMin,qMax)
336        FFTEraseMatrixButtonProc("")
337
338        rad = 100
339        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
340        DoBinnedSpheresCalcFFTPanel(num,qMin,qMax)
341        FFTEraseMatrixButtonProc("")
342
343        rad = 80
344        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
345        DoBinnedSpheresCalcFFTPanel(num,qMin,qMax)
346        FFTEraseMatrixButtonProc("")
347
348        rad = 50
349        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
350        DoBinnedSpheresCalcFFTPanel(num,qMin,qMax)
351        FFTEraseMatrixButtonProc("")
352       
353
354// type 3
355        fill = 3
356        rad = 150
357        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
358        fill = 2
359        rad = 100
360        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
361        fill = 1
362        rad = 50
363        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
364        DoBinnedSLDCalcFFTPanel(num,qMin,qMax)
365        FFTEraseMatrixButtonProc("")
366
367        fill = 3
368        rad = 120
369        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
370        fill = 2
371        rad = 100
372        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
373        fill = 1
374        rad = 50
375        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
376        DoBinnedSLDCalcFFTPanel(num,qMin,qMax)
377        FFTEraseMatrixButtonProc("")
378       
379        fill = 3
380        rad = 100
381        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
382        fill = 2
383        rad = 70
384        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
385        fill = 1
386        rad = 50
387        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
388        DoBinnedSLDCalcFFTPanel(num,qMin,qMax)
389        FFTEraseMatrixButtonProc("")
390
391        fill = 3
392        rad = 80
393        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
394        fill = 2
395        rad = 70
396        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
397        fill = 1
398        rad = 50
399        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
400        DoBinnedSLDCalcFFTPanel(num,qMin,qMax)
401        FFTEraseMatrixButtonProc("")
402
403        fill = 3
404        rad = 50
405        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
406        fill = 2
407        rad = 40
408        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
409        fill = 1
410        rad = 30
411        FillSphereRadius(mat,grid,rad,xc,yc,zc,fill)
412        DoBinnedSLDCalcFFTPanel(num,qMin,qMax)
413        FFTEraseMatrixButtonProc("")
414
415
416End
Note: See TracBrowser for help on using the repository browser.