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

Last change on this file since 1242 was 1242, checked in by srkline, 3 years ago

updating the IgorVersion? pragma to v7.0 for all files to be consistent.

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