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

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

a few changes to update the calculation of transmission values from the tables, and updating the correct units in a few places. Also updated how a transmission file is matched with a scattering file.
File size: 12.7 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
Note: See TracBrowser for help on using the repository browser.