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

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

lots of changes here:
many little fixes to clean up TODO items and marke them DONE

changed the handling of the panel "gap" to split the gap evenly. Q-calculations have been re-verified with this change.

re-named the list of "bin Type" values, and added a few more choices. Streamlined how the averaging and plotting works with this list so that it can be more easily modified as different combinations of binning are envisioned. This resulted in a lot of excess code being cut out and replaced with cleaner logic. This change has also been verified to work as intended.

Attenuation is now always calculated from the table. The table also by (NEW) definition has values for the white beam (one waelength) and graphite (multiple possible wavelengths) where the wavelengths are artificially scaled (*1000) or *1e6) so that the interpolations can be done internally without the need for multiple attenuator tables.

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 V_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.