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

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

additions to VCALC procedures to correctly account for panel motion (individual, not symmetric). Updated the plotting routines to all (mostly) pass through the same subroutines so that additional averaging modes will be easier to add.

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// - SANS -- called by CircSectAvg.ipf and RectAnnulAvg.ipf
57//
58// - VSANS -- called in VC_fDoBinning_QxQy2D(folderStr, binningType)
59//
60// DDet is the detector pixel resolution
61// apOff is the offset between the sample aperture and the sample position
62//
63//
64// INPUT:
65// inQ = q-value [1/A]
66// lambda = wavelength [A]
67// lambdaWidth = [dimensionless]
68// DDet = detector pixel resolution [cm]        **assumes square pixel
69// apOff = sample aperture to sample distance [cm]
70// S1 = source aperture diameter [mm]
71// S2 = sample aperture diameter [mm]
72// L1 = source to sample distance [m]
73// L2 = sample to detector distance [m]
74// BS = beam stop diameter [mm]
75// del_r = step size [mm] = binWidth*(mm/pixel)
76// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam
77//
78// OUPUT:
79// SigmaQ
80// QBar
81// fSubS
82//
83//
84Function V_getResolution(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,SigmaQ,QBar,fSubS)
85        Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses
86        Variable &SigmaQ, &QBar, &fSubS         //these are the output quantities at the input Q value
87       
88        //lots of calculation variables
89        Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g
90        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
91
92        //Constants
93        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
94        Variable g = 981.0                              //gravity acceleration [cm/s^2]
95
96
97        S1 *= 0.5*0.1                   //convert to radius and [cm]
98        S2 *= 0.5*0.1
99
100        L1 *= 100.0                     // [cm]
101        L1 -= apOff                             //correct the distance
102
103        L2 *= 100.0
104        L2 += apOff
105        del_r *= 0.1                            //width of annulus, convert mm to [cm]
106       
107        BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm]
108        // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture
109        Variable LB
110        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
111        BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax
112       
113        //Start resolution calculation
114        a2 = S1*L2/L1 + S2*(L1+L2)/L1
115        q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2)
116        lp = 1.0/( 1.0/L1 + 1.0/L2)
117
118        v_lambda = lambdaWidth^2/6.0
119       
120//      if(usingLenses==1)                      //SRK 2007
121        if(usingLenses != 0)                    //SRK 2008 allows for the possibility of different numbers of lenses in header
122                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term
123        else
124                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form
125        endif
126       
127        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
128        vz = vz_1 / lambda
129        yg = 0.5*g*L2*(L1+L2)/vz^2
130        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007
131
132        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
133        delta = 0.5*(BS - r0)^2/v_d
134
135        if (r0 < BS)
136                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
137        else
138                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
139        endif
140
141        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
142        if (fSubS <= 0.0)
143                fSubS = 1.e-10
144        endif
145        fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
146        fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
147
148        rmd = fr*r0
149        v_r1 = v_b + fv*v_d +v_g
150
151        rm = rmd + 0.5*v_r1/rmd
152        v_r = v_r1 - 0.5*(v_r1/rmd)^2
153        if (v_r < 0.0)
154                v_r = 0.0
155        endif
156        QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2))
157        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda)
158
159
160// more readable method for calculating the variance in Q
161// EXCEPT - this is calculated for Qo, NOT qBar
162// (otherwise, they are nearly equivalent, except for close to the beam stop)
163//      Variable kap,a_val,a_val_l2,m_h
164//      g = 981.0                               //gravity acceleration [cm/s^2]
165//      m_h     = 252.8                 // m/h [=] s/cm^2
166//     
167//      kap = 2*pi/lambda
168//      a_val = L2*(L1+L2)*g/2*(m_h)^2
169//      a_val_L2 = a_val/L2*1e-16               //convert 1/cm^2 to 1/A^2
170//
171//      sigmaQ = 0
172//      sigmaQ = 3*(S1/L1)^2
173//     
174//      if(usingLenses != 0)
175//              sigmaQ += 2*(S2/lp)^2*(lambdaWidth)^2   //2nd term w/ lenses
176//      else   
177//              sigmaQ += 2*(S2/lp)^2                                           //2nd term w/ no lenses
178//      endif
179//     
180//      sigmaQ += (del_r/L2)^2
181//      sigmaQ += 2*(r0/L2)^2*(lambdaWidth)^2
182//      sigmaQ += 4*(a_val_l2)^2*lambda^4*(lambdaWidth)^2
183//     
184//      sigmaQ *= kap^2/12
185//      sigmaQ = sqrt(sigmaQ)
186
187
188        Return (0)
189End
190
191
192//
193//**********************
194// 2D resolution function calculation - ***NOT*** in terms of X and Y
195// but written in terms of Parallel and perpendicular to the Q vector at each point
196//
197// -- it is more naturally written this way since the 2D function is an ellipse with its major
198// axis pointing in the direction of Q_parallel. Hence there is no way to properly define the
199// elliptical gaussian in terms of sigmaX and sigmaY
200//
201// For a full description of the gravity effect on the resolution, see:
202//
203//      "The effect of gravity on the resolution of small-angle neutron diffraction peaks"
204//      D.F.R Mildner, J.G. Barker & S.R. Kline J. Appl. Cryst. (2011). 44, 1127-1129.
205//      [ doi:10.1107/S0021889811033322 ]
206//
207//              2/17/12 SRK
208//              NOTE: the first 2/3 of this code is the 1D code, copied here just to have the beam stop
209//                              calculation here, if I decide to implement it. The real calculation is all at the
210//                              bottom and is quite compact
211//
212//
213//
214//
215// - 21 MAR 07 uses projected BS diameter on the detector
216// - APR 07 still need to add resolution with lenses. currently there is no flag in the
217//          raw data header to indicate the presence of lenses.
218//
219// - Aug 07 - added input to switch calculation based on lenses (==1 if in)
220//
221// passed values are read from RealsRead
222// except DDet and apOff, which are set from globals before passing
223//
224// phi is the azimuthal angle, CCW from +x axis
225// r_dist is the real-space distance from ctr of detector to QxQy pixel location
226//
227// MAR 2011 - removed the del_r terms, they don't apply since no bining is done to the 2D data
228//
229Function V_get2DResolution(inQ,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,SigmaQX,SigmaQY,fSubS)
230        Variable inQ, phi,lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses,r_dist
231        Variable &SigmaQX,&SigmaQY,&fSubS               //these are the output quantities at the input Q value
232       
233        //lots of calculation variables
234        Variable a2, lp, v_lambda, v_b, v_d, vz, yg, v_g
235        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
236
237        //Constants
238        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
239        Variable g = 981.0                              //gravity acceleration [cm/s^2]
240        Variable m_h    = 252.8                 // m/h [=] s/cm^2
241
242
243        S1 *= 0.5*0.1                   //convert to radius and [cm]
244        S2 *= 0.5*0.1
245
246        L1 *= 100.0                     // [cm]
247        L1 -= apOff                             //correct the distance
248
249        L2 *= 100.0
250        L2 += apOff
251        del_r *= 0.1                            //width of annulus, convert mm to [cm]
252       
253        BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm]
254        // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture
255        Variable LB
256        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
257        BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax
258       
259        //Start resolution calculation
260        a2 = S1*L2/L1 + S2*(L1+L2)/L1
261        lp = 1.0/( 1.0/L1 + 1.0/L2)
262
263        v_lambda = lambdaWidth^2/6.0
264       
265//      if(usingLenses==1)                      //SRK 2007
266        if(usingLenses != 0)                    //SRK 2008 allows for the possibility of different numbers of lenses in header
267                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term
268        else
269                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form
270        endif
271       
272        v_d = (DDet/2.3548)^2 + del_r^2/12.0
273        vz = vz_1 / lambda
274        yg = 0.5*g*L2*(L1+L2)/vz^2
275        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007
276
277        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
278        delta = 0.5*(BS - r0)^2/v_d
279
280        if (r0 < BS)
281                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
282        else
283                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
284        endif
285
286        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
287        if (fSubS <= 0.0)
288                fSubS = 1.e-10
289        endif
290//      fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
291//      fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
292//
293//      rmd = fr*r0
294//      v_r1 = v_b + fv*v_d +v_g
295//
296//      rm = rmd + 0.5*v_r1/rmd
297//      v_r = v_r1 - 0.5*(v_r1/rmd)^2
298//      if (v_r < 0.0)
299//              v_r = 0.0
300//      endif
301
302        Variable kap,a_val,a_val_L2,proj_DDet
303       
304        kap = 2*pi/lambda
305        a_val = L2*(L1+L2)*g/2*(m_h)^2
306        a_val_L2 = a_val/L2*1e-16               //convert 1/cm^2 to 1/A^2
307
308
309        // the detector pixel is square, so correct for phi
310        proj_DDet = DDet*cos(phi) + DDet*sin(phi)
311
312
313///////// OLD - don't use ---
314//in terms of Q_parallel ("x") and Q_perp ("y") - this works, since parallel is in the direction of Q and I
315// can calculate that from the QxQy (I just need the projection)
316//// for test case with no gravity, set a_val = 0
317//// note that gravity has no wavelength dependence. the lambda^4 cancels out.
318////
319////    a_val = 0
320////    a_val_l2 = 0
321//
322//     
323//      // this is really sigma_Q_parallel
324//      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)
325//      SigmaQX += inQ*inQ*v_lambda
326//
327//      //this is really sigma_Q_perpendicular
328//      proj_DDet = DDet*sin(phi) + DDet*cos(phi)               //not necessary, since DDet is the same in both X and Y directions
329//
330//      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)
331//     
332//      SigmaQX = sqrt(SigmaQX)
333//      SigmaQy = sqrt(SigmaQY)
334//     
335
336/////////////////////////////////////////////////
337/////   
338//      ////// this is all new, inclusion of gravity effect into the parallel component
339//       perpendicular component is purely geometric, no gravity component
340//
341// the shadow factor is calculated as above -so keep the above calculations, even though
342// most of them are redundant.
343//
344       
345////    //
346        Variable yg_d,acc,sdd,ssd,lambda0,DL_L,sig_l
347        Variable var_qlx,var_qly,var_ql,qx,qy,sig_perp,sig_para, sig_para_new
348       
349        G = 981.  //!   ACCELERATION OF GRAVITY, CM/SEC^2
350        acc = vz_1              //      3.956E5 //!     CONVERT WAVELENGTH TO VELOCITY CM/SEC
351        SDD = L2                //1317
352        SSD = L1                //1627          //cm
353        lambda0 = lambda                //              15
354        DL_L = lambdaWidth              //0.236
355        SIG_L = DL_L/sqrt(6)
356        YG_d = -0.5*G*SDD*(SSD+SDD)*(LAMBDA0/acc)^2
357/////   Print "DISTANCE BEAM FALLS DUE TO GRAVITY (CM) = ",YG
358//              Print "Gravity q* = ",-2*pi/lambda0*2*yg_d/sdd
359       
360        sig_perp = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2)
361        sig_perp = sqrt(sig_perp)
362       
363// TODO -- not needed???       
364//      FindQxQy(inQ,phi,qx,qy)
365
366
367// missing a factor of 2 here, and the form is different than the paper, so re-write   
368//      VAR_QLY = SIG_L^2 * (QY+4*PI*YG_d/(2*SDD*LAMBDA0))^2
369//      VAR_QLX = (SIG_L*QX)^2
370//      VAR_QL = VAR_QLY + VAR_QLX  //! WAVELENGTH CONTRIBUTION TO VARIANCE
371//      sig_para = (sig_perp^2 + VAR_QL)^0.5
372       
373        // r_dist is passed in, [=]cm
374        // from the paper
375        a_val = 0.5*G*SDD*(SSD+SDD)*m_h^2 * 1e-16               //units now are cm /(A^2)
376       
377        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)
378        sig_para_new = (sig_perp^2 + VAR_QL)^0.5
379       
380       
381///// return values PBR
382        SigmaQX = sig_para_new
383        SigmaQy = sig_perp
384       
385////   
386
387        Return (0)
388End
389
390
391
Note: See TracBrowser for help on using the repository browser.