source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/RectAnnulAvg.ipf @ 740

Last change on this file since 740 was 665, checked in by srkline, 13 years ago

Made preferences a common panel (moved to PlotUtilsMacro?.ipf and globals to root:Packages:NIST:) and added menu items for all packages. Many files had to be modified so that the preferences could be properly accessed

File Open dialog now is set to "All files" so that XML can be selected. I think that all open access that doesn't already have the full path go through this common function.

File size: 20.0 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
4
5//***********************
6// Vers. 1.2 092101
7//
8// functions to perform either ractangular averages (similar to sector averages)
9// or annular averages ( fixed Q, I(angle) )
10//
11// dispatched to this point by ExecuteProtocol()
12//
13//**************************
14
15////////////////////////////////////
16//
17//              For AVERAGE and for DRAWING
18//                      DRAWING routines only use a subset of the total list, since saving, naming, etc. don't apply
19//              (10) possible keywords, some numerical, some string values
20//              AVTYPE=string           string from set {Circular,Annular,Rectangular,Sector,2D_ASCII,PNG_Graphic}
21//              PHI=value                       azimuthal angle (-90,90)
22//              DPHI=value                      +/- angular range around phi for average
23//              WIDTH=value             total width of rectangular section, in pixels
24//              SIDE=string             string from set {left,right,both} **note NOT capitalized
25//              QCENTER=value           q-value (1/A) of center of annulus for annular average
26//              QDELTA=value            total width of annulus centered at QCENTER
27//              PLOT=string             string from set {Yes,No} = truth of generating plot of averaged data
28//              SAVE=string             string from set {Yes,No} = truth of saving averaged data to disk
29//              NAME=string             string from set {Auto,Manual} = Automatic name generation or Manual(dialog)
30//
31//////////////////////////////////
32
33
34//function to do average of a rectangular swath of the detector
35//a sector average seems to be more appropriate, but there may be some
36//utility in rectangular averages
37//the parameters in the global keyword-string must have already been set somewhere
38//either directly or from the protocol
39//
40// 2-D data in the folder must already be on a linear scale. The calling routine is
41//responsible for this -
42//writes out the averaged waves to the "type" data folder
43//data is not written to disk by this routine
44//
45Function RectangularAverageTo1D(type)
46        String type
47       
48        SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr
49       
50        //type is the data type to do the averaging on, and will be set as the current folder
51        //get the current displayed data (so the correct folder is used)
52        String destPath = "root:Packages:NIST:"+type
53        //
54        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist,dr,ddr
55        Variable lambda,trans
56        Wave reals = $(destPath + ":RealsRead")
57        WAVE/T textread = $(destPath + ":TextRead")
58        String fileStr = textread[3]
59       
60        // center of detector, for non-linear corrections
61        NVAR pixelsX = root:myGlobals:gNPixelsX
62        NVAR pixelsY = root:myGlobals:gNPixelsY
63       
64        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
65        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
66       
67        // beam center, in pixels
68        x0 = reals[16]
69        y0 = reals[17]
70        //detector calibration constants
71        sx = reals[10]          //mm/pixel (x)
72        sx3 = reals[11]         //nonlinear coeff
73        sy = reals[13]          //mm/pixel (y)
74        sy3 = reals[14]         //nonlinear coeff
75
76       
77        dtsize = 10*reals[20]           //det size in mm
78        dtdist = 1000*reals[18] // det distance in mm
79       
80        NVAR pixelsX = root:myGlobals:gNPixelsX
81        NVAR pixelsY = root:myGlobals:gNPixelsY
82       
83        NVAR binWidth=root:Packages:NIST:gBinWidth
84        dr = binWidth           // annulus width set by user, default is one
85        ddr = dr*sx             //step size, in mm (this is the value to pass to the resolution calculation, not dr 18NOV03)
86               
87        Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,dxbm,dybm,ii
88        Variable phi_rad,dphi_rad,phi_x,phi_y
89        Variable forward,mirror
90       
91        String side = StringByKey("SIDE",keyListStr,"=",";")
92//      Print "side = ",side
93
94        //convert from degrees to radians
95        phi_rad = (Pi/180)*NumberByKey("PHI",keyListStr,"=",";")
96        dphi_rad = (Pi/180)*NumberByKey("DPHI",keyListStr,"=",";")
97        //create cartesian values for unit vector in phi direction
98        phi_x = cos(phi_rad)
99        phi_y = sin(phi_rad)
100       
101        //get (total) width of band
102        Variable width = NumberByKey("WIDTH",keyListStr,"=",";")
103
104        /// data wave is data in the current folder which was set at the top of the function
105        Wave data=$(destPath + ":data")
106        //Check for the existence of the mask, if not, make one (local to this folder) that is null
107       
108        if(WaveExists($"root:Packages:NIST:MSK:data") == 0)
109                Print "There is no mask file loaded (WaveExists)- the data is not masked"
110                Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
111                WAVE mask = $(destPath + ":mask")
112                mask = 0
113        else
114                Wave mask=$"root:Packages:NIST:MSK:data"
115        Endif
116       
117        rcentr = 100            //pixels within rcentr of beam center are broken into 9 parts
118        // values for error if unable to estimate value
119        //large_num = 1e10
120        large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1
121        small_num = 1e-10
122       
123        // output wave are expected to exist (?) initialized to zero, what length?
124        // 300 points on VAX ---
125        Variable wavePts=500
126        Make/O/N=(wavePts) $(destPath + ":qval"),$(destPath + ":aveint")
127        Make/O/N=(wavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave")
128        Make/O/N=(wavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar")
129        WAVE qval = $(destPath + ":qval")
130        WAVE aveint = $(destPath + ":aveint")
131        WAVE ncells = $(destPath + ":ncells")
132        WAVE dsq = $(destPath + ":dsq")
133        WAVE sigave = $(destPath + ":sigave")
134        WAVE qbar = $(destPath + ":QBar")
135        WAVE sigmaq = $(destPath + ":SigmaQ")
136        WAVE fsubs = $(destPath + ":fSubS")
137
138        qval = 0
139        aveint = 0
140        ncells = 0
141        dsq = 0
142        sigave = 0
143        qbar = 0
144        sigmaq = 0
145        fsubs = 0
146
147        dtdis2 = dtdist^2
148        nq = 1
149        xoffst=0
150        //distance of beam center from detector center
151        dxbm = FX(x0,sx3,xcenter,sx)
152        dybm = FY(y0,sy3,ycenter,sy)
153               
154        //BEGIN AVERAGE **********
155        Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1
156        Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p,d_per,d_pll
157        Make/O/N=2 $(destPath + ":par")
158        WAVE par = $(destPath + ":par")
159       
160        // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
161        // loop index corresponds to FORTRAN (old code)
162        // and the IGOR array indices must be adjusted (-1) to the correct address
163        ii=1
164        do
165                xi = ii
166                dxi = FX(xi,sx3,xcenter,sx)
167                dx = dxi-dxbm           //dx and dy are in mm
168               
169                jj = 1
170                do
171                        data_pixel = data[ii-1][jj-1]           //assign to local variable
172                        yj = jj
173                        dyj = FY(yj,sy3,ycenter,sy)
174                        dy = dyj - dybm
175                        if(!(mask[ii][jj]))                     //masked pixels = 1, skip if masked (this way works...)
176                                dr2 = (dx^2 + dy^2)^(0.5)               //distance from beam center NOTE dr2 used here - dr used above
177                                if(dr2>rcentr)          //keep pixel whole
178                                        nd = 1
179                                        fd = 1
180                                else                            //break pixel into 9 equal parts
181                                        nd = 3
182                                        fd = 2
183                                endif
184                                nd2 = nd^2
185                                ll = 1          //"el-el" loop index
186                                do
187                                        dxx = dx + (ll - fd)*sx/3
188                                        kk = 1
189                                        do
190                                                dyy = dy + (kk - fd)*sy/3
191                                                //determine distance pixel is from beam center (d_pll)
192                                                //and distance off-line (d_per) and if in forward direction
193                                                par = 0                 //initialize the wave
194                                                forward = distance(dxx,dyy,phi_x,phi_y,par)
195                                                d_per = par[0]
196                                                d_pll = par[1]
197                                                //check whether pixel lies within width band
198                                                if(d_per <= (0.5*width*ddr))
199                                                        //check if pixel lies within allowed sector(s)
200                                                        if(cmpstr(side,"both")==0)              //both sectors
201                                                                        //increment
202                                                                        nq = IncrementPixel_Rec(data_pixel,ddr,d_pll,aveint,dsq,ncells,nq,nd2)
203                                                        else
204                                                                if(cmpstr(side,"right")==0)             //forward sector only
205                                                                        if(forward)
206                                                                                //increment
207                                                                                nq = IncrementPixel_Rec(data_pixel,ddr,d_pll,aveint,dsq,ncells,nq,nd2)
208                                                                        Endif
209                                                                else                    //mirror sector only
210                                                                        if(!forward)
211                                                                                //increment
212                                                                                nq = IncrementPixel_Rec(data_pixel,ddr,d_pll,aveint,dsq,ncells,nq,nd2)
213                                                                        Endif
214                                                                Endif
215                                                        Endif           //allowable sectors
216                                                Endif           //check if in band
217                                                kk+=1
218                                        while(kk<=nd)
219                                        ll += 1
220                                while(ll<=nd)
221                        Endif           //masked pixel check
222                        jj += 1
223                while (jj<=pixelsY)
224                ii += 1
225        while(ii<=pixelsX)              //end of the averaging
226               
227        //compute q-values and errors
228        Variable ntotal,rr,theta,avesq,aveisq,var
229       
230        lambda = reals[26]
231        ntotal = 0
232        kk = 1
233        do
234                rr = (2*kk-1)*ddr/2
235                theta = 0.5*atan(rr/dtdist)
236                qval[kk-1] = (4*Pi/lambda)*sin(theta)
237                if(ncells[kk-1] == 0)
238                        //no pixels in annuli, data unknown
239                        aveint[kk-1] = 0
240                        sigave[kk-1] = large_num
241                else
242                        if(ncells[kk-1] <= 1)
243                                //need more than one pixel to determine error
244                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
245                                sigave[kk-1] = large_num
246                        else
247                                //assume that the intensity in each pixel in annuli is normally
248                                // distributed about mean...
249                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
250                                avesq = aveint[kk-1]^2
251                                aveisq = dsq[kk-1]/ncells[kk-1]
252                                var = aveisq-avesq
253                                if(var<=0)
254                                        sigave[kk-1] = small_num
255                                else
256                                        sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1))
257                                endif
258                        endif
259                endif
260                ntotal += ncells[kk-1]
261                kk+=1
262        while(kk<=nq)
263       
264        //Print "NQ = ",nq
265        // data waves were defined as 200 points (=wavePts), but now have less than that (nq) points
266        // use DeletePoints to remove junk from end of waves
267        //WaveStats would be a more foolproof implementation, to get the # points in the wave
268        Variable startElement,numElements
269        startElement = nq
270        numElements = wavePts - startElement
271        DeletePoints startElement,numElements, qval,aveint,ncells,dsq,sigave
272       
273        //////////////end of VAX sector_ave()
274               
275        //angle dependent transmission correction
276        Variable uval,arg,cos_th
277        lambda = reals[26]
278        trans = reals[4]
279//
280//  The transmission correction is now done at the ADD step, in DetCorr()
281//             
282//      ////this section is the trans_correct() VAX routine
283//      if(trans<0.1)
284//              Print "***transmission is less than 0.1*** and is a significant correction"
285//      endif
286//      if(trans==0)
287//              Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
288//              trans = 1
289//      endif
290//      //optical thickness
291//      uval = -ln(trans)
292//      //apply correction to aveint[]
293//      //index from zero here, since only working with IGOR waves
294//      ii=0
295//      do
296//              theta = 2*asin(lambda*qval[ii]/(4*pi))
297//              cos_th = cos(theta)
298//              arg = (1-cos_th)/cos_th
299//              if((uval<0.01) || (cos_th>0.99))                //OR
300//                      //small arg, approx correction
301//                      aveint[ii] /= 1-0.5*uval*arg
302//              else
303//                      //large arg, exact correction
304//                      aveint[ii] /= (1-exp(-uval*arg))/(uval*arg)
305//              endif
306//              ii+=1
307//      while(ii<nq)
308//      //end of transmission/pathlength correction
309//     
310// ***************************************************************
311//
312// Do the extra 3 columns of resolution calculations starting here.
313//
314// ***************************************************************
315
316        Variable L2 = reals[18]
317        Variable BS = reals[21]
318        Variable S1 = reals[23]
319        Variable S2 = reals[24]
320        Variable L1 = reals[25]
321        lambda = reals[26]
322        Variable lambdaWidth = reals[27]
323        String detStr=textRead[9]
324       
325        Variable usingLenses = reals[28]                //new 2007
326
327        //Two parameters DDET and APOFF are instrument dependent.  Determine
328        //these from the instrument name in the header.
329        //From conversation with JB on 01.06.99 these are the current
330        //good values
331
332        Variable DDet
333        NVAR apOff = root:myGlobals:apOff               //in cm
334       
335//      DDet = DetectorPixelResolution(fileStr,detStr)          //needs detector type and beamline
336        //note that reading the detector pixel size from the header ASSUMES SQUARE PIXELS! - Jan2008
337        DDet = reals[10]/10                     // header value (X) is in mm, want cm here
338       
339        //Width of annulus used for the average is gotten from the
340        //input dialog before.  This also must be passed to the resol
341        //calculator. Currently the default is dr=1 so just keeping that.
342
343        //Go from 0 to nq doing the calc for all three values at
344        //every Q value
345
346        ii=0
347
348        Variable ret1,ret2,ret3
349        do
350                getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,usingLenses,ret1,ret2,ret3)
351                sigmaq[ii] = ret1
352                qbar[ii] = ret2
353                fsubs[ii] = ret3
354                ii+=1
355        while(ii<nq)
356
357        DeletePoints startElement,numElements, sigmaq, qbar, fsubs
358
359// ***************************************************************
360//
361// End of resolution calculations
362//
363// ***************************************************************
364
365        Avg_1D_Graph(aveint,qval,sigave)
366
367        //get rid of the default mask, if one was created (it is in the current folder)
368        //don't just kill "mask" since it might be pointing to the one in the MSK folder
369        Killwaves/Z $(destPath+":mask")
370       
371        KillWaves/Z $(destPath+":par")          //parameter wave used in function distance()
372       
373        //return to root folder (redundant)
374        SetDataFolder root:
375       
376        Return 0
377End
378
379//returns nq, new number of q-values
380//arrays aveint,dsq,ncells are also changed by this function
381//
382Function IncrementPixel_Rec(dataPixel,ddr,d_pll,aveint,dsq,ncells,nq,nd2)
383        Variable dataPixel,ddr,d_pll
384        Wave aveint,dsq,ncells
385        Variable nq,nd2
386       
387        Variable ir
388       
389        ir = trunc(abs(d_pll)/ddr)+1
390        if (ir>nq)
391                nq = ir         //resets maximum number of q-values
392        endif
393        aveint[ir-1] += dataPixel/nd2           //ir-1 must be used, since ir is physical
394        dsq[ir-1] += dataPixel*dataPixel/nd2
395        ncells[ir-1] += 1/nd2
396       
397        Return nq
398End
399
400//function determines disatnce in mm  that pixel is from line
401//intersecting cetner of detector and direction phi
402//at chosen azimuthal angle phi -> [cos(phi),sin(phi0] = [phi_x,phi_y]
403//distance is always positive
404//
405// distances are returned in  a wave
406// forward (truth) is the function return value
407//
408Function distance(dxx,dyy,phi_x,phi_y,par)
409        Variable dxx,dyy,phi_x,phi_y
410        Wave par                //par[0] = d_per
411                                        //par[1] = d_pll        , both are returned values
412       
413        Variable val,rr,dot_prod,forward,d_per,d_pll,dphi_pixel
414       
415        rr = sqrt(dxx^2 + dyy^2)
416        dot_prod = (dxx*phi_x + dyy*phi_y)/rr
417        if(dot_prod >= 0)
418                forward = 1
419        else
420                forward = 0
421        Endif
422        //? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
423        if(dot_prod > 1)
424                dot_prod =1
425        Endif
426        if(dot_prod < -1)
427                dot_prod = -1
428        Endif
429        dphi_pixel = acos(dot_prod)
430       
431        //distance (in mm) that pixel is from  line (perpendicular)
432        d_per = sin(dphi_pixel)*rr
433        //distance (in mm) that pixel projected onto line is from beam center (parallel)
434        d_pll = cos(dphi_pixel)*rr
435       
436        //assign to wave for return
437        par[0] = d_per
438        par[1] = d_pll
439       
440        return (forward)
441
442End
443
444//performs an average around an annulus of specified width, centered on a
445//specified q-value (Intensity vs. angle)
446//the parameters in the global keyword-string must have already been set somewhere
447//either directly or from the protocol
448//
449//the input (data in the "type" folder) must be on linear scale - the calling routine is
450//responsible for this
451//averaged data is written to the data folder and plotted. data is not written
452//to disk from this routine.
453//
454Function AnnularAverageTo1D(type)
455        String type
456       
457        SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr
458       
459        //type is the data type to do the averaging on, and will be set as the current folder
460        //get the current displayed data (so the correct folder is used)
461        String destPath = "root:Packages:NIST:"+type
462       
463        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist
464        Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,xbm,ybm,ii
465        Variable rc,delr,rlo,rhi,dphi,nphi,dr
466        Variable lambda,trans
467        Wave reals = $(destPath + ":RealsRead")
468
469        // center of detector, for non-linear corrections
470        NVAR pixelsX = root:myGlobals:gNPixelsX
471        NVAR pixelsY = root:myGlobals:gNPixelsY
472       
473        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
474        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
475       
476        // beam center, in pixels
477        x0 = reals[16]
478        y0 = reals[17]
479        //detector calibration constants
480        sx = reals[10]          //mm/pixel (x)
481        sx3 = reals[11]         //nonlinear coeff
482        sy = reals[13]          //mm/pixel (y)
483        sy3 = reals[14]         //nonlinear coeff
484       
485        dtsize = 10*reals[20]           //det size in mm
486        dtdist = 1000*reals[18] // det distance in mm
487        lambda = reals[26]
488       
489        Variable qc = NumberByKey("QCENTER",keyListStr,"=",";")
490        Variable nw = NumberByKey("QDELTA",keyListStr,"=",";")
491       
492        dr = 1                  //minimum annulus width, keep this fixed at one
493        NVAR numPhiSteps = root:Packages:NIST:gNPhiSteps
494        nphi = numPhiSteps              //number of anular sectors is set by users
495       
496        rc = 2*dtdist*asin(qc*lambda/4/Pi)              //in mm
497        delr = nw*sx/2
498        rlo = rc-delr
499        rhi = rc + delr
500        dphi = 360/nphi
501
502        /// data wave is data in the current folder which was set at the top of the function
503        Wave data=$(destPath + ":data")
504        //Check for the existence of the mask, if not, make one (local to this folder) that is null
505       
506        if(WaveExists($"root:Packages:NIST:MSK:data") == 0)
507                Print "There is no mask file loaded (WaveExists)- the data is not masked"
508                Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
509                WAVE mask = $(destPath + ":mask")
510                mask = 0
511        else
512                Wave mask=$"root:Packages:NIST:MSK:data"
513        Endif
514       
515        rcentr = 150            //pixels within rcentr of beam center are broken into 9 parts
516        // values for error if unable to estimate value
517        //large_num = 1e10
518        large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1
519        small_num = 1e-10
520       
521        // output wave are expected to exist (?) initialized to zero, what length?
522        // 300 points on VAX ---
523        Variable wavePts=500
524        Make/O/N=(wavePts) $(destPath + ":phival"),$(destPath + ":aveint")
525        Make/O/N=(wavePts) $(destPath + ":ncells"),$(destPath + ":sig"),$(destPath + ":sigave")
526        WAVE phival = $(destPath + ":phival")
527        WAVE aveint = $(destPath + ":aveint")
528        WAVE ncells = $(destPath + ":ncells")
529        WAVE sig = $(destPath + ":sig")
530        WAVE sigave = $(destPath + ":sigave")
531
532        phival = 0
533        aveint = 0
534        ncells = 0
535        sig = 0
536        sigave = 0
537
538        dtdis2 = dtdist^2
539        nq = 1
540        xoffst=0
541        //distance of beam center from detector center
542        xbm = FX(x0,sx3,xcenter,sx)
543        ybm = FY(y0,sy3,ycenter,sy)
544               
545        //BEGIN AVERAGE **********
546        Variable xi,xd,x,y,yd,yj,nd,fd,nd2,iphi,ntotal,var
547        Variable jj,data_pixel,xx,yy,ll,kk,rij,phiij,avesq,aveisq
548
549        // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
550        // loop index corresponds to FORTRAN (old code)
551        // and the IGOR array indices must be adjusted (-1) to the correct address
552        ntotal = 0
553        ii=1
554        do
555                xi = ii
556                xd = FX(xi,sx3,xcenter,sx)
557                x = xoffst + xd -xbm            //x and y are in mm
558               
559                jj = 1
560                do
561                        data_pixel = data[ii-1][jj-1]           //assign to local variable
562                        yj = jj
563                        yd = FY(yj,sy3,ycenter,sy)
564                        y = yd - ybm
565                        if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...)
566                                nd = 1
567                                fd = 1
568                                if( (abs(x) > rcentr) || (abs(y) > rcentr))     //break pixel into 9 equal parts
569                                        nd = 3
570                                        fd = 2
571                                Endif
572                                nd2 = nd^2
573                                ll = 1          //"el-el" loop index
574                                do
575                                        xx = x + (ll - fd)*sx/3
576                                        kk = 1
577                                        do
578                                                yy = y + (kk - fd)*sy/3
579                                                //test to see if center of pixel (i,j) lies in annulus
580                                                rij = sqrt(x*x + y*y)/dr + 1.001
581                                                //check whether pixel lies within width band
582                                                if((rij > rlo) && (rij < rhi))
583                                                        //in the annulus, do something
584                                                        if (yy >= 0)
585                                                                //phiij is in degrees
586                                                                phiij = atan2(yy,xx)*180/Pi             //0 to 180 deg
587                                                        else
588                                                                phiij = 360 + atan2(yy,xx)*180/Pi               //180 to 360 deg
589                                                        Endif
590                                                        if (phiij > (360-0.5*dphi))
591                                                                phiij -= 360
592                                                        Endif
593                                                        iphi = trunc(phiij/dphi + 1.501)
594                                                        aveint[iphi-1] += 9*data_pixel/nd2
595                                                        sig[iphi-1] += 9*data_pixel*data_pixel/nd2
596                                                        ncells[iphi-1] += 9/nd2
597                                                        ntotal += 9/nd2
598                                                Endif           //check if in annulus
599                                                kk+=1
600                                        while(kk<=nd)
601                                        ll += 1
602                                while(ll<=nd)
603                        Endif           //masked pixel check
604                        jj += 1
605                while (jj<=pixelsY)
606                ii += 1
607        while(ii<=pixelsX)              //end of the averaging
608               
609        //compute phi-values and errors
610       
611        ntotal /=9
612       
613        kk = 1
614        do
615                phival[kk-1] = dphi*(kk-1)
616                if(ncells[kk-1] != 0)
617                        aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
618                        avesq = aveint[kk-1]*aveint[kk-1]
619                        aveisq = sig[kk-1]/ncells[kk-1]
620                        var = aveisq - avesq
621                        if (var <=0 )
622                                sig[kk-1] = 0
623                                sigave[kk-1] = 0
624                                ncells[kk-1] /=9
625                        else
626                                if(ncells[kk-1] > 9)
627                                        sigave[kk-1] = sqrt(9*var/(ncells[kk-1]-9))
628                                        sig[kk-1] = sqrt( abs(aveint[kk-1])/(ncells[kk-1]/9) )
629                                        ncells[kk-1] /=9
630                                else
631                                        sig[kk-1] = 0
632                                        sigave[kk-1] = 0
633                                        ncells[kk-1] /=9
634                                Endif
635                        Endif
636                Endif
637                kk+=1
638        while(kk<=nphi)
639       
640        // data waves were defined as 200 points (=wavePts), but now have less than that (nphi) points
641        // use DeletePoints to remove junk from end of waves
642        Variable startElement,numElements
643        startElement = nphi
644        numElements = wavePts - startElement
645        DeletePoints startElement,numElements, phival,aveint,ncells,sig,sigave
646       
647        //////////////end of VAX Phibin.for
648               
649        //angle dependent transmission correction is not done in phiave
650        Ann_1D_Graph(aveint,phival,sigave)
651       
652        //get rid of the default mask, if one was created (it is in the current folder)
653        //don't just kill "mask" since it might be pointing to the one in the MSK folder
654        Killwaves/z $(destPath+":mask")
655               
656        //return to root folder (redundant)
657        SetDataFolder root:
658       
659        Return 0
660End
Note: See TracBrowser for help on using the repository browser.