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

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

multiple changes to improve the functionality of VCALC

more values are reported, and the IQ plot now accounts for a beam stop shadowing the low q region. Qmin and qmax values are reported for each panel. the beam intensity value is more realistic, with correct SSD values.

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