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

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

more bug fixes for the VCALC functionality

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