source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_KR_Parser.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: 13.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3
4
5
6
7/// seems to work - but what do I do about fractional positions? when converting to a matrix?
8//
9//
10
11Function KR_Load()
12        Variable I, J, K, L, PT //integer indices loops, num cylinders, include or exclude sphere in circle
13        Variable STH, SPH, CTH, CPH, FTR  //sine and cosines and deg-->rad conversion: x rotn theta & y rotn phi
14        Variable  XMID, YMID, ZMID, XOUT, YOUT, ZOUT  //cartesian positions used in various calculations
15        Variable RR,HH  //RR is limit of loops, GG used as end of read param files--exit=2, NUM of cylinder
16        Variable  P5  //spheres half diameter shift from grid points (avoids zeros)
17        Variable X0, Y0,Z0
18        Variable PI2
19        Variable ix,nptW
20
21
22        LoadWave /G /N 
23        Print S_filename
24        Print S_wavenames
25       
26        //Make / O /N=0 OutputPoints
27        //      wave out=OutputPoints
28        //      variable num=numpnts(out)
29       
30        KillWaves/Z xx,yy,zz,rri,hti,sbp,rotx,roty,sld,gg
31       
32        Rename wave0, xx
33        Rename wave1, yy
34        Rename wave2, zz
35        Rename wave3, RRI
36        Rename wave4, HTI
37        Rename wave5, SBP
38        Rename wave6, ROTX
39        Rename wave7, ROTY
40        Rename wave8, SLD
41        Rename wave9, GG
42       
43        //print  NUM,xx,yy,zz,rri,hti,sbp,rotx,roty,sld,gg
44
45       
46        wave gg = gg
47        variable nn =-1,npts,cyl
48        npts = numpnts(GG)
49       
50        for (i=0;i<=npts;i+=1)
51                if (gg[i]==2)
52                        cyl = i+1
53                        break
54                        print "gg[i],i=",gg,i
55                endif
56        endfor
57        print"cyl=",cyl
58       
59
60        wave xx=xx
61        wave yy=yy
62        wave zz=zz
63        wave rri=rri
64        wave hti=hti
65        wave sbp=sbp
66        wave rotx=rotx
67        wave roty=roty
68        wave sld=sld
69       
70        // SBP = diameter of the spheres
71        NVAR FFT_T = root:FFT_T
72        FFT_T = SBP[0]
73        //
74       
75        Make/O/D/N=0 xoutW,youtW,zoutW,sbpW,sldW
76       
77        PI2=pi*2
78        FTR=PI2/360
79        // print "ftr=", ftr
80       
81        nptW = 0
82       
83        for(l=0;l<(cyl);L+=1)   //only change from run4
84        //for each cylinder of loop use index NUM
85        //calculate x & y rotation cos and sin
86                STH=SIN(Rotx[L]*FTR)
87                SPH=sin(roty[L]*FTR)
88                CTH=cos(rotx[L]*FTR)
89                CPH=cos(roty[L]*FTR)
90                //print "sth",sth
91                //print"L=",L
92                P5=SBP[L]/2  //set sphere centers' half-diameter displacement from grid (avoids glitches)
93                // print "p5 & sbp[L]",p5,sbp[L]
94       
95                RR=(RRI[L]/SBP[L])//as an index, Igor truncates the number to an integer....does NOT round it
96                RR=RR+1 //rr is the loop limit for square around final circle
97                HH=(HTI[L]/(2*SBP[L]))  //as an index, Igor truncates the number to an integer....does NOT round it
98                for(k=-HH;k<HH;k+=1)  // should have +1 for HH to complete to k=HH?????
99                        for(i=-RR;i<RR;i+=1)  //should this have i<RR+1 or in above RR=RR+2????
100                                for(j=-RR;j<RR;J+=1)
101                                        x0=sbp*i+P5
102                                        y0=SBP*j+P5
103                                        z0=SBP*k+p5
104                                        if((((y0^2)/(RRI[L]^2))+((x0^2)/(RRI[L]^2)))<=1)
105                                                IX=-1
106                                        else
107                                                IX=0
108                                        endif
109                                        xmid=x0
110                                        ymid=y0*cth+z0*sth
111                                        zmid=-y0*sth+z0*cth
112                                        // end rotation about x begin rotn about y on rotated pts
113                                        //
114                                        xout=xmid*cph-zmid*sph
115                                        xout=xx[L]+xout
116                                        yout=ymid
117                                        yout=yy[L]+yout
118                                        zout=xmid*sph+zmid*cph
119                                        zout=zz[L]+zout
120
121                                        // now print to wave file the point or not depending on whether ix<0 or not
122
123                                        if (ix<0)
124                                        //write to wave file
125                                                InsertPoints nptW,1,xoutW,youtW,zoutW,sbpW,sldW
126                                                xoutW[nptW] = xout
127                                                youtW[nptW] = yout
128                                                zoutW[nptW] = zout
129                                                sbpW[nptW] = sbp[L]
130                                                sldW[nptW] = sld[L]
131                                               
132                                                nptW +=1
133                                       
134                                                //print  xout,yout,zout,sbp[L],sld[L]
135                                        //else
136                                                //continue
137                                        endif  //for write or not
138                                endfor  // for j
139                        endfor  //  for i
140                endfor  //for k
141        endfor // for L
142
143        // rescale to the sphere size
144        xoutW /= FFT_T
145        youtW /= FFT_T
146        zoutW /= FFT_T
147       
148        return(0) // end do loop cycle for cylinders
149end
150
151
152
153/// seems to work - but what do I do about fractional positions? when converting to a matrix?
154//
155// the wave "gg" has been dropped, since it's only used as a flag in an old file loader
156//
157// NOW - SBP is FORCED to the value of FFT_T - no matter what is in the file.
158//
159Function KR_MultiCylinder(xx,yy,zz,rri,hti,sbp,rotx,roty,sld)
160        Wave xx,yy,zz,rri,hti,sbp,rotx,roty,sld
161
162        Variable I, J, K, L, PT //integer indices loops, num cylinders, include or exclude sphere in circle
163        Variable STH, SPH, CTH, CPH, FTR  //sine and cosines and deg-->rad conversion: x rotn theta & y rotn phi
164        Variable  XMID, YMID, ZMID, XOUT, YOUT, ZOUT  //cartesian positions used in various calculations
165        Variable RR,HH  //RR is limit of loops, GG used as end of read param files--exit=2, NUM of cylinder
166        Variable  P5  //spheres half diameter shift from grid points (avoids zeros)
167        Variable X0, Y0,Z0
168        Variable PI2
169        Variable ix,nptW
170       
171        NVAR FFT_T = root:FFT_T
172//      FFT_T = sbp[0]
173        sbp[0] = FFT_T
174       
175        variable npts,cyl
176        npts = numpnts(xx)
177        cyl = npts
178               
179        Make/O/D/N=0 xoutW,youtW,zoutW,sbpW,sldW
180       
181        PI2=pi*2
182        FTR=PI2/360
183       
184        nptW = 0
185       
186        for(l=0;l<(cyl);L+=1)   //only change from run4
187        //for each cylinder of loop use index NUM
188        //calculate x & y rotation cos and sin
189                STH=SIN(Rotx[L]*FTR)
190                SPH=sin(roty[L]*FTR)
191                CTH=cos(rotx[L]*FTR)
192                CPH=cos(roty[L]*FTR)
193                //print "sth",sth
194                //print"L=",L
195                P5=SBP[L]/2  //set sphere centers' half-diameter displacement from grid (avoids glitches)
196                // print "p5 & sbp[L]",p5,sbp[L]
197       
198                RR=(RRI[L]/SBP[L])//as an index, Igor truncates the number to an integer....does NOT round it
199                RR=RR+1 //rr is the loop limit for square around final circle
200                HH=(HTI[L]/(2*SBP[L]))  //as an index, Igor truncates the number to an integer....does NOT round it
201                for(k=-HH;k<HH;k+=1)  // should have +1 for HH to complete to k=HH?????
202                        for(i=-RR;i<RR;i+=1)  //should this have i<RR+1 or in above RR=RR+2????
203                                for(j=-RR;j<RR;J+=1)
204                                        x0=sbp*i+P5
205                                        y0=SBP*j+P5
206                                        z0=SBP*k+p5
207                                        if((((y0^2)/(RRI[L]^2))+((x0^2)/(RRI[L]^2)))<=1)
208                                                IX=-1
209                                        else
210                                                IX=0
211                                        endif
212                                        xmid=x0
213                                        ymid=y0*cth+z0*sth
214                                        zmid=-y0*sth+z0*cth
215                                        // end rotation about x begin rotn about y on rotated pts
216                                        //
217                                        xout=xmid*cph-zmid*sph
218                                        xout=xx[L]+xout
219                                        yout=ymid
220                                        yout=yy[L]+yout
221                                        zout=xmid*sph+zmid*cph
222                                        zout=zz[L]+zout
223
224                                        // now print to wave file the point or not depending on whether ix<0 or not
225
226                                        if (ix<0)
227                                        //write to wave file
228                                                InsertPoints nptW,1,xoutW,youtW,zoutW,sbpW,sldW
229                                                xoutW[nptW] = xout
230                                                youtW[nptW] = yout
231                                                zoutW[nptW] = zout
232                                                sbpW[nptW] = sbp[L]
233                                                sldW[nptW] = sld[L]
234                                               
235                                                nptW +=1
236                                       
237                                                //print  xout,yout,zout,sbp[L],sld[L]
238                                        //else
239                                                //continue
240                                        endif  //for write or not
241                                endfor  // for j
242                        endfor  //  for i
243                endfor  //for k
244        endfor // for L
245
246
247        // rescale to the sphere size
248        xoutW /= FFT_T
249        youtW /= FFT_T
250        zoutW /= FFT_T
251       
252        return(0) // end do loop cycle for cylinders
253end
254
255// triplet to display as a scatter plot in Gizmo
256//
257// will overwrite the triplet
258//
259Function MakeTriplet(xoutW,youtW,zoutW)
260        Wave xoutW,youtW,zoutW
261       
262        KillWaves/Z triplet
263        concatenate/O {xoutW,youtW,zoutW},triplet
264end
265
266
267
268Proc KR_LoadAndFill()
269
270        KR_Load()
271        XYZV_FillMat(xoutW,youtW,ZoutW,sldW,1)                  //last 1 will erase the matrix
272        MakeTriplet(xoutW,youtW,zoutW)
273         
274        DoBinned_KR_FFTPanel()
275        Print "now display the gizmo, triplet or use one of the calculation methods"
276       
277End
278
279Proc KR_CalcFromInput()
280
281        KR_MultiCylinder(xCtr,yCtr,zCtr,rad,length,sphereDiam,rot_x,rot_y,SLD_sph)
282        // these are really just for display, or if the FFT of mat is wanted later.
283        XYZV_FillMat(xoutW,youtW,ZoutW,sldW,1)                  //last 1 will erase the matrix
284        MakeTriplet(xoutW,youtW,zoutW)
285
286// and the calculation. Assumes SLDs are all the same   
287        DoBinned_KR_FFTPanel(100,0.004,0.5)
288       
289End
290
291//called from the FFT method panel
292//
293// in this method, the distances are binned as by Otto Glatter, and has been partially XOPed
294//
295// if the number of bins is too high (say 100000), then using the non-integer XYZ will
296// be 2-3 times slower since there will be a lot more bins - then the loop over the q-values at the
297// very end will be what is significantly slower. If the number of bins is reduced to 10000 (as suggested
298// in Otto's book, p.160), then the two methods (types 12 and 2) give very similar timing, and the
299// results are indistinguishable.
300//
301Proc DoBinned_KR_FFTPanel(num,qMin,qMax)
302        Variable num=100,qmin=0.004,qmax=0.5
303       
304        Variable t1
305        String qStr="qval_KR",iStr="ival_KR"            //default wave names, always overwritten
306        Variable grid
307       
308        grid=root:FFT_T
309
310        Make/O/D/N=(num) $qStr,$iStr
311        $qStr = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))         
312       
313        Variable estTime,nx
314        String str = ""
315
316        nx = NonZeroValues(mat)
317       
318        estTime = EstimatedTime(nx,num,2)               // 0 =  XOP, 1 = no XOP, 2 = binned distances
319        sprintf str, "Estimated time for the calculation is %g seconds. Proceed?",estTime
320        DoAlert 1,str
321        if(V_Flag==1)           //yes, proceed
322                t1=ticks
323                fDoCalc($qStr,$iStr,grid,12,1)
324//              Printf "Elapsed AltiSpheres time = %g seconds\r\r",(ticks-t1)/60.15
325        Endif
326End
327
328Function fDoBinned_KR_FFTPanel(num,qMin,qMax)
329        Variable num,qmin,qmax
330       
331        Variable t1,multiSLD,mode
332        String qStr="qval_KR",iStr="ival_KR"            //default wave names, always overwritten
333       
334        NVAR grid=root:FFT_T
335        ControlInfo/W=MultiCyl check_0
336        multiSLD = V_Value
337        if(multiSLD)
338                mode=13
339        else
340                mode=12
341        endif
342
343        Make/O/D/N=(num) qval_KR,ival_KR
344        qval_KR = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))               
345       
346        Variable estTime,nx,tooLong
347        String str = ""
348        tooLong = 300
349       
350        nx = NonZeroValues(mat)
351       
352        estTime = EstimatedTime(nx,num,2)               // 0 =  XOP, 1 = no XOP, 2 = binned distances
353        if(estTime > tooLong)
354                sprintf str, "Estimated time for the calculation is %g seconds. Proceed?",estTime
355                DoAlert 1,str
356        endif
357        if(V_Flag==1 || estTime < tooLong)              //yes, proceed
358                t1=ticks
359                fDoCalc(qval_KR,ival_KR,grid,mode,1)
360//              Printf "Elapsed AltiSpheres time = %g seconds\r\r",(ticks-t1)/60.15
361        Endif
362End
363
364/////////////////////////////////////
365// for each cylinder:
366//
367// xx,yy,zz,rri,hti,sbp,rotx,roty,sld
368// xx,yy,zz = center of cylinder
369// rri,hti = radius, height (units??)
370// sbp = ??? -- I think this is the diameter of the primary sphere
371// rotx, rotx = rotation angles (in degrees, but defined as ??)
372// sld = SLD of cylinder
373//
374//
375//
376// Put this into a panel with the table and the data
377// and fields for all of the inputs
378Macro Setup_KR_MultiCylinder()
379       
380        Make/O/D/N=0 xx,yy,zz,rri,hti,sbp,rotx,roty,sld
381        Variable/G root:KR_Qmin = 0.004
382        Variable/G root:KR_Qmax = 0.4
383        Variable/G root:KR_Npt = 100
384
385        NewPanel /W=(241,44,1169,458)/N=MultiCyl/K=1 as "Multi-Cylinder"
386        ModifyPanel cbRGB=(49825,57306,65535)
387
388        Button button_0,pos={45,79},size={100,20},proc=KR_Show3DButtonProc,title="Show 3D"
389        Button button_1,pos={46,51},size={100,20},proc=KR_Plot1DButtonProc,title="Plot 1D"
390        Button button_2,pos={178,80},size={120,20},proc=KR_DoCalcButtonProc,title="Do Calculation"
391        ValDisplay valdisp_0,pos={339,16},size={80,13},title="FFT_T"
392        ValDisplay valdisp_0,limits={0,0,0},barmisc={0,1000},value= #"root:FFT_T"
393        SetVariable setvar_0,pos={339,40},size={140,15},title="Q min (A)"
394        SetVariable setvar_0,limits={0,10,0},value= KR_Qmin
395        SetVariable setvar_1,pos={339,65},size={140,15},title="Q max (A)"
396        SetVariable setvar_1,limits={0,10,0},value= KR_Qmax
397        SetVariable setvar_2,pos={339,90},size={140,15},title="Num Pts"
398        SetVariable setvar_2,limits={10,500,0},value= KR_Npt
399        CheckBox check_0,pos={599,93},size={59,14},title="Multi SLD",value= 0
400
401        Edit/W=(18,117,889,378)/HOST=#  xx,yy,zz,rri,hti,rotx,roty,sld
402        ModifyTable format(Point)=1,width(Point)=0
403        RenameWindow #,T0
404        SetActiveSubwindow ##
405
406End
407
408
409Function KR_Plot1DButtonProc(ba) : ButtonControl
410        STRUCT WMButtonAction &ba
411
412        switch( ba.eventCode )
413                case 2: // mouse up                     
414                        DoWindow/F KR_IQ
415                        if(V_flag==0)
416                                Execute "KR_IQ()"
417                        Endif
418                       
419                        break
420        endswitch
421
422        return 0
423End
424
425Function KR_Show3DButtonProc(ba) : ButtonControl
426        STRUCT WMButtonAction &ba
427
428        switch( ba.eventCode )
429                case 2: // mouse up
430                        DoWindow/F Gizmo_VoxelMat
431                        if(V_flag==0)
432                                Execute "Gizmo_VoxelMat()"
433                        endif
434       
435                        break
436        endswitch
437
438        return 0
439End
440
441Function KR_DoCalcButtonProc(ba) : ButtonControl
442        STRUCT WMButtonAction &ba
443
444        switch( ba.eventCode )
445                case 2: // mouse up
446                        Wave xx=root:xx
447                        if(numpnts(xx)==0)
448                                return(0)
449                        endif
450                        wave yy=yy
451                        wave zz=zz
452                        wave rri=rri
453                        wave hti=hti
454//                      wave sbp=sbp
455                        wave rotx=rotx
456                        wave roty=roty
457                        wave sld=sld
458       
459                        Duplicate/O xx, sbp
460                        NVAR FFT_T=root:FFT_T
461                        sbp = FFT_T
462                       
463                        // parse
464                        KR_MultiCylinder(xx,yy,zz,rri,hti,sbp,rotx,roty,sld)
465
466                        // these are really just for display, or if the FFT of mat is wanted later.
467                        WAVE xoutW=root:xoutW
468                        WAVE youtW=root:youtW
469                        WAVE zoutW=root:zoutW
470                        WAVE sldW=root:sldW
471       
472                        XYZV_FillMat(xoutW,youtW,ZoutW,sldW,1)                  //last 1 will erase the matrix
473                        MakeTriplet(xoutW,youtW,zoutW)
474               
475                // and the calculation. Assumes SLDs are all the same   
476                        NVAR qmin = root:KR_Qmin
477                        NVAR qmax = root:KR_Qmax
478                        NVAR npt = root:KR_Npt
479                 
480                        fDoBinned_KR_FFTPanel(npt,qmin,qmax)
481       
482                       
483                        break
484        endswitch
485
486        return 0
487End
488
489
490Proc KR_IQ() : Graph
491        PauseUpdate; Silent 1           // building window...
492        Display /W=(295,44,627,302) ival_KR vs qval_KR
493        DoWindow/C KR_IQ
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/N=text0/J "\\s(ival_KR) ival_KR"
502EndMacro
503
504/////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.