Ignore:
Timestamp:
Feb 1, 2019 2:25:12 PM (4 years ago)
Author:
srkline
Message:

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

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Instrument_Resolution.ipf

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