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

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

added a stub procedure file to hold the instrument resolution calculations. Currently none is functional, since even the SANS-like calculations are missing lots of physical dimensions, units, etc.

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