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

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

SA_Includes_v410 : now include Smear_2D

PeakGauss_2D, Sphere_2D : included threaded resolution smearing calculations for testing

DataSetHandling? : Included a quick and dirty batch converter for XML->6-col. See the top
of the file for the command to run

GaussUtils? : re-define the ResSemear_2D_AAOStruct. Relocated q-value and phi calculations from
RawWindowHook? to this file so they would be available for reduction and analysis

Smear_2D : now has a generic (non-threaded) smearing routine. Threading must be done in
individual functions since FUNCREF can't be passed to threads (plus a few other issues)

PlotUtils_2D : updated loader for new QxQy? columns. Fixes to Wrapper_2D to enable smeared fits

RawWindowHook? : removed Q-value calculation functions and moved these to GaussUtils?

WriteQIS : now writes out 8-columns for QxQy? data, defining the resolution
function in terms of directions parallel and perpendicular to Q. TEMPORARILY in the data
file an error in intensity is generated that is SQRT(I), being careful to
replace any NaN, inf, or zero with an average error value

MultiScatter_MonteCarlo_2D : 4-processor aware

NCNR_Utils : 2D resolution calculation is now in terms of parallel and perpendicular
rather than x and y. Gravity is included in the y-component

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