source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_KR_Parser.ipf @ 800

Last change on this file since 800 was 800, checked in by srkline, 12 years ago

Changed the 2D resolution calculation to include gravity. This seems to work, but will still require some testing. It requires no modifications to the smearing calculation, still using parallel and perpendicular directions.

Added a Fake Update feature to the RealTime? update. There are specific, separate instructions for how to set this up. Usef (possibly) for summer schools or demos.

Adjusted the 2D MonteCarlo? simulation to a corrected beam center. The Gauss Peak was not symmetric around the old beam center.

Corrected the AAO resolution smearing functions to work with cursors.
File size: 13.8 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        FFT_MakeMatrixButtonProc("")
386
387        NewPanel /W=(241,44,1169,458)/N=MultiCyl/K=1 as "Multi-Cylinder"
388        ModifyPanel cbRGB=(49825,57306,65535)
389
390        Button button_0,pos={45,80},size={100,20},proc=KR_Show3DButtonProc,title="Show 3D"
391        Button button_1,pos={46,51},size={100,20},proc=KR_Plot1DButtonProc,title="Plot 1D"
392        Button button_2,pos={178,80},size={120,20},proc=KR_DoCalcButtonProc,title="Do Calculation"
393        Button button_3,pos={600,50},size={120,20},proc=KR_DeleteRow,title="Delete Row(s)"
394        ValDisplay valdisp_0,pos={339,16},size={80,13},title="FFT_T"
395        ValDisplay valdisp_0,limits={0,0,0},barmisc={0,1000},value= #"root:FFT_T"
396        SetVariable setvar_0,pos={339,40},size={140,15},title="Q min (A)"
397        SetVariable setvar_0,limits={0,10,0},value= KR_Qmin
398        SetVariable setvar_1,pos={339,65},size={140,15},title="Q max (A)"
399        SetVariable setvar_1,limits={0,10,0},value= KR_Qmax
400        SetVariable setvar_2,pos={339,90},size={140,15},title="Num Pts"
401        SetVariable setvar_2,limits={10,500,0},value= KR_Npt
402        CheckBox check_0,pos={599,93},size={59,14},title="Multi SLD",value= 0
403
404        Edit/W=(18,117,889,378)/HOST=#  xx,yy,zz,rri,hti,rotx,roty,sld
405        ModifyTable format(Point)=1,width(Point)=0
406        RenameWindow #,T0
407        SetActiveSubwindow ##
408
409End
410
411
412Function KR_Plot1DButtonProc(ba) : ButtonControl
413        STRUCT WMButtonAction &ba
414
415        switch( ba.eventCode )
416                case 2: // mouse up                     
417                        DoWindow/F KR_IQ
418                        if(V_flag==0)
419                                Execute "KR_IQ()"
420                        Endif
421                       
422                        break
423        endswitch
424
425        return 0
426End
427
428Function KR_Show3DButtonProc(ba) : ButtonControl
429        STRUCT WMButtonAction &ba
430
431        switch( ba.eventCode )
432                case 2: // mouse up
433                        DoWindow/F Gizmo_VoxelMat
434                        if(V_flag==0)
435                                Execute "Gizmo_VoxelMat()"
436                        endif
437       
438                        break
439        endswitch
440
441        return 0
442End
443
444Function KR_DeleteRow(ba) : ButtonControl
445        STRUCT WMButtonAction &ba
446
447        switch( ba.eventCode )
448                case 2: // mouse up
449                       
450                        GetSelection table, MultiCyl#T0,3
451//                      Print V_flag, V_startRow, V_startCol, V_endRow, V_endCol
452                        DoAlert 1, "Do want to delete rows "+num2Str(V_StartRow)+" through "+num2str(V_endRow)+" ?"
453                        if(V_flag==1)                   
454                                DeletePoints V_StartRow,(V_endRow-V_StartRow+1),xx,yy,zz,rri,hti,sbp,rotx,roty,sld
455                        endif
456                       
457                        break
458        endswitch
459
460        return 0
461End
462
463Function KR_DoCalcButtonProc(ba) : ButtonControl
464        STRUCT WMButtonAction &ba
465
466        switch( ba.eventCode )
467                case 2: // mouse up
468                        Wave xx=root:xx
469                        if(numpnts(xx)==0)
470                                return(0)
471                        endif
472                        wave yy=yy
473                        wave zz=zz
474                        wave rri=rri
475                        wave hti=hti
476//                      wave sbp=sbp
477                        wave rotx=rotx
478                        wave roty=roty
479                        wave sld=sld
480       
481                        Duplicate/O xx, sbp
482                        NVAR FFT_T=root:FFT_T
483                        sbp = FFT_T
484                       
485                        // parse
486                        KR_MultiCylinder(xx,yy,zz,rri,hti,sbp,rotx,roty,sld)
487
488                        // these are really just for display, or if the FFT of mat is wanted later.
489                        WAVE xoutW=root:xoutW
490                        WAVE youtW=root:youtW
491                        WAVE zoutW=root:zoutW
492                        WAVE sldW=root:sldW
493       
494                        XYZV_FillMat(xoutW,youtW,ZoutW,sldW,1)                  //last 1 will erase the matrix
495                        MakeTriplet(xoutW,youtW,zoutW)
496               
497                // and the calculation. Assumes SLDs are all the same   
498                        NVAR qmin = root:KR_Qmin
499                        NVAR qmax = root:KR_Qmax
500                        NVAR npt = root:KR_Npt
501                 
502                        fDoBinned_KR_FFTPanel(npt,qmin,qmax)
503       
504                       
505                        break
506        endswitch
507
508        return 0
509End
510
511
512Proc KR_IQ() : Graph
513        PauseUpdate; Silent 1           // building window...
514        Display /W=(295,44,627,302) ival_KR vs qval_KR
515        DoWindow/C KR_IQ
516        ModifyGraph mode=4
517        ModifyGraph marker=19
518        ModifyGraph msize=2
519        ModifyGraph gaps=0
520        ModifyGraph grid=1
521        ModifyGraph log=1
522        ModifyGraph mirror=2
523        Legend/N=text0/J "\\s(ival_KR) ival_KR"
524EndMacro
525
526/////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.