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

Last change on this file since 925 was 925, checked in by srkline, 9 years ago

Adding changes to the real space modeling that allow calculation of the USANS (slit-smeared) intensity from an anisotropic structure that has been generated using a real-space construction. The detector plane from the FFT result is converted to USANS, on absolute scale. Will be extended as it gets some use, but it is functional and is found from the "3D Examples" Macros submenu. Documentation is on the way too.

File size: 36.1 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=4.00
3#pragma IgorVersion=6.1
4
5//
6// TO DO:
7//
8// - more intelligent beam stop masking
9// - constraints
10// - interactive masking
11//
12//
13// - OCT 2010   - added fDoBinning_QxQy2D(folderStr) that takes the QxQyQz data and bins it to I(Q)
14//                                      1D result must still be plotted manually, can be automated later.
15//                                      - for Scaled Image data (QxQy), the function fDoBinning_Scaled2D (in FFT_Cubes) could be
16//                                              modified and added to this file. It was meant for an FFT slice, but applies to
17//                                              2D model calculations (model_lin) 2D data files as well.
18
19// 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//
22
23
24//
25// changed June 2010 to read in the resolution information (now 8! columns)
26// -- subject to change --
27//
28// look for either the old-style 3-column (no resolution information) or the newer 8-column format
29Proc LoadQxQy()
30
31        SetDataFolder root:
32       
33        LoadWave/G/D/W/A
34        String fileName = S_fileName
35        String path = S_Path
36        Variable numCols = V_flag
37
38        String w0,w1,w2,w3,w4,w5,w6,w7
39        String n0,n1,n2,n3,n4,n5,n6,n7
40               
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)
61       
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
100
101
102// the very bad case where users don't supply the error data...
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 ,";" )
108
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)
113                w3 = CleanupName((S_fileName + "_iErr"),0)              //make a name for the error wave, to be generated here
114
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
146       
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       
161        endif           //3-columns
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
209       
210        endif           //4-columns     
211       
212       
213        /// do this for all 2D data, whether or not resolution information was read in
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:
227        KillWaves/Z $n0,$n1,$n2,$n3,$n4,$n5,$n6,$n7
228       
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}
288        AppendToGizmo freeAxesCue={0,0,0,1.5},name=freeAxesCue0
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
292        ModifyGizmo setDisplayList=3, object=freeAxesCue0
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
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)
316        Variable minQx=-0.1,maxQx=0.1,minQy=-0.1,maxQy=0.1
317        Variable numPix=64
318        String dataFolder="fake2DData"
319        Variable forceOverwrite = 0
320        Variable doPlot = 1
321
322        String baseStr=dataFolder
323       
324        if(DataFolderExists("root:"+baseStr) && forceOverwrite)
325                SetDataFolder $("root:"+baseStr)
326        else
327                if(DataFolderExists("root:"+baseStr))
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)
334                else
335                        NewDataFolder/S $("root:"+baseStr)
336                endif
337        endif
338       
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       
365        $(baseStr+"_qx") = tmpX[mod(p,numPix)]
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
374        if(doPlot)     
375                PlotQxQy(baseStr)               //this sets the data folder back to root:, but may not be executed
376        endif
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)
394        NVAR gIsLogScale=$(folderStr+"gIsLogScale")
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
463// object 3 is the axisCue
464        ModifyGizmo setDisplayList=4, object=surface1
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
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
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
531                        suffix = getModelSuffix(funcStr)
532                        if(stringmatch(funcStr, "smear*") == 1)
533                                suffix = "sm_"+suffix
534                        endif
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)
569                                w = w[p][q] == -inf ? NaN : w[p][q]             //remove the -inf for display
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:
653// - I have no error waves for the intensity (fixed 10/2010)
654// - There is no smeared model (coming soon after 10/2010)
655// - Cursors can't be used
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//
659Function FitWrapper2D(folderStr,funcStr,coefStr,useCursors,useEps,useConstr)
660        String folderStr,funcStr,coefStr
661        Variable useCursors,useEps,useConstr
662
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       
668        String suffix=getModelSuffix(funcStr)
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
687// note that the resolution waves may or may not exist, and may or may not be used in the fitting
688        String DF="root:"+folderStr+":"
689       
690        WAVE inten=$(DF+folderStr+"_i")
691        WAVE sw=$(DF+folderStr+"_iErr")
692        WAVE qx=$(DF+folderStr+"_qx")
693        WAVE qy=$(DF+folderStr+"_qy")
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")
698
699//just a dummy - I shouldn't need this
700        Duplicate/O qx resultW
701        resultW=0
702       
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
708        WAVE/Z s.qz = qz
709        WAVE/Z s.sQpl = sQpl
710        WAVE/Z s.sQpp = sQpp
711        WAVE/Z s.fs = shad
712       
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)
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
722       
723       
724        Duplicate/O inten inten_masked
725        inten_masked = (mask[p][q] == 0) ? NaN : inten[p][q]
726       
727
728        //for now, use res will always be 0 for 2D functions   
729        Variable useResol=0
730        if(stringmatch(funcStr, "Smear*"))              // if it's a smeared function, need a struct
731                useResol=1
732        endif
733
734        // can't use constraints defined as a single text wave for multivariate fits. See the curve fitting help file
735        // and "Contraint Matrix and Vector"
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)
739       
740//      if(useConstr)
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
766
767        endif
768
769        if(useCursors)
770                Print "Cursors not yet implemented"
771                useCursors = 0
772        endif   
773///// NO CURSORS for 2D waves
774        //if useCursors, and the data is USANS, need to feed a (reassigned) trimmed matrix to the fit
775        Variable pt1,pt2,newN
776        pt1 = 0
777        pt2 = numpnts(inten)-1
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
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
801        Variable t1=StopMSTimer(-2)
802        Variable tb = 0         //no textbox
803
804// /NTHR=1 means just one thread for the fit (since the function itself is threaded)
805// NTHR = 0 == "Auto" mode, using as many processors as are available (not appropriate here since the function itself is threaded?)
806
807        do
808       
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
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
827               
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               
834                if(useResol)            //res only
835                        Print "useRes only"
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
837                        break
838                endif
839                               
840/////   same as above, but all without useResol (no /STRC flag)
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
856               
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
859                        break
860                Endif
861               
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
865               
866        while(0)
867       
868        Print "elapsed 2D fit time  = ",(StopMSTimer(-2) - t1)/1e6," seconds = ",(StopMSTimer(-2) - t1)/1e6/60," minutes"
869
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=""
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
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)
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
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       
942        Variable killWhat=0             //kill nothing as default
943
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       
976                        KillWhat = 2            //kill just the table, leave the 2d visible for now
977                        KillTopGraphAndTable(killWhat)          // crude
978
979                        break
980        endswitch
981
982        return 0
983End
984
985// unused - all of this functionality has been added to the Wrapper Panel
986Window Plot_2D_Controls()
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
999
1000Function/S TopGizmoWindow()
1001        Return(StringFromList(0,WinList("*",";","WIN:4096")))
1002end
1003
1004////
1005//
1006// unused proc for testing
1007Proc Plot2DVsPointNumber(str)
1008        String str
1009        Prompt str,"Pick the data folder containing 2D data",popup,getAList(4)
1010
1011        PauseUpdate; Silent 1           // building window...
1012        String fldrSav0= GetDataFolder(1)
1013        SetDataFolder $("root:"+str)
1014        Display /W=(241,352,810,842) $(str+"_i")
1015        SetDataFolder fldrSav0
1016        ModifyGraph mode=3
1017        ModifyGraph marker=19
1018        ModifyGraph lSize=2
1019        ModifyGraph rgb($(str+"_i"))=(0,0,0)
1020        ModifyGraph msize=1
1021        ModifyGraph grid=1
1022        ModifyGraph log(left)=1
1023        ModifyGraph mirror=2
1024        if(exists("root:"+str+":sw"))
1025                ErrorBars/T=0 $(str+"_i") Y,wave=($("root:"+str+":sw"),$("root:"+str+":sw"))
1026        endif
1027        SetDataFolder fldrSav0
1028EndMacro
1029
1030// testing procedure, used when I had x,y swapped and was getting nonsensical results
1031//Macro CalculateChiSquared(str)
1032//      String str
1033//      Prompt str,"Pick the data folder containing 2D data",popup,getAList(4)
1034//
1035//
1036//      String fldrSav0= GetDataFolder(1)
1037//      SetDataFolder $("root:"+str)
1038//
1039//      Duplicate/O $(str+"_i") chi
1040//      chi = ((zwave_cyl2D_D - $(str+"_i"))/sw )^2
1041//     
1042//      chi = (mask == 1) ? chi : 0
1043//     
1044//      Print sum(chi,-inf,inf)
1045//     
1046//      SetDataFolder fldrSav0
1047//EndMacro
1048
1049
1050Function MakeBSMask(mask,rad)
1051        Wave mask
1052        Variable rad
1053
1054
1055// find the center based on the wave scaling
1056        Variable xCtr,yCtr,Qzero=0
1057        xCtr = (Qzero - DimOffset(Mask, 0))/DimDelta(Mask,0)
1058        yCtr = (Qzero - DimOffset(Mask, 1))/DimDelta(Mask,1)
1059
1060        Print xctr,yctr
1061
1062//      Variable center = sqrt(numpnts(mask))/2 -0.5
1063
1064        mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1
1065
1066        Variable xDim, yDim
1067        xDim = DimSize(mask,0)
1068        yDim = DimSize(mask,1)
1069        mask[][0] = 0
1070        mask[][1] = 0
1071        mask[][yDim-2] = 0
1072        mask[][yDim-1] = 0
1073        mask[0][] = 0
1074        mask[1][] = 0
1075        mask[xDim-2][] = 0
1076        mask[xDim-1][] = 0
1077       
1078        Redimension/N=(xDim*yDim) mask          //now 1D
1079       
1080       
1081End
1082
1083// This routine assumes that the 2D data was loaded with the NCNR loader, so that the
1084// data is in a data folder, and the extensions are known. A more generic form could
1085// be made too, if needed.
1086//
1087// X- need error on I(q) - now calculated both ways, from 2D error propagation, and from deviation from the mean
1088// X- need to set "proper" number of data points (delta and qMax?)
1089// X- need to remove points at high Q end
1090//
1091// -- like the routines in CircSectAve, start with 500 points, and trim after binning is done.
1092//      you'l end up with < 200 points.
1093//
1094// the results are in iBin_qxqy, qBin_qxqy, and eBin_qxqy, in the folder passed
1095//
1096
1097Proc BinQxQy_to_1D(folderStr)
1098        String folderStr
1099        Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
1100
1101        fDoBinning_QxQy2D(folderStr)
1102       
1103        SetDataFolder $("root:"+folderStr)
1104        Display iBin_qxqy vs qBin_qxqy
1105        ModifyGraph mirror=2,grid=1,log=1
1106        ModifyGraph mode=4,marker=19,msize=2
1107        ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)
1108        legend
1109       
1110        SetDataFolder root:
1111End
1112
1113
1114
1115//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
1116//      Wave inten,qx,qy,qz
1117
1118Function fDoBinning_QxQy2D(folderStr)
1119        String folderStr
1120
1121
1122        SetDataFolder $("root:"+folderStr)
1123       
1124//      WAVE inten = $("smeared_sf2D")
1125
1126        WAVE inten = $(folderStr + "_i")
1127        WAVE iErr = $(folderStr + "_iErr")
1128        WAVE qx = $(folderStr + "_qx")
1129        WAVE qy = $(folderStr + "_qy")
1130        WAVE qz = $(folderStr + "_qz")
1131       
1132        Variable xDim=numpnts(qx),yDim
1133        Variable ii,jj,delQ
1134        Variable qTot,nq,var,avesq,aveisq
1135        Variable binIndex,val
1136       
1137        nq = 500
1138       
1139        yDim = XDim
1140        Make/O/D/N=(nq) iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1141        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
1142        qBin_qxqy[] =  p*       delQ   
1143        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1144
1145        iBin_qxqy = 0
1146        iBin2_qxqy = 0
1147        eBin_qxqy = 0
1148        eBin2D_qxqy = 0
1149        nBin_qxqy = 0   //number of intensities added to each bin
1150       
1151        for(ii=0;ii<xDim;ii+=1)
1152                qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1153                binIndex = trunc(x2pnt(qBin_qxqy, qTot))
1154                val = inten[ii]
1155                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
1156                        iBin_qxqy[binIndex] += val
1157                        iBin2_qxqy[binIndex] += val*val
1158                        eBin2D_qxqy[binIndex] += iErr[ii]*iErr[ii]
1159                        nBin_qxqy[binIndex] += 1
1160                endif
1161        endfor
1162
1163//calculate errors, just like in CircSectAve.ipf
1164        for(ii=0;ii<nq;ii+=1)
1165                if(nBin_qxqy[ii] == 0)
1166                        //no pixels in annuli, data unknown
1167                        iBin_qxqy[ii] = 0
1168                        eBin_qxqy[ii] = 1
1169                        eBin2D_qxqy[ii] = NaN
1170                else
1171                        if(nBin_qxqy[ii] <= 1)
1172                                //need more than one pixel to determine error
1173                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1174                                eBin_qxqy[ii] = 1
1175                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1176                        else
1177                                //assume that the intensity in each pixel in annuli is normally
1178                                // distributed about mean...
1179                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1180                                avesq = iBin_qxqy[ii]^2
1181                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1182                                var = aveisq-avesq
1183                                if(var<=0)
1184                                        eBin_qxqy[ii] = 1e-6
1185                                else
1186                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1187                                endif
1188                                // and calculate as it is propagated pixel-by-pixel
1189                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1190                        endif
1191                endif
1192        endfor
1193       
1194        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1195       
1196        // find the last non-zero point, working backwards
1197        val=nq
1198        do
1199                val -= 1
1200        while(nBin_qxqy[val] == 0)
1201       
1202//      print val, nBin_qxqy[val]
1203        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1204       
1205        SetDataFolder root:
1206       
1207        return(0)
1208End
Note: See TracBrowser for help on using the repository browser.