Changeset 1120 for sans/Dev/trunk/NCNR_User_Procedures/Reduction
- Timestamp:
- Feb 1, 2019 3:47:27 PM (4 years ago)
- Location:
- sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Instrument_Resolution.ipf
r1119 r1120 709 709 710 710 711 //********************** 712 // Resolution calculation - used by the averaging routines 713 // to calculate the resolution function at each q-value 714 // - the return value is not used 715 // 716 // equivalent to John's routine on the VAX Q_SIGMA_AVE.FOR 717 // Incorporates eqn. 3-15 from J. Appl. Cryst. (1995) v. 28 p105-114 718 // 719 // - 21 MAR 07 uses projected BS diameter on the detector 720 // - APR 07 still need to add resolution with lenses. currently there is no flag in the 721 // raw data header to indicate the presence of lenses. 722 // 723 // - Aug 07 - added input to switch calculation based on lenses (==1 if in) 724 // 725 // - SANS -- called by CircSectAvg.ipf and RectAnnulAvg.ipf 726 // 727 // - VSANS -- called in VC_fDoBinning_QxQy2D(folderStr, binningType) 728 // 729 // DDet is the detector pixel resolution 730 // apOff is the offset between the sample aperture and the sample position 731 // 732 // 733 // INPUT: 734 // inQ = q-value [1/A] 735 // lambda = wavelength [A] 736 // lambdaWidth = [dimensionless] 737 // DDet = detector pixel resolution [cm] **assumes square pixel 738 // apOff = sample aperture to sample distance [cm] 739 // S1 = source aperture diameter [mm] 740 // S2 = sample aperture diameter [mm] 741 // L1 = source to sample distance [m] 742 // L2 = sample to detector distance [m] 743 // BS = beam stop diameter [mm] 744 // del_r = step size [mm] = binWidth*(mm/pixel) 745 // usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam 746 // 747 // OUPUT: 748 // SigmaQ 749 // QBar 750 // fSubS 751 // 752 // 753 Function V_getResolution_old(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,SigmaQ,QBar,fSubS) 754 Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses 755 Variable &SigmaQ, &QBar, &fSubS //these are the output quantities at the input Q value 756 757 //lots of calculation variables 758 Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g 759 Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r 760 761 //Constants 762 Variable vz_1 = 3.956e5 //velocity [cm/s] of 1 A neutron 763 Variable g = 981.0 //gravity acceleration [cm/s^2] 764 765 766 S1 *= 0.5*0.1 //convert to radius and [cm] 767 S2 *= 0.5*0.1 768 769 L1 *= 100.0 // [cm] 770 L1 -= apOff //correct the distance 771 772 L2 *= 100.0 773 L2 += apOff 774 del_r *= 0.1 //width of annulus, convert mm to [cm] 775 776 BS *= 0.5*0.1 //nominal BS diameter passed in, convert to radius and [cm] 777 // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture 778 Variable LB 779 LB = 20.1 + 1.61*BS //distance in cm from beamstop to anode plane (empirical) 780 BS = bs + bs*lb/(l2-lb) //adjusted diameter of shadow from parallax 781 782 //Start resolution calculation 783 a2 = S1*L2/L1 + S2*(L1+L2)/L1 784 q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2) 785 lp = 1.0/( 1.0/L1 + 1.0/L2) 786 787 v_lambda = lambdaWidth^2/6.0 788 789 // if(usingLenses==1) //SRK 2007 790 if(usingLenses != 0) //SRK 2008 allows for the possibility of different numbers of lenses in header 791 v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2 //correction to 2nd term 792 else 793 v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2 //original form 794 endif 795 796 v_d = (DDet/2.3548)^2 + del_r^2/12.0 //the 2.3548 is a conversion from FWHM->Gauss, see http://mathworld.wolfram.com/GaussianFunction.html 797 vz = vz_1 / lambda 798 yg = 0.5*g*L2*(L1+L2)/vz^2 799 v_g = 2.0*(2.0*yg^2*v_lambda) //factor of 2 correction, B. Hammouda, 2007 800 801 r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) )) 802 delta = 0.5*(BS - r0)^2/v_d 803 804 if (r0 < BS) 805 inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta)) 806 else 807 inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta)) 808 endif 809 810 fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) ) 811 if (fSubS <= 0.0) 812 fSubS = 1.e-10 813 endif 814 fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi)) 815 fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d 816 817 rmd = fr*r0 818 v_r1 = v_b + fv*v_d +v_g 819 820 rm = rmd + 0.5*v_r1/rmd 821 v_r = v_r1 - 0.5*(v_r1/rmd)^2 822 if (v_r < 0.0) 823 v_r = 0.0 824 endif 825 QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2)) 826 SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda) 827 828 829 // more readable method for calculating the variance in Q 830 // EXCEPT - this is calculated for Qo, NOT qBar 831 // (otherwise, they are nearly equivalent, except for close to the beam stop) 832 // Variable kap,a_val,a_val_l2,m_h 833 // g = 981.0 //gravity acceleration [cm/s^2] 834 // m_h = 252.8 // m/h [=] s/cm^2 835 // 836 // kap = 2*pi/lambda 837 // a_val = L2*(L1+L2)*g/2*(m_h)^2 838 // a_val_L2 = a_val/L2*1e-16 //convert 1/cm^2 to 1/A^2 839 // 840 // sigmaQ = 0 841 // sigmaQ = 3*(S1/L1)^2 842 // 843 // if(usingLenses != 0) 844 // sigmaQ += 2*(S2/lp)^2*(lambdaWidth)^2 //2nd term w/ lenses 845 // else 846 // sigmaQ += 2*(S2/lp)^2 //2nd term w/ no lenses 847 // endif 848 // 849 // sigmaQ += (del_r/L2)^2 850 // sigmaQ += 2*(r0/L2)^2*(lambdaWidth)^2 851 // sigmaQ += 4*(a_val_l2)^2*lambda^4*(lambdaWidth)^2 852 // 853 // sigmaQ *= kap^2/12 854 // sigmaQ = sqrt(sigmaQ) 855 856 857 Return (0) 858 End 859 860 711 712 713 -
sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Write_VSANS_QIS.ipf
r1119 r1120 518 518 // qval = CalcQval(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 519 519 Duplicate/O qTot,phi,r_dist 520 Variable pixSizeX,pixSizeY 520 Variable pixSizeX,pixSizeY,xctr,yctr 521 521 pixSizeX = V_getDet_x_pixel_size(type,detStr) 522 522 pixSizeY = V_getDet_y_pixel_size(type,detStr) 523 523 524 phi = V_FindPhi( pixSizeX*((p+1)-bCentX) , pixSizeY*((q+1)-bCentY)+(2)*yg_d) //(dx,dy+yg_d) 525 r_dist = sqrt( (pixSizeX*((p+1)-bCentX))^2 + (pixSizeY*((q+1)-bCentY)+(2)*yg_d)^2 ) //radial distance from ctr to pt 524 xctr = V_getDet_beam_center_x_pix(type,detStr) 525 yctr = V_getDet_beam_center_y_pix(type,detStr) 526 phi = V_FindPhi( pixSizeX*((p+1)-xctr) , pixSizeY*((q+1)-yctr)+(2)*yg_d) //(dx,dy+yg_d) 527 r_dist = sqrt( (pixSizeX*((p+1)-xctr))^2 + (pixSizeY*((q+1)-yctr)+(2)*yg_d)^2 ) //radial distance from ctr to pt 526 528 527 529 //make everything in 1D now … … 591 593 KillWaves/Z qx_val_s,qy_val_s,z_val_s,qz_val_s,SigmaQx_s,SigmaQy_s,fSubS_s,sw,sw_s 592 594 593 Killwaves/Z q x_val,qy_val,z_val,qval,qz_val,sigmaQx,SigmaQy,fSubS595 Killwaves/Z qval,sigmaQx,SigmaQy,fSubS,phi,r_dist 594 596 595 597 Print "QxQy_Export File written: ", V_GetFileNameFromPathNoSemi(detSavePath)
Note: See TracChangeset
for help on using the changeset viewer.