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

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

bug fix for #306:

now the 2D QxQy? loader switches correctly for 3-column or 8-column data. 8-column is the default output, 3-column is only for generic cases of Qx-Qy-I files.

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