source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Instrument_Resolution.ipf @ 1053

Last change on this file since 1053 was 1053, checked in by srkline, 5 years ago

a few quick bug fixes

File size: 20.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
4
5
6//
7// resolution calculations for VSANS, under a variety of collimation conditions
8//
9// Partially converted (July 2017)
10//
11//
12// -- still missing a lot of physical dimensions for the SANS (1D) case
13// let alone anything more complex
14//
15//
16// SANS-like (pinhole) conditions are largely copied from the SANS calcuations
17// and are the traditional extra three columns
18//
19// Other conditions, such as white beam, or narrow slit mode, will likely require some
20// format for the resolution information that is different than the three column format.
21// The USANS solution of a "flag" is clunky, and depends entirely on the analysis package to
22// know exactly what to do.
23//
24// the 2D SANS-like resolution calculation is also expected to be similar to SANS, but is
25// unverified at this point (July 2017). 2D errors are also unverified.
26// -- Most importantly for 2D VSANS data, there is no defined output format.
27//
28
29
30// TODO:
31// -- some of the input geometry is hidden in other locations:
32// Sample Aperture to Gate Valve (cm)  == /instrument/sample_aperture/distance
33// Sample [position] to Gate Valve (cm) = /instrument/sample_table/offset_distance
34//
35// -- the dimensions and the units for the beam stops are very odd, and what is written to the
36//   file is not what is noted in the GUI - so verify the units that I'm actually reading.
37//
38
39
40
41
42//**********************
43// Resolution calculation - used by the averaging routines
44// to calculate the resolution function at each q-value
45// - the return value is not used
46//
47// equivalent to John's routine on the VAX Q_SIGMA_AVE.FOR
48// Incorporates eqn. 3-15 from J. Appl. Cryst. (1995) v. 28 p105-114
49//
50// - 21 MAR 07 uses projected BS diameter on the detector
51// - APR 07 still need to add resolution with lenses. currently there is no flag in the
52//          raw data header to indicate the presence of lenses.
53//
54// - Aug 07 - added input to switch calculation based on lenses (==1 if in)
55//
56// - called by CircSectAvg.ipf and RectAnnulAvg.ipf
57//
58// passed values are read from RealsRead
59// except DDet and apOff, which are set from globals before passing
60//
61// DDet is the detector pixel resolution
62// apOff is the offset between the sample aperture and the sample position
63//
64//
65Function getResolution(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,SigmaQ,QBar,fSubS)
66        Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses
67        Variable &fSubS, &QBar, &SigmaQ         //these are the output quantities at the input Q value
68       
69        //lots of calculation variables
70        Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g
71        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
72
73        //Constants
74        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
75        Variable g = 981.0                              //gravity acceleration [cm/s^2]
76
77
78        S1 *= 0.5*0.1                   //convert to radius and [cm]
79        S2 *= 0.5*0.1
80
81        L1 *= 100.0                     // [cm]
82        L1 -= apOff                             //correct the distance
83
84        L2 *= 100.0
85        L2 += apOff
86        del_r *= 0.1                            //width of annulus, convert mm to [cm]
87       
88        BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm]
89        // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture
90        Variable LB
91        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
92        BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax
93       
94        //Start resolution calculation
95        a2 = S1*L2/L1 + S2*(L1+L2)/L1
96        q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2)
97        lp = 1.0/( 1.0/L1 + 1.0/L2)
98
99        v_lambda = lambdaWidth^2/6.0
100       
101//      if(usingLenses==1)                      //SRK 2007
102        if(usingLenses != 0)                    //SRK 2008 allows for the possibility of different numbers of lenses in header
103                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term
104        else
105                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form
106        endif
107       
108        v_d = (DDet/2.3548)^2 + del_r^2/12.0                    //the 2.3548 is a conversion from FWHM->Gauss, see http://mathworld.wolfram.com/GaussianFunction.html
109        vz = vz_1 / lambda
110        yg = 0.5*g*L2*(L1+L2)/vz^2
111        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007
112
113        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
114        delta = 0.5*(BS - r0)^2/v_d
115
116        if (r0 < BS)
117                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
118        else
119                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
120        endif
121
122        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
123        if (fSubS <= 0.0)
124                fSubS = 1.e-10
125        endif
126        fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
127        fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
128
129        rmd = fr*r0
130        v_r1 = v_b + fv*v_d +v_g
131
132        rm = rmd + 0.5*v_r1/rmd
133        v_r = v_r1 - 0.5*(v_r1/rmd)^2
134        if (v_r < 0.0)
135                v_r = 0.0
136        endif
137        QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2))
138        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda)
139
140
141// more readable method for calculating the variance in Q
142// EXCEPT - this is calculated for Qo, NOT qBar
143// (otherwise, they are nearly equivalent, except for close to the beam stop)
144//      Variable kap,a_val,a_val_l2,m_h
145//      g = 981.0                               //gravity acceleration [cm/s^2]
146//      m_h     = 252.8                 // m/h [=] s/cm^2
147//     
148//      kap = 2*pi/lambda
149//      a_val = L2*(L1+L2)*g/2*(m_h)^2
150//      a_val_L2 = a_val/L2*1e-16               //convert 1/cm^2 to 1/A^2
151//
152//      sigmaQ = 0
153//      sigmaQ = 3*(S1/L1)^2
154//     
155//      if(usingLenses != 0)
156//              sigmaQ += 2*(S2/lp)^2*(lambdaWidth)^2   //2nd term w/ lenses
157//      else   
158//              sigmaQ += 2*(S2/lp)^2                                           //2nd term w/ no lenses
159//      endif
160//     
161//      sigmaQ += (del_r/L2)^2
162//      sigmaQ += 2*(r0/L2)^2*(lambdaWidth)^2
163//      sigmaQ += 4*(a_val_l2)^2*lambda^4*(lambdaWidth)^2
164//     
165//      sigmaQ *= kap^2/12
166//      sigmaQ = sqrt(sigmaQ)
167
168
169        Return (0)
170End
171
172
173//
174//**********************
175// 2D resolution function calculation - ***NOT*** in terms of X and Y
176// but written in terms of Parallel and perpendicular to the Q vector at each point
177//
178// -- it is more naturally written this way since the 2D function is an ellipse with its major
179// axis pointing in the direction of Q_parallel. Hence there is no way to properly define the
180// elliptical gaussian in terms of sigmaX and sigmaY
181//
182// For a full description of the gravity effect on the resolution, see:
183//
184//      "The effect of gravity on the resolution of small-angle neutron diffraction peaks"
185//      D.F.R Mildner, J.G. Barker & S.R. Kline J. Appl. Cryst. (2011). 44, 1127-1129.
186//      [ doi:10.1107/S0021889811033322 ]
187//
188//              2/17/12 SRK
189//              NOTE: the first 2/3 of this code is the 1D code, copied here just to have the beam stop
190//                              calculation here, if I decide to implement it. The real calculation is all at the
191//                              bottom and is quite compact
192//
193//
194//
195//
196// - 21 MAR 07 uses projected BS diameter on the detector
197// - APR 07 still need to add resolution with lenses. currently there is no flag in the
198//          raw data header to indicate the presence of lenses.
199//
200// - Aug 07 - added input to switch calculation based on lenses (==1 if in)
201//
202// passed values are read from RealsRead
203// except DDet and apOff, which are set from globals before passing
204//
205// phi is the azimuthal angle, CCW from +x axis
206// r_dist is the real-space distance from ctr of detector to QxQy pixel location
207//
208// MAR 2011 - removed the del_r terms, they don't apply since no bining is done to the 2D data
209//
210Function get2DResolution(inQ,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,SigmaQX,SigmaQY,fSubS)
211        Variable inQ, phi,lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses,r_dist
212        Variable &SigmaQX,&SigmaQY,&fSubS               //these are the output quantities at the input Q value
213       
214        //lots of calculation variables
215        Variable a2, lp, v_lambda, v_b, v_d, vz, yg, v_g
216        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
217
218        //Constants
219        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
220        Variable g = 981.0                              //gravity acceleration [cm/s^2]
221        Variable m_h    = 252.8                 // m/h [=] s/cm^2
222
223
224        S1 *= 0.5*0.1                   //convert to radius and [cm]
225        S2 *= 0.5*0.1
226
227        L1 *= 100.0                     // [cm]
228        L1 -= apOff                             //correct the distance
229
230        L2 *= 100.0
231        L2 += apOff
232        del_r *= 0.1                            //width of annulus, convert mm to [cm]
233       
234        BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm]
235        // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture
236        Variable LB
237        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
238        BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax
239       
240        //Start resolution calculation
241        a2 = S1*L2/L1 + S2*(L1+L2)/L1
242        lp = 1.0/( 1.0/L1 + 1.0/L2)
243
244        v_lambda = lambdaWidth^2/6.0
245       
246//      if(usingLenses==1)                      //SRK 2007
247        if(usingLenses != 0)                    //SRK 2008 allows for the possibility of different numbers of lenses in header
248                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term
249        else
250                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form
251        endif
252       
253        v_d = (DDet/2.3548)^2 + del_r^2/12.0
254        vz = vz_1 / lambda
255        yg = 0.5*g*L2*(L1+L2)/vz^2
256        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007
257
258        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
259        delta = 0.5*(BS - r0)^2/v_d
260
261        if (r0 < BS)
262                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
263        else
264                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
265        endif
266
267        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
268        if (fSubS <= 0.0)
269                fSubS = 1.e-10
270        endif
271//      fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
272//      fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
273//
274//      rmd = fr*r0
275//      v_r1 = v_b + fv*v_d +v_g
276//
277//      rm = rmd + 0.5*v_r1/rmd
278//      v_r = v_r1 - 0.5*(v_r1/rmd)^2
279//      if (v_r < 0.0)
280//              v_r = 0.0
281//      endif
282
283        Variable kap,a_val,a_val_L2,proj_DDet
284       
285        kap = 2*pi/lambda
286        a_val = L2*(L1+L2)*g/2*(m_h)^2
287        a_val_L2 = a_val/L2*1e-16               //convert 1/cm^2 to 1/A^2
288
289
290        // the detector pixel is square, so correct for phi
291        proj_DDet = DDet*cos(phi) + DDet*sin(phi)
292
293
294///////// OLD - don't use ---
295//in terms of Q_parallel ("x") and Q_perp ("y") - this works, since parallel is in the direction of Q and I
296// can calculate that from the QxQy (I just need the projection)
297//// for test case with no gravity, set a_val = 0
298//// note that gravity has no wavelength dependence. the lambda^4 cancels out.
299////
300////    a_val = 0
301////    a_val_l2 = 0
302//
303//     
304//      // this is really sigma_Q_parallel
305//      SigmaQX = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2 + (sin(phi))^2*8*(a_val_L2)^2*lambda^4*lambdaWidth^2)
306//      SigmaQX += inQ*inQ*v_lambda
307//
308//      //this is really sigma_Q_perpendicular
309//      proj_DDet = DDet*sin(phi) + DDet*cos(phi)               //not necessary, since DDet is the same in both X and Y directions
310//
311//      SigmaQY = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2 + (cos(phi))^2*8*(a_val_L2)^2*lambda^4*lambdaWidth^2)
312//     
313//      SigmaQX = sqrt(SigmaQX)
314//      SigmaQy = sqrt(SigmaQY)
315//     
316
317/////////////////////////////////////////////////
318/////   
319//      ////// this is all new, inclusion of gravity effect into the parallel component
320//       perpendicular component is purely geometric, no gravity component
321//
322// the shadow factor is calculated as above -so keep the above calculations, even though
323// most of them are redundant.
324//
325       
326////    //
327        Variable yg_d,acc,sdd,ssd,lambda0,DL_L,sig_l
328        Variable var_qlx,var_qly,var_ql,qx,qy,sig_perp,sig_para, sig_para_new
329       
330        G = 981.  //!   ACCELERATION OF GRAVITY, CM/SEC^2
331        acc = vz_1              //      3.956E5 //!     CONVERT WAVELENGTH TO VELOCITY CM/SEC
332        SDD = L2                //1317
333        SSD = L1                //1627          //cm
334        lambda0 = lambda                //              15
335        DL_L = lambdaWidth              //0.236
336        SIG_L = DL_L/sqrt(6)
337        YG_d = -0.5*G*SDD*(SSD+SDD)*(LAMBDA0/acc)^2
338/////   Print "DISTANCE BEAM FALLS DUE TO GRAVITY (CM) = ",YG
339//              Print "Gravity q* = ",-2*pi/lambda0*2*yg_d/sdd
340       
341        sig_perp = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2)
342        sig_perp = sqrt(sig_perp)
343       
344// TODO -- not needed???       
345//      FindQxQy(inQ,phi,qx,qy)
346
347
348// missing a factor of 2 here, and the form is different than the paper, so re-write   
349//      VAR_QLY = SIG_L^2 * (QY+4*PI*YG_d/(2*SDD*LAMBDA0))^2
350//      VAR_QLX = (SIG_L*QX)^2
351//      VAR_QL = VAR_QLY + VAR_QLX  //! WAVELENGTH CONTRIBUTION TO VARIANCE
352//      sig_para = (sig_perp^2 + VAR_QL)^0.5
353       
354        // r_dist is passed in, [=]cm
355        // from the paper
356        a_val = 0.5*G*SDD*(SSD+SDD)*m_h^2 * 1e-16               //units now are cm /(A^2)
357       
358        var_QL = 1/6*(kap/SDD)^2*(DL_L)^2*(r_dist^2 - 4*r_dist*a_val*lambda0^2*sin(phi) + 4*a_val^2*lambda0^4)
359        sig_para_new = (sig_perp^2 + VAR_QL)^0.5
360       
361       
362///// return values PBR
363        SigmaQX = sig_para_new
364        SigmaQy = sig_perp
365       
366////   
367
368        Return (0)
369End
370
371
372
373
374
375
376////////Transmission
377//******************
378//lookup tables for attenuator transmissions
379//
380//
381// new calibration done June 2007, John Barker
382//
383Proc MakeNG3AttenTable()
384
385        NewDataFolder/O root:myGlobals:Attenuators
386        //do explicitly to avoid data folder problems, redundant, but it must work without fail
387        Variable num=10         //10 needed for tables after June 2007
388
389        Make/O/N=(num) root:myGlobals:Attenuators:ng3att0
390        Make/O/N=(num) root:myGlobals:Attenuators:ng3att1
391        Make/O/N=(num) root:myGlobals:Attenuators:ng3att2
392        Make/O/N=(num) root:myGlobals:Attenuators:ng3att3
393        Make/O/N=(num) root:myGlobals:Attenuators:ng3att4
394        Make/O/N=(num) root:myGlobals:Attenuators:ng3att5
395        Make/O/N=(num) root:myGlobals:Attenuators:ng3att6
396        Make/O/N=(num) root:myGlobals:Attenuators:ng3att7
397        Make/O/N=(num) root:myGlobals:Attenuators:ng3att8
398        Make/O/N=(num) root:myGlobals:Attenuators:ng3att9
399        Make/O/N=(num) root:myGlobals:Attenuators:ng3att10
400       
401        // and a wave for the errors at each attenuation factor
402        Make/O/N=(num) root:myGlobals:Attenuators:ng3att0_err
403        Make/O/N=(num) root:myGlobals:Attenuators:ng3att1_err
404        Make/O/N=(num) root:myGlobals:Attenuators:ng3att2_err
405        Make/O/N=(num) root:myGlobals:Attenuators:ng3att3_err
406        Make/O/N=(num) root:myGlobals:Attenuators:ng3att4_err
407        Make/O/N=(num) root:myGlobals:Attenuators:ng3att5_err
408        Make/O/N=(num) root:myGlobals:Attenuators:ng3att6_err
409        Make/O/N=(num) root:myGlobals:Attenuators:ng3att7_err
410        Make/O/N=(num) root:myGlobals:Attenuators:ng3att8_err
411        Make/O/N=(num) root:myGlobals:Attenuators:ng3att9_err
412        Make/O/N=(num) root:myGlobals:Attenuators:ng3att10_err
413       
414       
415        //each wave has 10 elements, the transmission of att# at the wavelengths
416        //lambda = 4,5,6,7,8,10,12,14,17,20 (4 A and 20 A are extrapolated values)
417        Make/O/N=(num) root:myGlobals:Attenuators:ng3lambda={4,5,6,7,8,10,12,14,17,20}
418       
419        // new calibration done June 2007, John Barker
420        root:myGlobals:Attenuators:ng3att0 = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }
421        root:myGlobals:Attenuators:ng3att1 = {0.444784,0.419,0.3935,0.3682,0.3492,0.3132,0.2936,0.2767,0.2477,0.22404}
422        root:myGlobals:Attenuators:ng3att2 = {0.207506,0.1848,0.1629,0.1447,0.1292,0.1056,0.09263,0.08171,0.06656,0.0546552}
423        root:myGlobals:Attenuators:ng3att3 = {0.092412,0.07746,0.06422,0.05379,0.04512,0.03321,0.02707,0.02237,0.01643,0.0121969}
424        root:myGlobals:Attenuators:ng3att4 = {0.0417722,0.03302,0.02567,0.02036,0.01604,0.01067,0.00812,0.006316,0.00419,0.00282411}
425        root:myGlobals:Attenuators:ng3att5 = {0.0187129,0.01397,0.01017,0.007591,0.005668,0.003377,0.002423,0.001771,0.001064,0.000651257}
426        root:myGlobals:Attenuators:ng3att6 = {0.00851048,0.005984,0.004104,0.002888,0.002029,0.001098,0.0007419,0.0005141,0.000272833,0.000150624}
427        root:myGlobals:Attenuators:ng3att7 = {0.00170757,0.001084,0.0006469,0.0004142,0.0002607,0.0001201,7.664e-05,4.06624e-05,1.77379e-05,7.30624e-06}
428        root:myGlobals:Attenuators:ng3att8 = {0.000320057,0.0001918,0.0001025,6.085e-05,3.681e-05,1.835e-05,6.74002e-06,3.25288e-06,1.15321e-06,3.98173e-07}
429        root:myGlobals:Attenuators:ng3att9 = {6.27682e-05,3.69e-05,1.908e-05,1.196e-05,8.738e-06,6.996e-06,6.2901e-07,2.60221e-07,7.49748e-08,2.08029e-08}
430        root:myGlobals:Attenuators:ng3att10 = {1.40323e-05,8.51e-06,5.161e-06,4.4e-06,4.273e-06,1.88799e-07,5.87021e-08,2.08169e-08,4.8744e-09,1.08687e-09}
431 
432  // percent errors as measured, May 2007 values
433  // zero error for zero attenuators, appropriate average values put in for unknown values
434        root:myGlobals:Attenuators:ng3att0_err = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
435        root:myGlobals:Attenuators:ng3att1_err = {0.15,0.142,0.154,0.183,0.221,0.328,0.136,0.13,0.163,0.15}
436        root:myGlobals:Attenuators:ng3att2_err = {0.25,0.257,0.285,0.223,0.271,0.405,0.212,0.223,0.227,0.25}
437        root:myGlobals:Attenuators:ng3att3_err = {0.3,0.295,0.329,0.263,0.323,0.495,0.307,0.28,0.277,0.3}
438        root:myGlobals:Attenuators:ng3att4_err = {0.35,0.331,0.374,0.303,0.379,0.598,0.367,0.322,0.33,0.35}
439        root:myGlobals:Attenuators:ng3att5_err = {0.4,0.365,0.418,0.355,0.454,0.745,0.411,0.367,0.485,0.4}
440        root:myGlobals:Attenuators:ng3att6_err = {0.45,0.406,0.473,0.385,0.498,0.838,0.454,0.49,0.5,0.5}
441        root:myGlobals:Attenuators:ng3att7_err = {0.6,0.554,0.692,0.425,0.562,0.991,0.715,0.8,0.8,0.8}
442        root:myGlobals:Attenuators:ng3att8_err = {0.7,0.705,0.927,0.503,0.691,1.27,1,1,1,1}
443        root:myGlobals:Attenuators:ng3att9_err = {1,0.862,1.172,0.799,1.104,1.891,1.5,1.5,1.5,1.5}
444        root:myGlobals:Attenuators:ng3att10_err = {1.5,1.054,1.435,1.354,1.742,2,2,2,2,2}
445 
446
447End
448
449
450
451
452
453//returns the transmission of the attenuator (at NG3) given the attenuator number
454//which must be an integer(to select the wave) and given the wavelength.
455//the wavelength may be any value between 4 and 20 (A), and is interpolated
456//between calibrated wavelengths for a given attenuator
457//
458// Mar 2010 - abs() added to attStr to account for ICE reporting -0.0001 as an attenuator position, which truncates to "-0"
459Function LookupAttenNG3(lambda,attenNo,atten_err)
460        Variable lambda, attenNo, &atten_err
461       
462        Variable trans
463        String attStr="root:myGlobals:Attenuators:ng3att"+num2str(trunc(abs(attenNo)))
464        String attErrWStr="root:myGlobals:Attenuators:ng3att"+num2str(trunc(abs(attenNo)))+"_err"
465        String lamStr = "root:myGlobals:Attenuators:ng3lambda"
466       
467        if(attenNo == 0)
468                return (1)              //no attenuation, return trans == 1
469        endif
470       
471        if( (lambda < 4) || (lambda > 20 ) )
472                Abort "Wavelength out of calibration range (4,20). You must manually enter the absolute parameters"
473        Endif
474       
475        if(!(WaveExists($attStr)) || !(WaveExists($lamStr)) || !(WaveExists($attErrWStr)))
476                Execute "MakeNG3AttenTable()"
477        Endif
478        //just in case creating the tables fails....
479        if(!(WaveExists($attStr)) || !(WaveExists($lamStr)) )
480                Abort "Attenuator lookup waves could not be found. You must manually enter the absolute parameters"
481        Endif
482       
483        //lookup the value by interpolating the wavelength
484        //the attenuator must always be an integer
485        Wave att = $attStr
486        Wave attErrW = $attErrWStr
487        Wave lam = $lamstr
488        trans = interp(lambda,lam,att)
489        atten_err = interp(lambda,lam,attErrW)
490
491// the error in the tables is % error. return the standard deviation instead
492        atten_err = trans*atten_err/100
493               
494//      Print "trans = ",trans
495//      Print "trans err = ",atten_err
496       
497        return trans
498End
499
500
501
502
503// a utility function so that I can get the values from the command line
504// since the atten_err is PBR
505//
506Function PrintAttenuation(instr,lam,attenNo)
507        String instr
508        Variable lam,attenNo
509       
510        Variable atten_err, attenFactor
511       
512        // 22 FEB 2013 - not sure what changed with the writeout of ICE data files... but ....
513        // to account for ICE occasionally writing out "3" as 2.9998, make sure I can construct
514        // a single digit -> string "3" to identify the proper wave in the lookup table
515       
516        attenNo = round(attenNo)
517       
518        strswitch(instr)
519                case "CGB":
520                case "NG3":
521                        attenFactor = LookupAttenNG3(lam,attenNo,atten_err)
522                        break
523                default:                                                       
524                        //return error?
525                        DoAlert 0, "No matching instrument -- PrintAttenuation"
526                        attenFactor=1
527        endswitch
528
529        Print "atten, err = ", attenFactor, atten_err
530       
531        return(0)
532End
533
534
535//
536//returns the proper attenuation factor based on the instrument (NG3, NG5, or NG7)
537//NG5 values are taken from the NG7 tables (there is very little difference in the
538//values, and NG5 attenuators have not been calibrated (as of 8/01)
539//
540// filestr is passed from TextRead[3] = the default directory
541// lam is passed from RealsRead[26]
542// AttenNo is passed from ReaslRead[3]
543//
544// Attenuation factor as defined here is <= 1
545//
546// HFIR can pass ("",1,attenuationFactor) and have this function simply
547// spit back the attenuationFactor (that was read into rw[3])
548//
549// called by Correct.ipf, ProtocolAsPanel.ipf, Transmission.ipf
550//
551//
552// as of March 2011, returns the error (one standard deviation) in the attenuation factor as the last parameter, by reference
553Function AttenuationFactor(fileStr,lam,attenNo,atten_err)
554        String fileStr
555        Variable lam,attenNo, &atten_err
556       
557        Variable attenFactor=1,loc
558        String instr=fileStr[1,3]       //filestr is "[NGnSANSn] " or "[NGnSANSnn]" (11 characters total)
559
560
561        // 22 FEB 2013 - not sure what changed with the writeout of ICE data files... but ....
562        // to account for ICE occasionally writing out "3" as 2.9998, make sure I can construct
563        // a single digit -> string "3" to identify the proper wave in the lookup table
564       
565        attenNo = round(attenNo)
566       
567               
568        strswitch(instr)
569                case "CGB":
570                case "NG3":
571                        attenFactor = LookupAttenNG3(lam,attenNo,atten_err)
572                        break
573                default:                                                       
574                        //return error?
575                        DoAlert 0, "No matching instrument -- PrintAttenuation"
576                        attenFactor=1
577        endswitch
578//      print "instr, lambda, attenNo,attenFactor = ",instr,lam,attenNo,attenFactor
579        return(attenFactor)
580End
581
582
Note: See TracBrowser for help on using the repository browser.