C	(C) Copr. 02/2011  Alexander Moroz
C	Wave-scattering.com
C	http://www.wave-scattering.com/
C	wavescattering@yahoo.com
C
C	Last update: February 2011
C
C	This is fsc2d, a computer code to generate data for 
C	the figures of my article
C
C	[1] A. Moroz, Electron mean-free path in metal coated nanowires,
C	              J. Opt. Soc. Am. B <b>28</b>(5), 1130-1138 (2011).
C         http://www.opticsinfobase.org/josab/abstract.cfm?URI=josab-28-5-1130
C =====================================================================
C
C	This code is made available on the following conditions:
C
C	- You can use and modify it, as long as acknowledgment is given
C	  to the author. 
C
C	- You can not sell it or use it in any way to get money out of it
C	  (save for the research grants you might get thanks to it)
C
C	- I'm making this code available to the best of my knowledge.
C	  It has been tested as throughly as possible. 
C            It work for Windows and Linux platforms.  However, it's
C	  your responsibility to check for its accuracy in the environment
C	  and size/shape ranges you use. 
C
C      DISCLAIMER: 
C
C      WHILE THE COMPUTER PROGRAMS HAVE BEEN TESTED FOR A VARIETY OF CASES,
C      IT IS NOT INCONCEIVABLE THAT THEY MAY CONTAIN UNDETECTED ERRORS. ALSO,
C      INPUT PARAMETERS CAN BE USED WHICH ARE OUTSIDE THE ENVELOPE OF
C      VALUES FOR WHICH RESULTS ARE COMPUTED ACCURATELY. FOR THIS REASON,
C      THE AUTHOR DISCLAIM ALL LIABILITY FOR ANY DAMAGES THAT MAY RESULT 
C      FROM THE USE OF THE PROGRAMS. 
C
C      I would highly appreciate informing me of any problems encountered 
C      with this code. Improvements, comments, additions, all welcome at 
C
C                       wavescattering@yahoo.com
C
C =====================================================================

	program fsc2d
C--------/---------/---------/---------/---------/---------/---------/--
C
C  THIS PROGRAM CALCULATES the effective mean free paths in a cylindrical
C  shell geometry and their asymptotics for various models of surface 
C  scattering
C
C  Output: 
C
C  efflngth2d.dat: r1 vs effective length
C  In columns: 
C              r1
C              xlbd  - billiard model 
C              xld   - diffusive model
C
C  efflasmpt2d.dat: r1 vs effective length
C  In columns: 
C              r1
C              xld    - diffusive model exact
C              xldas  - diffusive model leading asymptotic
C              xldas  - diffusive model 2-term asymptotic
C
C The code employs a combination  of the integration routines 
C
C    qromo.f 
C    polint.f
C    midpnt.f
C
C from Numerical Recipes.
C--------/---------/---------/---------/---------/---------/---------/--
	implicit none
	integer nout
	parameter(nout=30)
	integer istep,nstep
	real*8 pi,r1,r2,r1in,rsnm,rff,xd,xdel,xld,xldas,xldasl,
     &       xlbd,xl1,xl2,xup
	real*8 xufnct1,xufnct2

	external xufnct1,xufnct2,midpnt
	common/ttt/ r1,r2
	DATA PI/3.141592653589793d0/

	nstep=500
	r1in=1.d0   !0.d0
	r2=100.d0
	xdel=(r2-r1in)/dble(nstep)
	
	OPEN(UNIT=NOUT,FILE='efflngth2d.dat')
	rewind(NOUT)
	WRITE(NOUT,*)'rff vs effective length'
	write(nout,*)'In columns: rff,xlbd,xld'
	write(nout,*)


	OPEN(UNIT=NOUT+2,FILE='efflasmpt2d.dat')
	rewind(NOUT+2)
	WRITE(NOUT+2,*)'Thin-shell asymptotic diffusive model'
	WRITE(NOUT+2,*)'rff vs effective length'
	write(nout+2,*)'In columns: xd,xld,xldasl,xldas'
	write(nout+2,*)
*
c	write(6,*) 'log(10)=', log(10.d0)
c	pause

	do 10 istep=1,nstep-1

	r1=r1in + dble(istep)*xdel  !0.001d0

	if (r1.gt.r2) goto 10

	write(6,*) 'r1,r2=', r1,r2

	rff=r1/r2
	xd=r2-r1	

c Chernov (see his Eq. (2.6))

	xlbd=pi*xd/2.d0	
	write(6,*)'xlbd=',xlbd


c Diffusive

c integrate from 0 up to pi/2	

      call qromo(xufnct2,0.d0,pi/2.d0,xl2,midpnt)

c integrate from 0 up to arccos(r1/r2)

	xup=acos(rff)
	call qromo(xufnct1,0.d0,xup,xl1,midpnt)
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 Calls midpnt.f and polint.f
C--------/---------/---------/---------/---------/---------/---------/--

	xld=(r1*xl1+r2*xl2)/(r1+r2)
	xld=2.d0*xld/pi

	write(6,*)'xld=',xld
	write(nout,9998) r1,xlbd,xld

	if (r2.ge.r1in) then
	xldasl=(xd/pi)*log(r2/(2.d0*xd)+1.d0/4.d0) 
      xldas =xldasl + (2.d0+4.d0*log(2.d0))*xd/pi
	write(6,*)'xldas=',xldas
	write(nout+2,9997) xd,xld,xldasl,xldas
	end if

c	pause

 10	enddo

	close(nout)
	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 xufnct1(xtht)   !DOUBLEPRECISION FUNCTION 
	real*8 r1,r2,bt,x,xtht

	common/ttt/ r1,r2

	bt=(r1**2+r2**2)/(2.d0*r1*r2)
	x=cos(xtht)
	xufnct1=sqrt(r1*r2/2.d0)*(2.d0*bt*x-x**2-1.d0)/
     1 ( (x-(r1/r2)) * sqrt(bt-x) )
      return
	end
*
	doubleprecision function xufnct2(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

	xufnct2=x

      return
	end