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

Last change on this file since 857 was 857, checked in by srkline, 10 years ago

Added support for reading 4-column data into the 2D reader

Calculation of kappa for ABS now does a specific byte check for ensuring a "good" DIV file. hopefully this will eliminate errors where incorrect DIV files are selected and bad kappa values are generated - and the program thinks (incorrectly) that there is a valid DVI file present, when there really isn't one (and then there is no easy recovery from this)

Added CheckIfDIVData() functions to NCNR_Utils and to FACILITY_Utils, and to other (facility)_Utils

Changed the calls to writing BT5 files after adding so that it could be more easily scripted in the future.

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