source: sans/SANSReduction/branches/kline_29MAR07/Put in User Procedures/SANS_Reduction_v5.00/RectAnnulAvg.ipf @ 72

Last change on this file since 72 was 47, checked in by srkline, 16 years ago

bug fixes and performance tweaks

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