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

Last change on this file since 1119 was 1119, checked in by srkline, 4 years ago

added procedures to output QxQy_ASCII data. Each panel is output into its own file. Output format is the same as for 2D SANS data, including the 2D resolution function. However, reading the data back in with the SANS analysis macros currently does not redimension the data back to the matrix correctly as it assumes a square detector.

I will need to add the X and Y dimensions of each panel into the header, and then make use of these values when they are read in. Also, writing the QxQy? data is quick for the M and F panels ( 0.8 s) but is extremely slow for the back, High-res panel (120 s) since there are 1.1.E6 points there vs. 6144 pts. I'll need to find a way to speed this operation up.

File size: 28.0 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// folderStr = folder with the current reduction step
67// type = binning type (not the same as the detStr)
68// collimationStr = collimation type, to switch for lenses, etc.
69
70// READ/DERIVED within the function
71// lambda = wavelength [A]
72// lambdaWidth = [dimensionless]
73// DDet = detector pixel resolution [cm]        **assumes square pixel
74// apOff = sample aperture to sample distance [cm]
75// S1 = source aperture diameter [mm]
76// S2 = sample aperture diameter [mm]
77// L1 = source to sample distance [m]
78// L2 = sample to detector distance [m]
79// BS = beam stop diameter [mm]
80// del_r = step size [mm] = binWidth*(mm/pixel)
81// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam
82//
83// OUPUT:
84// SigmaQ
85// QBar
86// fSubS
87//
88//
89Function V_getResolution(inQ,folderStr,type,collimationStr,SigmaQ,QBar,fSubS)
90        Variable inQ
91        String folderStr,type,collimationStr
92        Variable &SigmaQ, &QBar, &fSubS         //these are the output quantities at the input Q value
93
94
95        Variable lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
96               
97        //lots of calculation variables
98        Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g
99        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
100
101        //Constants
102        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
103        Variable g = 981.0                              //gravity acceleration [cm/s^2]
104
105///////// get all of the values from the header
106// TODO: check the units of all of the inputs
107// lambda = wavelength [A]
108        lambda = V_getWavelength(folderStr)
109       
110// lambdaWidth = [dimensionless]
111        lambdaWidth = V_getWavelength_spread(folderStr)
112       
113// DDet = detector pixel resolution [cm]        **assumes square pixel
114        // V_getDet_pixel_fwhm_x(folderStr,detStr)
115        // V_getDet_pixel_fwhm_y(folderStr,detStr)
116//      DDet = 0.8              // TODO -- this is hard-wired
117
118        if(strlen(type) == 1)
119                // it's "B"
120                DDet = V_getDet_pixel_fwhm_x(folderStr,type)            // value is already in cm
121        else
122                DDet = V_getDet_pixel_fwhm_x(folderStr,type[0,1])               // value is already in cm
123        endif
124               
125// apOff = sample aperture to sample distance [cm]
126        apOff = 10              // TODO -- this is hard-wired
127
128
129       
130// S1 = source aperture diameter [mm]
131// may be either circle or rectangle
132        String s1_shape="",bs_shape=""
133        Variable width,height,equiv_S1,equiv_bs
134       
135       
136        s1_shape = V_getSourceAp_shape(folderStr)
137        if(cmpstr(s1_shape,"CIRCLE") == 0)
138                S1 = str2num(V_getSourceAp_size(folderStr))
139        else
140                S1 = V_getSourceAp_height(folderStr)            // TODO: need the width or at least an equivalent diameter
141        endif
142       
143       
144// S2 = sample aperture diameter [cm]
145// as of 3/2018, the "internal" sample aperture is not in use, only the external
146// TODO : verify the units on the Ap2 (external)
147// sample aperture 1(internal) is set to report "12.7 mm" as a STRING
148// sample aperture 2(external) reports the number typed in...
149//
150// so I'm trusting [cm] is in the file
151        S2 = V_getSampleAp2_size(folderStr)*10          // sample ap 1 or 2? 2 = the "external", convert to [mm]
152       
153// L1 = source to sample distance [m]
154        L1 = V_getSourceAp_distance(folderStr)/100
155
156// L2 = sample to detector distance [m]
157// take the first two characters of the "type" to get the correct distance.
158// if the type is say, MLRTB, then the implicit assumption in combining all four panels is that the resolution
159// is not an issue for the slightly different distances.
160        if(strlen(type) == 1)
161                // it's "B"
162                L2 = V_getDet_ActualDistance(folderStr,type)/100                //convert cm to m
163        else
164                L2 = V_getDet_ActualDistance(folderStr,type[0,1])/100           //convert cm to m
165        endif
166       
167// BS = beam stop diameter [mm]
168//TODO:? which BS is in? carr2, carr3, none?
169// -- need to check the detector, num_beamstops field, then description, then shape/size or shape/height and shape/width
170//
171// TODO: the values in the file are incorrect!!! BS = 1000 mm diameter!!!
172        BS = V_DeduceBeamstopDiameter(folderStr,type)           //returns diameter in [mm]
173//      BS = V_getBeamStopC2_size(folderStr)            // Units are [mm]
174//      BS = 25.4                       //TODO hard-wired value
175       
176//      bs_shape = V_getBeamStopC2_shape(folderStr)
177//      if(cmpstr(s1_shape,"CIRCLE") == 0)
178//              bs = V_getBeamStopC2_size(folderStr)
179//      else
180//              bs = V_getBeamStopC2_height(folderStr) 
181//      endif
182
183
184       
185// del_r = step size [mm] = binWidth*(mm/pixel)
186        del_r = 1*DDet*10               // TODO: this is probably not the correct value
187
188// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam
189        usingLenses = 0
190
191//if(cmpstr(type[0,1],"FL")==0)
192//      Print "(FL) Resolution lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses"
193//      Print lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
194//endif
195
196
197
198// TODO:
199// this is the point where I need to switch on the different collimation types (white beam, slit, Xtal, etc)
200// to calculate the correct resolution, or fill the waves with the correct "flags"
201//
202
203// For white beam data, the wavelength distribution can't be represented as a gaussian, but all of the other
204//  geometric corrections still apply. Passing zero for the lambdaWidth will return the geometry contribution,
205//  as long as the wavelength can be handled separately. It appears to be correct to do as a double integral,
206//  with the inner(lambda) calculated first, then the outer(geometry).
207//
208
209// possible values are:
210//
211// pinhole
212// pinhole_whiteBeam
213// convergingPinholes
214//
215// *slit data should be reduced using the slit routine, not here, proceed but warn
216// narrowSlit
217// narrowSlit_whiteBeam
218
219
220        if(cmpstr(collimationStr,"pinhole") == 0)
221                //nothing to change     
222        endif
223
224        if(cmpstr(collimationStr,"pinhole_whiteBeam") == 0)
225                //              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution.
226                // the white beam distribution will need to be flagged some other way
227                //
228                lambdaWidth = 0
229        endif
230
231        if(cmpstr(collimationStr,"convergingPinholes") == 0)
232
233                //              set usingLenses == 1 so that the Gaussian resolution calculation will be for a focus condition
234                usingLenses = 1
235        endif
236
237
238// should not end up here, except for odd testing cases
239        if(cmpstr(collimationStr,"narrowSlit") == 0)
240
241                Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???"
242        endif
243       
244// should not end up here, except for odd testing cases
245        if(cmpstr(collimationStr,"narrowSlit_whiteBeam") == 0)
246
247                //              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution.
248                // the white beam distribution will need to be flagged some other way
249                //
250                Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???"
251
252                lambdaWidth = 0
253        endif
254
255
256/////////////////////////////
257/////////////////////////////
258// do the calculation
259        S1 *= 0.5*0.1                   //convert to radius and [cm]
260        S2 *= 0.5*0.1
261
262        L1 *= 100.0                     // [cm]
263        L1 -= apOff                             //correct the distance
264
265        L2 *= 100.0
266        L2 += apOff
267        del_r *= 0.1                            //width of annulus, convert mm to [cm]
268       
269        BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm]
270        // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture
271        Variable LB
272        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
273        BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax
274       
275        //Start resolution calculation
276        a2 = S1*L2/L1 + S2*(L1+L2)/L1
277        q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2)
278        lp = 1.0/( 1.0/L1 + 1.0/L2)
279
280        v_lambda = lambdaWidth^2/6.0
281       
282//      if(usingLenses==1)                      //SRK 2007
283        if(usingLenses != 0)                    //SRK 2008 allows for the possibility of different numbers of lenses in header
284                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term
285        else
286                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form
287        endif
288       
289        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
290        vz = vz_1 / lambda
291        yg = 0.5*g*L2*(L1+L2)/vz^2
292        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007
293
294        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
295        delta = 0.5*(BS - r0)^2/v_d
296
297        if (r0 < BS)
298                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
299        else
300                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
301        endif
302
303        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
304        if (fSubS <= 0.0)
305                fSubS = 1.e-10
306        endif
307        fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
308        fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
309
310        rmd = fr*r0
311        v_r1 = v_b + fv*v_d +v_g
312
313        rm = rmd + 0.5*v_r1/rmd
314        v_r = v_r1 - 0.5*(v_r1/rmd)^2
315        if (v_r < 0.0)
316                v_r = 0.0
317        endif
318        QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2))
319        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda)
320
321
322// more readable method for calculating the variance in Q
323// EXCEPT - this is calculated for Qo, NOT qBar
324// (otherwise, they are nearly equivalent, except for close to the beam stop)
325//      Variable kap,a_val,a_val_l2,m_h
326//      g = 981.0                               //gravity acceleration [cm/s^2]
327//      m_h     = 252.8                 // m/h [=] s/cm^2
328//     
329//      kap = 2*pi/lambda
330//      a_val = L2*(L1+L2)*g/2*(m_h)^2
331//      a_val_L2 = a_val/L2*1e-16               //convert 1/cm^2 to 1/A^2
332//
333//      sigmaQ = 0
334//      sigmaQ = 3*(S1/L1)^2
335//     
336//      if(usingLenses != 0)
337//              sigmaQ += 2*(S2/lp)^2*(lambdaWidth)^2   //2nd term w/ lenses
338//      else   
339//              sigmaQ += 2*(S2/lp)^2                                           //2nd term w/ no lenses
340//      endif
341//     
342//      sigmaQ += (del_r/L2)^2
343//      sigmaQ += 2*(r0/L2)^2*(lambdaWidth)^2
344//      sigmaQ += 4*(a_val_l2)^2*lambda^4*(lambdaWidth)^2
345//     
346//      sigmaQ *= kap^2/12
347//      sigmaQ = sqrt(sigmaQ)
348
349
350        Return (0)
351End
352
353
354//
355//**********************
356// 2D resolution function calculation - ***NOT*** in terms of X and Y
357// but written in terms of Parallel and perpendicular to the Q vector at each point
358//
359// -- it is more naturally written this way since the 2D function is an ellipse with its major
360// axis pointing in the direction of Q_parallel. Hence there is no way to properly define the
361// elliptical gaussian in terms of sigmaX and sigmaY
362//
363// For a full description of the gravity effect on the resolution, see:
364//
365//      "The effect of gravity on the resolution of small-angle neutron diffraction peaks"
366//      D.F.R Mildner, J.G. Barker & S.R. Kline J. Appl. Cryst. (2011). 44, 1127-1129.
367//      [ doi:10.1107/S0021889811033322 ]
368//
369//              2/17/12 SRK
370//              NOTE: the first 2/3 of this code is the 1D code, copied here just to have the beam stop
371//                              calculation here, if I decide to implement it. The real calculation is all at the
372//                              bottom and is quite compact
373//
374//
375//
376//
377// - 21 MAR 07 uses projected BS diameter on the detector
378// - APR 07 still need to add resolution with lenses. currently there is no flag in the
379//          raw data header to indicate the presence of lenses.
380//
381// - Aug 07 - added input to switch calculation based on lenses (==1 if in)
382//
383// passed values are read from RealsRead
384// except DDet and apOff, which are set from globals before passing
385//
386// phi is the azimuthal angle, CCW from +x axis
387// r_dist is the real-space distance from ctr of detector to QxQy pixel location
388//
389// MAR 2011 - removed the del_r terms, they don't apply since no binning is done to the 2D data
390//
391Function V_get2DResolution(inQ,phi,r_dist,folderStr,type,collimationStr,SigmaQX,SigmaQY,fSubS)
392        Variable inQ,phi,r_dist
393        String folderStr,type,collimationStr
394        Variable &SigmaQX,&SigmaQY,&fSubS               //these are the output quantities at the input Q value
395//      Variable SigmaQX,SigmaQY,fSubS          //these are the output quantities at the input Q value
396
397       
398        Variable lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses
399
400        //      phi = FindPhi( pixSize*((p+1)-xctr) , pixSize*((q+1)-yctr)+(2)*yg_d)            //(dx,dy+yg_d)
401        //      r_dist = sqrt(  (pixSize*((p+1)-xctr))^2 +  (pixSize*((q+1)-yctr)+(2)*yg_d)^2 )         //radial distance from ctr to pt
402
403///////// get all of the values from the header
404// TODO: check the units of all of the inputs
405// lambda = wavelength [A]
406        lambda = V_getWavelength(folderStr)
407       
408// lambdaWidth = [dimensionless]
409        lambdaWidth = V_getWavelength_spread(folderStr)
410       
411// DDet = detector pixel resolution [cm]        **assumes square pixel
412        // V_getDet_pixel_fwhm_x(folderStr,detStr)
413        // V_getDet_pixel_fwhm_y(folderStr,detStr)
414//      DDet = 0.8              // TODO -- this is hard-wired
415
416        if(strlen(type) == 1)
417                // it's "B"
418                DDet = V_getDet_pixel_fwhm_x(folderStr,type)            // value is already in cm
419        else
420                DDet = V_getDet_pixel_fwhm_x(folderStr,type[0,1])               // value is already in cm
421        endif
422               
423// apOff = sample aperture to sample distance [cm]
424        apOff = 10              // TODO -- this is hard-wired
425
426
427// S1 = source aperture diameter [mm]
428// may be either circle or rectangle
429        String s1_shape="",bs_shape=""
430        Variable width,height,equiv_S1,equiv_bs
431       
432        s1_shape = V_getSourceAp_shape(folderStr)
433        if(cmpstr(s1_shape,"CIRCLE") == 0)
434                S1 = str2num(V_getSourceAp_size(folderStr))
435        else
436                S1 = V_getSourceAp_height(folderStr)            // TODO: need the width or at least an equivalent diameter
437        endif
438       
439       
440// S2 = sample aperture diameter [cm]
441// as of 3/2018, the "internal" sample aperture is not in use, only the external
442// TODO : verify the units on the Ap2 (external)
443// sample aperture 1(internal) is set to report "12.7 mm" as a STRING
444// sample aperture 2(external) reports the number typed in...
445//
446// so I'm trusting [cm] is in the file
447        S2 = V_getSampleAp2_size(folderStr)*10          // sample ap 1 or 2? 2 = the "external", convert to [mm]
448       
449// L1 = source to sample distance [m]
450        L1 = V_getSourceAp_distance(folderStr)/100
451
452// L2 = sample to detector distance [m]
453// take the first two characters of the "type" to get the correct distance.
454// if the type is say, MLRTB, then the implicit assumption in combining all four panels is that the resolution
455// is not an issue for the slightly different distances.
456        if(strlen(type) == 1)
457                // it's "B"
458                L2 = V_getDet_ActualDistance(folderStr,type)/100                //convert cm to m
459        else
460                L2 = V_getDet_ActualDistance(folderStr,type[0,1])/100           //convert cm to m
461        endif
462       
463// BS = beam stop diameter [mm]
464//TODO:? which BS is in? carr2, carr3, none?
465// -- need to check the detector, num_beamstops field, then description, then shape/size or shape/height and shape/width
466//
467// TODO: the values in the file are incorrect!!! BS = 1000 mm diameter!!!
468        BS = V_DeduceBeamstopDiameter(folderStr,type)           //returns diameter in [mm]
469//      BS = V_getBeamStopC2_size(folderStr)            // Units are [mm]
470//      BS = 25.4                       //TODO hard-wired value
471       
472//      bs_shape = V_getBeamStopC2_shape(folderStr)
473//      if(cmpstr(s1_shape,"CIRCLE") == 0)
474//              bs = V_getBeamStopC2_size(folderStr)
475//      else
476//              bs = V_getBeamStopC2_height(folderStr) 
477//      endif
478
479
480       
481// del_r = step size [mm] = binWidth*(mm/pixel)
482        del_r = 1*DDet*10               // TODO: this is probably not the correct value
483
484// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam
485        usingLenses = 0
486
487//if(cmpstr(type[0,1],"FL")==0)
488//      Print "(FL) Resolution lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses"
489//      Print lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
490//endif
491
492
493
494// TODO:
495// this is the point where I need to switch on the different collimation types (white beam, slit, Xtal, etc)
496// to calculate the correct resolution, or fill the waves with the correct "flags"
497//
498
499// For white beam data, the wavelength distribution can't be represented as a gaussian, but all of the other
500//  geometric corrections still apply. Passing zero for the lambdaWidth will return the geometry contribution,
501//  as long as the wavelength can be handled separately. It appears to be correct to do as a double integral,
502//  with the inner(lambda) calculated first, then the outer(geometry).
503//
504
505// possible values are:
506//
507// pinhole
508// pinhole_whiteBeam
509// convergingPinholes
510//
511// *slit data should be reduced using the slit routine, not here, proceed but warn
512// narrowSlit
513// narrowSlit_whiteBeam
514
515
516        if(cmpstr(collimationStr,"pinhole") == 0)
517                //nothing to change     
518        endif
519
520        if(cmpstr(collimationStr,"pinhole_whiteBeam") == 0)
521                //              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution.
522                // the white beam distribution will need to be flagged some other way
523                //
524                lambdaWidth = 0
525        endif
526
527        if(cmpstr(collimationStr,"convergingPinholes") == 0)
528
529                //              set usingLenses == 1 so that the Gaussian resolution calculation will be for a focus condition
530                usingLenses = 1
531        endif
532
533
534// should not end up here, except for odd testing cases
535        if(cmpstr(collimationStr,"narrowSlit") == 0)
536
537                Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???"
538        endif
539       
540// should not end up here, except for odd testing cases
541        if(cmpstr(collimationStr,"narrowSlit_whiteBeam") == 0)
542
543                //              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution.
544                // the white beam distribution will need to be flagged some other way
545                //
546                Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???"
547
548                lambdaWidth = 0
549        endif
550
551
552       
553        //lots of calculation variables
554        Variable a2, lp, v_lambda, v_b, v_d, vz, yg, v_g
555        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
556
557        //Constants
558        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
559        Variable g = 981.0                              //gravity acceleration [cm/s^2]
560        Variable m_h    = 252.8                 // m/h [=] s/cm^2
561
562
563        S1 *= 0.5*0.1                   //convert to radius and [cm]
564        S2 *= 0.5*0.1
565
566        L1 *= 100.0                     // [cm]
567        L1 -= apOff                             //correct the distance
568
569        L2 *= 100.0
570        L2 += apOff
571        del_r *= 0.1                            //width of annulus, convert mm to [cm]
572       
573        BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm]
574        // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture
575        Variable LB
576        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
577        BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax
578       
579        //Start resolution calculation
580        a2 = S1*L2/L1 + S2*(L1+L2)/L1
581        lp = 1.0/( 1.0/L1 + 1.0/L2)
582
583        v_lambda = lambdaWidth^2/6.0
584       
585//      if(usingLenses==1)                      //SRK 2007
586        if(usingLenses != 0)                    //SRK 2008 allows for the possibility of different numbers of lenses in header
587                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term
588        else
589                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form
590        endif
591       
592        v_d = (DDet/2.3548)^2 + del_r^2/12.0
593        vz = vz_1 / lambda
594        yg = 0.5*g*L2*(L1+L2)/vz^2
595        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007
596
597        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
598        delta = 0.5*(BS - r0)^2/v_d
599
600        if (r0 < BS)
601                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
602        else
603                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
604        endif
605
606        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
607        if (fSubS <= 0.0)
608                fSubS = 1.e-10
609        endif
610//      fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
611//      fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
612//
613//      rmd = fr*r0
614//      v_r1 = v_b + fv*v_d +v_g
615//
616//      rm = rmd + 0.5*v_r1/rmd
617//      v_r = v_r1 - 0.5*(v_r1/rmd)^2
618//      if (v_r < 0.0)
619//              v_r = 0.0
620//      endif
621
622        Variable kap,a_val,a_val_L2,proj_DDet
623       
624        kap = 2*pi/lambda
625        a_val = L2*(L1+L2)*g/2*(m_h)^2
626        a_val_L2 = a_val/L2*1e-16               //convert 1/cm^2 to 1/A^2
627
628
629        // the detector pixel is square, so correct for phi
630        proj_DDet = DDet*cos(phi) + DDet*sin(phi)
631
632
633///////// OLD - don't use ---
634//in terms of Q_parallel ("x") and Q_perp ("y") - this works, since parallel is in the direction of Q and I
635// can calculate that from the QxQy (I just need the projection)
636//// for test case with no gravity, set a_val = 0
637//// note that gravity has no wavelength dependence. the lambda^4 cancels out.
638////
639////    a_val = 0
640////    a_val_l2 = 0
641//
642//     
643//      // this is really sigma_Q_parallel
644//      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)
645//      SigmaQX += inQ*inQ*v_lambda
646//
647//      //this is really sigma_Q_perpendicular
648//      proj_DDet = DDet*sin(phi) + DDet*cos(phi)               //not necessary, since DDet is the same in both X and Y directions
649//
650//      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)
651//     
652//      SigmaQX = sqrt(SigmaQX)
653//      SigmaQy = sqrt(SigmaQY)
654//     
655
656/////////////////////////////////////////////////
657/////   
658//      ////// this is all new, inclusion of gravity effect into the parallel component
659//       perpendicular component is purely geometric, no gravity component
660//
661// the shadow factor is calculated as above -so keep the above calculations, even though
662// most of them are redundant.
663//
664       
665////    //
666        Variable yg_d,acc,sdd,ssd,lambda0,DL_L,sig_l
667        Variable var_qlx,var_qly,var_ql,qx,qy,sig_perp,sig_para, sig_para_new
668       
669        G = 981.  //!   ACCELERATION OF GRAVITY, CM/SEC^2
670        acc = vz_1              //      3.956E5 //!     CONVERT WAVELENGTH TO VELOCITY CM/SEC
671        SDD = L2                //1317
672        SSD = L1                //1627          //cm
673        lambda0 = lambda                //              15
674        DL_L = lambdaWidth              //0.236
675        SIG_L = DL_L/sqrt(6)
676        YG_d = -0.5*G*SDD*(SSD+SDD)*(LAMBDA0/acc)^2
677/////   Print "DISTANCE BEAM FALLS DUE TO GRAVITY (CM) = ",YG
678//              Print "Gravity q* = ",-2*pi/lambda0*2*yg_d/sdd
679       
680        sig_perp = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2)
681        sig_perp = sqrt(sig_perp)
682       
683// TODO -- not needed???       
684//      FindQxQy(inQ,phi,qx,qy)
685
686
687// missing a factor of 2 here, and the form is different than the paper, so re-write   
688//      VAR_QLY = SIG_L^2 * (QY+4*PI*YG_d/(2*SDD*LAMBDA0))^2
689//      VAR_QLX = (SIG_L*QX)^2
690//      VAR_QL = VAR_QLY + VAR_QLX  //! WAVELENGTH CONTRIBUTION TO VARIANCE
691//      sig_para = (sig_perp^2 + VAR_QL)^0.5
692       
693        // r_dist is passed in, [=]cm
694        // from the paper
695        a_val = 0.5*G*SDD*(SSD+SDD)*m_h^2 * 1e-16               //units now are cm /(A^2)
696       
697        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)
698        sig_para_new = (sig_perp^2 + VAR_QL)^0.5
699       
700       
701///// return values PBR
702        SigmaQX = sig_para_new
703        SigmaQy = sig_perp
704       
705////   
706
707        Return (0)
708End
709
710
711//**********************
712// Resolution calculation - used by the averaging routines
713// to calculate the resolution function at each q-value
714// - the return value is not used
715//
716// equivalent to John's routine on the VAX Q_SIGMA_AVE.FOR
717// Incorporates eqn. 3-15 from J. Appl. Cryst. (1995) v. 28 p105-114
718//
719// - 21 MAR 07 uses projected BS diameter on the detector
720// - APR 07 still need to add resolution with lenses. currently there is no flag in the
721//          raw data header to indicate the presence of lenses.
722//
723// - Aug 07 - added input to switch calculation based on lenses (==1 if in)
724//
725// - SANS -- called by CircSectAvg.ipf and RectAnnulAvg.ipf
726//
727// - VSANS -- called in VC_fDoBinning_QxQy2D(folderStr, binningType)
728//
729// DDet is the detector pixel resolution
730// apOff is the offset between the sample aperture and the sample position
731//
732//
733// INPUT:
734// inQ = q-value [1/A]
735// lambda = wavelength [A]
736// lambdaWidth = [dimensionless]
737// DDet = detector pixel resolution [cm]        **assumes square pixel
738// apOff = sample aperture to sample distance [cm]
739// S1 = source aperture diameter [mm]
740// S2 = sample aperture diameter [mm]
741// L1 = source to sample distance [m]
742// L2 = sample to detector distance [m]
743// BS = beam stop diameter [mm]
744// del_r = step size [mm] = binWidth*(mm/pixel)
745// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam
746//
747// OUPUT:
748// SigmaQ
749// QBar
750// fSubS
751//
752//
753Function V_getResolution_old(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,SigmaQ,QBar,fSubS)
754        Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses
755        Variable &SigmaQ, &QBar, &fSubS         //these are the output quantities at the input Q value
756       
757        //lots of calculation variables
758        Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g
759        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
760
761        //Constants
762        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
763        Variable g = 981.0                              //gravity acceleration [cm/s^2]
764
765
766        S1 *= 0.5*0.1                   //convert to radius and [cm]
767        S2 *= 0.5*0.1
768
769        L1 *= 100.0                     // [cm]
770        L1 -= apOff                             //correct the distance
771
772        L2 *= 100.0
773        L2 += apOff
774        del_r *= 0.1                            //width of annulus, convert mm to [cm]
775       
776        BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm]
777        // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture
778        Variable LB
779        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
780        BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax
781       
782        //Start resolution calculation
783        a2 = S1*L2/L1 + S2*(L1+L2)/L1
784        q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2)
785        lp = 1.0/( 1.0/L1 + 1.0/L2)
786
787        v_lambda = lambdaWidth^2/6.0
788       
789//      if(usingLenses==1)                      //SRK 2007
790        if(usingLenses != 0)                    //SRK 2008 allows for the possibility of different numbers of lenses in header
791                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term
792        else
793                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form
794        endif
795       
796        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
797        vz = vz_1 / lambda
798        yg = 0.5*g*L2*(L1+L2)/vz^2
799        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007
800
801        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
802        delta = 0.5*(BS - r0)^2/v_d
803
804        if (r0 < BS)
805                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
806        else
807                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
808        endif
809
810        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
811        if (fSubS <= 0.0)
812                fSubS = 1.e-10
813        endif
814        fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
815        fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
816
817        rmd = fr*r0
818        v_r1 = v_b + fv*v_d +v_g
819
820        rm = rmd + 0.5*v_r1/rmd
821        v_r = v_r1 - 0.5*(v_r1/rmd)^2
822        if (v_r < 0.0)
823                v_r = 0.0
824        endif
825        QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2))
826        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda)
827
828
829// more readable method for calculating the variance in Q
830// EXCEPT - this is calculated for Qo, NOT qBar
831// (otherwise, they are nearly equivalent, except for close to the beam stop)
832//      Variable kap,a_val,a_val_l2,m_h
833//      g = 981.0                               //gravity acceleration [cm/s^2]
834//      m_h     = 252.8                 // m/h [=] s/cm^2
835//     
836//      kap = 2*pi/lambda
837//      a_val = L2*(L1+L2)*g/2*(m_h)^2
838//      a_val_L2 = a_val/L2*1e-16               //convert 1/cm^2 to 1/A^2
839//
840//      sigmaQ = 0
841//      sigmaQ = 3*(S1/L1)^2
842//     
843//      if(usingLenses != 0)
844//              sigmaQ += 2*(S2/lp)^2*(lambdaWidth)^2   //2nd term w/ lenses
845//      else   
846//              sigmaQ += 2*(S2/lp)^2                                           //2nd term w/ no lenses
847//      endif
848//     
849//      sigmaQ += (del_r/L2)^2
850//      sigmaQ += 2*(r0/L2)^2*(lambdaWidth)^2
851//      sigmaQ += 4*(a_val_l2)^2*lambda^4*(lambdaWidth)^2
852//     
853//      sigmaQ *= kap^2/12
854//      sigmaQ = sqrt(sigmaQ)
855
856
857        Return (0)
858End
859
860
Note: See TracBrowser for help on using the repository browser.