source: sans/SANSReduction/trunk/Put in User Procedures/SANS_Reduction_v5.00/CircSectAve.ipf @ 40

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

Initial checkin of V5.00 SANSReduction files

File size: 16.0 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=4.0
4
5//*********************************************
6//              For AVERAGE and for DRAWING
7//                      DRAWING routines only use a subset of the total list, since saving, naming, etc. don't apply
8//              (10) possible keywords, some numerical, some string values
9//              AVTYPE=string           string from set {Circular,Annular,Rectangular,Sector,2D_ASCII,QxQy_ASCII,PNG_Graphic}
10//              PHI=value                       azimuthal angle (-90,90)
11//              DPHI=value                      +/- angular range around phi for average
12//              WIDTH=value             total width of rectangular section, in pixels
13//              SIDE=string             string from set {left,right,both} **note NOT capitalized
14//              QCENTER=value           q-value (1/A) of center of annulus for annular average
15//              QDELTA=value            total width of annulus centered at QCENTER
16//              PLOT=string             string from set {Yes,No} = truth of generating plot of averaged data
17//              SAVE=string             string from set {Yes,No} = truth of saving averaged data to disk
18//              NAME=string             string from set {Auto,Manual} = Automatic name generation or Manual(dialog)
19//***********************************************
20
21
22// this fuction also does sector averaging
23//the parameters in the global keyword-string must have already been set somewhere
24//either directly, from the protocol, or from the Average_Panel
25//** the keyword-list has already been "pre-parsed" to send only Circular or Sector
26//averages to this routine. Rectangualr or annular averages get done elsewhere
27// TYPE parameter determines which data folder to work from
28//
29//annnulus (step) size is currently fixed at 1 (variable dr, below)
30Function CircularAverageTo1D(type)
31        String type
32       
33        SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr          //this is the list that has it all
34        Variable isCircular = 0
35       
36        if( cmpstr("Circular",StringByKey("AVTYPE",keyListStr,"=",";")) ==0)
37                isCircular = 1          //set a switch for later
38        Endif
39       
40        //type is the data type to do the averaging on, and will be set as the current folder
41        //get the current displayed data (so the correct folder is used)
42        String destPath = "root:"+type
43       
44        //
45        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist,dr,ddr
46        Variable lambda,trans
47        WAVE reals = $(destPath + ":RealsRead")
48        WAVE/T textread = $(destPath + ":TextRead")
49        String fileStr = textread[3]
50       
51        // center of detector, for non-linear corrections
52        NVAR pixelsX = root:myGlobals:gNPixelsX
53        NVAR pixelsY = root:myGlobals:gNPixelsY
54       
55        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
56        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
57       
58        // beam center, in pixels
59        x0 = reals[16]
60        y0 = reals[17]
61        //detector calibration constants
62        sx = reals[10]          //mm/pixel (x)
63        sx3 = reals[11]         //nonlinear coeff
64        sy = reals[13]          //mm/pixel (y)
65        sy3 = reals[14]         //nonlinear coeff
66       
67        dtsize = 10*reals[20]           //det size in mm
68        dtdist = 1000*reals[18]         // det distance in mm
69       
70        NVAR binWidth=root:myGlobals:gBinWidth
71       
72        dr = binWidth           // ***********annulus width set by user, default is one***********
73        ddr = dr*sx             //step size, in mm (this value should be passed to the resolution calculation, not dr 18NOV03)
74               
75        Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,dxbm,dybm,ii
76        Variable phi_rad,dphi_rad,phi_x,phi_y
77        Variable forward,mirror
78       
79        String side = StringByKey("SIDE",keyListStr,"=",";")
80//      Print "side = ",side
81       
82        if(!isCircular)         //must be sector avg (rectangular not sent to this function)
83                //convert from degrees to radians
84                phi_rad = (Pi/180)*NumberByKey("PHI",keyListStr,"=",";")
85                dphi_rad = (Pi/180)*NumberByKey("DPHI",keyListStr,"=",";")
86                //create cartesian values for unit vector in phi direction
87                phi_x = cos(phi_rad)
88                phi_y = sin(phi_rad)
89        Endif
90       
91        /// data wave is data in the current folder which was set at the top of the function
92        WAVE data=$(destPath + ":data")
93        //Check for the existence of the mask, if not, make one (local to this folder) that is null
94       
95        if(WaveExists($"root:MSK:data") == 0)
96                Print "There is no mask file loaded (WaveExists)- the data is not masked"
97                Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
98                Wave mask = $(destPath + ":mask")
99                mask = 0
100        else
101                Wave mask=$"root:MSK:data"
102        Endif
103       
104        //
105        //pixels within rcentr of beam center are broken into 9 parts (units of mm)
106        rcentr = 100
107        // values for error if unable to estimate value
108        //large_num = 1e10
109        large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1
110        small_num = 1e-10
111       
112        // output wave are expected to exist (?) initialized to zero, what length?
113        // 200 points on VAX --- use 300 here, or more if SAXS data is used with 1024x1024 detector (1000 pts seems good)
114        Variable defWavePts=500
115        Make/O/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint")
116        Make/O/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave")
117        Make/O/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar")
118
119        WAVE qval = $(destPath + ":qval")
120        WAVE aveint = $(destPath + ":aveint")
121        WAVE ncells = $(destPath + ":ncells")
122        WAVE dsq = $(destPath + ":dsq")
123        WAVE sigave = $(destPath + ":sigave")
124        WAVE qbar = $(destPath + ":QBar")
125        WAVE sigmaq = $(destPath + ":SigmaQ")
126        WAVE fsubs = $(destPath + ":fSubS")
127       
128        qval = 0
129        aveint = 0
130        ncells = 0
131        dsq = 0
132        sigave = 0
133        qbar = 0
134        sigmaq = 0
135        fsubs = 0
136
137        dtdis2 = dtdist^2
138        nq = 1
139        xoffst=0
140        //distance of beam center from detector center
141        dxbm = FX(x0,sx3,xcenter,sx)
142        dybm = FY(y0,sy3,ycenter,sy)
143               
144        //BEGIN AVERAGE **********
145        Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1
146        Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p
147       
148        // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
149        // loop index corresponds to FORTRAN (old code)
150        // and the IGOR array indices must be adjusted (-1) to the correct address
151        ii=1
152        do
153                xi = ii
154                dxi = FX(xi,sx3,xcenter,sx)
155                dx = dxi-dxbm           //dx and dy are in mm
156               
157                jj = 1
158                do
159                        data_pixel = data[ii-1][jj-1]           //assign to local variable
160                        yj = jj
161                        dyj = FY(yj,sy3,ycenter,sy)
162                        dy = dyj - dybm
163                        if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...)
164                                dr2 = (dx^2 + dy^2)^(0.5)               //distance from beam center NOTE dr2 used here - dr used above
165                                if(dr2>rcentr)          //keep pixel whole
166                                        nd = 1
167                                        fd = 1
168                                else                            //break pixel into 9 equal parts
169                                        nd = 3
170                                        fd = 2
171                                endif
172                                nd2 = nd^2
173                                ll = 1          //"el-el" loop index
174                                do
175                                        dxx = dx + (ll - fd)*sx/3
176                                        kk = 1
177                                        do
178                                                dyy = dy + (kk - fd)*sy/3
179                                                if(isCircular)
180                                                        //circular average, use all pixels
181                                                        //(increment)
182                                                        nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
183                                                else
184                                                        //a sector average - determine azimuthal angle
185                                                        dphi_p = dphi_pixel(dxx,dyy,phi_x,phi_y)
186                                                        if(dphi_p < dphi_rad)
187                                                                forward = 1                     //within forward sector
188                                                        else
189                                                                forward = 0
190                                                        Endif
191                                                        if((Pi - dphi_p) < dphi_rad)
192                                                                mirror = 1              //within mirror sector
193                                                        else
194                                                                mirror = 0
195                                                        Endif
196                                                        //check if pixel lies within allowed sector(s)
197                                                        if(cmpstr(side,"both")==0)              //both sectors
198                                                                if ( mirror || forward)
199                                                                        //increment
200                                                                        nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
201                                                                Endif
202                                                        else
203                                                                if(cmpstr(side,"right")==0)             //forward sector only
204                                                                        if(forward)
205                                                                                //increment
206                                                                                nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
207                                                                        Endif
208                                                                else                    //mirror sector only
209                                                                        if(mirror)
210                                                                                //increment
211                                                                                nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
212                                                                        Endif
213                                                                Endif
214                                                        Endif           //allowable sectors
215                                                Endif   //circular or sector check
216                                                kk+=1
217                                        while(kk<=nd)
218                                        ll += 1
219                                while(ll<=nd)
220                        Endif           //masked pixel check
221                        jj += 1
222                while (jj<=pixelsY)
223                ii += 1
224        while(ii<=pixelsX)              //end of the averaging
225               
226        //compute q-values and errors
227        Variable ntotal,rr,theta,avesq,aveisq,var
228       
229        lambda = reals[26]
230        ntotal = 0
231        kk = 1
232        do
233                rr = (2*kk-1)*ddr/2
234                theta = 0.5*atan(rr/dtdist)
235                qval[kk-1] = (4*Pi/lambda)*sin(theta)
236                if(ncells[kk-1] == 0)
237                        //no pixels in annuli, data unknown
238                        aveint[kk-1] = 0
239                        sigave[kk-1] = large_num
240                else
241                        if(ncells[kk-1] <= 1)
242                                //need more than one pixel to determine error
243                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
244                                sigave[kk-1] = large_num
245                        else
246                                //assume that the intensity in each pixel in annuli is normally
247                                // distributed about mean...
248                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
249                                avesq = aveint[kk-1]^2
250                                aveisq = dsq[kk-1]/ncells[kk-1]
251                                var = aveisq-avesq
252                                if(var<=0)
253                                        sigave[kk-1] = small_num
254                                else
255                                        sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1))
256                                endif
257                        endif
258                endif
259                ntotal += ncells[kk-1]
260                kk+=1
261        while(kk<=nq)
262       
263        //Print "NQ = ",nq
264        // data waves were defined as 300 points (=defWavePts), but now have less than that (nq) points
265        // use DeletePoints to remove junk from end of waves
266        //WaveStats would be a more foolproof implementation, to get the # points in the wave
267        Variable startElement,numElements
268        startElement = nq
269        numElements = defWavePts - startElement
270        DeletePoints startElement,numElements, qval,aveint,ncells,dsq,sigave
271       
272        //////////////end of VAX sector_ave()
273               
274        //angle dependent transmission correction
275        Variable uval,arg,cos_th
276        lambda = reals[26]
277        trans = reals[4]
278       
279        ////this section is the trans_correct() VAX routine
280        if(trans<0.1)
281                Print "***transmission is less than 0.1*** and is a significant correction"
282        endif
283        if(trans==0)
284                Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
285                trans = 1
286        endif
287        //optical thickness
288        uval = -ln(trans)               //use natural logarithm
289        //apply correction to aveint[]
290        //index from zero here, since only working with IGOR waves
291        ii=0
292        do
293                theta = 2*asin(lambda*qval[ii]/(4*pi))
294                cos_th = cos(theta)
295                arg = (1-cos_th)/cos_th
296                if((uval<0.01) || (cos_th>0.99))                //OR
297                        //small arg, approx correction
298                        aveint[ii] /= 1-0.5*uval*arg
299                else
300                        //large arg, exact correction
301                        aveint[ii] *= uval*arg/(1-exp(-uval*arg))
302                endif
303                ii+=1
304        while(ii<nq)
305        //end of transmission/pathlength correction
306
307// ***************************************************************
308//
309// Do the extra 3 columns of resolution calculations starting here.
310//
311// ***************************************************************
312
313        Variable L2 = reals[18]
314        Variable BS = reals[21]
315        Variable S1 = reals[23]
316        Variable S2 = reals[24]
317        Variable L1 = reals[25]
318        lambda = reals[26]
319        Variable lambdaWidth = reals[27]
320        String detStr=textRead[9]
321
322        //Two parameters DDET and APOFF are instrument dependent.  Determine
323        //these from the instrument name in the header.
324        //From conversation with JB on 01.06.99 these are the current
325        //good values
326
327        Variable DDet, apOff=0.0
328        DDet = DetectorPixelResolution(fileStr,detStr)          //needs detector type and beamline
329       
330        //print fileStr
331//      do
332//        if(strsearch(fileStr, "NG5", 0) != -1)
333//          //string was found, it's an NG5 file
334//          DDet = 1.0
335//          break
336//        Endif
337//        if(strsearch(fileStr, "NG3", 0) != -1)
338//          //string was found, it's an NG3 file
339//          DDet = 0.5
340//          break
341//        Endif
342//        if(strsearch(fileStr, "NG7", 0) != -1)
343//          //string was found, it's an NG7 file
344//          DDet = 0.5
345//          break
346//        Endif
347//      while(0)
348
349        //Width of annulus used for the average is gotten from the
350        //input dialog before.  This also must be passed to the resolution
351        //calculator. Currently the default is dr=1 so just keeping that.
352
353        //Go from 0 to nq doing the calc for all three values at
354        //every Q value
355
356        ii=0
357//      String res_string="root:myGlobals:Res_vals"
358//      Make/O/D/N=3 $res_string
359//      Wave res_wave=$res_string
360        Variable ret1,ret2,ret3
361        do
362                getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,ret1,ret2,ret3)
363                sigmaq[ii] = ret1       //res_wave[0]
364                qbar[ii] = ret2         //res_wave[1]
365                fsubs[ii] = ret3                //res_wave[2]
366                ii+=1
367        while(ii<nq)
368        DeletePoints startElement,numElements, sigmaq, qbar, fsubs
369
370// End of resolution calculations
371// ***************************************************************
372       
373        //Plot the data in the Plot_1d window
374        Avg_1D_Graph(aveint,qval,sigave)
375//      DoWindow/F Plot_1d
376//      If(V_flag == 1)         //the graph window already exists
377//              //kill the old graph and make a new one
378//              //easier than adjusting the old one
379//              DoWindow/K Plot_1d
380//      Endif
381//      //make a completely new graph
382//      Display /W=(412,51,727,302)/K=1 aveint vs qval
383//      ModifyGraph log=1
384//      ModifyGraph mode=3,marker=19,msize=1,rgb=(0,0,0)
385//      ErrorBars aveint Y,wave=(sigave,sigave)
386//      Label left "Counts";DelayUpdate
387//      Label bottom "q\\U"
388//      SVAR angst = root:myGlobals:gAngstStr
389//      Label bottom "q ("+angst+"\\S-1\\M)"
390//      DoWindow/C Plot_1d
391       
392        //get rid of the default mask, if one was created (it is in the current folder)
393        //don't just kill "mask" since it might be pointing to the one in the MSK folder
394        Killwaves/z $(destPath+":mask")
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(dataPixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
406        Variable dataPixel,ddr,dxx,dyy
407        Wave aveint,dsq,ncells
408        Variable nq,nd2
409       
410        Variable ir
411       
412        ir = trunc(sqrt(dxx*dxx+dyy*dyy)/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 azimuthal angle dphi that a vector connecting
424//center of detector to pixel makes with respect to vector
425//at chosen azimuthal angle phi -> [cos(phi),sin(phi0] = [phi_x,phi_y]
426//dphi is always positive, varying from 0 to Pi
427//
428Function dphi_pixel(dxx,dyy,phi_x,phi_y)
429        Variable dxx,dyy,phi_x,phi_y
430       
431        Variable val,rr,dot_prod
432       
433        rr = sqrt(dxx^2 + dyy^2)
434        dot_prod = (dxx*phi_x + dyy*phi_y)/rr
435        //? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
436        if(dot_prod > 1)
437                dot_prod =1
438        Endif
439        if(dot_prod < -1)
440                dot_prod = -1
441        Endif
442       
443        val = acos(dot_prod)
444       
445        return val
446
447End
448
449//calculates the x distance from the center of the detector, w/nonlinear corrections
450//
451Function FX(xx,sx3,xcenter,sx)         
452        Variable xx,sx3,xcenter,sx
453       
454        Variable retval
455       
456        retval = sx3*tan((xx-xcenter)*sx/sx3)
457        Return retval
458End
459
460//calculates the y distance from the center of the detector, w/nonlinear corrections
461//
462Function FY(yy,sy3,ycenter,sy)         
463        Variable yy,sy3,ycenter,sy
464       
465        Variable retval
466       
467        retval = sy3*tan((yy-ycenter)*sy/sy3)
468        Return retval
469End
470
471//old function not called anymore, now "ave" button calls routine from AvgGraphics.ipf
472//to get input from panel rather than large prompt for missing parameters
473Function Ave_button(button0) : ButtonControl
474        String button0
475
476        // the button on the graph will average the currently displayed data
477        SVAR type=root:myGlobals:gDataDisplayType
478       
479        //Check for logscale data in "type" folder
480        SetDataFolder "root:"+type              //use the full path, so it will always work
481        String dest = "root:" + type
482       
483        NVAR isLogScale = $(dest + ":gIsLogScale")
484        if(isLogScale == 1)
485                //data is logscale, convert it back and reset the global
486                Duplicate/O $(dest + ":linear_data") $(dest + ":data")
487//              WAVE vlegend=$(dest + ":vlegend")
488        //  Make the color table linear scale
489//              vlegend = y
490                Variable/G $(dest + ":gIsLogScale") = 0         //copy to keep with the current data folder
491                SetDataFolder root:
492                //rename the button to reflect "isLin" - the displayed name must have been isLog
493                Button bisLog,title="isLin",rename=bisLin
494               
495        Endif
496
497        //set data folder back to root
498        SetDataFolder root:
499       
500        //do the average - ask the user for what type of average
501        //ask the user for averaging paramters
502        Execute "GetAvgInfo()"
503       
504        //dispatch to correct averaging routine
505        //if you want to save the files, see Panel_DoAverageButtonProc(ctrlName) function
506        //for making a fake protocol (needed to write out data)
507        SVAR tempStr = root:myGlobals:Protocols:gAvgInfoStr
508        String choice = StringByKey("AVTYPE",tempStr,"=",";")
509        if(cmpstr("Rectangular",choice)==0)
510                //dispatch to rectangular average
511                RectangularAverageTo1D(type)
512        else
513                if(cmpstr("Annular",choice)==0)
514                        AnnularAverageTo1D(type)
515                else
516                        //circular or sector
517                        CircularAverageTo1D(type)
518                Endif
519        Endif
520       
521        Return 0
522End
523//above is unused function
Note: See TracBrowser for help on using the repository browser.