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

Last change on this file since 794 was 794, checked in by srkline, 11 years ago

Lots of changes:
-2D resolution smearing
-2D error propagation

1) 2D resolution smearing has been corrected to use sigma (perp) correctly
rather than phi. This shows up in the quadrature loop in all of the 2D models
and in the Smear_2D "generic" function.

2) 1D resolutionn smearing is now corrected to account for the integration
range of +/- 3 sigma (99.73% of distribution). Now the result is divided by
0.9973 to rescale it to the correct value.

3) Smeared models are now AAO to improve speed and to allow easier use with
functions that are inherently AAO. No changes are needed, since the call is
behind the scenes, replacing Smear_Model_N() with Smear_Model_N_AAO().

4) in PlotUtils_2D added functions to re-bin the QxQy? data into a 1D format
BinQxQy_to_1D(). This also re-bins the errors in two ways, adding the per-pixel
errors in quadrature, or the deviation from the mean of the intensity. Manually
editing the intensity declaration allows 2D->1D binning of smeared models.

5) Per-pixel propagation of errors has been added through the entire data
reduction sequence. Data errors are generated on loading using Poisson
statistics (specifically tailored for accuracy at low counts) and then is
propagated through each manipulation of the data using standard error
propagation. The error matrix is linear_data_error. As a by-product, all of
the math operations on data are explicitly done on linear_data, to avoid
any potential mistakes of log/linear scaling. Results of this propagation
largely match J. Barker's /ERROR routines from the VAX, with some differences
at low data count values (as expected) and at higher count values near the
beam stop (somewhat unexpected). This per-pixel error is ouput in the QxQy_ASCII
data files. NO CHANGE has been made to the 1D data, which uses the deviation from
the mean as the error - since this is correct.

6) Added tables for the uncertainty in attenuator transmission (from JGB)

7) added two more REAL values to the VAX header to store information
necessary for error propagation. These are couting error that are part of
the transmission error and of the absolute scaling error. These add Read/Write?
functions in NCNR_DataReadWrite

The transmission error (one standard deviation) is written at byte 396 (4 bytes)

RealsRead?[41]

The Box Sum error (one standard deviation) is written at byte 400 (4 bytes)

RealsRead?[42]

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