source: sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotUtils2D_v40.ipf @ 762

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

Fixing bugs associated with 2D fitting, making sure that smeared and unsmeared 2D models both fit correctly and report errors to the FitManager?.

Added "_v40" extension to the Raspberry models. These models still need to be added to the Procedure list

File size: 32.4 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=4.00
3#pragma IgorVersion=6.1
4
5//
6// TO DO:
7//
8// - more intelligent beam stop masking
9// - constraints
10// - interactive masking
11//
12//
13// - OCT 2010   - added fDoBinning_QxQy2D(folderStr) that takes the QxQyQz data and bins it to I(Q)
14//                                      1D result must still be plotted manually, can be automated later.
15//                                      - for Scaled Image data (QxQy), the function fDoBinning_Scaled2D (in FFT_Cubes) could be
16//                                              modified and added to this file. It was meant for an FFT slice, but applies to
17//                                              2D model calculations (model_lin) 2D data files as well.
18
19
20//
21// changed June 2010 to read in the resolution information (now 8! columns)
22// -- subject to change --
23//
24// look for either the old-style 3-column (no resolution information) or the newer 8-column format
25Proc LoadQxQy()
26
27        LoadWave/G/D/W/A
28        String fileName = S_fileName
29        String path = S_Path
30        Variable numCols = V_flag
31
32        String w0,w1,w2,w3,w4,w5,w6,w7
33        String n0,n1,n2,n3,n4,n5,n6,n7
34               
35        if(numCols == 8)
36                // put the names of the 8 loaded waves into local names
37                n0 = StringFromList(0, S_waveNames ,";" )
38                n1 = StringFromList(1, S_waveNames ,";" )
39                n2 = StringFromList(2, S_waveNames ,";" )
40                n3 = StringFromList(3, S_waveNames ,";" )
41                n4 = StringFromList(4, S_waveNames ,";" )
42                n5 = StringFromList(5, S_waveNames ,";" )
43                n6 = StringFromList(6, S_waveNames ,";" )
44                n7 = StringFromList(7, S_waveNames ,";" )
45               
46                //remove the semicolon AND period from file names
47                w0 = CleanupName((S_fileName + "_qx"),0)
48                w1 = CleanupName((S_fileName + "_qy"),0)
49                w2 = CleanupName((S_fileName + "_i"),0)
50                w3 = CleanupName((S_fileName + "_iErr"),0)
51                w4 = CleanupName((S_fileName + "_qz"),0)
52                w5 = CleanupName((S_fileName + "_sQpl"),0)
53                w6 = CleanupName((S_fileName + "_sQpp"),0)
54                w7 = CleanupName((S_fileName + "_fs"),0)
55       
56                String baseStr=w1[0,strlen(w1)-4]
57                if(DataFolderExists("root:"+baseStr))
58                                DoAlert 1,"The file "+S_filename+" has already been loaded. Do you want to load the new data file, overwriting the data in memory?"
59                                if(V_flag==2)   //user selected No, don't load the data
60                                        SetDataFolder root:
61                                        KillWaves/Z $n0,$n1,$n2,$n3,$n4,$n5,$n6,$n7             // kill the default waveX that were loaded
62                                        return          //quits the macro
63                                endif
64                                SetDataFolder $("root:"+baseStr)
65                else
66                        NewDataFolder/S $("root:"+baseStr)
67                endif
68               
69                //read in the 18 lines of header (18th line starts w/ ASCII... 19th line is blank)
70                Make/O/T/N=18 header
71                Variable refnum,ii
72                string tmpStr=""
73                Open/R refNum  as (path+filename)
74                ii=0
75                do
76                        tmpStr = ""
77                        FReadLine refNum, tmpStr
78                        header[ii] = tmpStr
79                        ii+=1
80                while(ii < 18)         
81                Close refnum           
82               
83                ////overwrite the existing data, if it exists
84                Duplicate/O $("root:"+n0), $w0
85                Duplicate/O $("root:"+n1), $w1
86                Duplicate/O $("root:"+n2), $w2
87                Duplicate/O $("root:"+n3), $w3
88                Duplicate/O $("root:"+n4), $w4
89                Duplicate/O $("root:"+n5), $w5
90                Duplicate/O $("root:"+n6), $w6
91                Duplicate/O $("root:"+n7), $w7
92       
93        endif           //8-columns
94       
95        if(numCols == 3)
96                // put the names of the 3 loaded waves into local names
97                n0 = StringFromList(0, S_waveNames ,";" )
98                n1 = StringFromList(1, S_waveNames ,";" )
99                n2 = StringFromList(2, S_waveNames ,";" )
100
101                //remove the semicolon AND period from file names
102                w0 = CleanupName((S_fileName + "_qx"),0)
103                w1 = CleanupName((S_fileName + "_qy"),0)
104                w2 = CleanupName((S_fileName + "_i"),0)
105                w3 = CleanupName((S_fileName + "_iErr"),0)              //make a name for the error wave, to be generated here
106
107                String baseStr=w1[0,strlen(w1)-4]
108                if(DataFolderExists("root:"+baseStr))
109                                DoAlert 1,"The file "+S_filename+" has already been loaded. Do you want to load the new data file, overwriting the data in memory?"
110                                if(V_flag==2)   //user selected No, don't load the data
111                                        SetDataFolder root:
112                                        KillWaves/Z $n0,$n1,$n2         // kill the default waveX that were loaded
113                                        return          //quits the macro
114                                endif
115                                SetDataFolder $("root:"+baseStr)
116                else
117                        NewDataFolder/S $("root:"+baseStr)
118                endif
119               
120                //read in the 18 lines of header (18th line starts w/ ASCII... 19th line is blank)
121                Make/O/T/N=18 header
122                Variable refnum,ii
123                string tmpStr=""
124                Open/R refNum  as (path+filename)
125                ii=0
126                do
127                        tmpStr = ""
128                        FReadLine refNum, tmpStr
129                        header[ii] = tmpStr
130                        ii+=1
131                while(ii < 18)         
132                Close refnum           
133               
134                ////overwrite the existing data, if it exists
135                Duplicate/O $("root:"+n0), $w0
136                Duplicate/O $("root:"+n1), $w1
137                Duplicate/O $("root:"+n2), $w2
138       
139
140
141                // generate my own error wave for I(qx,qy). This is exactly the same procedure that is used in the QxQy_Export function
142                Duplicate/O $("root:"+n0), $w3
143                $w3 = sqrt($w2)         //assumes Poisson statistics for each cell (counter)
144                //      sw = 0.05*sw            // uniform 5% error? tends to favor the low intensity too strongly
145                // get rid of the "bad" errorsby replacing the NaN, Inf, and zero with V_avg
146                // THIS IS EXTREMEMLY IMPORTANT - if this is not done, there are some "bad" values in the
147                // error wave (things that are not numbers) - and this wrecks the smeared model fitting.
148                // It appears to have no effect on the unsmeared model.
149                WaveStats/Q $w3
150                $w3 = numtype($w3[p]) == 0 ? $w3[p] : V_avg
151                $w3 = $w3[p] != 0 ? $w3[p] : V_avg
152       
153        endif           //3-columns
154       
155       
156        /// do this for all 2D data, whether or not resolution information was read in
157       
158        Variable/G gIsLogScale = 0
159       
160        Variable num=numpnts($w0)
161        // assume that the Q-grid is "uniform enough" for DISPLAY ONLY
162        // use the 3 original waves for all of the fitting...
163        ConvertQxQy2Mat($w0,$w1,$w2,baseStr+"_mat")
164        Duplicate/O $(baseStr+"_mat"),$(baseStr+"_lin")                 //keep a linear-scaled version of the data
165       
166        PlotQxQy(baseStr)               //this sets the data folder back to root:!!
167
168        //clean up             
169        SetDataFolder root:
170        KillWaves/Z $n0,$n1,$n2,$n3,$n4,$n5,$n6,$n7
171       
172EndMacro
173
174//does not seem to need to be flipped at all from the standard QxQy output
175//
176// plots the data, as surface0, YellowHot color table
177//
178Proc PlotQxQy(str)
179        String str
180        Prompt str,"Pick the data folder containing 2D data",popup,getAList(4)
181
182        PauseUpdate; Silent 1   // Building Gizmo 6 window...
183
184        // Do nothing if the Gizmo XOP is not available.
185        if(exists("NewGizmo")!=4)
186                DoAlert 0, "Gizmo XOP must be installed"
187                return
188        endif
189
190        String data = "root:"+str+":"+str+"_mat"
191       
192        NewGizmo/N=$str/T=str /W=(169,44,495,359)
193        ModifyGizmo startRecMacro
194        AppendToGizmo Surface=$data,name=surface0
195
196        // for constant color (data is a blue color, lighter on the bottom)
197//      ModifyGizmo ModifyObject=surface0 property={ surfaceColorType,1}
198//      ModifyGizmo ModifyObject=surface0 property={ srcMode,0}
199//      ModifyGizmo ModifyObject=surface0 property={ frontColor,0.749996,0.811093,1,1}
200//      ModifyGizmo ModifyObject=surface0 property={ backColor,0.250019,0.433326,1,1}
201        //e constant color
202       
203        //for color table       
204        ModifyGizmo ModifyObject=surface0 property={ srcMode,0}
205//      ModifyGizmo ModifyObject=surface0 property={ surfaceCTab,Blue}
206        ModifyGizmo ModifyObject=surface0 property={ surfaceCTab,ColdWarm}
207        ModifyGizmo ModifyObject=surface0 property={ SurfaceCTABScaling,4}
208        //e color table
209       
210        AppendToGizmo Axes=boxAxes,name=axes0
211        ModifyGizmo ModifyObject=axes0,property={0,axisRange,-1,-1,-1,1,-1,-1}
212        ModifyGizmo ModifyObject=axes0,property={1,axisRange,-1,-1,-1,-1,1,-1}
213        ModifyGizmo ModifyObject=axes0,property={2,axisRange,-1,-1,-1,-1,-1,1}
214        ModifyGizmo ModifyObject=axes0,property={3,axisRange,-1,1,-1,-1,1,1}
215        ModifyGizmo ModifyObject=axes0,property={4,axisRange,1,1,-1,1,1,1}
216        ModifyGizmo ModifyObject=axes0,property={5,axisRange,1,-1,-1,1,-1,1}
217        ModifyGizmo ModifyObject=axes0,property={6,axisRange,-1,-1,1,-1,1,1}
218        ModifyGizmo ModifyObject=axes0,property={7,axisRange,1,-1,1,1,1,1}
219        ModifyGizmo ModifyObject=axes0,property={8,axisRange,1,-1,-1,1,1,-1}
220        ModifyGizmo ModifyObject=axes0,property={9,axisRange,-1,1,-1,1,1,-1}
221        ModifyGizmo ModifyObject=axes0,property={10,axisRange,-1,1,1,1,1,1}
222        ModifyGizmo ModifyObject=axes0,property={11,axisRange,-1,-1,1,1,-1,1}
223        ModifyGizmo ModifyObject=axes0,property={-1,axisScalingMode,1}
224        ModifyGizmo ModifyObject=axes0,property={-1,axisColor,0,0,0,1}
225        ModifyGizmo ModifyObject=axes0,property={0,ticks,3}
226        ModifyGizmo ModifyObject=axes0,property={1,ticks,3}
227        ModifyGizmo ModifyObject=axes0,property={2,ticks,3}
228        ModifyGizmo ModifyObject=axes0,property={0,fontScaleFactor,1.5}
229        ModifyGizmo ModifyObject=axes0,property={1,fontScaleFactor,1.5}
230        ModifyGizmo ModifyObject=axes0,property={2,fontScaleFactor,1.5}
231        AppendToGizmo freeAxesCue={0,0,0,1.5},name=freeAxesCue0
232        ModifyGizmo setDisplayList=0, opName=clearColor0, operation=clearColor, data={1,0.917,0.75,1}
233        ModifyGizmo setDisplayList=1, object=surface0
234        ModifyGizmo setDisplayList=2, object=axes0
235        ModifyGizmo setDisplayList=3, object=freeAxesCue0
236        ModifyGizmo SETQUATERNION={0.521287,-0.097088,-0.138769,0.836408}
237        ModifyGizmo autoscaling=1
238        ModifyGizmo currentGroupObject=""
239        ModifyGizmo compile
240
241        ModifyGizmo NamedHook={GizmoContours,WMGizmoContoursNamedHook}
242        ModifyGizmo endRecMacro
243       
244// don't bother with the flat image plot right now
245       
246//      Display $yw vs $xw
247//      modifygraph log=0
248//      ModifyGraph mode=3,marker=16,zColor($yw)={$zw,*,*,YellowHot,0}
249//      ModifyGraph standoff=0
250//      ModifyGraph width={Aspect,1}
251//      ModifyGraph lowTrip=0.001
252
253End
254
255Proc FakeQxQy(minQx,maxQx,minQy,maxQy,numPix,dataFolder)
256        Variable minQx=-0.1,maxQx=0.1,minQy=-0.1,maxQy=0.1
257        Variable numPix=64
258        String dataFolder="fake2DData"
259
260        String baseStr=dataFolder
261       
262
263        if(DataFolderExists("root:"+baseStr))
264                        DoAlert 1,"Fake data "+baseStr+" has already been created. Do you want to overwrite this fake data set?"
265                        if(V_flag==2)   //user selected No, don't load the data
266                                SetDataFolder root:
267                                return          //quits the macro
268                        endif
269                        SetDataFolder $("root:"+baseStr)
270        else
271                NewDataFolder/S $("root:"+baseStr)
272        endif
273
274        Make/O/D/N=(numPix*numPix) $(baseStr+"_qx"),$(baseStr+"_qy"),$(baseStr+"_i")
275       
276        $(baseStr+"_i") = 1
277       
278        Make/O/D/N=(numPix) tmpX,tmpY
279       
280        variable delX,delY
281        delX = (maxQx - minQx)/numPix
282        delY = (maxQy - minQy)/numPix
283        tmpX = minQx + x*delX
284        tmpY = minQy + x*delY
285       
286        // get rid of Q=0 values in the waves (may be a bad actor in model calculations)
287        //
288        tmpX = (tmpX == 0) ? 1e-6 : tmpX
289        tmpY = (tmpY == 0) ? 1e-6 : tmpY
290       
291        // X wave varies more rapidly
292        // Y wave is blocks of values
293        Variable ii=0
294        do
295                $(baseStr+"_qy")[ii*numPix,(ii+1)*numPix] = tmpY[ii]
296                ii+=1
297        while(ii<numPix)
298       
299        $(baseStr+"_qx") = tmpX[mod(p,numPix)]
300       
301        Variable/G gIsLogScale = 0
302       
303        // assume that the Q-grid is "uniform enough" for DISPLAY ONLY
304        // use the 3 original waves for all of the fitting...
305        ConvertQxQy2Mat($(baseStr+"_qx"),$(baseStr+"_qy"),$(baseStr+"_i"),baseStr+"_mat")
306        Duplicate/O $(baseStr+"_mat"),$(baseStr+"_lin")                 //keep a linear-scaled version of the data
307       
308        PlotQxQy(baseStr)               //this sets the data folder back to root:!!
309
310        //clean up             
311        SetDataFolder root:
312       
313EndMacro
314
315// this assumes that:
316// --QxQy data was written out in the format specified by the Igor macros, that is the x varies most rapidly
317// --the matrix is square!
318//
319// probably some other stuff...
320//
321Function UpdateQxQy2Mat(Qx,Qy,inten,linMat,mat)
322        Wave Qx,Qy,inten,linMat,mat
323       
324        Variable xrows=DimSize(mat, 0 )                 //assumes square matrix!!
325       
326        String folderStr=GetWavesDataFolder(Qx,1)
327        NVAR gIsLogScale=$(folderStr+"gIsLogScale")
328       
329        linMat = inten[q*xrows+p]
330       
331        if(gIsLogScale)
332                mat = log(linMat)
333        else
334                mat = linMat
335        endif
336       
337        return(0)
338End
339
340// this assumes that:
341// --QxQy data was written out in the format specified by the Igor macros, that is the x varies most rapidly
342// --the matrix is square!
343//
344// probably some other stuff...
345//
346Function ConvertQxQy2Mat(Qx,Qy,inten,matStr)
347        Wave Qx,Qy,inten
348        String matStr
349       
350        String folderStr=GetWavesDataFolder(Qx,1)
351       
352        Variable num=sqrt(numpnts(Qx))          //assumes square matrix, Qx = num x num points long
353        Make/O/D/N=(num,num) $(folderStr + matStr)
354        Wave mat=$matStr
355       
356        WaveStats/Q Qx
357        SetScale/I x, V_min, V_max, "", mat
358        WaveStats/Q Qy
359        SetScale/I y, V_min, V_max, "", mat
360       
361        Variable xrows=num
362       
363        mat = inten[q*xrows+p]
364       
365        return(0)
366End
367
368
369//str is the full path to the surface to append
370Proc AppendSurfaceToGizmo(str)
371        String str
372       
373        PauseUpdate; Silent 1   // Building Gizmo 6 window...
374
375        AppendToGizmo/Z Surface=$str,name=surface1
376       
377        // for a constant color (model is darker on top, lighter on the bottom)
378        //need these two lines, plus a color
379        ModifyGizmo ModifyObject=surface1 property={ surfaceColorType,1}
380        ModifyGizmo ModifyObject=surface1 property={ srcMode,0}
381        //green
382        ModifyGizmo ModifyObject=surface1 property={ frontColor,0.528923,0.882353,0.321584,1}
383        ModifyGizmo ModifyObject=surface1 property={ backColor,0.300221,0.6,1.5259e-05,1}
384        //red
385//      ModifyGizmo ModifyObject=surface1 property={ frontColor,1,0.749996,0.749996,1}
386//      ModifyGizmo ModifyObject=surface1 property={ backColor,1,0.250019,0.250019,1}
387
388        ModifyGizmo ModifyObject=surface1 property={ SurfaceCTABScaling,4}
389       
390        // for a color table (maybe not a good choice for vizualization if data uses a color table)
391//      ModifyGizmo ModifyObject=surface1 property={ srcMode,0}
392//      ModifyGizmo ModifyObject=surface1 property={ surfaceCTab,Red}
393       
394//      ModifyGizmo setDisplayList=0, object=surface0
395//      ModifyGizmo setDisplayList=1, object=axes0
396// object 3 is the axisCue
397        ModifyGizmo setDisplayList=4, object=surface1
398        ModifyGizmo SETQUATERNION={0.565517,-0.103105,-0.139134,0.806350}
399//      ModifyGizmo autoscaling=1
400//      ModifyGizmo currentGroupObject=""
401        ModifyGizmo compile
402
403//      ModifyGizmo endRecMacro
404End
405
406// would be nice, but I can't get this to work...
407//
408//Macro AdjustColorTables()
409//
410//      ModifyGizmo ModifyObject=surface1 property={ srcMode,0}
411//      ModifyGizmo ModifyObject=surface1 property={ surfaceCTab,Rainbow}
412//      ModifyGizmo setDisplayList=3, object=surface1
413////    ModifyGizmo SETQUATERNION={0.565517,-0.103105,-0.139134,0.806350}
414//      ModifyGizmo autoscaling=1
415//      ModifyGizmo currentGroupObject=""
416//      ModifyGizmo compile
417//     
418//end
419
420Function LogToggle2DButtonProc(ba) : ButtonControl
421        STRUCT WMButtonAction &ba
422
423        switch( ba.eventCode )
424                case 2: // mouse up
425                        // click code here
426                        String winStr=TopGizmoWindow()
427                        if(strlen(winstr)==0)
428                                Print "no gizmo window"
429                                break
430                        endif
431                        // the winStr is also the data folder
432                        // toggle everything in the data folder
433                        ToggleFolderScaling(winStr)             
434                        break
435        endswitch
436
437        return 0
438End
439
440Function Plot2DButtonProc(ba) : ButtonControl
441        STRUCT WMButtonAction &ba
442
443        switch( ba.eventCode )
444                case 2: // mouse up
445                        // click code here
446                        Execute "PlotQxQy()"
447                        break
448        endswitch
449
450        return 0
451End
452
453Function Append2DModelButtonProc(ba) : ButtonControl
454        STRUCT WMButtonAction &ba
455
456        switch( ba.eventCode )
457                case 2: // mouse up
458                        // click code here
459                        //str is the full path to the surface to append
460                        String DF=TopGizmoWindow()
461                        String str,funcStr,suffix
462                        ControlInfo/W=WrapperPanel popup_1
463                        funcStr = S_Value
464                        suffix = getModelSuffix(funcStr)
465                        if(stringmatch(funcStr, "smear*") == 1)
466                                suffix = "sm_"+suffix
467                        endif
468                        str = "root:"+DF+":"+suffix+"_mat"
469                       
470                        Execute "AppendSurfaceToGizmo(\""+str+"\")"
471//                      Print str
472                        break
473        endswitch
474
475        return 0
476End
477
478// toggle the scaling of every matrix in the folder
479Function ToggleFolderScaling(DF)
480        String DF
481
482        // look for waves DF+"_mat" (the data)
483        // and models prefix + "_mat"
484        SetDataFolder $("root:"+DF)
485        //check the global to see the state of the data
486        NVAR gIsLogScale = gIsLogScale
487       
488        String matrixList=WaveList("*_mat",";",""),item
489        Variable ii=0,num=itemsinlist(matrixlist,";"),len
490       
491        for(ii=0;ii<num;ii+=1)
492                        item = StringFromList(ii,matrixList,";")
493                        len = strlen(item)
494                        Wave w = $item
495                        Wave linW = $(item[0,len-5]+"_lin")
496                        if(gisLogScale)
497                                //make linear
498                                w=linW
499                        else   
500                                //make log scale
501                                w=log(linW)
502                                w = w[p][q] == -inf ? NaN : w[p][q]             //remove the -inf for display
503                        endif
504        endfor 
505       
506        // toggle the global
507        gIsLogScale = !gIsLogScale
508       
509        SetDataFolder root:
510        return(0)
511End
512
513
514Function Load2DDataButtonProc(ba) : ButtonControl
515        STRUCT WMButtonAction &ba
516
517        switch( ba.eventCode )
518                case 2: // mouse up
519                        // click code here
520                        Execute "LoadQxQy()"
521                        break
522        endswitch
523
524        return 0
525End
526
527
528// dispatch the fit
529Function Do2DFitButtonProc(ba) : ButtonControl
530        STRUCT WMButtonAction &ba
531
532        String folderStr,funcStr,coefStr
533        Variable useCursors,useEps,useConstr
534       
535        switch( ba.eventCode )
536                case 2: // mouse up
537                        ControlInfo/W=WrapperPanel popup_0
538                        folderStr=S_Value
539                       
540                        ControlInfo/W=WrapperPanel popup_1
541                        funcStr=S_Value
542                       
543                        ControlInfo/W=WrapperPanel popup_2
544                        coefStr=S_Value
545                       
546                        ControlInfo/W=WrapperPanel check_0
547                        useCursors=V_Value
548                        //
549                        // NO CURSORS for 2D waves - force to zero
550                        //
551                        useCursors = 0
552                        //
553                        ControlInfo/W=WrapperPanel check_1
554                        useEps=V_Value
555                        ControlInfo/W=WrapperPanel check_2
556                        useConstr=V_Value
557                       
558                        if(!CheckFunctionAndCoef(funcStr,coefStr))
559                                DoAlert 0,"The coefficients and function type do not match. Please correct the selections in the popup menus."
560                                break
561                        endif
562                       
563                        FitWrapper2D(folderStr,funcStr,coefStr,useCursors,useEps,useConstr)
564                       
565                        //      DoUpdate (does not work!)
566                        //?? why do I need to force an update ??
567                        if(!exists("root:"+folderStr+":"+coefStr))
568                                Wave w=$coefStr
569                        else
570                                Wave w=$("root:"+folderStr+":"+coefStr) //smeared coefs in data folder
571                        endif
572                        w[0] += 1e-6
573                        w[0] -= 1e-6
574       
575                        break
576        endswitch
577
578        return 0
579End
580
581// wrapper to do the desired fit, pretty much a clone of FitWrapper, the 1D version
582//
583// folderStr is the data folder for the desired data set
584//
585// Currently the limitations are:
586// - I have no error waves for the intensity (fixed 10/2010)
587// - There is no smeared model (coming soon after 10/2010)
588// - Cursors can't be used
589// - the report works OK, but I have little control over the graphics
590// - the mask is generated here with a default radius of 8 pixels around the beam center
591//
592Function FitWrapper2D(folderStr,funcStr,coefStr,useCursors,useEps,useConstr)
593        String folderStr,funcStr,coefStr
594        Variable useCursors,useEps,useConstr
595
596        //These only make sense for the 1D fits, but put them here so keep the look of the dispatching the same
597        Variable useResiduals, useTextBox
598        useResiduals = 0
599        useTextBox = 0
600       
601        String suffix=getModelSuffix(funcStr)
602       
603        SetDataFolder $("root:"+folderStr)
604        if(!exists(coefStr))
605                // must be unsmeared model, work in the root folder
606                SetDataFolder root:                             
607                if(!exists(coefStr))            //this should be fine if the coef filter is working, but check anyhow
608                        DoAlert 0,"the coefficient and data sets do not match"
609                        return 0
610                endif
611        endif
612               
613        WAVE cw=$(coefStr)     
614        Wave hold=$("Hold_"+suffix)
615        Wave/T lolim=$("LoLim_"+suffix)
616        Wave/T hilim=$("HiLim_"+suffix)
617        Wave eps=$("epsilon_"+suffix)
618       
619// fill a struct instance whether I need one or not
620// note that the resolution waves may or may not exist, and may or may not be used in the fitting
621        String DF="root:"+folderStr+":"
622       
623        WAVE inten=$(DF+folderStr+"_i")
624        WAVE sw=$(DF+folderStr+"_iErr")
625        WAVE qx=$(DF+folderStr+"_qx")
626        WAVE qy=$(DF+folderStr+"_qy")
627        WAVE/Z qz=$(DF+folderStr+"_qz")
628        WAVE/Z sQpl=$(DF+folderStr+"_sQpl")
629        WAVE/Z sQpp=$(DF+folderStr+"_sQpp")
630        WAVE/Z shad=$(DF+folderStr+"_fs")
631
632//just a dummy - I shouldn't need this
633        Duplicate/O qx resultW
634        resultW=0
635       
636        STRUCT ResSmear_2D_AAOStruct s
637        WAVE s.coefW = cw       
638        WAVE s.zw = resultW     
639        WAVE s.xw[0] = qx
640        WAVE s.xw[1] = qy
641        WAVE/Z s.qz = qz
642        WAVE/Z s.sQpl = sQpl
643        WAVE/Z s.sQpp = sQpp
644        WAVE/Z s.fs = shad
645       
646
647        // generate my own mask wave - as a matrix first, then redimension to N*N vector
648        // default mask is two pixels all the way around, (0 is excluded, 1 is included)
649        WAVE DataMat=$(DF+folderStr+"_lin")
650        if(exists(DF+"mask") == 0)
651                Duplicate/O dataMat mask
652                Variable bsRadius=8             //pixels?
653                MakeBSMask(mask,bsRadius)
654        Endif
655       
656       
657        Duplicate/O inten inten_masked
658        inten_masked = (mask[p][q] == 0) ? NaN : inten[p][q]
659       
660
661        //for now, use res will always be 0 for 2D functions   
662        Variable useResol=0
663        if(stringmatch(funcStr, "Smear*"))              // if it's a smeared function, need a struct
664                useResol=1
665        endif
666
667        // can't use constraints in this way for multivariate fits. See the curve fitting help file
668        // and "Contraint Matrix and Vector"
669        if(useConstr)
670                Print "Constraints not yet implemented"
671                useConstr = 0
672        endif   
673        WAVE/Z constr=constr            //will be a null reference
674       
675//      // do not construct constraints for any of the coefficients that are being held
676//      // -- this will generate an "unknown error" from the curve fitting
677//      Make/O/T/N=0 constr
678//      if(useConstr)
679//              String constraintExpression
680//              Variable i, nPnts=DimSize(lolim, 0),nextRow=0
681//              for (i=0; i < nPnts; i += 1)
682//                      if (strlen(lolim[i]) > 0 && hold[i] == 0)
683//                              InsertPoints nextRow, 1, constr
684//                              sprintf constraintExpression, "K%d > %s", i, lolim[i]
685//                              constr[nextRow] = constraintExpression
686//                              nextRow += 1
687//                      endif
688//                      if (strlen(hilim[i]) > 0 && hold[i] == 0)
689//                              InsertPoints nextRow, 1, constr
690//                              sprintf constraintExpression, "K%d < %s", i, hilim[i]
691//                              constr[nextRow] = constraintExpression
692//                              nextRow += 1
693//                      endif
694//              endfor
695//      endif
696
697        if(useCursors)
698                Print "Cursors not yet implemented"
699                useCursors = 0
700        endif   
701///// NO CURSORS for 2D waves
702        //if useCursors, and the data is USANS, need to feed a (reassigned) trimmed matrix to the fit
703        Variable pt1,pt2,newN
704        pt1 = 0
705        pt2 = numpnts(inten)-1
706//      if(useCursors && (dimsize(resW,1) > 4) )
707//              if(pcsr(A) > pcsr(B))
708//                      pt1 = pcsr(B)
709//                      pt2 = pcsr(A)
710//              else
711//                      pt1 = pcsr(A)
712//                      pt2 = pcsr(B)
713//              endif
714//              newN = pt2 - pt1 + 1            // +1 includes both cursors in the fit
715//              Make/O/D/N=(newN,newN) $(DF+"crsrResW")
716//              WAVE crsrResW = $(DF+"crsrResW")
717//              crsrResW = resW[p+pt1][q+pt1]
718//              //assign to the struct
719//              WAVE fs.resW =  crsrResW               
720//      endif
721
722// create these variables so that FuncFit will set them on exit
723        Variable/G V_FitError=0                         //0=no err, 1=error,(2^1+2^0)=3=singular matrix
724        Variable/G V_FitQuitReason=0            //0=ok,1=maxiter,2=user stop,3=no chisq decrease
725       
726// don't use the auto-destination with no flag, it doesn't appear to work correctly
727// dispatch the fit
728        //      FuncFit/H="11110111111"/NTHR=0 Cylinder2D_D :cyl2d_c_txt:coef_Cyl2D_D  :cyl2d_c_txt:cyl2d_c_txt_i /X={:cyl2d_c_txt:cyl2d_c_txt_qy,:cyl2d_c_txt:cyl2d_c_txt_qx} /W=:cyl2d_c_txt:sw /I=1 /M=:cyl2d_c_txt:mask /D
729        Variable t1=StopMSTimer(-2)
730        Variable tb = 0         //no textbox
731
732// /NTHR=1 means just one thread for the fit (since the function itself is threaded)
733// NTHR = 0 == "Auto" mode, using as many processors as are available (not appropriate here since the function itself is threaded?)
734
735        do
736       
737                        // now useCursors, useEps, and useConstr are all handled w/ /NWOK, just like FitWrapper
738
739
740//              if(useResol && useResiduals && useTextBox)              //do it all
741//                      FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /R /NWOK
742//                      break
743//              endif
744//             
745//              if(useResol && useResiduals)            //res + resid
746//                      FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /R /NWOK
747//                      break
748//              endif
749//
750//             
751//              if(useResol && useTextBox)              //res + text
752//                      FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /NWOK
753//                      break
754//              endif
755               
756                if(useResol)            //res only
757                        Print "useRes only"
758                        FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /NWOK
759                        break
760                endif
761                       
762                               
763/////   same as above, but all without useResol (no /STRC flag)
764                if(useResiduals && useTextBox)          //resid+ text
765                        FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /R /NWOK
766                        break
767                endif
768               
769                if(useResiduals)                //resid
770                        FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /R /NWOK
771                        break
772                endif
773
774               
775                if(useTextBox)          //text
776                        FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /NWOK
777                        break
778                endif
779               
780                //just a plain vanilla fit
781
782                FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /NWOK
783               
784        while(0)
785       
786        Print "elapsed 2D fit time  = ",(StopMSTimer(-2) - t1)/1e6," s = ",(StopMSTimer(-2) - t1)/1e6/60," min"
787
788        // append the fit
789        // need to manage duplicate copies
790        // Don't plot the full curve if cursors were used (set fitYw to NaN on entry...)
791//      String traces=TraceNameList("", ";", 1 )                //"" as first parameter == look on the target graph
792//      if(strsearch(traces,"FitYw",0) == -1)
793//              AppendToGraph FitYw vs xw
794//      else
795//              RemoveFromGraph FitYw
796//              AppendToGraph FitYw vs xw
797//      endif
798//      ModifyGraph lsize(FitYw)=2,rgb(FitYw)=(0,0,0)
799       
800//      DoUpdate                //force update of table and graph with fitted values (why doesn't this work? - the table still does not update)
801       
802        // report the results (to the panel?)
803        print "V_chisq = ",V_chisq
804        print cw
805        WAVE w_sigma
806        print w_sigma
807        String resultStr=""
808       
809        if(waveexists(W_sigma))
810                //append it to the table, if it's not already there
811                CheckDisplayed/W=WrapperPanel#T0 W_sigma
812                if(V_flag==0)
813                        //not there, append it
814                        AppendtoTable/W=wrapperPanel#T0 W_sigma
815                else
816                        //remove it, and put it back on to make sure it's the right one (do I need to do this?)
817                        // -- not really, since any switch of the function menu takes W_Sigma off
818                endif
819        endif
820               
821        //now re-write the results
822        sprintf resultStr,"Chi^2 = %g  Sqrt(X^2/N) = %g",V_chisq,sqrt(V_chisq/V_Npnts)
823        resultStr = PadString(resultStr,63,0x20)
824        GroupBox grpBox_2 title=resultStr
825        ControlUpdate/W=WrapperPanel grpBox_2
826        sprintf resultStr,"FitErr = %s : FitQuit = %s",W_ErrorMessage(V_FitError),W_QuitMessage(V_FitQuitReason)
827        resultStr = PadString(resultStr,63,0x20)
828        GroupBox grpBox_3 title=resultStr
829        ControlUpdate/W=WrapperPanel grpBox_3
830       
831        Variable yesSave=0,yesReport=0
832        ControlInfo/W=WrapperPanel check_4
833        yesReport = V_Value
834        ControlInfo/W=WrapperPanel check_5
835        yesSave = V_Value
836       
837       
838        if(yesReport)
839                String parStr = getFunctionParams(funcStr)
840//              String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+suffix, "", "TEXT:1," )              //old way, but doesn't work in 2D folders
841                String topGraph= TopGizmoWindow()       //this is the topmost Gizmo (XOP) window
842               
843                DoUpdate                //force an update of the graph before making a copy of it for the report
844       
845                W_GenerateReport(funcStr,folderStr,$parStr,cw,yesSave,V_chisq,W_sigma,V_npnts,V_FitError,V_FitQuitReason,V_startRow,V_endRow,topGraph)
846        endif
847       
848        SetDataFolder root:
849        return(0)
850End
851
852
853// right now, there are only unsmeared models, but I want experimental points
854// for the QxQy of the calculation - so get the folder string
855Function Plot2DFunctionButtonProc(ba) : ButtonControl
856        STRUCT WMButtonAction &ba
857
858        String folderStr,funcStr,coefStr,cmdStr=""
859       
860        Variable killWhat=0             //kill nothing as default
861
862        switch( ba.eventCode )
863                case 2: // mouse up
864                        ControlInfo/W=WrapperPanel popup_0
865                        folderStr=S_Value
866                       
867                        ControlInfo/W=WrapperPanel popup_1
868                        funcStr=S_Value
869                       
870                        // check for smeared or smeared function
871//                      if(stringmatch(funcStr, "Smear*" )==1)
872                                //it's a smeared model
873                                // check for the special case of RPA that has an extra parameter
874//                              if(strsearch(funcStr, "RPAForm", 0 ,0) == -1)
875                                        sprintf cmdStr, "Plot%s(\"%s\")",funcStr,folderStr              //not RPA
876//                              else
877//                                      sprintf cmdStr, "Plot%s(\"%s\",)",funcStr,folderStr             //yes RPA, leave a comma for input
878//                              endif
879//                      else
880//                              // it's not,                   
881//                              sprintf cmdStr, "Plot%s()",funcStr
882//                      endif
883                       
884                        //Print cmdStr
885                        Execute cmdStr
886                       
887                        //pop the function menu to set the proper coefficients
888                        DoWindow/F WrapperPanel
889                        STRUCT WMPopupAction pa
890                        pa.popStr = funcStr
891                        pa.eventcode = 2
892                        Function_PopMenuProc(pa)
893       
894                        KillWhat = 2            //kill just the table, leave the 2d visible for now
895                        KillTopGraphAndTable(killWhat)          // crude
896
897                        break
898        endswitch
899
900        return 0
901End
902
903// unused - all of this functionality has been added to the Wrapper Panel
904Window Plot_2D_Controls()
905        PauseUpdate; Silent 1           // building window...
906        NewPanel /W=(950,44,1250,244)
907        Button button0,pos={172,37},size={70,20},proc=LogToggle2DButtonProc,title="Log/Lin"
908        Button button1,pos={146,9},size={100,20},proc=Plot2DButtonProc,title="Plot 2D Data"
909        Button button2,pos={164,118},size={120,20},proc=Append2DModelButtonProc,title="Append Model"
910        PopupMenu popup_1,pos={9,84},size={132,20},proc=Function_PopMenuProc,title="Function"
911        PopupMenu popup_1,mode=1,popvalue="Cylinder2D",value= #"W_FunctionPopupList()"
912        Button button3,pos={11,9},size={100,20},proc=Load2DDataButtonProc,title="Load 2D Data"
913        Button button4,pos={164,161},size={120,20},proc=Do2DFitButtonProc,title="Do 2D Fit"
914        Button button5,pos={164,85},size={120,20},proc=Plot2DFunctionButtonProc,title="Plot 2D Function"
915EndMacro
916
917
918Function/S TopGizmoWindow()
919        Return(StringFromList(0,WinList("*",";","WIN:4096")))
920end
921
922////
923//
924// unused proc for testing
925Proc Plot2DVsPointNumber(str)
926        String str
927        Prompt str,"Pick the data folder containing 2D data",popup,getAList(4)
928
929        PauseUpdate; Silent 1           // building window...
930        String fldrSav0= GetDataFolder(1)
931        SetDataFolder $("root:"+str)
932        Display /W=(241,352,810,842) $(str+"_i")
933        SetDataFolder fldrSav0
934        ModifyGraph mode=3
935        ModifyGraph marker=19
936        ModifyGraph lSize=2
937        ModifyGraph rgb($(str+"_i"))=(0,0,0)
938        ModifyGraph msize=1
939        ModifyGraph grid=1
940        ModifyGraph log(left)=1
941        ModifyGraph mirror=2
942        if(exists("root:"+str+":sw"))
943                ErrorBars/T=0 $(str+"_i") Y,wave=($("root:"+str+":sw"),$("root:"+str+":sw"))
944        endif
945        SetDataFolder fldrSav0
946EndMacro
947
948// testing procedure, used when I had x,y swapped and was getting nonsensical results
949//Macro CalculateChiSquared(str)
950//      String str
951//      Prompt str,"Pick the data folder containing 2D data",popup,getAList(4)
952//
953//
954//      String fldrSav0= GetDataFolder(1)
955//      SetDataFolder $("root:"+str)
956//
957//      Duplicate/O $(str+"_i") chi
958//      chi = ((zwave_cyl2D_D - $(str+"_i"))/sw )^2
959//     
960//      chi = (mask == 1) ? chi : 0
961//     
962//      Print sum(chi,-inf,inf)
963//     
964//      SetDataFolder fldrSav0
965//EndMacro
966
967
968Function MakeBSMask(mask,rad)
969        Wave mask
970        Variable rad
971
972
973// find the center based on the wave scaling
974        Variable xCtr,yCtr,Qzero=0
975        xCtr = (Qzero - DimOffset(Mask, 0))/DimDelta(Mask,0)
976        yCtr = (Qzero - DimOffset(Mask, 1))/DimDelta(Mask,1)
977
978        Print xctr,yctr
979
980//      Variable center = sqrt(numpnts(mask))/2 -0.5
981
982        mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1
983
984        Variable xDim, yDim
985        xDim = DimSize(mask,0)
986        yDim = DimSize(mask,1)
987        mask[][0] = 0
988        mask[][1] = 0
989        mask[][yDim-2] = 0
990        mask[][yDim-1] = 0
991        mask[0][] = 0
992        mask[1][] = 0
993        mask[xDim-2][] = 0
994        mask[xDim-1][] = 0
995       
996        Redimension/N=(xDim*yDim) mask          //now 1D
997       
998       
999End
1000
1001// This routine assumes that the 2D data was loaded with the NCNR loader, so that the
1002// data is in a data folder, and the extensions are known. A more generic form could
1003// be made too, if needed.
1004//
1005// X- need error on I(q)
1006// -- need to set "proper" number of data points (delta and qMax?)
1007// X- need to remove points at high Q end
1008//
1009// -- like the routines in CircSectAve, start with 500 points, and trim after binning is done.
1010//      you'l end up with < 200 points.
1011//
1012// the results are in iBin_qxqy, qBin_qxqy, and eBin_qxqy, in the folder passed
1013//
1014//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
1015Function fDoBinning_QxQy2D(folderStr)
1016        String folderStr
1017
1018//      Wave inten,qx,qy,qz
1019
1020        SetDataFolder $("root:"+folderStr)
1021       
1022        WAVE inten = $(folderStr + "_i")
1023        WAVE qx = $(folderStr + "_qx")
1024        WAVE qy = $(folderStr + "_qy")
1025        WAVE qz = $(folderStr + "_qz")
1026       
1027        Variable xDim=numpnts(qx),yDim
1028        Variable ii,jj,delQ
1029        Variable qTot,nq,var,avesq,aveisq
1030        Variable binIndex,val
1031       
1032        nq = 500
1033       
1034        yDim = XDim
1035        Make/O/D/N=(nq) iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy
1036        delQ = abs(sqrt(qx[2]^2+qy[2]^2+qz[2]^2) - sqrt(qx[1]^2+qy[1]^2+qz[1]^2))               //use bins of 1 pixel width
1037        qBin_qxqy[] =  p*       delQ   
1038        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1039
1040        iBin_qxqy = 0
1041        iBin2_qxqy = 0
1042        eBin_qxqy = 0
1043        nBin_qxqy = 0   //number of intensities added to each bin
1044       
1045        for(ii=0;ii<xDim;ii+=1)
1046                qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1047                binIndex = trunc(x2pnt(qBin_qxqy, qTot))
1048                val = inten[ii]
1049                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
1050                        iBin_qxqy[binIndex] += val
1051                        iBin2_qxqy[binIndex] += val*val
1052                        nBin_qxqy[binIndex] += 1
1053                endif
1054        endfor
1055
1056//calculate errors, just like in CircSectAve.ipf
1057        for(ii=0;ii<nq;ii+=1)
1058                if(nBin_qxqy[ii] == 0)
1059                        //no pixels in annuli, data unknown
1060                        iBin_qxqy[ii] = 0
1061                        eBin_qxqy[ii] = 1
1062                else
1063                        if(nBin_qxqy[ii] <= 1)
1064                                //need more than one pixel to determine error
1065                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1066                                eBin_qxqy[ii] = 1
1067                        else
1068                                //assume that the intensity in each pixel in annuli is normally
1069                                // distributed about mean...
1070                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1071                                avesq = iBin_qxqy[ii]^2
1072                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1073                                var = aveisq-avesq
1074                                if(var<=0)
1075                                        eBin_qxqy[ii] = 1e-6
1076                                else
1077                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1078                                endif
1079                        endif
1080                endif
1081        endfor
1082       
1083        // find the last non-zero point, working backwards
1084        val=nq
1085        do
1086                val -= 1
1087        while(nBin_qxqy[val] == 0)
1088       
1089//      print val, nBin_qxqy[val]
1090        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy
1091       
1092        SetDataFolder root:
1093       
1094        return(0)
1095End
Note: See TracBrowser for help on using the repository browser.