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 function 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. Rectangular 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) |
---|
30 | Function 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:Packages:NIST: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 | 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 300 points (=defWavePts), 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 = defWavePts - 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 | // |
---|
282 | // The transmission correction is now done at the ADD step, in DetCorr() |
---|
283 | // |
---|
284 | // ////this section is the trans_correct() VAX routine |
---|
285 | // if(trans<0.1) |
---|
286 | // Print "***transmission is less than 0.1*** and is a significant correction" |
---|
287 | // endif |
---|
288 | // if(trans==0) |
---|
289 | // Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation" |
---|
290 | // trans = 1 |
---|
291 | // endif |
---|
292 | // //optical thickness |
---|
293 | // uval = -ln(trans) //use natural logarithm |
---|
294 | // //apply correction to aveint[] |
---|
295 | // //index from zero here, since only working with IGOR waves |
---|
296 | // ii=0 |
---|
297 | // do |
---|
298 | // theta = 2*asin(lambda*qval[ii]/(4*pi)) |
---|
299 | // cos_th = cos(theta) |
---|
300 | // arg = (1-cos_th)/cos_th |
---|
301 | // if((uval<0.01) || (cos_th>0.99)) //OR |
---|
302 | // //small arg, approx correction |
---|
303 | // aveint[ii] /= 1-0.5*uval*arg |
---|
304 | // else |
---|
305 | // //large arg, exact correction |
---|
306 | // aveint[ii] /= (1-exp(-uval*arg))/(uval*arg) |
---|
307 | // endif |
---|
308 | // ii+=1 |
---|
309 | // while(ii<nq) |
---|
310 | // //end of transmission/pathlength correction |
---|
311 | |
---|
312 | // *************************************************************** |
---|
313 | // |
---|
314 | // Do the extra 3 columns of resolution calculations starting here. |
---|
315 | // |
---|
316 | // *************************************************************** |
---|
317 | |
---|
318 | Variable L2 = reals[18] |
---|
319 | Variable BS = reals[21] |
---|
320 | Variable S1 = reals[23] |
---|
321 | Variable S2 = reals[24] |
---|
322 | Variable L1 = reals[25] |
---|
323 | lambda = reals[26] |
---|
324 | Variable lambdaWidth = reals[27] |
---|
325 | String detStr=textRead[9] |
---|
326 | |
---|
327 | Variable usingLenses = reals[28] //new 2007 |
---|
328 | |
---|
329 | //Two parameters DDET and APOFF are instrument dependent. Determine |
---|
330 | //these from the instrument name in the header. |
---|
331 | //From conversation with JB on 01.06.99 these are the current |
---|
332 | //good values |
---|
333 | |
---|
334 | Variable DDet |
---|
335 | NVAR apOff = root:myGlobals:apOff //in cm |
---|
336 | |
---|
337 | // DDet = DetectorPixelResolution(fileStr,detStr) //needs detector type and beamline |
---|
338 | //note that reading the detector pixel size from the header ASSUMES SQUARE PIXELS! - Jan2008 |
---|
339 | DDet = reals[10]/10 // header value (X) is in mm, want cm here |
---|
340 | |
---|
341 | |
---|
342 | //Width of annulus used for the average is gotten from the |
---|
343 | //input dialog before. This also must be passed to the resolution |
---|
344 | //calculator. Currently the default is dr=1 so just keeping that. |
---|
345 | |
---|
346 | //Go from 0 to nq doing the calc for all three values at |
---|
347 | //every Q value |
---|
348 | |
---|
349 | ii=0 |
---|
350 | |
---|
351 | Variable ret1,ret2,ret3 |
---|
352 | do |
---|
353 | getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,usingLenses,ret1,ret2,ret3) |
---|
354 | sigmaq[ii] = ret1 |
---|
355 | qbar[ii] = ret2 |
---|
356 | fsubs[ii] = ret3 |
---|
357 | ii+=1 |
---|
358 | while(ii<nq) |
---|
359 | DeletePoints startElement,numElements, sigmaq, qbar, fsubs |
---|
360 | |
---|
361 | // End of resolution calculations |
---|
362 | // *************************************************************** |
---|
363 | |
---|
364 | //Plot the data in the Plot_1d window |
---|
365 | Avg_1D_Graph(aveint,qval,sigave) |
---|
366 | |
---|
367 | //get rid of the default mask, if one was created (it is in the current folder) |
---|
368 | //don't just kill "mask" since it might be pointing to the one in the MSK folder |
---|
369 | Killwaves/Z $(destPath+":mask") |
---|
370 | |
---|
371 | //return to root folder (redundant) |
---|
372 | SetDataFolder root: |
---|
373 | |
---|
374 | Return 0 |
---|
375 | End |
---|
376 | |
---|
377 | //returns nq, new number of q-values |
---|
378 | //arrays aveint,dsq,ncells are also changed by this function |
---|
379 | // |
---|
380 | Function IncrementPixel(dataPixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2) |
---|
381 | Variable dataPixel,ddr,dxx,dyy |
---|
382 | Wave aveint,dsq,ncells |
---|
383 | Variable nq,nd2 |
---|
384 | |
---|
385 | Variable ir |
---|
386 | |
---|
387 | ir = trunc(sqrt(dxx*dxx+dyy*dyy)/ddr)+1 |
---|
388 | if (ir>nq) |
---|
389 | nq = ir //resets maximum number of q-values |
---|
390 | endif |
---|
391 | aveint[ir-1] += dataPixel/nd2 //ir-1 must be used, since ir is physical |
---|
392 | dsq[ir-1] += dataPixel*dataPixel/nd2 |
---|
393 | ncells[ir-1] += 1/nd2 |
---|
394 | |
---|
395 | Return nq |
---|
396 | End |
---|
397 | |
---|
398 | //function determines azimuthal angle dphi that a vector connecting |
---|
399 | //center of detector to pixel makes with respect to vector |
---|
400 | //at chosen azimuthal angle phi -> [cos(phi),sin(phi)] = [phi_x,phi_y] |
---|
401 | //dphi is always positive, varying from 0 to Pi |
---|
402 | // |
---|
403 | Function dphi_pixel(dxx,dyy,phi_x,phi_y) |
---|
404 | Variable dxx,dyy,phi_x,phi_y |
---|
405 | |
---|
406 | Variable val,rr,dot_prod |
---|
407 | |
---|
408 | rr = sqrt(dxx^2 + dyy^2) |
---|
409 | dot_prod = (dxx*phi_x + dyy*phi_y)/rr |
---|
410 | //? correct for roundoff error? - is this necessary in IGOR, w/ double precision? |
---|
411 | if(dot_prod > 1) |
---|
412 | dot_prod =1 |
---|
413 | Endif |
---|
414 | if(dot_prod < -1) |
---|
415 | dot_prod = -1 |
---|
416 | Endif |
---|
417 | |
---|
418 | val = acos(dot_prod) |
---|
419 | |
---|
420 | return val |
---|
421 | |
---|
422 | End |
---|
423 | |
---|
424 | //calculates the x distance from the center of the detector, w/nonlinear corrections |
---|
425 | // |
---|
426 | Function FX(xx,sx3,xcenter,sx) |
---|
427 | Variable xx,sx3,xcenter,sx |
---|
428 | |
---|
429 | Variable retval |
---|
430 | |
---|
431 | retval = sx3*tan((xx-xcenter)*sx/sx3) |
---|
432 | Return retval |
---|
433 | End |
---|
434 | |
---|
435 | //calculates the y distance from the center of the detector, w/nonlinear corrections |
---|
436 | // |
---|
437 | Function FY(yy,sy3,ycenter,sy) |
---|
438 | Variable yy,sy3,ycenter,sy |
---|
439 | |
---|
440 | Variable retval |
---|
441 | |
---|
442 | retval = sy3*tan((yy-ycenter)*sy/sy3) |
---|
443 | Return retval |
---|
444 | End |
---|
445 | |
---|
446 | //old function not called anymore, now "ave" button calls routine from AvgGraphics.ipf |
---|
447 | //to get input from panel rather than large prompt for missing parameters |
---|
448 | Function Ave_button(button0) : ButtonControl |
---|
449 | String button0 |
---|
450 | |
---|
451 | // the button on the graph will average the currently displayed data |
---|
452 | SVAR type=root:myGlobals:gDataDisplayType |
---|
453 | |
---|
454 | //Check for logscale data in "type" folder |
---|
455 | SetDataFolder "root:Packages:NIST:"+type //use the full path, so it will always work |
---|
456 | String dest = "root:Packages:NIST:" + type |
---|
457 | |
---|
458 | NVAR isLogScale = $(dest + ":gIsLogScale") |
---|
459 | if(isLogScale == 1) |
---|
460 | //data is logscale, convert it back and reset the global |
---|
461 | Duplicate/O $(dest + ":linear_data") $(dest + ":data") |
---|
462 | // WAVE vlegend=$(dest + ":vlegend") |
---|
463 | // Make the color table linear scale |
---|
464 | // vlegend = y |
---|
465 | Variable/G $(dest + ":gIsLogScale") = 0 //copy to keep with the current data folder |
---|
466 | SetDataFolder root: |
---|
467 | //rename the button to reflect "isLin" - the displayed name must have been isLog |
---|
468 | Button bisLog,title="isLin",rename=bisLin |
---|
469 | Endif |
---|
470 | |
---|
471 | //set data folder back to root |
---|
472 | SetDataFolder root: |
---|
473 | |
---|
474 | //do the average - ask the user for what type of average |
---|
475 | //ask the user for averaging paramters |
---|
476 | Execute "GetAvgInfo()" |
---|
477 | |
---|
478 | //dispatch to correct averaging routine |
---|
479 | //if you want to save the files, see Panel_DoAverageButtonProc(ctrlName) function |
---|
480 | //for making a fake protocol (needed to write out data) |
---|
481 | SVAR tempStr = root:myGlobals:Protocols:gAvgInfoStr |
---|
482 | String choice = StringByKey("AVTYPE",tempStr,"=",";") |
---|
483 | if(cmpstr("Rectangular",choice)==0) |
---|
484 | //dispatch to rectangular average |
---|
485 | RectangularAverageTo1D(type) |
---|
486 | else |
---|
487 | if(cmpstr("Annular",choice)==0) |
---|
488 | AnnularAverageTo1D(type) |
---|
489 | else |
---|
490 | //circular or sector |
---|
491 | CircularAverageTo1D(type) |
---|
492 | Endif |
---|
493 | Endif |
---|
494 | |
---|
495 | Return 0 |
---|
496 | End |
---|
497 | |
---|
498 | |
---|
499 | |
---|
500 | // -- seems to work, now I need to give it a name, add it to the list, and |
---|
501 | // make sure I've thought of all of the cases - then the average can be passed as case "Sector_PlusMinus" |
---|
502 | // and run through the normal average and writing routines. |
---|
503 | // |
---|
504 | // |
---|
505 | // -- depending on what value PHI has - it's [-90,90] "left" and "right" may not be what |
---|
506 | // you expect. so sorting the concatenated values may be necessary (always) |
---|
507 | // |
---|
508 | // -- need documentation of the definition of PHI, left, and right so that it can make better sense |
---|
509 | // which quadrants of the detector become "negative" depending on the choice of phi. may need a |
---|
510 | // switch after a little thinking. |
---|
511 | // |
---|
512 | // may want a variation of this where both sides are done, in separate files. but I think that's already |
---|
513 | // called a "sector" average. save them. load them. plot them. |
---|
514 | // |
---|
515 | // |
---|
516 | Function Sector_PlusMinus1D(type) |
---|
517 | String type |
---|
518 | |
---|
519 | // do the left side (-) |
---|
520 | // then hold that data in tmp_ waves |
---|
521 | // then do the right (+) |
---|
522 | // then concatenate the data |
---|
523 | |
---|
524 | // the button on the pink panel copies the two strings so they're the same |
---|
525 | SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr //this is the list that has it all |
---|
526 | |
---|
527 | String oldStr = "" |
---|
528 | String CurPath="root:myGlobals:Plot_1D:" |
---|
529 | String destPath = "root:Packages:NIST:"+type+":" |
---|
530 | |
---|
531 | oldStr = StringByKey("SIDE",keyListStr,"=",";") |
---|
532 | |
---|
533 | // do the left first, and call it negative |
---|
534 | keyListStr = ReplaceStringByKey("SIDE",keyListStr,"left","=",";") |
---|
535 | |
---|
536 | CircularAverageTo1D(type) |
---|
537 | |
---|
538 | WAVE qval = $(destPath + "qval") |
---|
539 | WAVE aveint = $(destPath + "aveint") |
---|
540 | WAVE sigave = $(destPath + "sigave") |
---|
541 | WAVE qbar = $(destPath + "QBar") |
---|
542 | WAVE sigmaq = $(destPath + "SigmaQ") |
---|
543 | WAVE fsubs = $(destPath + "fSubS") |
---|
544 | |
---|
545 | // copy the average, set the q's negative |
---|
546 | qval *= -1 |
---|
547 | Duplicate/O qval $(destPath+"tmp_q") |
---|
548 | Duplicate/O aveint $(destPath+"tmp_i") |
---|
549 | Duplicate/O sigave $(destPath+"tmp_s") |
---|
550 | Duplicate/O qbar $(destPath+"tmp_qb") |
---|
551 | Duplicate/O sigmaq $(destPath+"tmp_sq") |
---|
552 | Duplicate/O fsubs $(destPath+"tmp_fs") |
---|
553 | |
---|
554 | |
---|
555 | // do the right side |
---|
556 | keyListStr = ReplaceStringByKey("SIDE",keyListStr,"right","=",";") |
---|
557 | |
---|
558 | CircularAverageTo1D(type) |
---|
559 | |
---|
560 | // concatenate |
---|
561 | WAVE tmp_q = $(destPath + "tmp_q") |
---|
562 | WAVE tmp_i = $(destPath + "tmp_i") |
---|
563 | WAVE tmp_s = $(destPath + "tmp_s") |
---|
564 | WAVE tmp_qb = $(destPath + "tmp_qb") |
---|
565 | WAVE tmp_sq = $(destPath + "tmp_sq") |
---|
566 | WAVE tmp_fs = $(destPath + "tmp_fs") |
---|
567 | |
---|
568 | SetDataFolder destPath //to get the concatenation in the right folder |
---|
569 | Concatenate/NP/O {tmp_q,qval},tmp_cat |
---|
570 | Duplicate/O tmp_cat qval |
---|
571 | Concatenate/NP/O {tmp_i,aveint},tmp_cat |
---|
572 | Duplicate/O tmp_cat aveint |
---|
573 | Concatenate/NP/O {tmp_s,sigave},tmp_cat |
---|
574 | Duplicate/O tmp_cat sigave |
---|
575 | Concatenate/NP/O {tmp_qb,qbar},tmp_cat |
---|
576 | Duplicate/O tmp_cat qbar |
---|
577 | Concatenate/NP/O {tmp_sq,sigmaq},tmp_cat |
---|
578 | Duplicate/O tmp_cat sigmaq |
---|
579 | Concatenate/NP/O {tmp_fs,fsubs},tmp_cat |
---|
580 | Duplicate/O tmp_cat fsubs |
---|
581 | |
---|
582 | // then sort |
---|
583 | Sort qval, qval,aveint,sigave,qbar,sigmaq,fsubs |
---|
584 | |
---|
585 | // move these to the Plot_1D folder for plotting |
---|
586 | Duplicate/O qval $(curPath+"xAxisWave") |
---|
587 | Duplicate/O aveint $(curPath+"yAxisWave") |
---|
588 | Duplicate/O sigave $(curPath+"yErrWave") |
---|
589 | |
---|
590 | keyListStr = ReplaceStringByKey("SIDE",keyListStr,oldStr,"=",";") |
---|
591 | |
---|
592 | DoUpdate/W=Plot_1d |
---|
593 | |
---|
594 | // clean up |
---|
595 | KillWaves/Z tmp_q,tmp_i,tmp_s,tmp_qb,tmp_sq,tmp_fs,tmp_cat |
---|
596 | |
---|
597 | SetDataFolder root: |
---|
598 | |
---|
599 | return(0) |
---|
600 | End |
---|