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

Last change on this file since 940 was 940, checked in by srkline, 9 years ago

Updated the reduction Read/Write? and corresponding factility stubs to accomodate detector dead time written to the VAX header. It is currently not written to the header, but may be with NICE (hopefully).

With the move of NG3 SANS to CGB(upper), the NG3 designation in the account name [NGxSANSxx] has been replaced with CGB. Folders on charlotte and radio button on SASCALC are labeled "NGB30" to be more obvious which instrument it is.

FFT routines have minor typos cleaned up.
File size: 15.7 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        sbp = FFT_T
175       
176        variable npts,cyl
177        npts = numpnts(xx)
178        cyl = npts
179               
180        Make/O/D/N=0 xoutW,youtW,zoutW,sbpW,sldW
181       
182        PI2=pi*2
183        FTR=PI2/360
184       
185        nptW = 0
186       
187        for(l=0;l<(cyl);L+=1)   //only change from run4
188        //for each cylinder of loop use index NUM
189        //calculate x & y rotation cos and sin
190                STH=SIN(Rotx[L]*FTR)
191                SPH=sin(roty[L]*FTR)
192                CTH=cos(rotx[L]*FTR)
193                CPH=cos(roty[L]*FTR)
194                //print "sth",sth
195                //print"L=",L
196                P5=SBP[L]/2  //set sphere centers' half-diameter displacement from grid (avoids glitches)
197                // print "p5 & sbp[L]",p5,sbp[L]
198       
199                RR=(RRI[L]/SBP[L])//as an index, Igor truncates the number to an integer....does NOT round it
200                RR=RR+1 //rr is the loop limit for square around final circle
201                HH=(HTI[L]/(2*SBP[L]))  //as an index, Igor truncates the number to an integer....does NOT round it
202                for(k=-HH;k<HH;k+=1)  // should have +1 for HH to complete to k=HH?????
203                        for(i=-RR;i<RR;i+=1)  //should this have i<RR+1 or in above RR=RR+2????
204                                for(j=-RR;j<RR;J+=1)
205                                        x0=sbp*i+P5
206                                        y0=SBP*j+P5
207                                        z0=SBP*k+p5
208                                        if((((y0^2)/(RRI[L]^2))+((x0^2)/(RRI[L]^2)))<=1)
209                                                IX=-1
210                                        else
211                                                IX=0
212                                        endif
213                                        xmid=x0
214                                        ymid=y0*cth+z0*sth
215                                        zmid=-y0*sth+z0*cth
216                                        // end rotation about x begin rotn about y on rotated pts
217                                        //
218                                        xout=xmid*cph-zmid*sph
219                                        xout=xx[L]+xout
220                                        yout=ymid
221                                        yout=yy[L]+yout
222                                        zout=xmid*sph+zmid*cph
223                                        zout=zz[L]+zout
224
225                                        // now print to wave file the point or not depending on whether ix<0 or not
226
227                                        if (ix<0)
228                                        //write to wave file
229                                                InsertPoints nptW,1,xoutW,youtW,zoutW,sbpW,sldW
230                                                xoutW[nptW] = xout
231                                                youtW[nptW] = yout
232                                                zoutW[nptW] = zout
233                                                sbpW[nptW] = sbp[L]
234                                                sldW[nptW] = sld[L]
235                                               
236                                                nptW +=1
237                                       
238                                                //print  xout,yout,zout,sbp[L],sld[L]
239                                        //else
240                                                //continue
241                                        endif  //for write or not
242                                endfor  // for j
243                        endfor  //  for i
244                endfor  //for k
245        endfor // for L
246
247
248        // rescale to the sphere size
249        xoutW /= FFT_T
250        youtW /= FFT_T
251        zoutW /= FFT_T
252       
253        return(0) // end do loop cycle for cylinders
254end
255
256// triplet to display as a scatter plot in Gizmo
257//
258// will overwrite the triplet
259//
260Function MakeTriplet(xoutW,youtW,zoutW)
261        Wave xoutW,youtW,zoutW
262       
263        KillWaves/Z triplet
264        concatenate/O {xoutW,youtW,zoutW},triplet
265end
266
267
268
269Proc KR_LoadAndFill()
270
271        KR_Load()
272        XYZV_FillMat(xoutW,youtW,ZoutW,sldW,1)                  //last 1 will erase the matrix
273        MakeTriplet(xoutW,youtW,zoutW)
274         
275        DoBinned_KR_FFTPanel()
276        Print "now display the gizmo, triplet or use one of the calculation methods"
277       
278End
279
280Proc KR_CalcFromInput()
281
282        KR_MultiCylinder(xCtr,yCtr,zCtr,rad,length,sphereDiam,rot_x,rot_y,SLD_sph)
283        // these are really just for display, or if the FFT of mat is wanted later.
284        XYZV_FillMat(xoutW,youtW,ZoutW,sldW,1)                  //last 1 will erase the matrix
285        MakeTriplet(xoutW,youtW,zoutW)
286
287// and the calculation. Assumes SLDs are all the same   
288        DoBinned_KR_FFTPanel(100,0.004,0.5)
289       
290End
291
292//called from the FFT method panel
293//
294// in this method, the distances are binned as by Otto Glatter, and has been partially XOPed
295//
296// if the number of bins is too high (say 100000), then using the non-integer XYZ will
297// be 2-3 times slower since there will be a lot more bins - then the loop over the q-values at the
298// very end will be what is significantly slower. If the number of bins is reduced to 10000 (as suggested
299// in Otto's book, p.160), then the two methods (types 12 and 2) give very similar timing, and the
300// results are indistinguishable.
301//
302Proc DoBinned_KR_FFTPanel(num,qMin,qMax)
303        Variable num=100,qmin=0.004,qmax=0.5
304       
305        Variable t1
306        String qStr="qval_KR",iStr="ival_KR"            //default wave names, always overwritten
307        Variable grid
308       
309        grid=root:FFT_T
310
311        Make/O/D/N=(num) $qStr,$iStr
312        $qStr = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))         
313       
314        Variable estTime,nx
315        String str = ""
316
317        nx = NonZeroValues(mat)
318       
319        estTime = EstimatedTime(nx,num,2)               // 0 =  XOP, 1 = no XOP, 2 = binned distances
320        sprintf str, "Estimated time for the calculation is %g seconds. Proceed?",estTime
321        DoAlert 1,str
322        if(V_Flag==1)           //yes, proceed
323                t1=ticks
324                fDoCalc($qStr,$iStr,grid,12,1)
325//              Printf "Elapsed AltiSpheres time = %g seconds\r\r",(ticks-t1)/60.15
326        Endif
327End
328
329Function fDoBinned_KR_FFTPanel(num,qMin,qMax)
330        Variable num,qmin,qmax
331       
332        Variable t1,multiSLD,mode
333        String qStr="qval_KR",iStr="ival_KR"            //default wave names, always overwritten
334       
335        NVAR grid=root:FFT_T
336        ControlInfo/W=MultiCyl check_0
337        multiSLD = V_Value
338        if(multiSLD)
339                mode=13
340        else
341                mode=12
342        endif
343
344        Make/O/D/N=(num) qval_KR,ival_KR
345        qval_KR = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))               
346       
347        Variable estTime,nx,tooLong
348        String str = ""
349        tooLong = 300
350       
351        nx = NonZeroValues(mat)
352       
353        estTime = EstimatedTime(nx,num,2)               // 0 =  XOP, 1 = no XOP, 2 = binned distances
354        if(estTime > tooLong)
355                sprintf str, "Estimated time for the calculation is %g seconds. Proceed?",estTime
356                DoAlert 1,str
357        endif
358        if(V_Flag==1 || estTime < tooLong)              //yes, proceed
359                t1=ticks
360                fDoCalc(qval_KR,ival_KR,grid,mode,1)
361//              Printf "Elapsed AltiSpheres time = %g seconds\r\r",(ticks-t1)/60.15
362        Endif
363End
364
365/////////////////////////////////////
366// for each cylinder:
367//
368// xx,yy,zz,rri,hti,sbp,rotx,roty,sld
369// xx,yy,zz = center of cylinder
370// rri,hti = radius, height (units??)
371// sbp = ??? -- I think this is the diameter of the primary sphere
372// rotx, rotx = rotation angles (in degrees, but defined as ??)
373// sld = SLD of cylinder
374//
375//
376//
377// Put this into a panel with the table and the data
378// and fields for all of the inputs
379Macro Setup_KR_MultiCylinder()
380       
381        Make/O/D/N=0 xx,yy,zz,rri,hti,sbp,rotx,roty,sld
382        Variable/G root:KR_Qmin = 0.004
383        Variable/G root:KR_Qmax = 0.4
384        Variable/G root:KR_Npt = 100
385
386        FFT_MakeMatrixButtonProc("")
387
388        NewPanel /W=(241,44,1169,458)/N=MultiCyl/K=1 as "Multi-Cylinder"
389        ModifyPanel cbRGB=(49825,57306,65535)
390
391        Button button_0,pos={45,80},size={100,20},proc=KR_Show3DButtonProc,title="Show 3D"
392        Button button_1,pos={46,51},size={100,20},proc=KR_Plot1DButtonProc,title="Plot 1D"
393        Button button_2,pos={178,50},size={150,20},proc=KR_GenerateButtonProc,title="Generate Structure"
394        Button button_4,pos={178,80},size={120,20},proc=KR_DoCalcButtonProc,title="Do Calculation"
395        Button button_3,pos={600,60},size={120,20},proc=KR_DeleteRow,title="Delete Row(s)"
396        Button button_5,pos={600,10},size={120,20},proc=KR_SaveTable,title="Save Table"
397        Button button_6,pos={600,35},size={120,20},proc=KR_ImportTable,title="Import Table"
398        ValDisplay valdisp_0,pos={339,16},size={80,13},title="FFT_T"
399        ValDisplay valdisp_0,limits={0,0,0},barmisc={0,1000},value= #"root:FFT_T"
400        SetVariable setvar_0,pos={339,40},size={140,15},title="Q min (A)"
401        SetVariable setvar_0,limits={0,10,0},value= KR_Qmin
402        SetVariable setvar_1,pos={339,65},size={140,15},title="Q max (A)"
403        SetVariable setvar_1,limits={0,10,0},value= KR_Qmax
404        SetVariable setvar_2,pos={339,90},size={140,15},title="Num Pts"
405        SetVariable setvar_2,limits={10,500,0},value= KR_Npt
406        CheckBox check_0,pos={599,93},size={59,14},title="Multi SLD",value= 0
407
408        Edit/W=(18,117,889,378)/HOST=#  xx,yy,zz,rri,hti,rotx,roty,sld
409        ModifyTable format(Point)=1,width(Point)=0
410        RenameWindow #,T0
411        SetActiveSubwindow ##
412
413End
414
415
416Function KR_Plot1DButtonProc(ba) : ButtonControl
417        STRUCT WMButtonAction &ba
418
419        switch( ba.eventCode )
420                case 2: // mouse up                     
421                        DoWindow/F KR_IQ
422                        if(V_flag==0)
423                                Execute "KR_IQ()"
424                        Endif
425                       
426                        break
427        endswitch
428
429        return 0
430End
431
432Function KR_Show3DButtonProc(ba) : ButtonControl
433        STRUCT WMButtonAction &ba
434
435        switch( ba.eventCode )
436                case 2: // mouse up
437                        DoWindow/F Gizmo_VoxelMat
438                        if(V_flag==0)
439                                Execute "Gizmo_VoxelMat()"
440                        endif
441       
442                        break
443        endswitch
444
445        return 0
446End
447
448Function KR_DeleteRow(ba) : ButtonControl
449        STRUCT WMButtonAction &ba
450
451        switch( ba.eventCode )
452                case 2: // mouse up
453                       
454                        GetSelection table, MultiCyl#T0,3
455//                      Print V_flag, V_startRow, V_startCol, V_endRow, V_endCol
456                        DoAlert 1, "Do want to delete rows "+num2Str(V_StartRow)+" through "+num2str(V_endRow)+" ?"
457                        if(V_flag==1)                   
458                                DeletePoints V_StartRow,(V_endRow-V_StartRow+1),xx,yy,zz,rri,hti,sbp,rotx,roty,sld
459                        endif
460                       
461                        break
462        endswitch
463
464        return 0
465End
466
467Function KR_SaveTable(ba) : ButtonControl
468        STRUCT WMButtonAction &ba
469
470        switch( ba.eventCode )
471                case 2: // mouse up
472                       
473//                      xx,yy,zz,rri,hti,rotx,roty,sld
474                        Save/T/P=home/I xx,yy,zz,rri,hti,rotx,roty,sld as "SavedCyl.itx"
475                       
476                        break
477        endswitch
478
479        return 0
480End
481
482Function KR_ImportTable(ba) : ButtonControl
483        STRUCT WMButtonAction &ba
484
485        switch( ba.eventCode )
486                case 2: // mouse up
487                       
488                        LoadWave/T/O/P=home                     
489                        break
490        endswitch
491
492        return 0
493End
494
495//just generates the structure, no calculation
496Function KR_GenerateButtonProc(ba) : ButtonControl
497        STRUCT WMButtonAction &ba
498
499        switch( ba.eventCode )
500                case 2: // mouse up
501                        Wave xx=root:xx
502                        if(numpnts(xx)==0)
503                                return(0)
504                        endif
505                        wave yy=yy
506                        wave zz=zz
507                        wave rri=rri
508                        wave hti=hti
509//                      wave sbp=sbp
510                        wave rotx=rotx
511                        wave roty=roty
512                        wave sld=sld
513       
514                        Duplicate/O xx, sbp
515                        NVAR FFT_T=root:FFT_T
516                        sbp = FFT_T
517                       
518                        // parse
519                        KR_MultiCylinder(xx,yy,zz,rri,hti,sbp,rotx,roty,sld)
520
521                        // these are really just for display, or if the FFT of mat is wanted later.
522                        WAVE xoutW=root:xoutW
523                        WAVE youtW=root:youtW
524                        WAVE zoutW=root:zoutW
525                        WAVE sldW=root:sldW
526       
527                        XYZV_FillMat(xoutW,youtW,ZoutW,sldW,1)                  //last 1 will erase the matrix
528//                      MakeTriplet(xoutW,youtW,zoutW)
529//             
530//              // and the calculation. Assumes SLDs are all the same   
531//                      NVAR qmin = root:KR_Qmin
532//                      NVAR qmax = root:KR_Qmax
533//                      NVAR npt = root:KR_Npt
534//               
535//                      fDoBinned_KR_FFTPanel(npt,qmin,qmax)
536//     
537
538                        //force a redraw (re-coloring) of the gizmo window
539                        FFTMakeGizmoButtonProc("")
540                       
541                        break
542        endswitch
543
544
545       
546        return 0
547End
548
549
550
551Function KR_DoCalcButtonProc(ba) : ButtonControl
552        STRUCT WMButtonAction &ba
553
554        switch( ba.eventCode )
555                case 2: // mouse up
556                        Wave xx=root:xx
557                        if(numpnts(xx)==0)
558                                return(0)
559                        endif
560                        wave yy=yy
561                        wave zz=zz
562                        wave rri=rri
563                        wave hti=hti
564//                      wave sbp=sbp
565                        wave rotx=rotx
566                        wave roty=roty
567                        wave sld=sld
568       
569                        Duplicate/O xx, sbp
570                        NVAR FFT_T=root:FFT_T
571                        sbp = FFT_T
572                       
573                        // parse
574                        KR_MultiCylinder(xx,yy,zz,rri,hti,sbp,rotx,roty,sld)
575
576                        // these are really just for display, or if the FFT of mat is wanted later.
577                        WAVE xoutW=root:xoutW
578                        WAVE youtW=root:youtW
579                        WAVE zoutW=root:zoutW
580                        WAVE sldW=root:sldW
581       
582                        XYZV_FillMat(xoutW,youtW,ZoutW,sldW,1)                  //last 1 will erase the matrix
583                        MakeTriplet(xoutW,youtW,zoutW)
584               
585                // and the calculation. Assumes SLDs are all the same   
586                        NVAR qmin = root:KR_Qmin
587                        NVAR qmax = root:KR_Qmax
588                        NVAR npt = root:KR_Npt
589                 
590                        fDoBinned_KR_FFTPanel(npt,qmin,qmax)
591       
592                       
593                        break
594        endswitch
595
596        return 0
597End
598
599
600Proc KR_IQ() : Graph
601        PauseUpdate; Silent 1           // building window...
602        Display /W=(295,44,627,302) ival_KR vs qval_KR
603        DoWindow/C KR_IQ
604        ModifyGraph mode=4
605        ModifyGraph marker=19
606        ModifyGraph msize=2
607        ModifyGraph gaps=0
608        ModifyGraph grid=1
609        ModifyGraph log=1
610        ModifyGraph mirror=2
611        Legend/N=text0/J "\\s(ival_KR) ival_KR"
612EndMacro
613
614/////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.