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

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

added fDoBinning_QxQy2D(folderStr) that takes the QxQyQz? data and bins it to I(Q) (PlotUtils2D)

added comments to Sphere_2D describing how the 2D smearing calculation is different than the 1D calculation

the trac ticket page button on the USANS panel now points to the correct function

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