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

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

Adding new model functions to the picker list

Adding resolution-smeared functions to the 2D analysis models

Constraints are now functional during a 2D fit

vesib.cor example data is now 6-column

File size: 33.1 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 defined as a single text wave for multivariate fits. See the curve fitting help file
668        // and "Contraint Matrix and Vector"
669        //
670        // -- generate a constraint text wave, then the /C flag automatically generates a constraint matrix and vector
671        // -- then use the proper logic to dispatch (can't use /NWOK anymore)
672       
673//      if(useConstr)
674//              Print "Constraints not yet implemented"
675//              useConstr = 0
676//      endif   
677//      WAVE/Z constr=constr            //will be a null reference
678       
679        // do not construct constraints for any of the coefficients that are being held
680        // -- this will generate an "unknown error" from the curve fitting
681        Make/O/T/N=0 constr
682        if(useConstr)
683                String constraintExpression
684                Variable i, nPnts=DimSize(lolim, 0),nextRow=0
685                for (i=0; i < nPnts; i += 1)
686                        if (strlen(lolim[i]) > 0 && hold[i] == 0)
687                                InsertPoints nextRow, 1, constr
688                                sprintf constraintExpression, "K%d > %s", i, lolim[i]
689                                constr[nextRow] = constraintExpression
690                                nextRow += 1
691                        endif
692                        if (strlen(hilim[i]) > 0 && hold[i] == 0)
693                                InsertPoints nextRow, 1, constr
694                                sprintf constraintExpression, "K%d < %s", i, hilim[i]
695                                constr[nextRow] = constraintExpression
696                                nextRow += 1
697                        endif
698                endfor
699
700        endif
701
702        if(useCursors)
703                Print "Cursors not yet implemented"
704                useCursors = 0
705        endif   
706///// NO CURSORS for 2D waves
707        //if useCursors, and the data is USANS, need to feed a (reassigned) trimmed matrix to the fit
708        Variable pt1,pt2,newN
709        pt1 = 0
710        pt2 = numpnts(inten)-1
711//      if(useCursors && (dimsize(resW,1) > 4) )
712//              if(pcsr(A) > pcsr(B))
713//                      pt1 = pcsr(B)
714//                      pt2 = pcsr(A)
715//              else
716//                      pt1 = pcsr(A)
717//                      pt2 = pcsr(B)
718//              endif
719//              newN = pt2 - pt1 + 1            // +1 includes both cursors in the fit
720//              Make/O/D/N=(newN,newN) $(DF+"crsrResW")
721//              WAVE crsrResW = $(DF+"crsrResW")
722//              crsrResW = resW[p+pt1][q+pt1]
723//              //assign to the struct
724//              WAVE fs.resW =  crsrResW               
725//      endif
726
727// create these variables so that FuncFit will set them on exit
728        Variable/G V_FitError=0                         //0=no err, 1=error,(2^1+2^0)=3=singular matrix
729        Variable/G V_FitQuitReason=0            //0=ok,1=maxiter,2=user stop,3=no chisq decrease
730       
731// don't use the auto-destination with no flag, it doesn't appear to work correctly
732// dispatch the fit
733        //      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
734        Variable t1=StopMSTimer(-2)
735        Variable tb = 0         //no textbox
736
737// /NTHR=1 means just one thread for the fit (since the function itself is threaded)
738// NTHR = 0 == "Auto" mode, using as many processors as are available (not appropriate here since the function itself is threaded?)
739
740        do
741       
742                        // now useEps, and useConstr are all handled w/ /NWOK, just like FitWrapper
743                        // useCursors needs to have the /C flag in the command for the constraint matrix and vector to be auto-generated
744
745//              if(useResol && useResiduals && useTextBox)              //do it all
746//                      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
747//                      break
748//              endif
749//             
750//              if(useResol && useResiduals)            //res + resid
751//                      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
752//                      break
753//              endif
754//
755//             
756//              if(useResol && useTextBox)              //res + text
757//                      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
758//                      break
759//              endif
760               
761                if(useResol && useConstr)               //res  and constraints
762                        Print "useRes only"
763                        FuncFit/C/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
764                        break
765                endif
766               
767                if(useResol)            //res only
768                        Print "useRes only"
769                        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
770                        break
771                endif
772                               
773/////   same as above, but all without useResol (no /STRC flag)
774//              if(useResiduals && useTextBox)          //resid+ text
775//                      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
776//                      break
777//              endif
778//             
779//              if(useResiduals)                //resid
780//                      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
781//                      break
782//              endif
783//
784//             
785//              if(useTextBox)          //text
786//                      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
787//                      break
788//              endif
789               
790                if(useConstr)
791                        FuncFit/C/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /NWOK
792                        break
793                Endif
794               
795                //just a plain vanilla fit
796
797                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
798               
799        while(0)
800       
801        Print "elapsed 2D fit time  = ",(StopMSTimer(-2) - t1)/1e6," seconds = ",(StopMSTimer(-2) - t1)/1e6/60," minutes"
802
803        // append the fit
804        // need to manage duplicate copies
805        // Don't plot the full curve if cursors were used (set fitYw to NaN on entry...)
806//      String traces=TraceNameList("", ";", 1 )                //"" as first parameter == look on the target graph
807//      if(strsearch(traces,"FitYw",0) == -1)
808//              AppendToGraph FitYw vs xw
809//      else
810//              RemoveFromGraph FitYw
811//              AppendToGraph FitYw vs xw
812//      endif
813//      ModifyGraph lsize(FitYw)=2,rgb(FitYw)=(0,0,0)
814       
815//      DoUpdate                //force update of table and graph with fitted values (why doesn't this work? - the table still does not update)
816       
817        // report the results (to the panel?)
818        print "V_chisq = ",V_chisq
819        print cw
820        WAVE w_sigma
821        print w_sigma
822        String resultStr=""
823       
824        if(waveexists(W_sigma))
825                //append it to the table, if it's not already there
826                CheckDisplayed/W=WrapperPanel#T0 W_sigma
827                if(V_flag==0)
828                        //not there, append it
829                        AppendtoTable/W=wrapperPanel#T0 W_sigma
830                else
831                        //remove it, and put it back on to make sure it's the right one (do I need to do this?)
832                        // -- not really, since any switch of the function menu takes W_Sigma off
833                endif
834        endif
835               
836        //now re-write the results
837        sprintf resultStr,"Chi^2 = %g  Sqrt(X^2/N) = %g",V_chisq,sqrt(V_chisq/V_Npnts)
838        resultStr = PadString(resultStr,63,0x20)
839        GroupBox grpBox_2 title=resultStr
840        ControlUpdate/W=WrapperPanel grpBox_2
841        sprintf resultStr,"FitErr = %s : FitQuit = %s",W_ErrorMessage(V_FitError),W_QuitMessage(V_FitQuitReason)
842        resultStr = PadString(resultStr,63,0x20)
843        GroupBox grpBox_3 title=resultStr
844        ControlUpdate/W=WrapperPanel grpBox_3
845       
846        Variable yesSave=0,yesReport=0
847        ControlInfo/W=WrapperPanel check_4
848        yesReport = V_Value
849        ControlInfo/W=WrapperPanel check_5
850        yesSave = V_Value
851       
852       
853        if(yesReport)
854                String parStr = getFunctionParams(funcStr)
855//              String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+suffix, "", "TEXT:1," )              //old way, but doesn't work in 2D folders
856                String topGraph= TopGizmoWindow()       //this is the topmost Gizmo (XOP) window
857               
858                DoUpdate                //force an update of the graph before making a copy of it for the report
859       
860                W_GenerateReport(funcStr,folderStr,$parStr,cw,yesSave,V_chisq,W_sigma,V_npnts,V_FitError,V_FitQuitReason,V_startRow,V_endRow,topGraph)
861        endif
862       
863        SetDataFolder root:
864        return(0)
865End
866
867
868// right now, there are only unsmeared models, but I want experimental points
869// for the QxQy of the calculation - so get the folder string
870Function Plot2DFunctionButtonProc(ba) : ButtonControl
871        STRUCT WMButtonAction &ba
872
873        String folderStr,funcStr,coefStr,cmdStr=""
874       
875        Variable killWhat=0             //kill nothing as default
876
877        switch( ba.eventCode )
878                case 2: // mouse up
879                        ControlInfo/W=WrapperPanel popup_0
880                        folderStr=S_Value
881                       
882                        ControlInfo/W=WrapperPanel popup_1
883                        funcStr=S_Value
884                       
885                        // check for smeared or smeared function
886//                      if(stringmatch(funcStr, "Smear*" )==1)
887                                //it's a smeared model
888                                // check for the special case of RPA that has an extra parameter
889//                              if(strsearch(funcStr, "RPAForm", 0 ,0) == -1)
890                                        sprintf cmdStr, "Plot%s(\"%s\")",funcStr,folderStr              //not RPA
891//                              else
892//                                      sprintf cmdStr, "Plot%s(\"%s\",)",funcStr,folderStr             //yes RPA, leave a comma for input
893//                              endif
894//                      else
895//                              // it's not,                   
896//                              sprintf cmdStr, "Plot%s()",funcStr
897//                      endif
898                       
899                        //Print cmdStr
900                        Execute cmdStr
901                       
902                        //pop the function menu to set the proper coefficients
903                        DoWindow/F WrapperPanel
904                        STRUCT WMPopupAction pa
905                        pa.popStr = funcStr
906                        pa.eventcode = 2
907                        Function_PopMenuProc(pa)
908       
909                        KillWhat = 2            //kill just the table, leave the 2d visible for now
910                        KillTopGraphAndTable(killWhat)          // crude
911
912                        break
913        endswitch
914
915        return 0
916End
917
918// unused - all of this functionality has been added to the Wrapper Panel
919Window Plot_2D_Controls()
920        PauseUpdate; Silent 1           // building window...
921        NewPanel /W=(950,44,1250,244)
922        Button button0,pos={172,37},size={70,20},proc=LogToggle2DButtonProc,title="Log/Lin"
923        Button button1,pos={146,9},size={100,20},proc=Plot2DButtonProc,title="Plot 2D Data"
924        Button button2,pos={164,118},size={120,20},proc=Append2DModelButtonProc,title="Append Model"
925        PopupMenu popup_1,pos={9,84},size={132,20},proc=Function_PopMenuProc,title="Function"
926        PopupMenu popup_1,mode=1,popvalue="Cylinder2D",value= #"W_FunctionPopupList()"
927        Button button3,pos={11,9},size={100,20},proc=Load2DDataButtonProc,title="Load 2D Data"
928        Button button4,pos={164,161},size={120,20},proc=Do2DFitButtonProc,title="Do 2D Fit"
929        Button button5,pos={164,85},size={120,20},proc=Plot2DFunctionButtonProc,title="Plot 2D Function"
930EndMacro
931
932
933Function/S TopGizmoWindow()
934        Return(StringFromList(0,WinList("*",";","WIN:4096")))
935end
936
937////
938//
939// unused proc for testing
940Proc Plot2DVsPointNumber(str)
941        String str
942        Prompt str,"Pick the data folder containing 2D data",popup,getAList(4)
943
944        PauseUpdate; Silent 1           // building window...
945        String fldrSav0= GetDataFolder(1)
946        SetDataFolder $("root:"+str)
947        Display /W=(241,352,810,842) $(str+"_i")
948        SetDataFolder fldrSav0
949        ModifyGraph mode=3
950        ModifyGraph marker=19
951        ModifyGraph lSize=2
952        ModifyGraph rgb($(str+"_i"))=(0,0,0)
953        ModifyGraph msize=1
954        ModifyGraph grid=1
955        ModifyGraph log(left)=1
956        ModifyGraph mirror=2
957        if(exists("root:"+str+":sw"))
958                ErrorBars/T=0 $(str+"_i") Y,wave=($("root:"+str+":sw"),$("root:"+str+":sw"))
959        endif
960        SetDataFolder fldrSav0
961EndMacro
962
963// testing procedure, used when I had x,y swapped and was getting nonsensical results
964//Macro CalculateChiSquared(str)
965//      String str
966//      Prompt str,"Pick the data folder containing 2D data",popup,getAList(4)
967//
968//
969//      String fldrSav0= GetDataFolder(1)
970//      SetDataFolder $("root:"+str)
971//
972//      Duplicate/O $(str+"_i") chi
973//      chi = ((zwave_cyl2D_D - $(str+"_i"))/sw )^2
974//     
975//      chi = (mask == 1) ? chi : 0
976//     
977//      Print sum(chi,-inf,inf)
978//     
979//      SetDataFolder fldrSav0
980//EndMacro
981
982
983Function MakeBSMask(mask,rad)
984        Wave mask
985        Variable rad
986
987
988// find the center based on the wave scaling
989        Variable xCtr,yCtr,Qzero=0
990        xCtr = (Qzero - DimOffset(Mask, 0))/DimDelta(Mask,0)
991        yCtr = (Qzero - DimOffset(Mask, 1))/DimDelta(Mask,1)
992
993        Print xctr,yctr
994
995//      Variable center = sqrt(numpnts(mask))/2 -0.5
996
997        mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1
998
999        Variable xDim, yDim
1000        xDim = DimSize(mask,0)
1001        yDim = DimSize(mask,1)
1002        mask[][0] = 0
1003        mask[][1] = 0
1004        mask[][yDim-2] = 0
1005        mask[][yDim-1] = 0
1006        mask[0][] = 0
1007        mask[1][] = 0
1008        mask[xDim-2][] = 0
1009        mask[xDim-1][] = 0
1010       
1011        Redimension/N=(xDim*yDim) mask          //now 1D
1012       
1013       
1014End
1015
1016// This routine assumes that the 2D data was loaded with the NCNR loader, so that the
1017// data is in a data folder, and the extensions are known. A more generic form could
1018// be made too, if needed.
1019//
1020// X- need error on I(q)
1021// -- need to set "proper" number of data points (delta and qMax?)
1022// X- need to remove points at high Q end
1023//
1024// -- like the routines in CircSectAve, start with 500 points, and trim after binning is done.
1025//      you'l end up with < 200 points.
1026//
1027// the results are in iBin_qxqy, qBin_qxqy, and eBin_qxqy, in the folder passed
1028//
1029//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
1030Function fDoBinning_QxQy2D(folderStr)
1031        String folderStr
1032
1033//      Wave inten,qx,qy,qz
1034
1035        SetDataFolder $("root:"+folderStr)
1036       
1037        WAVE inten = $(folderStr + "_i")
1038        WAVE qx = $(folderStr + "_qx")
1039        WAVE qy = $(folderStr + "_qy")
1040        WAVE qz = $(folderStr + "_qz")
1041       
1042        Variable xDim=numpnts(qx),yDim
1043        Variable ii,jj,delQ
1044        Variable qTot,nq,var,avesq,aveisq
1045        Variable binIndex,val
1046       
1047        nq = 500
1048       
1049        yDim = XDim
1050        Make/O/D/N=(nq) iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy
1051        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
1052        qBin_qxqy[] =  p*       delQ   
1053        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1054
1055        iBin_qxqy = 0
1056        iBin2_qxqy = 0
1057        eBin_qxqy = 0
1058        nBin_qxqy = 0   //number of intensities added to each bin
1059       
1060        for(ii=0;ii<xDim;ii+=1)
1061                qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1062                binIndex = trunc(x2pnt(qBin_qxqy, qTot))
1063                val = inten[ii]
1064                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
1065                        iBin_qxqy[binIndex] += val
1066                        iBin2_qxqy[binIndex] += val*val
1067                        nBin_qxqy[binIndex] += 1
1068                endif
1069        endfor
1070
1071//calculate errors, just like in CircSectAve.ipf
1072        for(ii=0;ii<nq;ii+=1)
1073                if(nBin_qxqy[ii] == 0)
1074                        //no pixels in annuli, data unknown
1075                        iBin_qxqy[ii] = 0
1076                        eBin_qxqy[ii] = 1
1077                else
1078                        if(nBin_qxqy[ii] <= 1)
1079                                //need more than one pixel to determine error
1080                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1081                                eBin_qxqy[ii] = 1
1082                        else
1083                                //assume that the intensity in each pixel in annuli is normally
1084                                // distributed about mean...
1085                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1086                                avesq = iBin_qxqy[ii]^2
1087                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1088                                var = aveisq-avesq
1089                                if(var<=0)
1090                                        eBin_qxqy[ii] = 1e-6
1091                                else
1092                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1093                                endif
1094                        endif
1095                endif
1096        endfor
1097       
1098        // find the last non-zero point, working backwards
1099        val=nq
1100        do
1101                val -= 1
1102        while(nBin_qxqy[val] == 0)
1103       
1104//      print val, nBin_qxqy[val]
1105        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy
1106       
1107        SetDataFolder root:
1108       
1109        return(0)
1110End
Note: See TracBrowser for help on using the repository browser.