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

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

Added method to RealSpace? modeling to allow drawing a single cylidner with given center and arbitrary orientation.

Updated the Simulation help file to note that the simulation will work with a real space structure as the "model" input. In this way, there is no limitation on what can be simulated through SASCALC.
File size: 27.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3
4Proc testAddRotCyl()
5
6        Variable ii,zval,rotx,roty,fill,rad
7       
8        FFT_T=root:FFT_T
9       
10        rotx=0
11        roty=90
12        zval = 0
13        rad=20
14        fill=20
15       
16        ii=rad/FFT_T+1
17//      print ii
18        do
19                zval=ii
20                FFTDrawRotCylinder("mat",rad,200,64,64,zval,rotx,roty,fill)
21                ii+=(2*rad)/FFT_T
22                rotx += 12
23        while(ii<(127-(rad/FFT_T)))
24       
25//      print ii
26End
27
28
29// draws a single cylinder with the specified center and rotation. does not clear the matrix
30//
31// x-rotation is done first, then the y-rotation
32//
33Proc FFTDrawRotCylinder(matStr,rad,len,xc,yc,zc,xrot,yrot,fill)
34        String matStr="mat"
35        Variable rad=25,len=300,xc=50,yc=50,zc=50,xrot=10,yrot=90,fill=10
36        Prompt matStr,"the wave"                //,popup,WaveList("*",";","")
37        Prompt rad,"enter real radius (A)"
38        Prompt len,"enter length (A)"
39        Prompt xc,"enter the X-center"
40        Prompt yc,"enter the Y-center"
41        Prompt zc,"enter the Z-center"
42        Prompt xrot,"x-rotation (degrees)"
43        Prompt yrot,"y-rotation (degrees)"
44        Prompt fill,"fill SLD value"   
45       
46        Make/O/D/N=1 xx,yy,zz,rri,hti,rotx,roty,sld
47       
48        xx[0] = xc
49        yy[0] = yc
50        zz[0] = zc
51        rri[0] = rad
52        hti[0] = len
53        rotx[0] = xrot
54        roty[0] = yrot
55        sld[0] = fill   
56
57        Duplicate/O xx, sbp
58       
59        // parse
60//      KR_MultiCylinder(xx,yy,zz,rri,hti,sbp,rotx,roty,sld)
61        KR_MultiCylinder_Units(xx,yy,zz,rri,hti,sbp,rotx,roty,sld)
62
63        XYZV_FillMat_Centered(xoutW,youtW,ZoutW,xc,yc,zc,rad,len,sldW,0)                        //last 1 will erase the matrix, 0 retains matrix
64
65        //force a redraw (re-coloring) of the gizmo window
66        FFTMakeGizmoButtonProc("")
67       
68End
69
70
71// -- now this will put the cylinders properly at the specified centers, in the units of "mat"
72//
73/// seems to work - but what do I do about fractional positions? when converting to a matrix?
74//
75// the wave "gg" has been dropped, since it's only used as a flag in an old file loader
76//
77// NOW - SBP is FORCED to the value of FFT_T - no matter what is in the file.
78//
79Function KR_MultiCylinder_Units(xx,yy,zz,rri,hti,sbp,rotx,roty,sld)
80        Wave xx,yy,zz,rri,hti,sbp,rotx,roty,sld
81
82        Variable I, J, K, L, PT //integer indices loops, num cylinders, include or exclude sphere in circle
83        Variable STH, SPH, CTH, CPH, FTR  //sine and cosines and deg-->rad conversion: x rotn theta & y rotn phi
84        Variable  XMID, YMID, ZMID, XOUT, YOUT, ZOUT  //cartesian positions used in various calculations
85        Variable RR,HH  //RR is limit of loops, GG used as end of read param files--exit=2, NUM of cylinder
86        Variable  P5  //spheres half diameter shift from grid points (avoids zeros)
87        Variable X0, Y0,Z0
88        Variable PI2
89        Variable ix,nptW
90       
91        NVAR FFT_T = root:FFT_T
92//      FFT_T = sbp[0]
93//      sbp[0] = FFT_T
94        sbp = FFT_T
95       
96        variable npts,cyl
97        npts = numpnts(xx)
98        cyl = npts
99               
100        Make/O/D/N=0 xoutW,youtW,zoutW,sbpW,sldW
101       
102        PI2=pi*2
103        FTR=PI2/360
104       
105        nptW = 0
106       
107        for(l=0;l<(cyl);L+=1)   //only change from run4
108        //for each cylinder of loop use index NUM
109        //calculate x & y rotation cos and sin
110                STH=SIN(Rotx[L]*FTR)
111                SPH=sin(roty[L]*FTR)
112                CTH=cos(rotx[L]*FTR)
113                CPH=cos(roty[L]*FTR)
114                //print "sth",sth
115                //print"L=",L
116                P5=SBP[L]/2  //set sphere centers' half-diameter displacement from grid (avoids glitches)
117                // print "p5 & sbp[L]",p5,sbp[L]
118       
119                RR=(RRI[L]/SBP[L])//as an index, Igor truncates the number to an integer....does NOT round it
120                RR=RR+1 //rr is the loop limit for square around final circle
121                HH=(HTI[L]/(2*SBP[L]))  //as an index, Igor truncates the number to an integer....does NOT round it
122                for(k=-HH;k<HH;k+=1)  // should have +1 for HH to complete to k=HH?????
123                        for(i=-RR;i<RR;i+=1)  //should this have i<RR+1 or in above RR=RR+2????
124                                for(j=-RR;j<RR;J+=1)
125                                        x0=sbp[L]*i+P5
126                                        y0=SBP[L]*j+P5
127                                        z0=SBP[L]*k+p5
128                                        if((((y0^2)/(RRI[L]^2))+((x0^2)/(RRI[L]^2)))<=1)
129                                                IX=-1
130                                        else
131                                                IX=0
132                                        endif
133                                        xmid=x0
134                                        ymid=y0*cth+z0*sth
135                                        zmid=-y0*sth+z0*cth
136                                        // end rotation about x begin rotn about y on rotated pts
137                                        //
138                                        xout=xmid*cph-zmid*sph
139                                        xout=xx[L]+xout/SBP[L]
140                                        yout=ymid
141                                        yout=yy[L]+yout/SBP[L]
142                                        zout=xmid*sph+zmid*cph
143                                        zout=zz[L]+zout/SBP[L]
144
145                                        // now print to wave file the point or not depending on whether ix<0 or not
146
147                                        if (ix<0)
148                                        //write to wave file
149                                                InsertPoints nptW,1,xoutW,youtW,zoutW,sbpW,sldW
150                                                xoutW[nptW] = xout
151                                                youtW[nptW] = yout
152                                                zoutW[nptW] = zout
153                                                sbpW[nptW] = sbp[L]
154                                                sldW[nptW] = sld[L]
155                                               
156                                                nptW +=1
157                                       
158                                                //print  xout,yout,zout,sbp[L],sld[L]
159                                        //else
160                                                //continue
161                                        endif  //for write or not
162                                endfor  // for j
163                        endfor  //  for i
164                endfor  //for k
165        endfor // for L
166
167
168        // rescale to the sphere size
169//      xoutW /= FFT_T
170//      youtW /= FFT_T
171//      zoutW /= FFT_T
172
173        xoutW = trunc(xoutW)
174        youtW = trunc(youtW)
175        zoutW = trunc(zoutW)
176       
177        return(0) // end do loop cycle for cylinders
178end
179
180
181
182
183
184///
185/// -- replaced by KR_MultiCylinder_Units()
186///
187/// seems to work - but what do I do about fractional positions? when converting to a matrix?
188//
189// the wave "gg" has been dropped, since it's only used as a flag in an old file loader
190//
191// NOW - SBP is FORCED to the value of FFT_T - no matter what is in the file.
192//
193Function KR_MultiCylinder(xx,yy,zz,rri,hti,sbp,rotx,roty,sld)
194        Wave xx,yy,zz,rri,hti,sbp,rotx,roty,sld
195
196        Variable I, J, K, L, PT //integer indices loops, num cylinders, include or exclude sphere in circle
197        Variable STH, SPH, CTH, CPH, FTR  //sine and cosines and deg-->rad conversion: x rotn theta & y rotn phi
198        Variable  XMID, YMID, ZMID, XOUT, YOUT, ZOUT  //cartesian positions used in various calculations
199        Variable RR,HH  //RR is limit of loops, GG used as end of read param files--exit=2, NUM of cylinder
200        Variable  P5  //spheres half diameter shift from grid points (avoids zeros)
201        Variable X0, Y0,Z0
202        Variable PI2
203        Variable ix,nptW
204       
205        NVAR FFT_T = root:FFT_T
206//      FFT_T = sbp[0]
207//      sbp[0] = FFT_T
208        sbp = FFT_T
209       
210        variable npts,cyl
211        npts = numpnts(xx)
212        cyl = npts
213               
214        Make/O/D/N=0 xoutW,youtW,zoutW,sbpW,sldW
215       
216        PI2=pi*2
217        FTR=PI2/360
218       
219        nptW = 0
220       
221        for(l=0;l<(cyl);L+=1)   //only change from run4
222        //for each cylinder of loop use index NUM
223        //calculate x & y rotation cos and sin
224                STH=SIN(Rotx[L]*FTR)
225                SPH=sin(roty[L]*FTR)
226                CTH=cos(rotx[L]*FTR)
227                CPH=cos(roty[L]*FTR)
228                //print "sth",sth
229                //print"L=",L
230                P5=SBP[L]/2  //set sphere centers' half-diameter displacement from grid (avoids glitches)
231                // print "p5 & sbp[L]",p5,sbp[L]
232       
233                RR=(RRI[L]/SBP[L])//as an index, Igor truncates the number to an integer....does NOT round it
234                RR=RR+1 //rr is the loop limit for square around final circle
235                HH=(HTI[L]/(2*SBP[L]))  //as an index, Igor truncates the number to an integer....does NOT round it
236                for(k=-HH;k<HH;k+=1)  // should have +1 for HH to complete to k=HH?????
237                        for(i=-RR;i<RR;i+=1)  //should this have i<RR+1 or in above RR=RR+2????
238                                for(j=-RR;j<RR;J+=1)
239                                        x0=sbp*i+P5
240                                        y0=SBP*j+P5
241                                        z0=SBP*k+p5
242                                        if((((y0^2)/(RRI[L]^2))+((x0^2)/(RRI[L]^2)))<=1)
243                                                IX=-1
244                                        else
245                                                IX=0
246                                        endif
247                                        xmid=x0
248                                        ymid=y0*cth+z0*sth
249                                        zmid=-y0*sth+z0*cth
250                                        // end rotation about x begin rotn about y on rotated pts
251                                        //
252                                        xout=xmid*cph-zmid*sph
253                                        xout=xx[L]+xout
254                                        yout=ymid
255                                        yout=yy[L]+yout
256                                        zout=xmid*sph+zmid*cph
257                                        zout=zz[L]+zout
258
259                                        // now print to wave file the point or not depending on whether ix<0 or not
260
261                                        if (ix<0)
262                                        //write to wave file
263                                                InsertPoints nptW,1,xoutW,youtW,zoutW,sbpW,sldW
264                                                xoutW[nptW] = xout
265                                                youtW[nptW] = yout
266                                                zoutW[nptW] = zout
267                                                sbpW[nptW] = sbp[L]
268                                                sldW[nptW] = sld[L]
269                                               
270                                                nptW +=1
271                                       
272                                                //print  xout,yout,zout,sbp[L],sld[L]
273                                        //else
274                                                //continue
275                                        endif  //for write or not
276                                endfor  // for j
277                        endfor  //  for i
278                endfor  //for k
279        endfor // for L
280
281
282        // rescale to the sphere size
283        xoutW /= FFT_T
284        youtW /= FFT_T
285        zoutW /= FFT_T
286       
287        return(0) // end do loop cycle for cylinders
288end
289
290
291
292
293
294
295// triplet to display as a scatter plot in Gizmo
296//
297// will overwrite the triplet
298//
299Function MakeTriplet(xoutW,youtW,zoutW)
300        Wave xoutW,youtW,zoutW
301       
302        KillWaves/Z triplet
303        concatenate/O {xoutW,youtW,zoutW},triplet
304end
305
306
307
308Proc KR_LoadAndFill()
309
310        KR_Load()
311        XYZV_FillMat(xoutW,youtW,ZoutW,sldW,1)                  //last 1 will erase the matrix
312        MakeTriplet(xoutW,youtW,zoutW)
313         
314        DoBinned_KR_FFTPanel()
315        Print "now display the gizmo, triplet or use one of the calculation methods"
316       
317End
318
319Proc KR_CalcFromInput()
320
321        KR_MultiCylinder(xCtr,yCtr,zCtr,rad,length,sphereDiam,rot_x,rot_y,SLD_sph)
322        // these are really just for display, or if the FFT of mat is wanted later.
323        XYZV_FillMat(xoutW,youtW,ZoutW,sldW,1)                  //last 1 will erase the matrix
324        MakeTriplet(xoutW,youtW,zoutW)
325
326// and the calculation. Assumes SLDs are all the same   
327        DoBinned_KR_FFTPanel(100,0.004,0.5)
328       
329End
330
331//called from the FFT method panel
332//
333// in this method, the distances are binned as by Otto Glatter, and has been partially XOPed
334//
335// if the number of bins is too high (say 100000), then using the non-integer XYZ will
336// be 2-3 times slower since there will be a lot more bins - then the loop over the q-values at the
337// very end will be what is significantly slower. If the number of bins is reduced to 10000 (as suggested
338// in Otto's book, p.160), then the two methods (types 12 and 2) give very similar timing, and the
339// results are indistinguishable.
340//
341Proc DoBinned_KR_FFTPanel(num,qMin,qMax)
342        Variable num=100,qmin=0.004,qmax=0.5
343       
344        Variable t1
345        String qStr="qval_KR",iStr="ival_KR"            //default wave names, always overwritten
346        Variable grid
347       
348        grid=root:FFT_T
349
350        Make/O/D/N=(num) $qStr,$iStr
351        $qStr = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))         
352       
353        Variable estTime,nx
354        String str = ""
355
356        nx = NonZeroValues(mat)
357       
358        estTime = EstimatedTime(nx,num,2)               // 0 =  XOP, 1 = no XOP, 2 = binned distances
359        sprintf str, "Estimated time for the calculation is %g seconds. Proceed?",estTime
360        DoAlert 1,str
361        if(V_Flag==1)           //yes, proceed
362                t1=ticks
363                fDoCalc($qStr,$iStr,grid,12,1)
364//              Printf "Elapsed AltiSpheres time = %g seconds\r\r",(ticks-t1)/60.15
365        Endif
366End
367
368Function fDoBinned_KR_FFTPanel(num,qMin,qMax)
369        Variable num,qmin,qmax
370       
371        Variable t1,multiSLD,mode
372        String qStr="qval_KR",iStr="ival_KR"            //default wave names, always overwritten
373       
374        NVAR grid=root:FFT_T
375        ControlInfo/W=MultiCyl check_0
376        multiSLD = V_Value
377        if(multiSLD)
378                mode=13
379        else
380                mode=12
381        endif
382
383        Make/O/D/N=(num) qval_KR,ival_KR
384        qval_KR = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))               
385       
386        Variable estTime,nx,tooLong
387        String str = ""
388        tooLong = 300
389       
390        nx = NonZeroValues(mat)
391       
392        estTime = EstimatedTime(nx,num,2)               // 0 =  XOP, 1 = no XOP, 2 = binned distances
393        if(estTime > tooLong)
394                sprintf str, "Estimated time for the calculation is %g seconds. Proceed?",estTime
395                DoAlert 1,str
396        endif
397        if(V_Flag==1 || estTime < tooLong)              //yes, proceed
398                t1=ticks
399                fDoCalc(qval_KR,ival_KR,grid,mode,1)
400//              Printf "Elapsed AltiSpheres time = %g seconds\r\r",(ticks-t1)/60.15
401        Endif
402End
403
404/////////////////////////////////////
405// for each cylinder:
406//
407// xx,yy,zz,rri,hti,sbp,rotx,roty,sld
408// xx,yy,zz = center of cylinder
409// rri,hti = radius, height (units??)
410// sbp = ??? -- I think this is the diameter of the primary sphere
411// rotx, rotx = rotation angles (in degrees, but defined as ??)
412// sld = SLD of cylinder
413//
414// Put this into a panel with the table and the data
415// and fields for all of the inputs
416Macro Setup_KR_MultiCylinder()
417       
418        Make/O/D/N=0 xx,yy,zz,rri,hti,sbp,rotx,roty,sld
419        Variable/G root:KR_Qmin = 0.004
420        Variable/G root:KR_Qmax = 0.4
421        Variable/G root:KR_Npt = 100
422
423        FFT_MakeMatrixButtonProc("")
424
425        NewPanel /W=(241,44,1169,458)/N=MultiCyl/K=1 as "Multi-Cylinder"
426        ModifyPanel cbRGB=(49825,57306,65535)
427
428        Button button_0,pos={45,80},size={100,20},proc=KR_Show3DButtonProc,title="Show 3D"
429        Button button_1,pos={46,51},size={100,20},proc=KR_Plot1DButtonProc,title="Plot 1D"
430        Button button_2,pos={178,50},size={150,20},proc=KR_GenerateButtonProc,title="Generate Structure"
431        Button button_4,pos={178,80},size={120,20},proc=KR_DoCalcButtonProc,title="Do Calculation"
432        Button button_3,pos={600,60},size={120,20},proc=KR_DeleteRow,title="Delete Row(s)"
433        Button button_5,pos={600,10},size={120,20},proc=KR_SaveTable,title="Save Table"
434        Button button_6,pos={600,35},size={120,20},proc=KR_ImportTable,title="Import Table"
435        ValDisplay valdisp_0,pos={339,16},size={80,13},title="FFT_T"
436        ValDisplay valdisp_0,limits={0,0,0},barmisc={0,1000},value= #"root:FFT_T"
437        SetVariable setvar_0,pos={339,40},size={140,15},title="Q min (A)"
438        SetVariable setvar_0,limits={0,10,0},value= KR_Qmin
439        SetVariable setvar_1,pos={339,65},size={140,15},title="Q max (A)"
440        SetVariable setvar_1,limits={0,10,0},value= KR_Qmax
441        SetVariable setvar_2,pos={339,90},size={140,15},title="Num Pts"
442        SetVariable setvar_2,limits={10,500,0},value= KR_Npt
443        CheckBox check_0,pos={599,93},size={59,14},title="Multi SLD",value= 0
444
445        Edit/W=(18,117,889,378)/HOST=#  xx,yy,zz,rri,hti,rotx,roty,sld
446        ModifyTable format(Point)=1,width(Point)=0
447        RenameWindow #,T0
448        SetActiveSubwindow ##
449
450End
451
452
453
454
455Function KR_Plot1DButtonProc(ba) : ButtonControl
456        STRUCT WMButtonAction &ba
457
458        switch( ba.eventCode )
459                case 2: // mouse up                     
460                        DoWindow/F KR_IQ
461                        if(V_flag==0)
462                                Execute "KR_IQ()"
463                        Endif
464                       
465                        break
466        endswitch
467
468        return 0
469End
470
471Function KR_Show3DButtonProc(ba) : ButtonControl
472        STRUCT WMButtonAction &ba
473
474        switch( ba.eventCode )
475                case 2: // mouse up
476                        DoWindow/F Gizmo_VoxelMat
477                        if(V_flag==0)
478                                Execute "Gizmo_VoxelMat()"
479                        endif
480       
481                        break
482        endswitch
483
484        return 0
485End
486
487Function KR_DeleteRow(ba) : ButtonControl
488        STRUCT WMButtonAction &ba
489
490        switch( ba.eventCode )
491                case 2: // mouse up
492                       
493                        GetSelection table, MultiCyl#T0,3
494//                      Print V_flag, V_startRow, V_startCol, V_endRow, V_endCol
495                        DoAlert 1, "Do want to delete rows "+num2Str(V_StartRow)+" through "+num2str(V_endRow)+" ?"
496                        if(V_flag==1)                   
497                                DeletePoints V_StartRow,(V_endRow-V_StartRow+1),xx,yy,zz,rri,hti,sbp,rotx,roty,sld
498                        endif
499                       
500                        break
501        endswitch
502
503        return 0
504End
505
506Function KR_SaveTable(ba) : ButtonControl
507        STRUCT WMButtonAction &ba
508
509        switch( ba.eventCode )
510                case 2: // mouse up
511                       
512//                      xx,yy,zz,rri,hti,rotx,roty,sld
513                        Save/T/P=home/I xx,yy,zz,rri,hti,rotx,roty,sld as "SavedCyl.itx"
514                       
515                        break
516        endswitch
517
518        return 0
519End
520
521Function KR_ImportTable(ba) : ButtonControl
522        STRUCT WMButtonAction &ba
523
524        switch( ba.eventCode )
525                case 2: // mouse up
526                       
527                        LoadWave/T/O/P=home                     
528                        break
529        endswitch
530
531        return 0
532End
533
534//just generates the structure, no calculation
535Function KR_GenerateButtonProc(ba) : ButtonControl
536        STRUCT WMButtonAction &ba
537
538        switch( ba.eventCode )
539                case 2: // mouse up
540                        Wave xx=root:xx
541                        if(numpnts(xx)==0)
542                                return(0)
543                        endif
544                        wave yy=yy
545                        wave zz=zz
546                        wave rri=rri
547                        wave hti=hti
548//                      wave sbp=sbp
549                        wave rotx=rotx
550                        wave roty=roty
551                        wave sld=sld
552       
553                        Duplicate/O xx, sbp
554                        NVAR FFT_T=root:FFT_T
555                        sbp = FFT_T
556                       
557                        // parse
558                        KR_MultiCylinder(xx,yy,zz,rri,hti,sbp,rotx,roty,sld)
559
560                        // these are really just for display, or if the FFT of mat is wanted later.
561                        WAVE xoutW=root:xoutW
562                        WAVE youtW=root:youtW
563                        WAVE zoutW=root:zoutW
564                        WAVE sldW=root:sldW
565       
566                        XYZV_FillMat(xoutW,youtW,ZoutW,sldW,1)                  //last 1 will erase the matrix
567//                      MakeTriplet(xoutW,youtW,zoutW)
568//             
569//              // and the calculation. Assumes SLDs are all the same   
570//                      NVAR qmin = root:KR_Qmin
571//                      NVAR qmax = root:KR_Qmax
572//                      NVAR npt = root:KR_Npt
573//               
574//                      fDoBinned_KR_FFTPanel(npt,qmin,qmax)
575//     
576
577                        //force a redraw (re-coloring) of the gizmo window
578                        FFTMakeGizmoButtonProc("")
579                       
580                        break
581        endswitch
582
583
584       
585        return 0
586End
587
588
589
590Function KR_DoCalcButtonProc(ba) : ButtonControl
591        STRUCT WMButtonAction &ba
592
593        switch( ba.eventCode )
594                case 2: // mouse up
595                        Wave xx=root:xx
596                        if(numpnts(xx)==0)
597                                return(0)
598                        endif
599                        wave yy=yy
600                        wave zz=zz
601                        wave rri=rri
602                        wave hti=hti
603//                      wave sbp=sbp
604                        wave rotx=rotx
605                        wave roty=roty
606                        wave sld=sld
607       
608                        Duplicate/O xx, sbp
609                        NVAR FFT_T=root:FFT_T
610                        sbp = FFT_T
611                       
612                        // parse
613                        KR_MultiCylinder(xx,yy,zz,rri,hti,sbp,rotx,roty,sld)
614
615                        // these are really just for display, or if the FFT of mat is wanted later.
616                        WAVE xoutW=root:xoutW
617                        WAVE youtW=root:youtW
618                        WAVE zoutW=root:zoutW
619                        WAVE sldW=root:sldW
620       
621                        XYZV_FillMat(xoutW,youtW,ZoutW,sldW,1)                  //last 1 will erase the matrix
622                        MakeTriplet(xoutW,youtW,zoutW)
623               
624                // and the calculation. Assumes SLDs are all the same   
625                        NVAR qmin = root:KR_Qmin
626                        NVAR qmax = root:KR_Qmax
627                        NVAR npt = root:KR_Npt
628                 
629                        fDoBinned_KR_FFTPanel(npt,qmin,qmax)
630       
631                       
632                        break
633        endswitch
634
635        return 0
636End
637
638
639Proc KR_IQ() : Graph
640        PauseUpdate; Silent 1           // building window...
641        Display /W=(295,44,627,302) ival_KR vs qval_KR
642        DoWindow/C KR_IQ
643        ModifyGraph mode=4
644        ModifyGraph marker=19
645        ModifyGraph msize=2
646        ModifyGraph gaps=0
647        ModifyGraph grid=1
648        ModifyGraph log=1
649        ModifyGraph mirror=2
650        Legend/N=text0/J "\\s(ival_KR) ival_KR"
651EndMacro
652
653/////////////////////////////////////
654//
655//
656// for the "manual" fitting
657//
658
659Proc Vary_One_Cyl_Param(waveStr,row,percent,numStep)
660// pick the wave and row and %
661        String waveStr="xx"
662        Variable row=1
663        Variable percent = 105
664        Variable numStep = 10
665        Prompt waveStr,"wave",popup,"xx;yy;zz;rri;hti;rotx;roty;sld;"
666
667        print waveStr
668       
669// dispatch to calculation
670        MultiCyl_Loop($waveStr,row,percent,numStep)
671       
672// plot the chi2_map
673        DoWindow/F MultiCyl_ChiMap
674        if(V_flag==0)
675                MultiCyl_ChiMap()
676        endif
677
678end
679
680Function MultiCyl_Loop(w,row,percent,numStep)
681        Wave w
682        Variable row,percent,numStep
683       
684        Variable loLim,hiLim,ii,minIndex,minChiSq
685        String folderStr
686       
687        Make/O/D/N=(numStep) chi2_map,testVals
688        Wave chi2_map=chi2_map
689        Wave testVals=testVals
690       
691        loLim = w[row] - percent*w[row]/100
692        hiLim = w[row] + percent*w[row]/100
693        testVals = loLim + x*(hiLim-loLim)/numStep
694
695//              the experimental data
696        ControlInfo/W=WrapperPanel popup_0
697        folderStr=S_Value
698        // wave references for the data (to pass)
699        String DF="root:"+folderStr+":"
700       
701        WAVE yw=$(DF+folderStr+"_i")
702        WAVE xw=$(DF+folderStr+"_q")
703        WAVE sw=$(DF+folderStr+"_s")
704
705        duplicate/o yw, interpCalc,chi2_data
706        Wave chi2_data=chi2_data
707        Wave interpCalc=interpCalc
708               
709        STRUCT WMButtonAction ba
710        ba.eventCode = 2
711       
712// loop
713        for(ii=0;ii<numStep;ii+=1)
714//      set the value
715                w[row] = testVals[ii]
716//              generate the structure
717//
718                KR_GenerateButtonProc(ba)
719
720//              do the calculation
721                KR_DoCalcButtonProc(ba)
722
723                WAVE ival_KR=ival_KR
724                WAVE qval_KR=qval_KR
725               
726                interpCalc = interp(xw, qval_KR, ival_KR )
727               
728//              calculate chi-squared
729                chi2_data = (yw-interpCalc)^2
730                chi2_data /= sw^2
731       
732
733                Wavestats/Q chi2_data
734                chi2_map[ii] = V_avg * V_npnts
735// end loop
736        endfor
737
738// find the best chi squared
739        WaveStats/Q chi2_map
740// reset the value to the best
741        minIndex = V_minRowLoc
742        w[row] = testVals[minIndex]
743       
744        minChiSq = chi2_map[minIndex]
745        print "Minimum chi2 = ",minChiSq
746
747
748// and then recalculate at the best solution
749        KR_GenerateButtonProc(ba)
750
751//              do the calculation
752        KR_DoCalcButtonProc(ba)
753       
754        return(0)
755End
756
757Proc MultiCyl_ChiMap()
758        PauseUpdate; Silent 1           // building window...
759        Display /W=(35,44,466,414) chi2_map vs testVals
760        DoWindow/C MultiCyl_ChiMap
761        ModifyGraph mode=4
762        ModifyGraph marker=19
763        ModifyGraph msize=2
764        Label left "chi-squared"
765        Label bottom "test values"
766end
767
768//Function testKRPar()
769//     
770//      Variable row, col
771//      String wStr
772//     
773//      getParamFromKRSetup(row,col,wStr)
774//      Print row, col, wStr
775//     
776//      wStr = StringFromList(0, wStr)          // some wave "xx.d"
777//      wStr = StringFromList(0, wStr, ".")  // removes the ".d"
778//      Wave w=$wStr
779//      print w[row]
780//     
781//      Variable numStep,loLim,hiLim,percent
782//      numStep = 25
783//      percent = 50
784//     
785//      loLim = w[row] - percent*w[row]/100
786//      hiLim = w[row] + percent*w[row]/100
787//     
788//     
789//      Make/O/D/N=(numStep) testKRVals
790//      testKRVals = loLim + x*(hiLim-loLim)/numStep   
791//
792//      print testKRvals
793//      return(0)
794//     
795//End
796//
797//Function getParamFromKRSetup(row,col,wStr)
798//      Variable &row,&col
799//      String &wStr
800//     
801//      Variable parNum
802//     
803//      GetSelection table, MultiCyl#T0, 3
804//      row = V_startRow
805//      col = V_startCol
806//      Print S_Selection
807//      wStr = S_Selection
808//     
809//     
810//      return(0)
811//End
812
813Proc Vary_Two_Cyl_Param(waveStr,row,waveStr2,row2,percent,percent2,numStep)
814// pick the wave and row and %
815        String waveStr="xx"
816        Variable row=1
817        String waveStr2="rri"
818        Variable row2=0
819        Variable percent = 105
820        Variable percent2 = 50
821        Variable numStep=5
822        Prompt waveStr,"wave",popup,"xx;yy;zz;rri;hti;rotx;roty;sld;"
823        Prompt waveStr2,"wave2",popup,"xx;yy;zz;rri;hti;rotx;roty;sld;"
824       
825// dispatch to calculation
826        MultiCyl_Loop_2D($waveStr,row,$waveStr2,row2,percent,percent2,numStep)
827       
828// plot the chi2_map
829        DoWindow/F MultiCyl_ChiMap_2D
830        if(V_flag==0)
831                MultiCyl_ChiMap_2D()
832        else
833                //V_min*1.01 = the 1% neighborhood around the solution
834                WaveStats/Q chi2_Map_2D
835                ModifyImage/W=MultiCyl_ChiMap_2D chi2_Map_2D ctab= {(V_min*1.01),*,ColdWarm,0}
836                ModifyImage/W=MultiCyl_ChiMap_2D chi2_Map_2D minRGB=(0,65535,0),maxRGB=(0,65535,0)
837        endif
838
839end
840
841
842Function MultiCyl_Loop_2D(w,row,w2,row2,percent,percent2,numStep)
843        Wave w
844        Variable row
845        Wave w2
846        Variable row2,percent,percent2,numStep
847       
848        Variable loLim,hiLim,ii,jj,minIndex,minChiSq
849        String folderStr
850       
851        Make/O/D/N=(numStep,numStep) chi2_Map_2D
852        Make/O/D/N=(numStep) testVals,testVals2
853        Wave chi2_Map_2D=chi2_Map_2D
854        Wave testVals=testVals
855        Wave testVals2=testVals2
856       
857        testVals = 0
858        testVals2 = 0
859        chi2_Map_2D = 0
860       
861       
862        loLim = w[row] - percent*w[row]/100
863        hiLim = w[row] + percent*w[row]/100
864        testVals = loLim + x*(hiLim-loLim)/(numStep-1)
865//      Print lolim,hilim
866
867        SetScale/I x LoLim,HiLim,"", chi2_Map_2D
868
869        loLim = w2[row2] - percent2*w2[row2]/100
870        hiLim = w2[row2] + percent2*w2[row2]/100
871        testVals2 = loLim + x*(hiLim-loLim)/(numStep-1)
872//      Print lolim,hilim
873
874        SetScale/I y LoLim,HiLim,"", chi2_Map_2D
875
876
877//              the experimental data
878        ControlInfo/W=WrapperPanel popup_0
879        folderStr=S_Value
880        // wave references for the data (to pass)
881        String DF="root:"+folderStr+":"
882       
883        WAVE yw=$(DF+folderStr+"_i")
884        WAVE xw=$(DF+folderStr+"_q")
885        WAVE sw=$(DF+folderStr+"_s")
886
887        duplicate/o yw, interpCalc,chi2_data
888        Wave chi2_data=chi2_data
889        Wave interpCalc=interpCalc
890       
891        STRUCT WMButtonAction ba
892        ba.eventCode = 2
893       
894// double loop
895        for(ii=0;ii<numStep;ii+=1)
896                Print "                 Outer Loop Index = ",ii," out of ",numStep
897                //set the value from the outer loop
898                w[row] = testVals[ii]
899
900                for(jj=0;jj<numStep;jj+=1)
901               
902        //      set the inner value
903                        w2[row2] = testVals2[jj]
904        //              generate the structure
905        //
906                        KR_GenerateButtonProc(ba)
907       
908        //              do the calculation
909                        KR_DoCalcButtonProc(ba)
910       
911                        WAVE ival_KR=ival_KR
912                        WAVE qval_KR=qval_KR
913                       
914                        interpCalc = interp(xw, qval_KR, ival_KR )
915                       
916        //              calculate chi-squared
917                        chi2_data = (yw-interpCalc)^2
918                        chi2_data /= sw^2
919               
920                        Wavestats/Q chi2_data
921                        chi2_Map_2D[ii][jj] = V_avg * V_npnts
922                       
923                endfor
924        endfor
925
926
927// find the best chi squared
928        WaveStats/Q chi2_Map_2D
929// reset the value to the best
930        w[row] = V_MinRowLoc
931        w2[row2] = V_MinColLoc 
932       
933        minChiSq = V_Min
934        print "Minimum chi2 = ",minChiSq
935
936// and then recalculate at the best solution
937        KR_GenerateButtonProc(ba)
938
939//              do the calculation
940        KR_DoCalcButtonProc(ba)
941       
942        return(0)
943End
944
945Proc MultiCyl_ChiMap_2D()
946        PauseUpdate; Silent 1           // building window...
947        Display /W=(35,44,466,414)
948        AppendImage chi2_Map_2D
949        DoWindow/C MultiCyl_ChiMap_2D
950        ModifyImage chi2_Map_2D ctab= {*,*,ColdWarm,0}
951       
952        //V_min*1.01 = the 1% neighborhood around the solution
953        WaveStats/Q chi2_Map_2D
954        ModifyImage/W=MultiCyl_ChiMap_2D chi2_Map_2D ctab= {(V_min*1.01),*,ColdWarm,0}
955        ModifyImage/W=MultiCyl_ChiMap_2D chi2_Map_2D minRGB=(0,65535,0),maxRGB=(0,65535,0)
956       
957        Label bottom "test values"
958        Label left "test values"
959end
960
961
962
963//
964///// seems to work - but what do I do about fractional positions? when converting to a matrix?
965////
966////
967//
968//Function KR_Load()
969//      Variable I, J, K, L, PT //integer indices loops, num cylinders, include or exclude sphere in circle
970//      Variable STH, SPH, CTH, CPH, FTR  //sine and cosines and deg-->rad conversion: x rotn theta & y rotn phi
971//      Variable  XMID, YMID, ZMID, XOUT, YOUT, ZOUT  //cartesian positions used in various calculations
972//      Variable RR,HH  //RR is limit of loops, GG used as end of read param files--exit=2, NUM of cylinder
973//      Variable  P5  //spheres half diameter shift from grid points (avoids zeros)
974//      Variable X0, Y0,Z0
975//      Variable PI2
976//      Variable ix,nptW
977//
978//
979//      LoadWave /G /N 
980//      Print S_filename
981//      Print S_wavenames
982//     
983//      //Make / O /N=0 OutputPoints
984//      //      wave out=OutputPoints
985//      //      variable num=numpnts(out)
986//     
987//      KillWaves/Z xx,yy,zz,rri,hti,sbp,rotx,roty,sld,gg
988//     
989//      Rename wave0, xx
990//      Rename wave1, yy
991//      Rename wave2, zz
992//      Rename wave3, RRI
993//      Rename wave4, HTI
994//      Rename wave5, SBP
995//      Rename wave6, ROTX
996//      Rename wave7, ROTY
997//      Rename wave8, SLD
998//      Rename wave9, GG
999//     
1000//      //print  NUM,xx,yy,zz,rri,hti,sbp,rotx,roty,sld,gg
1001//
1002//     
1003//      wave gg = gg
1004//      variable nn =-1,npts,cyl
1005//      npts = numpnts(GG)
1006//     
1007//      for (i=0;i<=npts;i+=1)
1008//              if (gg[i]==2)
1009//                      cyl = i+1
1010//                      break
1011//                      print "gg[i],i=",gg,i
1012//              endif
1013//      endfor
1014//      print"cyl=",cyl
1015//     
1016//
1017//      wave xx=xx
1018//      wave yy=yy
1019//      wave zz=zz
1020//      wave rri=rri
1021//      wave hti=hti
1022//      wave sbp=sbp
1023//      wave rotx=rotx
1024//      wave roty=roty
1025//      wave sld=sld
1026//     
1027//      // SBP = diameter of the spheres
1028//      NVAR FFT_T = root:FFT_T
1029//      FFT_T = SBP[0]
1030//      //
1031//     
1032//      Make/O/D/N=0 xoutW,youtW,zoutW,sbpW,sldW
1033//     
1034//      PI2=pi*2
1035//      FTR=PI2/360
1036//      // print "ftr=", ftr
1037//     
1038//      nptW = 0
1039//     
1040//      for(l=0;l<(cyl);L+=1)   //only change from run4
1041//      //for each cylinder of loop use index NUM
1042//      //calculate x & y rotation cos and sin
1043//              STH=SIN(Rotx[L]*FTR)
1044//              SPH=sin(roty[L]*FTR)
1045//              CTH=cos(rotx[L]*FTR)
1046//              CPH=cos(roty[L]*FTR)
1047//              //print "sth",sth
1048//              //print"L=",L
1049//              P5=SBP[L]/2  //set sphere centers' half-diameter displacement from grid (avoids glitches)
1050//              // print "p5 & sbp[L]",p5,sbp[L]
1051//     
1052//              RR=(RRI[L]/SBP[L])//as an index, Igor truncates the number to an integer....does NOT round it
1053//              RR=RR+1 //rr is the loop limit for square around final circle
1054//              HH=(HTI[L]/(2*SBP[L]))  //as an index, Igor truncates the number to an integer....does NOT round it
1055//              for(k=-HH;k<HH;k+=1)  // should have +1 for HH to complete to k=HH?????
1056//                      for(i=-RR;i<RR;i+=1)  //should this have i<RR+1 or in above RR=RR+2????
1057//                              for(j=-RR;j<RR;J+=1)
1058//                                      x0=sbp*i+P5
1059//                                      y0=SBP*j+P5
1060//                                      z0=SBP*k+p5
1061//                                      if((((y0^2)/(RRI[L]^2))+((x0^2)/(RRI[L]^2)))<=1)
1062//                                              IX=-1
1063//                                      else
1064//                                              IX=0
1065//                                      endif
1066//                                      xmid=x0
1067//                                      ymid=y0*cth+z0*sth
1068//                                      zmid=-y0*sth+z0*cth
1069//                                      // end rotation about x begin rotn about y on rotated pts
1070//                                      //
1071//                                      xout=xmid*cph-zmid*sph
1072//                                      xout=xx[L]+xout
1073//                                      yout=ymid
1074//                                      yout=yy[L]+yout
1075//                                      zout=xmid*sph+zmid*cph
1076//                                      zout=zz[L]+zout
1077//
1078//                                      // now print to wave file the point or not depending on whether ix<0 or not
1079//
1080//                                      if (ix<0)
1081//                                      //write to wave file
1082//                                              InsertPoints nptW,1,xoutW,youtW,zoutW,sbpW,sldW
1083//                                              xoutW[nptW] = xout
1084//                                              youtW[nptW] = yout
1085//                                              zoutW[nptW] = zout
1086//                                              sbpW[nptW] = sbp[L]
1087//                                              sldW[nptW] = sld[L]
1088//                                             
1089//                                              nptW +=1
1090//                                     
1091//                                              //print  xout,yout,zout,sbp[L],sld[L]
1092//                                      //else
1093//                                              //continue
1094//                                      endif  //for write or not
1095//                              endfor  // for j
1096//                      endfor  //  for i
1097//              endfor  //for k
1098//      endfor // for L
1099//
1100//      // rescale to the sphere size
1101//      xoutW /= FFT_T
1102//      youtW /= FFT_T
1103//      zoutW /= FFT_T
1104//     
1105//      return(0) // end do loop cycle for cylinders
1106//end
Note: See TracBrowser for help on using the repository browser.