program fsc
C--------/---------/---------/---------/---------/---------/---------/--
C
C  THIS PROGRAM CALCULATES the effective mean free paths in a spherical
C  shell geometry and their asymptotics for various models of surface 
C  scattering
C
C  Output: 
C
C  efflength.dat: r1 vs effective length
C  In columns: 
C              r1
C              xlavr  - shell thickness
C              xlqh   - Granqvist and Hunderi 
C              xlxu   - Xu 
C              xlkp   - Kachan and Ponyavina 
C              xlchrn - billiard model 
C
C  efflp.dat: r1 vs effective length
C  In columns: 
C              r1
C              xld    - diffusive model
C              xli    - isotropic model
C              xlchrn - billiard model
C
C  efflasmpt.dat: r1 vs effective length
C  In columns: 
C              r1
C              xld    - diffusive model
C              xli    - isotropic model
C              xlchrn - billiard model
C              2*xd
C
C Provided that you are in possession of the integration 
C routine qrombt.f from Numerical Recipes, you can comment
C out the lines beginning with "cq" and view how numerical integration
C of xlxu compares against analytical result.
C--------/---------/---------/---------/---------/---------/---------/--
	implicit none
	integer nout
	parameter(nout=30)
	integer istep,nstep
	real*8 pi,r1,r2,rsnm,rff,xdel,xlavr,xlkp,xlqh,xlxu,xld,xli,
     &       xlchrn,xd,xn,xdn,xx,xx1,xl1,xl2,xl1d,xl2d,
     &       xlxuas,xldas,xlias,xlchas,xl1as,xl2as
	real*8 xufnct

	external xufnct
	common/ttt/ r1,r2
	DATA PI/3.141592653589793d0/

	nstep=250
	r2=100.d0
	xdel=r2/dble(nstep)
	
	OPEN(UNIT=NOUT,FILE='efflength.dat')
	rewind(NOUT)
	WRITE(NOUT,*)'rff vs effective length'
	write(nout,*)'In columns: rff,xlavr,xlqh,xlxu,xlkp,xlchrn'
	write(nout,*)

	OPEN(UNIT=NOUT+1,FILE='efflp.dat')
	rewind(NOUT+1)
	WRITE(NOUT+1,*)'rff vs effective length'
	write(nout+1,*)'In columns: rff,xld,xli,xlchrn'
	write(nout+1,*)

	OPEN(UNIT=NOUT+2,FILE='efflasmpt.dat')
	rewind(NOUT+2)
	WRITE(NOUT+2,*)'Thin-shell asymptotic for basic scattering models'
	WRITE(NOUT+2,*)'rff vs effective length'
	write(nout+2,*)'In columns: rff,xldas,xlias,xlchas'
	write(nout+2,*)


*
c	write(6,*) 'log(10)=', log(10.d0)
c	pause

	do 10 istep=0,nstep

	r1=0.001d0+ dble(istep)*xdel

	if (r1.gt.r2) goto 10

	write(6,*) 'r1,r2=', r1,r2
	rff=r1/r2
	
c Averitt
	xd=r2-r1

	write(6,*) 'xlavr=',xd

c Granqvist&Hunderi

	xlqh=(xd*(r2**2-r1**2))**(1.d0/3.d0)
	
	write(6,*)'xlqh=',xlqh	

c Xu	- integrate from 0 up to pi/2	

cq      call qrombt(xufnct,0.d0,pi/2.d0,xlxu)
C--------/---------/---------/---------/---------/---------/---------/--
C Returns the integral xlxu of the function xufnct from 0.d0 to pi/2.d0. 
C Integration is performed by Romberg's method of order 2K, where, 
C e.g.,K=2 is Simpson's rule.
C Internal QROMB Parameters: 
C EPS  ...  the fractional accuracy desired, as determined 
C           by the extrapolation error estimate; 
C JMAX ...  limits the total number of steps
C K    ...  the number of points used in the extrapolation.
C
C              Includes trapzd.f and polint.f
C--------/---------/---------/---------/---------/---------/---------/--

cq	write(6,*)'Numerical xlxu=',xlxu

	xlxu=xd*(2.d0*r2+r1)/(2.d0*r2)-
     &  ((r2**2-r1**2)/(4.d0*r2))*log(xd/xx1)

	write(6,*)'Analytical value of xlxu=',xlxu

	xlxuas=3*xd/2.d0+(xd/2.d0)*log(2.d0*r2/xd-1.d0)

	write(6,*)'xlxuas=',xlxuas

c Kachan&Ponyavina:
	if(rff.ne.1.d0) then
	xlkp=(1.d0-rff**2)*(1.d0-rff)*log((1.d0-rff)/(1.d0+rff))/
     & (4.d0*(1.d0+rff**2))
	xlkp=1.0/(1.d0+rff**2)-rff/2.d0 -xlkp
	xlkp=r2*xlkp

	else

	xlkp=0.d0

	end if

	write(6,*)'xlkp=',xlkp

c Chernov (see his Eq. (2.6))

	xn=(r2**3-r1**3)/3.d0	
      xdn=r2**2+r1**2
	xlchrn=4.d0*xn/xdn	

	write(6,*)'xlchrn=',xlchrn

	xlchas=2.d0*xd-xd**3/(3.d0*r2**2)

	write(6,*)'xlchas=',xlchas

	write(nout,9999) r1,xd,xlqh,xlxu,xlkp,xlchrn

c Diffusive
	xx=r1**2+r2**2
	xx1=r1+r2
	xn=4.d0*r2**3-2.d0*r1*xx-xx1**2*xd*log(xd/xx1)
	xld=xn/(4.d0*xx)

	write(6,*)'xld=',xld

c	xldas=(xd/2.d0)*log(r2/xd)
c	xldas=xd*log(r2/xd)

	xldas=xd +(xd/2.d0)*log(2.d0*r2/xd-1.d0)

	write(6,*)'xldas=',xldas

c Isotropic

	xl2=5.d0*r2**2*(1.d0-(r1/r2)**2)**(3.d0/2.d0)/3.d0 
     &  -(2.d0*r1+r2)*xd**2/(3.d0*r2)

	xl2d=xd*(2.d0*r2+r1)/(2.d0*r2)-
     &  ((r2**2-r1**2)/(4.d0*r2))*log(xd/xx1)

	xl2=xl2/xl2d

	xl2as=20.d0*sqrt(2.d0*r2*xd)/(3.d0*(3.d0+log(-1.d0+2.d0*r2/xd)))

	write(6,*)'xl2=',xl2
	write(6,*)'xl2as=',xl2as 
c	write(6,*)'log(10.d0)=',log(10.d0) 
	 
	xl1=2.d0*(r2**2-r1**2)**(3.d0/2.d0)/(3.d0*r1)
     &  -(2.d0*r2+r1)*xd**2/(3.d0*r1)

	xl1d=xd/2.d0-((r2**2-r1**2)/(4.d0*r1))*log(xd/xx1)

	xl1=xl1/xl1d

	xl1as=8.d0*sqrt(2.d0*r2*xd)/(3.d0*(1.d0+log(-1.d0+2.d0*r2/xd)))

	write(6,*)'xl1=',xl1
	write(6,*)'xl1as=',xl1as 
	 
	xli=(r2**2*xl2+r1**2*xl1)/(r2**2+r1**2)

	write(6,*)'xli=',xli

	xlias=14.d0*sqrt(2.d0*r2*xd)/
     &           (3.d0*(3.d0+log(-1.d0+2.d0*r2/xd)))

	write(6,*)'xlias=',xlias

c Writing

	xn=xd/r2

	write(nout+1,9997) r1,xld,xli,xlchrn,2.*xd

	if (rff.ge.0.74) write(nout+2,9997) xd,xld,xli,xlchrn,2.*xd

c	pause

 10	enddo

	close(nout)
	close(nout+1)
	close(nout+2)

 9997 format(F8.4,1X,4F10.4)
 9998 format(F8.4,1X,3F10.4)
 9999 format(F8.4,1X,5F10.4)


	end

*
	doubleprecision function xufnct(xtht)   !DOUBLEPRECISION FUNCTION 
	real*8 r1,r2,x,xtht

	common/ttt/ r1,r2

	if (sin(xtht).le.(r1/r2)) then

	x=r2*cos(xtht)-sqrt(r1**2-r2**2*sin(xtht)**2)

	else

	x=2.d0*r2*cos(xtht)

	endif

	xufnct=x*sin(xtht)
      return
	end