source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/CircSectAve.ipf @ 575

Last change on this file since 575 was 570, checked in by srkline, 13 years ago

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

Change (2):
As a side benefit of requiring 6.1, we can use the MultiThread? keyword to thread any model function we like. The speed benefit is only noticeable on functions that require at least one integration and at least 100 points (resolution smearing is NOT threaded, too many threadSafe issues, too little benefit). I have chosen to use the MultiThread? only on the XOP assignment. In the Igor code there are too many functions that are not explicitly declared threadsafe, making for a mess.

File size: 15.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
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:Packages:NIST:"+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:Packages:NIST: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:Packages:NIST:MSK:data"
102        Endif
103       
104        //
105        //pixels within rcentr of beam center are broken into 9 parts (units of mm)
106        rcentr = 100            //original
107//      rcentr = 0
108        // values for error if unable to estimate value
109        //large_num = 1e10
110        large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1
111        small_num = 1e-10
112       
113        // output wave are expected to exist (?) initialized to zero, what length?
114        // 200 points on VAX --- use 300 here, or more if SAXS data is used with 1024x1024 detector (1000 pts seems good)
115        Variable defWavePts=500
116        Make/O/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint")
117        Make/O/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave")
118        Make/O/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar")
119
120        WAVE qval = $(destPath + ":qval")
121        WAVE aveint = $(destPath + ":aveint")
122        WAVE ncells = $(destPath + ":ncells")
123        WAVE dsq = $(destPath + ":dsq")
124        WAVE sigave = $(destPath + ":sigave")
125        WAVE qbar = $(destPath + ":QBar")
126        WAVE sigmaq = $(destPath + ":SigmaQ")
127        WAVE fsubs = $(destPath + ":fSubS")
128       
129        qval = 0
130        aveint = 0
131        ncells = 0
132        dsq = 0
133        sigave = 0
134        qbar = 0
135        sigmaq = 0
136        fsubs = 0
137
138        dtdis2 = dtdist^2
139        nq = 1
140        xoffst=0
141        //distance of beam center from detector center
142        dxbm = FX(x0,sx3,xcenter,sx)
143        dybm = FY(y0,sy3,ycenter,sy)
144               
145        //BEGIN AVERAGE **********
146        Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1
147        Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p
148       
149        // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
150        // loop index corresponds to FORTRAN (old code)
151        // and the IGOR array indices must be adjusted (-1) to the correct address
152        ii=1
153        do
154                xi = ii
155                dxi = FX(xi,sx3,xcenter,sx)
156                dx = dxi-dxbm           //dx and dy are in mm
157               
158                jj = 1
159                do
160                        data_pixel = data[ii-1][jj-1]           //assign to local variable
161                        yj = jj
162                        dyj = FY(yj,sy3,ycenter,sy)
163                        dy = dyj - dybm
164                        if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...)
165                                dr2 = (dx^2 + dy^2)^(0.5)               //distance from beam center NOTE dr2 used here - dr used above
166                                if(dr2>rcentr)          //keep pixel whole
167                                        nd = 1
168                                        fd = 1
169                                else                            //break pixel into 9 equal parts
170                                        nd = 3
171                                        fd = 2
172                                endif
173                                nd2 = nd^2
174                                ll = 1          //"el-el" loop index
175                                do
176                                        dxx = dx + (ll - fd)*sx/3
177                                        kk = 1
178                                        do
179                                                dyy = dy + (kk - fd)*sy/3
180                                                if(isCircular)
181                                                        //circular average, use all pixels
182                                                        //(increment)
183                                                        nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
184                                                else
185                                                        //a sector average - determine azimuthal angle
186                                                        dphi_p = dphi_pixel(dxx,dyy,phi_x,phi_y)
187                                                        if(dphi_p < dphi_rad)
188                                                                forward = 1                     //within forward sector
189                                                        else
190                                                                forward = 0
191                                                        Endif
192                                                        if((Pi - dphi_p) < dphi_rad)
193                                                                mirror = 1              //within mirror sector
194                                                        else
195                                                                mirror = 0
196                                                        Endif
197                                                        //check if pixel lies within allowed sector(s)
198                                                        if(cmpstr(side,"both")==0)              //both sectors
199                                                                if ( mirror || forward)
200                                                                        //increment
201                                                                        nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
202                                                                Endif
203                                                        else
204                                                                if(cmpstr(side,"right")==0)             //forward sector only
205                                                                        if(forward)
206                                                                                //increment
207                                                                                nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
208                                                                        Endif
209                                                                else                    //mirror sector only
210                                                                        if(mirror)
211                                                                                //increment
212                                                                                nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
213                                                                        Endif
214                                                                Endif
215                                                        Endif           //allowable sectors
216                                                Endif   //circular or sector check
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 300 points (=defWavePts), 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 = defWavePts - 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//
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)               //use natural logarithm
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       
341        //Width of annulus used for the average is gotten from the
342        //input dialog before.  This also must be passed to the resolution
343        //calculator. Currently the default is dr=1 so just keeping that.
344
345        //Go from 0 to nq doing the calc for all three values at
346        //every Q value
347
348        ii=0
349
350        Variable ret1,ret2,ret3
351        do
352                getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,usingLenses,ret1,ret2,ret3)
353                sigmaq[ii] = ret1       
354                qbar[ii] = ret2
355                fsubs[ii] = ret3       
356                ii+=1
357        while(ii<nq)
358        DeletePoints startElement,numElements, sigmaq, qbar, fsubs
359
360// End of resolution calculations
361// ***************************************************************
362       
363        //Plot the data in the Plot_1d window
364        Avg_1D_Graph(aveint,qval,sigave)
365
366        //get rid of the default mask, if one was created (it is in the current folder)
367        //don't just kill "mask" since it might be pointing to the one in the MSK folder
368        Killwaves/Z $(destPath+":mask")
369       
370        //return to root folder (redundant)
371        SetDataFolder root:
372       
373        Return 0
374End
375
376//returns nq, new number of q-values
377//arrays aveint,dsq,ncells are also changed by this function
378//
379Function IncrementPixel(dataPixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
380        Variable dataPixel,ddr,dxx,dyy
381        Wave aveint,dsq,ncells
382        Variable nq,nd2
383       
384        Variable ir
385       
386        ir = trunc(sqrt(dxx*dxx+dyy*dyy)/ddr)+1
387        if (ir>nq)
388                nq = ir         //resets maximum number of q-values
389        endif
390        aveint[ir-1] += dataPixel/nd2           //ir-1 must be used, since ir is physical
391        dsq[ir-1] += dataPixel*dataPixel/nd2
392        ncells[ir-1] += 1/nd2
393       
394        Return nq
395End
396
397//function determines azimuthal angle dphi that a vector connecting
398//center of detector to pixel makes with respect to vector
399//at chosen azimuthal angle phi -> [cos(phi),sin(phi)] = [phi_x,phi_y]
400//dphi is always positive, varying from 0 to Pi
401//
402Function dphi_pixel(dxx,dyy,phi_x,phi_y)
403        Variable dxx,dyy,phi_x,phi_y
404       
405        Variable val,rr,dot_prod
406       
407        rr = sqrt(dxx^2 + dyy^2)
408        dot_prod = (dxx*phi_x + dyy*phi_y)/rr
409        //? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
410        if(dot_prod > 1)
411                dot_prod =1
412        Endif
413        if(dot_prod < -1)
414                dot_prod = -1
415        Endif
416       
417        val = acos(dot_prod)
418       
419        return val
420
421End
422
423//calculates the x distance from the center of the detector, w/nonlinear corrections
424//
425Function FX(xx,sx3,xcenter,sx)         
426        Variable xx,sx3,xcenter,sx
427       
428        Variable retval
429       
430        retval = sx3*tan((xx-xcenter)*sx/sx3)
431        Return retval
432End
433
434//calculates the y distance from the center of the detector, w/nonlinear corrections
435//
436Function FY(yy,sy3,ycenter,sy)         
437        Variable yy,sy3,ycenter,sy
438       
439        Variable retval
440       
441        retval = sy3*tan((yy-ycenter)*sy/sy3)
442        Return retval
443End
444
445//old function not called anymore, now "ave" button calls routine from AvgGraphics.ipf
446//to get input from panel rather than large prompt for missing parameters
447Function Ave_button(button0) : ButtonControl
448        String button0
449
450        // the button on the graph will average the currently displayed data
451        SVAR type=root:myGlobals:gDataDisplayType
452       
453        //Check for logscale data in "type" folder
454        SetDataFolder "root:Packages:NIST:"+type                //use the full path, so it will always work
455        String dest = "root:Packages:NIST:" + type
456       
457        NVAR isLogScale = $(dest + ":gIsLogScale")
458        if(isLogScale == 1)
459                //data is logscale, convert it back and reset the global
460                Duplicate/O $(dest + ":linear_data") $(dest + ":data")
461//              WAVE vlegend=$(dest + ":vlegend")
462        //  Make the color table linear scale
463//              vlegend = y
464                Variable/G $(dest + ":gIsLogScale") = 0         //copy to keep with the current data folder
465                SetDataFolder root:
466                //rename the button to reflect "isLin" - the displayed name must have been isLog
467                Button bisLog,title="isLin",rename=bisLin
468        Endif
469
470        //set data folder back to root
471        SetDataFolder root:
472       
473        //do the average - ask the user for what type of average
474        //ask the user for averaging paramters
475        Execute "GetAvgInfo()"
476       
477        //dispatch to correct averaging routine
478        //if you want to save the files, see Panel_DoAverageButtonProc(ctrlName) function
479        //for making a fake protocol (needed to write out data)
480        SVAR tempStr = root:myGlobals:Protocols:gAvgInfoStr
481        String choice = StringByKey("AVTYPE",tempStr,"=",";")
482        if(cmpstr("Rectangular",choice)==0)
483                //dispatch to rectangular average
484                RectangularAverageTo1D(type)
485        else
486                if(cmpstr("Annular",choice)==0)
487                        AnnularAverageTo1D(type)
488                else
489                        //circular or sector
490                        CircularAverageTo1D(type)
491                Endif
492        Endif
493       
494        Return 0
495End
Note: See TracBrowser for help on using the repository browser.