263 lines
6.7 KiB
Fortran
263 lines
6.7 KiB
Fortran
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
|