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

Last change on this file since 961 was 961, checked in by srkline, 8 years ago

removed the function "distance" from RectAnnulAvg? to avoid potential conflicts with the same-named field in Nexus files.The function only occured locally in that one file.
File size: 20.1 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 = s_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                                        sigave[kk-1] = large_num                //if there are zero counts, make error large_num to be a warning flag
256                                else
257                                        sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1))
258                                endif
259                        endif
260                endif
261                ntotal += ncells[kk-1]
262                kk+=1
263        while(kk<=nq)
264       
265        //Print "NQ = ",nq
266        // data waves were defined as 200 points (=wavePts), but now have less than that (nq) points
267        // use DeletePoints to remove junk from end of waves
268        //WaveStats would be a more foolproof implementation, to get the # points in the wave
269        Variable startElement,numElements
270        startElement = nq
271        numElements = wavePts - startElement
272        DeletePoints startElement,numElements, qval,aveint,ncells,dsq,sigave
273       
274        //////////////end of VAX sector_ave()
275               
276        //angle dependent transmission correction
277        Variable uval,arg,cos_th
278        lambda = reals[26]
279        trans = reals[4]
280//
281//  The transmission correction is now done at the ADD step, in DetCorr()
282//             
283//      ////this section is the trans_correct() VAX routine
284//      if(trans<0.1)
285//              Print "***transmission is less than 0.1*** and is a significant correction"
286//      endif
287//      if(trans==0)
288//              Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
289//              trans = 1
290//      endif
291//      //optical thickness
292//      uval = -ln(trans)
293//      //apply correction to aveint[]
294//      //index from zero here, since only working with IGOR waves
295//      ii=0
296//      do
297//              theta = 2*asin(lambda*qval[ii]/(4*pi))
298//              cos_th = cos(theta)
299//              arg = (1-cos_th)/cos_th
300//              if((uval<0.01) || (cos_th>0.99))                //OR
301//                      //small arg, approx correction
302//                      aveint[ii] /= 1-0.5*uval*arg
303//              else
304//                      //large arg, exact correction
305//                      aveint[ii] /= (1-exp(-uval*arg))/(uval*arg)
306//              endif
307//              ii+=1
308//      while(ii<nq)
309//      //end of transmission/pathlength correction
310//     
311// ***************************************************************
312//
313// Do the extra 3 columns of resolution calculations starting here.
314//
315// ***************************************************************
316
317        Variable L2 = reals[18]
318        Variable BS = reals[21]
319        Variable S1 = reals[23]
320        Variable S2 = reals[24]
321        Variable L1 = reals[25]
322        lambda = reals[26]
323        Variable lambdaWidth = reals[27]
324        String detStr=textRead[9]
325       
326        Variable usingLenses = reals[28]                //new 2007
327
328        //Two parameters DDET and APOFF are instrument dependent.  Determine
329        //these from the instrument name in the header.
330        //From conversation with JB on 01.06.99 these are the current
331        //good values
332
333        Variable DDet
334        NVAR apOff = root:myGlobals:apOff               //in cm
335       
336//      DDet = DetectorPixelResolution(fileStr,detStr)          //needs detector type and beamline
337        //note that reading the detector pixel size from the header ASSUMES SQUARE PIXELS! - Jan2008
338        DDet = reals[10]/10                     // header value (X) is in mm, want cm here
339       
340        //Width of annulus used for the average is gotten from the
341        //input dialog before.  This also must be passed to the resol
342        //calculator. Currently the default is dr=1 so just keeping that.
343
344        //Go from 0 to nq doing the calc for all three values at
345        //every Q value
346
347        ii=0
348
349        Variable ret1,ret2,ret3
350        do
351                getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,usingLenses,ret1,ret2,ret3)
352                sigmaq[ii] = ret1
353                qbar[ii] = ret2
354                fsubs[ii] = ret3
355                ii+=1
356        while(ii<nq)
357
358        DeletePoints startElement,numElements, sigmaq, qbar, fsubs
359
360// ***************************************************************
361//
362// End of resolution calculations
363//
364// ***************************************************************
365
366        Avg_1D_Graph(aveint,qval,sigave)
367
368        //get rid of the default mask, if one was created (it is in the current folder)
369        //don't just kill "mask" since it might be pointing to the one in the MSK folder
370        Killwaves/Z $(destPath+":mask")
371       
372        KillWaves/Z $(destPath+":par")          //parameter wave used in function distance()
373       
374        //return to root folder (redundant)
375        SetDataFolder root:
376       
377        Return 0
378End
379
380//returns nq, new number of q-values
381//arrays aveint,dsq,ncells are also changed by this function
382//
383Function IncrementPixel_Rec(dataPixel,ddr,d_pll,aveint,dsq,ncells,nq,nd2)
384        Variable dataPixel,ddr,d_pll
385        Wave aveint,dsq,ncells
386        Variable nq,nd2
387       
388        Variable ir
389       
390        ir = trunc(abs(d_pll)/ddr)+1
391        if (ir>nq)
392                nq = ir         //resets maximum number of q-values
393        endif
394        aveint[ir-1] += dataPixel/nd2           //ir-1 must be used, since ir is physical
395        dsq[ir-1] += dataPixel*dataPixel/nd2
396        ncells[ir-1] += 1/nd2
397       
398        Return nq
399End
400
401//function determines disatnce in mm  that pixel is from line
402//intersecting cetner of detector and direction phi
403//at chosen azimuthal angle phi -> [cos(phi),sin(phi0] = [phi_x,phi_y]
404//distance is always positive
405//
406// distances are returned in  a wave
407// forward (truth) is the function return value
408//
409Function s_distance(dxx,dyy,phi_x,phi_y,par)
410        Variable dxx,dyy,phi_x,phi_y
411        Wave par                //par[0] = d_per
412                                        //par[1] = d_pll        , both are returned values
413       
414        Variable val,rr,dot_prod,forward,d_per,d_pll,dphi_pixel
415       
416        rr = sqrt(dxx^2 + dyy^2)
417        dot_prod = (dxx*phi_x + dyy*phi_y)/rr
418        if(dot_prod >= 0)
419                forward = 1
420        else
421                forward = 0
422        Endif
423        //? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
424        if(dot_prod > 1)
425                dot_prod =1
426        Endif
427        if(dot_prod < -1)
428                dot_prod = -1
429        Endif
430        dphi_pixel = acos(dot_prod)
431       
432        //distance (in mm) that pixel is from  line (perpendicular)
433        d_per = sin(dphi_pixel)*rr
434        //distance (in mm) that pixel projected onto line is from beam center (parallel)
435        d_pll = cos(dphi_pixel)*rr
436       
437        //assign to wave for return
438        par[0] = d_per
439        par[1] = d_pll
440       
441        return (forward)
442
443End
444
445//performs an average around an annulus of specified width, centered on a
446//specified q-value (Intensity vs. angle)
447//the parameters in the global keyword-string must have already been set somewhere
448//either directly or from the protocol
449//
450//the input (data in the "type" folder) must be on linear scale - the calling routine is
451//responsible for this
452//averaged data is written to the data folder and plotted. data is not written
453//to disk from this routine.
454//
455Function AnnularAverageTo1D(type)
456        String type
457       
458        SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr
459       
460        //type is the data type to do the averaging on, and will be set as the current folder
461        //get the current displayed data (so the correct folder is used)
462        String destPath = "root:Packages:NIST:"+type
463       
464        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist
465        Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,xbm,ybm,ii
466        Variable rc,delr,rlo,rhi,dphi,nphi,dr
467        Variable lambda,trans
468        Wave reals = $(destPath + ":RealsRead")
469
470        // center of detector, for non-linear corrections
471        NVAR pixelsX = root:myGlobals:gNPixelsX
472        NVAR pixelsY = root:myGlobals:gNPixelsY
473       
474        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
475        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
476       
477        // beam center, in pixels
478        x0 = reals[16]
479        y0 = reals[17]
480        //detector calibration constants
481        sx = reals[10]          //mm/pixel (x)
482        sx3 = reals[11]         //nonlinear coeff
483        sy = reals[13]          //mm/pixel (y)
484        sy3 = reals[14]         //nonlinear coeff
485       
486        dtsize = 10*reals[20]           //det size in mm
487        dtdist = 1000*reals[18] // det distance in mm
488        lambda = reals[26]
489       
490        Variable qc = NumberByKey("QCENTER",keyListStr,"=",";")
491        Variable nw = NumberByKey("QDELTA",keyListStr,"=",";")
492       
493        dr = 1                  //minimum annulus width, keep this fixed at one
494        NVAR numPhiSteps = root:Packages:NIST:gNPhiSteps
495        nphi = numPhiSteps              //number of anular sectors is set by users
496       
497        rc = 2*dtdist*asin(qc*lambda/4/Pi)              //in mm
498        delr = nw*sx/2
499        rlo = rc-delr
500        rhi = rc + delr
501        dphi = 360/nphi
502
503        /// data wave is data in the current folder which was set at the top of the function
504        Wave data=$(destPath + ":data")
505        //Check for the existence of the mask, if not, make one (local to this folder) that is null
506       
507        if(WaveExists($"root:Packages:NIST:MSK:data") == 0)
508                Print "There is no mask file loaded (WaveExists)- the data is not masked"
509                Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
510                WAVE mask = $(destPath + ":mask")
511                mask = 0
512        else
513                Wave mask=$"root:Packages:NIST:MSK:data"
514        Endif
515       
516        rcentr = 150            //pixels within rcentr of beam center are broken into 9 parts
517        // values for error if unable to estimate value
518        //large_num = 1e10
519        large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1
520        small_num = 1e-10
521       
522        // output wave are expected to exist (?) initialized to zero, what length?
523        // 300 points on VAX ---
524        Variable wavePts=500
525        Make/O/N=(wavePts) $(destPath + ":phival"),$(destPath + ":aveint")
526        Make/O/N=(wavePts) $(destPath + ":ncells"),$(destPath + ":sig"),$(destPath + ":sigave")
527        WAVE phival = $(destPath + ":phival")
528        WAVE aveint = $(destPath + ":aveint")
529        WAVE ncells = $(destPath + ":ncells")
530        WAVE sig = $(destPath + ":sig")
531        WAVE sigave = $(destPath + ":sigave")
532
533        phival = 0
534        aveint = 0
535        ncells = 0
536        sig = 0
537        sigave = 0
538
539        dtdis2 = dtdist^2
540        nq = 1
541        xoffst=0
542        //distance of beam center from detector center
543        xbm = FX(x0,sx3,xcenter,sx)
544        ybm = FY(y0,sy3,ycenter,sy)
545               
546        //BEGIN AVERAGE **********
547        Variable xi,xd,x,y,yd,yj,nd,fd,nd2,iphi,ntotal,var
548        Variable jj,data_pixel,xx,yy,ll,kk,rij,phiij,avesq,aveisq
549
550        // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
551        // loop index corresponds to FORTRAN (old code)
552        // and the IGOR array indices must be adjusted (-1) to the correct address
553        ntotal = 0
554        ii=1
555        do
556                xi = ii
557                xd = FX(xi,sx3,xcenter,sx)
558                x = xoffst + xd -xbm            //x and y are in mm
559               
560                jj = 1
561                do
562                        data_pixel = data[ii-1][jj-1]           //assign to local variable
563                        yj = jj
564                        yd = FY(yj,sy3,ycenter,sy)
565                        y = yd - ybm
566                        if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...)
567                                nd = 1
568                                fd = 1
569                                if( (abs(x) > rcentr) || (abs(y) > rcentr))     //break pixel into 9 equal parts
570                                        nd = 3
571                                        fd = 2
572                                Endif
573                                nd2 = nd^2
574                                ll = 1          //"el-el" loop index
575                                do
576                                        xx = x + (ll - fd)*sx/3
577                                        kk = 1
578                                        do
579                                                yy = y + (kk - fd)*sy/3
580                                                //test to see if center of pixel (i,j) lies in annulus
581                                                rij = sqrt(x*x + y*y)/dr + 1.001
582                                                //check whether pixel lies within width band
583                                                if((rij > rlo) && (rij < rhi))
584                                                        //in the annulus, do something
585                                                        if (yy >= 0)
586                                                                //phiij is in degrees
587                                                                phiij = atan2(yy,xx)*180/Pi             //0 to 180 deg
588                                                        else
589                                                                phiij = 360 + atan2(yy,xx)*180/Pi               //180 to 360 deg
590                                                        Endif
591                                                        if (phiij > (360-0.5*dphi))
592                                                                phiij -= 360
593                                                        Endif
594                                                        iphi = trunc(phiij/dphi + 1.501)
595                                                        aveint[iphi-1] += 9*data_pixel/nd2
596                                                        sig[iphi-1] += 9*data_pixel*data_pixel/nd2
597                                                        ncells[iphi-1] += 9/nd2
598                                                        ntotal += 9/nd2
599                                                Endif           //check if in annulus
600                                                kk+=1
601                                        while(kk<=nd)
602                                        ll += 1
603                                while(ll<=nd)
604                        Endif           //masked pixel check
605                        jj += 1
606                while (jj<=pixelsY)
607                ii += 1
608        while(ii<=pixelsX)              //end of the averaging
609               
610        //compute phi-values and errors
611       
612        ntotal /=9
613       
614        kk = 1
615        do
616                phival[kk-1] = dphi*(kk-1)
617                if(ncells[kk-1] != 0)
618                        aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
619                        avesq = aveint[kk-1]*aveint[kk-1]
620                        aveisq = sig[kk-1]/ncells[kk-1]
621                        var = aveisq - avesq
622                        if (var <=0 )
623                                sig[kk-1] = 0
624                                sigave[kk-1] = 0
625                                ncells[kk-1] /=9
626                        else
627                                if(ncells[kk-1] > 9)
628                                        sigave[kk-1] = sqrt(9*var/(ncells[kk-1]-9))
629                                        sig[kk-1] = sqrt( abs(aveint[kk-1])/(ncells[kk-1]/9) )
630                                        ncells[kk-1] /=9
631                                else
632                                        sig[kk-1] = 0
633                                        sigave[kk-1] = 0
634                                        ncells[kk-1] /=9
635                                Endif
636                        Endif
637                Endif
638                kk+=1
639        while(kk<=nphi)
640       
641        // data waves were defined as 200 points (=wavePts), but now have less than that (nphi) points
642        // use DeletePoints to remove junk from end of waves
643        Variable startElement,numElements
644        startElement = nphi
645        numElements = wavePts - startElement
646        DeletePoints startElement,numElements, phival,aveint,ncells,sig,sigave
647       
648        //////////////end of VAX Phibin.for
649               
650        //angle dependent transmission correction is not done in phiave
651        Ann_1D_Graph(aveint,phival,sigave)
652       
653        //get rid of the default mask, if one was created (it is in the current folder)
654        //don't just kill "mask" since it might be pointing to the one in the MSK folder
655        Killwaves/z $(destPath+":mask")
656               
657        //return to root folder (redundant)
658        SetDataFolder root:
659       
660        Return 0
661End
Note: See TracBrowser for help on using the repository browser.