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 | // |
---|
45 | Function 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 |
---|
378 | End |
---|
379 | |
---|
380 | //returns nq, new number of q-values |
---|
381 | //arrays aveint,dsq,ncells are also changed by this function |
---|
382 | // |
---|
383 | Function 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 |
---|
399 | End |
---|
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 | // |
---|
409 | Function 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 | |
---|
443 | End |
---|
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 | // |
---|
455 | Function 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 |
---|
661 | End |
---|