subroutine fit_abskor(muRinp, rinp) real muRinp, rinp call fit_abskor2(muRinp, rinp, 0.0) end subroutine fit_abskor2(muRinp, rinp, radc) c absorption correction for neutron powder diffraction c cylindrical geometry c implicit none include 'fit.inc' real muRinp, rinp, radc real muR,ri,ra character*64 value real cvt_ratio if (muRinp .eq. 0) then write(isyswr,'(x,a)') 'Absorption correction' write(isyswr,'(x,2a)') '-- for hollow cyclinders, the muR' 1 ,' of a full cylinder has to be given' muR=0 call sym_get_real('muR', muR) if (muR .eq. 0) then write(isyswr,'(x,a,$)') 'muR: ' else write(isyswr,'(x,a,f8.4,a,$)') 'muR (default:', muR, '): ' endif read(isysrd,'(F20.0)',end=999,err=999) muR write(isyswr,'(x,a,$)') 'r/R (default: 0.0 / full cyclinder): ' read(isysrd,'(A)',end=999,err=999) value ri=cvt_ratio(value, -1.0) if (ri .lt. 0) goto 999 write(isyswr,'(/3(x,a/),x,a,$)') 'Radial collimator correction' 1 ,'on DMC, in general, no correction is needed' 1 ,'on HRPT, coll. is 7 or 14 mm' 1 ,'sample radius/collimation fwhm (default: 0 = no corr.): ' read(isysrd,'(A)',end=999,err=999) value ra=cvt_ratio(value, -1.0) if (ra .lt. 0.0) goto 999 else muR=muRinp ri=rinp ra=radc endif call sym_put_real('muR', abs(muR)) if (muR .gt. 0.0) then write(isyswr,'(x,a,f8.4)') 'Correct for absorption, muR=',muR call abscorn(muR, ri, ra, 2, nxmax+1-nxmin 1 , xval(nxmin), yval(nxmin), sig(nxmin)) elseif (muR .lt. 0.0) then write(isyswr,'(x,a,f8.4)') 1 'Inverse absorption correction, muR=',-muR call abscorn(-muR, ri, ra, -2, nxmax+1-nxmin 1 , xval(nxmin), yval(nxmin), sig(nxmin)) endif if (ri .eq. 0) then write(isyswr,'(x,a)') 'for a full cyclinder' else write(isyswr,'(x,a,f8.4)') 'for a hollow cyclinder, r/R=',ri endif if (ra .ne. 0) then write(isyswr,'(x,a,f5.1,a)') 'correction for radial collimator ' 1 ,ra,' ra/rc' endif return 999 print *,'input error' end subroutine abscorn(muR, ri, ra, mode, npt, twoth, yy, sig) ! calculate absorption correction of hollow cylinders real muR ! absorption coefficient real ri ! (inner radius)/(outer radius) real ra ! radial collimator 1/fwhm integer mode ! mode = 0 transmission is returned in YY ! mode = 1 correction of YY ! mode =-1 inverse correction of YY ! mode = 2 correction of YY and SIG ! mode =-2 inverse correction of YY and SIG integer npt ! number of data points real twoth(npt) ! two_theta values real yy(npt) ! returned transmission or Y to correct real sig(*) ! sigma to correct integer nmax parameter (nmax=60) double precision th1, th2, dt, tc0, tcn real t, tnorm real tcalc(0:nmax) real slope(0:nmax), x integer i, j, n external spline, abscor1 real spline double precision abscor1 if (abs(mode) .gt. 3) stop 'illegal mode' if (mode .eq. 0) then do i=1,npt yy(i)=1.0 enddo else tnorm=abscor1(dble(muR), dble(ri), dble(ra), 0.0D0) endif if (npt .le. nmax) then do i=1,npt t=abscor1(dble(muR), dble(ri), dble(ra), dble(twoth(i)*0.5)) if (mode .gt. 0) t=1/t yy(i)=yy(i)*t if (abs(mode) .eq. 2) then sig(i)=sig(i)*t endif enddo else th1=90. th2=0. do i=1,npt th1=min(th1, max(0.0D0,dble(twoth(i)*0.5-90/nmax))) th2=max(th2, min(90.D0,dble(twoth(i)*0.5+90/nmax))) enddo n=nint(nmax*min(1.D0,(th2-th1)/90)) dt=(th2-th1)/n tc0=abscor1(dble(muR), dble(ri), dble(ra), th1) tcalc(0)=tc0 tcn=abscor1(dble(muR), dble(ri), dble(ra), th2) tcalc(n)=tcn do i=1,n-1 tcalc(i)=abscor1(dble(muR), dble(ri), dble(ra), th1+dble(i*dt)) enddo slope(0)=-10.*(tc0-abscor1(dble(muR),dble(ri),dble(ra),th1+dt*0.1)) slope(n)= 10.*(tcn-abscor1(dble(muR),dble(ri),dble(ra),th2-dt*0.1)) do i=1,n-1 slope(i)=(tcalc(i+1)-tcalc(i-1))*0.5 enddo do i=1,npt x=(twoth(i)*0.5-th1)/dt j=min(n-1,int(x)) t=spline(x-j, tcalc(j), slope(j)) if (mode .gt. 0) t=1/t yy(i)=yy(i)*t if (abs(mode) .eq. 2) then sig(i)=sig(i)*t endif enddo endif if (mode .ne. 0) then if (mode .lt. 0) tnorm=1/tnorm do i=1,npt yy(i)=yy(i)*tnorm if (abs(mode) .eq. 2) then sig(i)=sig(i)*tnorm endif enddo endif end real function spline(x, y, s) real x, y(2), s(2) real dy dy=y(2)-y(1) spline=y(1)+x*((x-1)*(s(2)*x+s(1)*(x-1))+dy*x*(3-2*x)) end double precision function abscor1(muR, ri, ra, th) double precision th, mur, ri, ra double precision sum, sums, p, w, r, dphi double precision dr, phi, ri2, pa integer nr, nf, j, k, ip double precision path double precision cosd,x cosd(x)=cos(x/180.*3.14159265) nr=100 nf=200 ri2=ri*ri dr=(1-ri)/nr ! simpson integration over r sums=0 do j=0,nr if (mod(j,2) .eq. 1) then w=4./3. else if (j .ne. 0 .and. j .ne. nr) then w=2./3. else w=1./3. endif r=min(1.D0,ri+j*dr) ! simple integration over phi (simpson would not help, as the integration goes over a circle) sum=0 dphi=360./nf do k=1,nf phi=(k-0.5)*dphi pa=path(phi, th, r, ri2) ip=nint(pa*10) p=mur*pa if (ra .eq. 0) then sum=sum+exp(-p) else sum=sum+exp(-p)*max(0D0,1D0-r*abs(cosd(phi+th))*ra) endif enddo sums=sums+w*sum*r*dphi enddo abscor1=sums/(1.+ri)/dble(nr*180) end double precision function path(phi, th, r, ri2) ! ! calculate path length through a hollow cylinder sample of an ! outer radius ro=1 and an inner radius ri ! ! phi: angle (scatterer,sample center,y-direction) ! r: distance scatterer - sample center ! th: scattering angle theta (deg) ! ri2 = ri*ri ! the incoming beam has direction -theta, the outgoing beam direction +theta double precision phi, th, r, ri2 double precision r1, r2, p1, p2, ph double precision sind, cosd, x cosd(x)=cos(x/180.*3.14159265) sind(x)=sin(x/180.*3.14159265) ph=mod(phi,360.D0) if (ph .gt. 180.D0) ph=360.D0-ph ! problem is symmetric (y-axis is mirror axis) r1=(r*cosd(th-ph))**2 ! square of distance of incoming beam from center p1=sqrt(1.0D0-r1)-r*sind(th-ph) ! path from entry point if (r1 .lt. ri2 .and. ph .gt. th) then ! incoming beam cuts inner cylinder before scatterer p1=p1-2.D0*sqrt(ri2-r1) ! subtract path through inner cylinder endif r2=(r*cosd(th+ph))**2 ! square of distance of outgoing beam from center p2=sqrt(1-r2)-r*sind(th+ph) ! path to exit point if (r2 .lt. ri2 .and. ph .gt. 180.D0-th) then ! outgoing beam cuts inner cylinder after scatterer p2=p2-2.D0*sqrt(ri2-r2) ! subtract path through inner cylinder endif path=p1+p2 end