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

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

Moved rescaling panel from Wrapper.ipf to PlotUtils?.ipf, a more natural location.

Added "Manual Optimization" utility to the SANS Models/1D operations menu. this is a simple panel that allows users to "optimize" the fit in 1 or 2 directions. Very instructive to see whether you're near a minimum, and what the chi2 surface looks like around the minimum.

Did a similar Manual Optimization for the Real-space MultiClyinder? calculations from Ken Rubinson. simpler interface here.

All of the "CGB" fixes are present here, including the calculation of the number of guides.

File size: 22.5 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// Put this into a panel with the table and the data
376// and fields for all of the inputs
377Macro Setup_KR_MultiCylinder()
378       
379        Make/O/D/N=0 xx,yy,zz,rri,hti,sbp,rotx,roty,sld
380        Variable/G root:KR_Qmin = 0.004
381        Variable/G root:KR_Qmax = 0.4
382        Variable/G root:KR_Npt = 100
383
384        FFT_MakeMatrixButtonProc("")
385
386        NewPanel /W=(241,44,1169,458)/N=MultiCyl/K=1 as "Multi-Cylinder"
387        ModifyPanel cbRGB=(49825,57306,65535)
388
389        Button button_0,pos={45,80},size={100,20},proc=KR_Show3DButtonProc,title="Show 3D"
390        Button button_1,pos={46,51},size={100,20},proc=KR_Plot1DButtonProc,title="Plot 1D"
391        Button button_2,pos={178,50},size={150,20},proc=KR_GenerateButtonProc,title="Generate Structure"
392        Button button_4,pos={178,80},size={120,20},proc=KR_DoCalcButtonProc,title="Do Calculation"
393        Button button_3,pos={600,60},size={120,20},proc=KR_DeleteRow,title="Delete Row(s)"
394        Button button_5,pos={600,10},size={120,20},proc=KR_SaveTable,title="Save Table"
395        Button button_6,pos={600,35},size={120,20},proc=KR_ImportTable,title="Import Table"
396        ValDisplay valdisp_0,pos={339,16},size={80,13},title="FFT_T"
397        ValDisplay valdisp_0,limits={0,0,0},barmisc={0,1000},value= #"root:FFT_T"
398        SetVariable setvar_0,pos={339,40},size={140,15},title="Q min (A)"
399        SetVariable setvar_0,limits={0,10,0},value= KR_Qmin
400        SetVariable setvar_1,pos={339,65},size={140,15},title="Q max (A)"
401        SetVariable setvar_1,limits={0,10,0},value= KR_Qmax
402        SetVariable setvar_2,pos={339,90},size={140,15},title="Num Pts"
403        SetVariable setvar_2,limits={10,500,0},value= KR_Npt
404        CheckBox check_0,pos={599,93},size={59,14},title="Multi SLD",value= 0
405
406        Edit/W=(18,117,889,378)/HOST=#  xx,yy,zz,rri,hti,rotx,roty,sld
407        ModifyTable format(Point)=1,width(Point)=0
408        RenameWindow #,T0
409        SetActiveSubwindow ##
410
411End
412
413
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/////////////////////////////////////
615//
616//
617// for the "manual" fitting
618//
619
620Proc Vary_One_Cyl_Param(waveStr,row,percent,numStep)
621// pick the wave and row and %
622        String waveStr="xx"
623        Variable row=1
624        Variable percent = 105
625        Variable numStep = 10
626        Prompt waveStr,"wave",popup,"xx;yy;zz;rri;hti;rotx;roty;sld;"
627
628        print waveStr
629       
630// dispatch to calculation
631        MultiCyl_Loop($waveStr,row,percent,numStep)
632       
633// plot the chi2_map
634        DoWindow/F MultiCyl_ChiMap
635        if(V_flag==0)
636                MultiCyl_ChiMap()
637        endif
638
639end
640
641Function MultiCyl_Loop(w,row,percent,numStep)
642        Wave w
643        Variable row,percent,numStep
644       
645        Variable loLim,hiLim,ii,minIndex,minChiSq
646        String folderStr
647       
648        Make/O/D/N=(numStep) chi2_map,testVals
649        Wave chi2_map=chi2_map
650        Wave testVals=testVals
651       
652        loLim = w[row] - percent*w[row]/100
653        hiLim = w[row] + percent*w[row]/100
654        testVals = loLim + x*(hiLim-loLim)/numStep
655
656//              the experimental data
657        ControlInfo/W=WrapperPanel popup_0
658        folderStr=S_Value
659        // wave references for the data (to pass)
660        String DF="root:"+folderStr+":"
661       
662        WAVE yw=$(DF+folderStr+"_i")
663        WAVE xw=$(DF+folderStr+"_q")
664        WAVE sw=$(DF+folderStr+"_s")
665
666        duplicate/o yw, interpCalc,chi2_data
667        Wave chi2_data=chi2_data
668        Wave interpCalc=interpCalc
669               
670        STRUCT WMButtonAction ba
671        ba.eventCode = 2
672       
673// loop
674        for(ii=0;ii<numStep;ii+=1)
675//      set the value
676                w[row] = testVals[ii]
677//              generate the structure
678//
679                KR_GenerateButtonProc(ba)
680
681//              do the calculation
682                KR_DoCalcButtonProc(ba)
683
684                WAVE ival_KR=ival_KR
685                WAVE qval_KR=qval_KR
686               
687                interpCalc = interp(xw, qval_KR, ival_KR )
688               
689//              calculate chi-squared
690                chi2_data = (yw-interpCalc)^2
691                chi2_data /= sw^2
692       
693
694                Wavestats/Q chi2_data
695                chi2_map[ii] = V_avg * V_npnts
696// end loop
697        endfor
698
699// find the best chi squared
700        WaveStats/Q chi2_map
701// reset the value to the best
702        minIndex = V_minRowLoc
703        w[row] = testVals[minIndex]
704       
705        minChiSq = chi2_map[minIndex]
706        print "Minimum chi2 = ",minChiSq
707
708
709// and then recalculate at the best solution
710        KR_GenerateButtonProc(ba)
711
712//              do the calculation
713        KR_DoCalcButtonProc(ba)
714       
715        return(0)
716End
717
718Proc MultiCyl_ChiMap()
719        PauseUpdate; Silent 1           // building window...
720        Display /W=(35,44,466,414) chi2_map vs testVals
721        DoWindow/C MultiCyl_ChiMap
722        ModifyGraph mode=4
723        ModifyGraph marker=19
724        ModifyGraph msize=2
725        Label left "chi-squared"
726        Label bottom "test values"
727end
728
729//Function testKRPar()
730//     
731//      Variable row, col
732//      String wStr
733//     
734//      getParamFromKRSetup(row,col,wStr)
735//      Print row, col, wStr
736//     
737//      wStr = StringFromList(0, wStr)          // some wave "xx.d"
738//      wStr = StringFromList(0, wStr, ".")  // removes the ".d"
739//      Wave w=$wStr
740//      print w[row]
741//     
742//      Variable numStep,loLim,hiLim,percent
743//      numStep = 25
744//      percent = 50
745//     
746//      loLim = w[row] - percent*w[row]/100
747//      hiLim = w[row] + percent*w[row]/100
748//     
749//     
750//      Make/O/D/N=(numStep) testKRVals
751//      testKRVals = loLim + x*(hiLim-loLim)/numStep   
752//
753//      print testKRvals
754//      return(0)
755//     
756//End
757//
758//Function getParamFromKRSetup(row,col,wStr)
759//      Variable &row,&col
760//      String &wStr
761//     
762//      Variable parNum
763//     
764//      GetSelection table, MultiCyl#T0, 3
765//      row = V_startRow
766//      col = V_startCol
767//      Print S_Selection
768//      wStr = S_Selection
769//     
770//     
771//      return(0)
772//End
773
774Proc Vary_Two_Cyl_Param(waveStr,row,waveStr2,row2,percent,percent2,numStep)
775// pick the wave and row and %
776        String waveStr="xx"
777        Variable row=1
778        String waveStr2="rri"
779        Variable row2=0
780        Variable percent = 105
781        Variable percent2 = 50
782        Variable numStep=5
783        Prompt waveStr,"wave",popup,"xx;yy;zz;rri;hti;rotx;roty;sld;"
784        Prompt waveStr2,"wave2",popup,"xx;yy;zz;rri;hti;rotx;roty;sld;"
785       
786// dispatch to calculation
787        MultiCyl_Loop_2D($waveStr,row,$waveStr2,row2,percent,percent2,numStep)
788       
789// plot the chi2_map
790        DoWindow/F MultiCyl_ChiMap_2D
791        if(V_flag==0)
792                MultiCyl_ChiMap_2D()
793        else
794                //V_min*1.01 = the 1% neighborhood around the solution
795                WaveStats/Q chi2_Map_2D
796                ModifyImage/W=MultiCyl_ChiMap_2D chi2_Map_2D ctab= {(V_min*1.01),*,ColdWarm,0}
797                ModifyImage/W=MultiCyl_ChiMap_2D chi2_Map_2D minRGB=(0,65535,0),maxRGB=(0,65535,0)
798        endif
799
800end
801
802
803Function MultiCyl_Loop_2D(w,row,w2,row2,percent,percent2,numStep)
804        Wave w
805        Variable row
806        Wave w2
807        Variable row2,percent,percent2,numStep
808       
809        Variable loLim,hiLim,ii,jj,minIndex,minChiSq
810        String folderStr
811       
812        Make/O/D/N=(numStep,numStep) chi2_Map_2D
813        Make/O/D/N=(numStep) testVals,testVals2
814        Wave chi2_Map_2D=chi2_Map_2D
815        Wave testVals=testVals
816        Wave testVals2=testVals2
817       
818        testVals = 0
819        testVals2 = 0
820        chi2_Map_2D = 0
821       
822       
823        loLim = w[row] - percent*w[row]/100
824        hiLim = w[row] + percent*w[row]/100
825        testVals = loLim + x*(hiLim-loLim)/(numStep-1)
826//      Print lolim,hilim
827
828        SetScale/I x LoLim,HiLim,"", chi2_Map_2D
829
830        loLim = w2[row2] - percent2*w2[row2]/100
831        hiLim = w2[row2] + percent2*w2[row2]/100
832        testVals2 = loLim + x*(hiLim-loLim)/(numStep-1)
833//      Print lolim,hilim
834
835        SetScale/I y LoLim,HiLim,"", chi2_Map_2D
836
837
838//              the experimental data
839        ControlInfo/W=WrapperPanel popup_0
840        folderStr=S_Value
841        // wave references for the data (to pass)
842        String DF="root:"+folderStr+":"
843       
844        WAVE yw=$(DF+folderStr+"_i")
845        WAVE xw=$(DF+folderStr+"_q")
846        WAVE sw=$(DF+folderStr+"_s")
847
848        duplicate/o yw, interpCalc,chi2_data
849        Wave chi2_data=chi2_data
850        Wave interpCalc=interpCalc
851       
852        STRUCT WMButtonAction ba
853        ba.eventCode = 2
854       
855// double loop
856        for(ii=0;ii<numStep;ii+=1)
857                Print "                 Outer Loop Index = ",ii," out of ",numStep
858                //set the value from the outer loop
859                w[row] = testVals[ii]
860
861                for(jj=0;jj<numStep;jj+=1)
862               
863        //      set the inner value
864                        w2[row2] = testVals2[jj]
865        //              generate the structure
866        //
867                        KR_GenerateButtonProc(ba)
868       
869        //              do the calculation
870                        KR_DoCalcButtonProc(ba)
871       
872                        WAVE ival_KR=ival_KR
873                        WAVE qval_KR=qval_KR
874                       
875                        interpCalc = interp(xw, qval_KR, ival_KR )
876                       
877        //              calculate chi-squared
878                        chi2_data = (yw-interpCalc)^2
879                        chi2_data /= sw^2
880               
881                        Wavestats/Q chi2_data
882                        chi2_Map_2D[ii][jj] = V_avg * V_npnts
883                       
884                endfor
885        endfor
886
887
888// find the best chi squared
889        WaveStats/Q chi2_Map_2D
890// reset the value to the best
891        w[row] = V_MinRowLoc
892        w2[row2] = V_MinColLoc 
893       
894        minChiSq = V_Min
895        print "Minimum chi2 = ",minChiSq
896
897// and then recalculate at the best solution
898        KR_GenerateButtonProc(ba)
899
900//              do the calculation
901        KR_DoCalcButtonProc(ba)
902       
903        return(0)
904End
905
906Proc MultiCyl_ChiMap_2D()
907        PauseUpdate; Silent 1           // building window...
908        Display /W=(35,44,466,414)
909        AppendImage chi2_Map_2D
910        DoWindow/C MultiCyl_ChiMap_2D
911        ModifyImage chi2_Map_2D ctab= {*,*,ColdWarm,0}
912       
913        //V_min*1.01 = the 1% neighborhood around the solution
914        WaveStats/Q chi2_Map_2D
915        ModifyImage/W=MultiCyl_ChiMap_2D chi2_Map_2D ctab= {(V_min*1.01),*,ColdWarm,0}
916        ModifyImage/W=MultiCyl_ChiMap_2D chi2_Map_2D minRGB=(0,65535,0),maxRGB=(0,65535,0)
917       
918        Label bottom "test values"
919        Label left "test values"
920end
Note: See TracBrowser for help on using the repository browser.