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

Last change on this file was 1001, checked in by srkline, 6 years ago

A number of simple fixes to make the macros compatible with Igor 7, and in some places back-compatible with Igor 6. There were only a few instances where the IgorVersion? was checked, so it does not caues a huge disruption.

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