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

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

many minor changes after real VSANS data collected.

additional procedures added to allow easy correction of the incorrect header information from NICE.

Most notable addition is the pinhole resolution added to the calculation and the I(q) output. White beam is also treated (incorrectly) as a gaussian distrivution, but the results of smeared fitting look to be quite good.

Trimming and sorting routines are now (pinhole) resolution aware.

File identification routines have been updated to use the proper definitions of "purpose" and "intent". Both fields are now in the catalog, to allow for better sorting.

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