Title                  :RECOIL.ABS
Keywords               :Sputtering,Defects,Displacement Damage,EELS,XEDS,HVEM
Computer               :DEC VAX 11/785
Operating System       :VAX/VMS
Programming Language   :VAX Fortran
Hardware Requirements  :None
Author(s)              :Charles R. Bradley
Correspondence Address :Argonne Nat. Lab, Electron Microscopy Center,  Bldg  212
                       :Materials Science Division, Argonne, Illinois 60439, USA
Abstract:
	This program calculates differential Mott cross-sections of elemental
nuclei in collisions with relativistic electrons. The cross-sections are
output as functions of the nuclear recoil angle and are proportional to the
probability that a nucleus will recoil in a given direction.  These
cross-sections are useful in the study of displacement damage, electron
transmission sputtering, electron beam mixing, and recoil implantation.  For
medium and high voltage transmission electron microscopes now available,
having clean vacuum systems and high brightness sources, modification of
sample composition due to sputtering and/or displacement by the probe will
become an important consideration during microanalysis or imaging.
Title                  :RECOIL.DOC
Keywords               :Sputtering,Defects,Displacement Damage,EELS,XEDS,HVEM
Computer               :DEC VAX 11/785
Operating System       :VAX/VMS
Programming Language   :VAX Fortran V4.0
Hardware Requirements  :None
Author(s)              :Charles R. Bradley
Correspondence Address :Argonne Nat. Lab, Electron Microscopy Center,  Bldg  212
                       :Materials Science Division, Argonne, Illinois 60439, USA
Description            :Outline/Description and the purpose of the code.

Version:
	1.0

Purpose:
	This program calculates differential Mott cross-sections of elemental
nuclei in collisions with relativistic electrons. The cross-sections are
output as functions of the nuclear recoil angle and are proportional to the
probability that a nucleus will recoil in a given direction.  These
cross-sections are useful in the study of displacement damage, electron
transmission sputtering, electron beam mixing, and recoil implantation.  For
medium and high voltage transmission electron microscopes now available,
having clean vacuum systems and high brightness sources, modification of
sample composition due to sputtering and/or displacement by the probe will
become an important consideration during microanalysis or imaging.


References: 
	McKinley & Feshbach
	Phys Rev 74,1759(1948)

	Feshbach
	Phys Rev 88,295(1952)

	Doggett & Spencer
	Phys Rev 103,1597(1956)

	O. S. Oen
	CROSS-SECTIONS FOR ATOMIC DISPLACEMENTS
	IN SOLIDS BY FAST ELECTRONS
	ORNL Report # TM-4897,(1973)

	T. J. I'A Bromwich
	AN INTRODUCTION TO THE THEORY OF
	INFINITE SERIES, 2nd ed.
	Macmillan 1955

	Francis Scheid
	NUMERICAL ANALYSIS
	Schaum's Outline Series
	McGraw Hill Book Co. 1968

	Abromowitz and Stegun
	HANDBOOK OF MATHEMATICAL FUNCTIONS
	Dover, New York 1972

	James W. Corbett
	ELECTRON RADIATION DAMAGE IN SEMICONDUCTORS AND METALS
	in: Solid State Physics Advances in Research and Applications
	Supplement 7, Ed. Frederick Seitz and David Turnbull
	Academic Press, New York 1966

	D. Cherns, F. J. Minter, and R. S. Nelson
	Nuclear Instruments and Methods 132,369(1976)

Compilation Procedure: 
	Using VAX-Fortran compiler
	1. Rename file to RECOIL.FOR
	2. At DCL prompt type FORTRAN/EXTEND/NOF77 RECOIL
	
Linking Procedure: 
	At DCL prompt type LINK RECOIL

Input/Output:
	UNIT 5	input	TT:
	UNIT 7	output	TT:
	UNIT 2	output	(to be assigned)

Test Data: 
	This routine is designed to read input from the keyboard

General Comments:
	Subroutines and function subprograms required are included in the
	listing.  Amoung those included are:
		D which returns coefficients in the legendre
		  expansions which express the cross-section.
		EULER which preforms an euler transformation of an
		  alternating series and returns the sum of the series.
		LGDRE which returns legendre polynomials.
		LNGAM which returns the complex value of the log of
		  a gamma function with complex arguments
		SIGMR which returns the ratio of the Mott cross-section
		  to the Rutherford Cross-section.

Sample run:

	prompt>assign test.out for002
	%DCL-I-SUPERSEDE, previous value of FOR002 has been superseded
	prompt>run recoil
         Atomic Number of target: 79
           Atomic Mass of target: 196.967
Surface Binding Energy of target: 3.820
       Accelerating Voltage (kV): 400.

	One or two minutes will be required before the prompt is
	displayed again. Output is directed to the file named TEST.OUT

Sample Output: 
	The following is output generated by the above sample input
	which corresponds to Au with the assumption that the surface
	binding energy is equal to the sublimation energy.

    Atomic number of target =        79.
 Atomic mass of target(amu) =    196.967
 Surface Binding Energy(eV) =      3.820
   Accelerating Voltage(kV) =    400.000
      Maximum  Recoil Angle =     38.300
    
 Recoil                      Mott
 Angle                    Cross-Section
                         
  
      1.000		     277.851
      2.000		     279.282
      3.000		     281.670
      4.000		     285.020
      5.000		     289.339
      6.000		     294.634
      7.000		     300.917
      8.000		     308.200
      9.000		     316.498
     10.000		     325.827
     11.000		     336.207
     12.000		     347.659
     13.000		     360.207
     14.000		     373.878
     15.000		     388.701
     16.000		     404.708
     17.000		     421.935
     18.000		     440.420
     19.000		     460.206
     20.000		     481.339
     21.000		     503.869
     22.000		     527.851
     23.000		     553.345
     24.000		     580.414
     25.000		     609.130
     26.000		     639.569
     27.000		     671.815
     28.000		     705.958
     29.000		     742.098
     30.000		     780.343
     31.000		     820.811
     32.000		     863.631
     33.000		     908.944
     34.000		     956.905
     35.000		    1007.686
     36.000		    1061.472
     37.000		    1118.471
     38.000		    1178.909
     38.300		    1197.719
Title                  :RECOIL.SRC
Keywords               :Sputtering,Defects,Displacement Damage,EELS,XEDS,HVEM
Computer               :DEC VAX 11/785
Operating System       :VAX/VMS
Programming Language   :VAX Fortran
Hardware Requirements  :None
Author(s)              :Charles R. Bradley
Correspondence Address :Argonne Nat. Lab, Electron Microscopy Center,  Bldg  212
                       :Materials Science Division, Argonne, Illinois 60439, USA
Code:
*     ..................................................................
*
*        PROGRAM recoil
*
*	Version		Date		Initials	Reason
*-----------------------------------------------------------------------
*	 1.0		05/10/88	CRB		Creation by modification
*							of dMOTTcs version 1.0
*
*        Purpose
*           Compute the differential Mott cross-section
*	    as a function of the nuclear recoil angle
*
*        Subroutines and function subprograms required
*	    are included in this listing
*
*	I/O UNITS used;
*	 UNIT=2	(to be assigned)	output
*	 UNIT=5	TT:			input
*	 UNIT=7	TT:			output
*
*	 REFERENCES:
*	    McKinley & Feshbach
*	    Phys Rev 74,1759(1948)
*
*	    Feshbach
*	    Phys Rev 88,295(1952)
*
*	    Doggett & Spencer
*	    Phys Rev 103,1597(1956)
*
*	    O. S. Oen
*	    CROSS-SECTIONS FOR ATOMIC DISPLACEMENTS
*	    IN SOLIDS BY FAST ELECTRONS
*	    ORNL Report # TM-4897,(1973)
*
*	    G. H. Kinchin & R. S. Pease
*	    Repts. Progr. Phys. 18,1 (1955)
*
*	    T. J. I'A Bromwich
*	    AN INTRODUCTION TO THE THEORY OF
*	    INFINITE SERIES, 2nd ed.
*	    Macmillan 1955
*
*	    Francis Scheid
*	    NUMERICAL ANALYSIS
*	    Schaum's Outline Series
*	    McGraw Hill Book Co. 1968
*
*           ABROMOWITZ AND STEGUN
*           HANDBOOK OF MATHEMATICAL FUNCTIONS
*           DOVER, NEW YORK 1972
*
*     ..................................................................

	program recoil

	COMMON /CONST/elmass,amu,evpKeV,fsc,CLrad,hbar,c,pi,one
	COMMON /REPEAT/AVold

	integer*2 ke	! utility counter

	real*8 zero
	real*8 xulog
	real*8 rth	! Rutherford C-S
	real*8 Mottcs	! Mott Cross-Section
	real*8 CS	! cross-section in barns (E-24 cm**2)
	real*16 Z	! atomic number of target
	real*16 AV	! accelerating voltage
	real*16 Ts	! sputtering energy
	real*16 Th	! threshold incident KE that delives Ts
	real*16 A	! Atomic mass of sputter target
	real*16 Em	! Maximum energy transfered
	real*16 beta	! Relativistic beta
	real*16 alpha	! Z*fine structure constant
	real*16 theta	! scattering angle
	real*16 psi	! nuclear recoil angle
	real*16 psilim	! Max nuclear recoil angle
	real*16 ratio	! ratio of cross-sections
	real*16 jacobi	! jacobian of transformation from theta to psi
	real*16 util
*	declaration for common block REPEAT
	real*16 AVold	! repeated calculation at the same energy?

*	declarations for common block const - constants
	real*16 elmass	! rest mass of electron (KeV)
	real*16 amu	! atomic mass unit (KeV)
	real*16 evpKeV	! ev/KeV conversion
	real*16 fsc	! fine structure constant
	real*16 CLrad	! classical radius of the electron
			! in E-12 cm -- e**2/mc**2
	real*16 hbar	! h bar quantum mechanical const.(eV s)
	real*16 c	! speed of light(cm/s)
	real*16 pi	! constant
	real*16 one	! constant

*-------INITIALIZE CONSTANTS IN COMMON------------------------
	AVold = -1.0		! new calc

	elmass = 510.976	! rest mass of electron (KeV)
	amu = 931141.		! atomic mass unit (KeV)
	evpKeV = 10**3		! ev/KeV conversion
	fsc = 1./137.0370	! fine structure constant
	CLrad = 0.281785	! classical radius of the electron
				! in E-12 cm -- e**2/mc**2
	hbar = 6.5817E-16	! h bar (eV s)
	c = 2.997930E10		! speed of light (cm/s)
	one = 1.0
	PI = 4. * qatan(one)	! calculate pi for latter use

*----------------INPUT Z,A,Ed ---- from DEVICE 1 ------
*	Z  = atomic number
*	A = atomic mass
*	Ts = sputtering or displacement energy (ev)

	write(7,10)
10	format('$','         Atomic Number of target: ')
	read(5,FMT=*) Z		! input atomic number
	write(7,11)
11	format('$','           Atomic Mass of target: ')
	read(5,FMT=*) A		! input atomic mass
	write(7,12)
12	format('$','Surface Binding Energy of target: ')
	read(5,FMT=*) Ts		! input atomic number
	write(7,20)
20	format('$','       Accelerating Voltage (kV): ')
	read(5,FMT=*) AV	! input accel. voltage

*----------------------initial calc ------------

	alpha = Z * fsc
	util = Ts*A*amu/(evpkev*2.*AV*(AV+2.*elmass))
	if(util.lt.1.)then
	  psilim = qacosd(sqrt(Ts*A*amu/(evpkev*2.*AV*(AV+2.*elmass))))
	  else
	  endif

*---------OUTPUT TO DEVICE 2 ------------------------
	write(2,40) Z
40	format(' ','   Atomic number of target = ',F10.0)
	write(2,41) A
41	format(' ','Atomic mass of target(amu) = ',F10.3)
	write(2,42) Ts
42	format(' ','Surface Binding Energy(eV) = ',F10.3)
	write(2,50) AV
50	format(' ','  Accelerating Voltage(kV) = ',F10.3)
	write(2,500) psilim
500	format(' ','     Maximum  Recoil Angle = ',F10.3)
	if(util.gt.1.)then
	 Th = -1.*elmass + qsqrt(elmass**2 + A*amu*Ts/(2.*evpKeV))
	 write(2,510) Th
510	 format(' ','Accelerating voltage is below the threshold of ',F10.3,'kV')
	 goto 3000
	 else
	 endif

	write(2,55)
55	format(' ','   ')
	write(2,60)
60	format(' ','Recoil                      Mott')
	write(2,65)
65	format(' ','Angle                    Cross-Section')
	write(2,67)
67	format(' ','                        ')

*-----------------------------------------


	write(2,80)
80	format(' ',' ')

	psi = 0.000
	psiinc = 1.000	! increment in psi used for calculation of table values

	do 2000 while(psi.lt.psilim)
	psi = psi + psiinc
	if(psi.gt.psilim) then
	  psi = psilim
	  else
	  endif

	theta = 2.0*qasind(qcosd(psi))
	beta = sqrt(1. - 1./(1.+ AV/elmass)**2)
c	Rth = (Z*Clrad)**2*(1.-beta**2)/(4.*(beta*qsind(theta/2.))**4)	! Rutherford CS
	Rth = (Z*Clrad)**2*(1.-beta**2)/(4.*(beta*qcosd(psi))**4)	! Rutherford CS

	call sigmr(Z,AV,theta,ratio)
	psi = qacosd(qsind(theta/2.))	! update psi since sigmr may modify
	jacobi = 4.0*qcosd(psi)
	MottCS = Rth*ratio*jacobi

	write(2,95) psi,mottcs
95	format(' ',F10.3,'		',F12.3)

2000	continue

3000	continue
	call exit
	end
*
*     ..................................................................
*
*        function D (k,Z,AV)
*
*        PURPOSE
*	    Calculate precursors of the coefficients
*	    of the Mott cross-section
*
*        USAGE
*           x = D(k,Z,AV)
*
*        DESCRIPTION OF PARAMETERS
*	    k	- order of the coefficient
*	    Z	- Atomic number of sputter target
*	    AV	- KE of incident electrons
*
*        REMARKS
*	    real*16
*
*        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
*	    LNGAM - natural logarithm of complex gamma fcn
*
*        METHOD
*	    see -
*	    McKinley & Feshbach
*	    Phys Rev 74,1759(1948)
*     ..................................................................
*
	function D(k,Z,AV)

	COMMON /CONST/elmass,amu,evpKeV,fsc,CLrad,hbar,c,pi,one

	integer*2 k	! order of the coefficient

	real*16 Z	! Atomic number of sputter target
	real*16 q	! parameter
	real*16 AV	! KE of incident electrons
	real*16 alpha	! fine structure constant * Z

*	DECLARATIONS FOR COMMON BLOCK CONST - CONSTANTS
	real*16 elmass	! rest mass of electron (KeV)
	real*16 amu	! atomic mass unit (KeV)
	real*16 evpKeV	! ev/KeV conversion
	real*16 fsc	! fine structure constant
	real*16 CLrad	! classical radius of the electron
			! in E-12 cm -- e**2/mc**2
	real*16 hbar	! h bar quantum mechanical const.(eV s)
	real*16 c	! speed of light(cm/s)
	real*16 pi	! constant
	real*16 one	! constant

	real*16 beta	! relativistic beta
	real*16 zero	! constant

	complex*16 D	! coefficient
	complex*16 rho	! parameter
	complex*16 ratioA,ratioB	! ratio of gamma functions in coef
	complex*16 zz		! utility variable
	complex*16 Im		! utility variable
	complex*16 lngam	! logarithm of gamma function

	EXTERNAL lngam		! ln(gamma) function routine

*	Initialize values
	alpha = Z * fsc
	zero = 0.0
	Im = dcmplx(zero,one)	! imaginary unit
	beta = sqrt(1. - 1./(1.+ AV/elmass)**2)

	q = alpha / beta
	zz = qfloat(k)**2 - alpha**2 + zero*Im
	rho = cdsqrt(zz)

	zz = qfloat(k) - Im*q
	ratioA = lngam(zz) - lngam(dconjg(zz))
	ratioA = cdexp(ratioA)
	ratioA = ratioA/dconjg(zz)
	ratioA = ratioA * cdexp(-1.*pi*Im*k)

	zz = rho - Im*q
	ratioB = lngam(zz)
	zz = rho + Im*q
	ratioB = ratioB - lngam(zz)
	ratioB = cdexp(ratioB)
	ratioB = ratioB/zz
	ratioB = ratioB * cdexp (-1.*pi*rho*Im)

	D = ratioA - ratioB

	return
	end

*        SUBROUTINE EULER
*
*        PURPOSE
*	    Accelerate convergence of oscillating series y(i)
*	    and return the sum of that series sumy
*
*        USAGE
*           CALL EULER (Y,NTERMS,SUMY)
*
*        DESCRIPTION OF PARAMETERS
*	    y	-	  input array - series to sum
*		    	   sequence to sum, sequence has form:
*			   y = y(1) - y(2) + y(3) - y(4) + ...
*	    nterms	-  number of terms in y
*	    sumy	-  result
*
*        REMARKS
*
*        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
*
*        METHOD
*	see T. J. I'A Bromwich
*	    AN INTRODUCTION TO THE THEORY OF
*	    INFINITE SERIES, 2nd ed.
*	    Macmillan 1955
*	    pg 63
*	and Francis Scheid
*	    NUMERICAL ANALYSIS
*	    Schaum's Outline Series
*	    McGraw Hill Book Co. 1968
*	    pg 24,39,75,160
*     ..................................................................
*
	subroutine euler(y,nterms,sumy)

	parameter (mxterm = 120)	! max number of terms in y(i)

	integer*2 maxdif		! max number of finite differences
	integer*2 nterms		! number of terms to use in sums
	integer*2 istart		! start calculating differences at
					! this term in the series

	real*16 cor			! correction to estimate of sumy
	real*16 y(mxterm)		! input array - series to sum
	real*16 sumy			! estimated sum
	real*16 B(mxterm,0:mxterm)	! binomial coefficient
	real*16 sigma			! utility sum
	real*16 delta(0:mxterm)		! finite difference of order ()

*	Initialize constants

	istart = 10			! start at y(10)

	maxdif = nterms - istart -1	! max number of differences
					! that can be generated
	if(maxdif.lt.2) then		! quit
	  type*,'too few points SUBROUTINE EULER'
	  stop
	 else
	 endif

*	sum the first istart terms

	sumy = 0.
	do 70 i=1,istart
	sumy = sumy + y(i)*(-1.)**(i-1)
70	continue

*	Generate Binomial Coefficients

*	               ( k )        k!
*	   B(k,n)  =   (   )  =  ---------
*                      ( n )      n!(k-n)!

*	using the recursion relation:

*	   ( k+1 )     (  k  )     ( k )
*	   (     )  =  (     )  +  (   )
*	   ( n+1 )     ( n+1 )     ( n )

	B(1,1) = 1.
	do 100 k=1,maxdif
100	B(k,0) = 1.

	do 200 k=2,maxdif
	do 200 n=1,k
	B(k,n) = B(k-1,n) + B(k-1,n-1)
200	continue

*	Generate finite differences of the series to be summed
*	using the relations:

*	   delta = 1-E  , delta**2 = (1-E)**2
*	   where Ey(n) = y(n+1)

*	see Scheid pg 39 & 160

	do 310 k=1,maxdif
	do 300 i=0,k
	 if(i.eq.0) then
	  sigma = 0.
	  else
	  endif

	 sigma = sigma + (-1)**(i)*B(k,i)*y(k-i+istart+1)

300	continue
310	delta(k) = sigma

*	Now sum Euler Transformed Series

	delta(0) = y(istart+1)	! first term in series
	do 400 i=0,maxdif
	cor = (-1)**(i+istart)*delta(i)/(2.**(i+1))
400	 sumy = sumy + cor

	return
	end

*
*     ..................................................................
*
*        SUBROUTINE LGDRE
*
*        PURPOSE
*	    Calculate legendre polynomials
*
*        USAGE
*           CALL lgdre (k,x,P)
*
*        DESCRIPTION OF PARAMETERS
*	    K	   - MAX ORDER OF POLYNOMIAL TO BE GENERATED
*	    X	   - ARGUMENT
*	    P	   - ARRAY OF LEGENDRE POLYNOMIALS RETURNED
*
*        REMARKS
*	    real*16
*
*        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
*
*        METHOD
*	    Use the recursion relation of
*	    Handbook of Mathematical Functions
*	    Dover, New York 1972
*	    Abromowitz & Stegun -
*	    8.5.3 with mu=0
*	    (n+1)P sub n+1 (x) = (2n+1)xP sub n (x) -nP sub n-1 (x)
*
*     ..................................................................
*
	subroutine lgdre(k,x,P)

	parameter(maxord = 120)

	real*16	P(0:maxord)	! values of legendre functions
				! up to k th order
	real*16 x		! argument of the polynomial
	integer*2 k		! order of the polynomial

*	check ranges

	if(k.lt.0.or.k.gt.maxord) then
	  type*,'Legendre Polynomial order must be >= 0'
	  type*,'and <=',maxord
	  else if(k.eq.0) then
	    P(0) = 1.0
	  else if(k.eq.1) then
	    P(1) = x
	  else
	    P(0) = 1.0
	    P(1) = x
	    do 100 m=1,k-1
	     P(m+1) = (2*m+1)*x*P(m) - m*P(m-1)
100	     P(m+1) = P(m+1)/(m+1)
	  endif

	return
	end
*
*     ..................................................................
*
*        function lngam(q)
*
*        PURPOSE
*	    Calculate the natural logarithm
*	    of the complex gamma function
*
*        USAGE
*           x = lngam(q)
*
*        DESCRIPTION OF PARAMETERS
*	    q	   - complex*16 input value
*
*        REMARKS
*	    complex*16 result
*
*        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
*
*        METHOD
*	    Use the asymptotic approximation of
*	    Abromowitz & Stegun -
*	    6.1.41
*	    STERLING'S APPROXIMATION
*	    and
*	    the recursion relation
*	    Abromowitz & Stegun-
*	    6.1.16
*     ..................................................................
*
	function lngam(q)

	complex*16 q		! input value
	complex*16 z		! input value + inc
	complex*16 lngam	! value of the ln(gamma) function
	complex*16 pi		! value of constant pi
	complex*16 corr		! correction for integer increment
				! added to accelerate convergence
	real*16 tangt		! tangent pi/4
	real*16 pix		! real part of pi
	real*16 piy		! imaginary part of pi
	integer*2 flag		! increment added flag recursion
				! correction required
	integer*2 inc		! convergence acceleration increment

*	initialize variables

	tangt = 1.
	pix = 4.*qatan(tangt)
	piy = 0.
	pi = dcmplx(pix,piy)
	inc = 10

	if(dreal(q).le.10.) then
	  z = q + inc		! add integer to accelerate
				! convergence

*	  calculate ln gamma using A&S 6.1.41

	lngam = (z - 0.5)*cdlog(z) - z + 0.5 * cdlog(2*pi) + 1./(12.*z)
     @   - 1./(360 * z**3) + 1./(1260. * z**5) - 1./(1680. * z**7)

*	  Since convergence has been accelerated by addition of
*	  an integer to the real part of q then use
*	  the log of recursion relation A&S 6.1.16 to extract
*	  the correct value

	  corr = 0.
	  do 100 k=1,inc
100	    corr = corr + cdlog((k-1)+q)
	  lngam = lngam - corr

	 else

*	  calculate ln gamma using A&S 6.1.41

	lngam = (q - 0.5)*cdlog(q) - q + 0.5 * cdlog(2*pi) + 1./(12.*q)
     @   - 1./(360 * q**3) + 1./(1260. * q**5) - 1./(1680. * q**7)

	 endif

	return
	end

*
*
*     ..................................................................
*
*        SUBROUTINE SIGMR
*
*        PURPOSE
*	    Calculate the ratio of the differential Mott Cross-section
*	    to the differential Rutherford Cross-section
*
*        USAGE
*           CALL SIGMR (Z,AV,THETA,RATIO)
*
*        DESCRIPTION OF PARAMETERS
*	    Z	   - atomic number of target
*	    AV	   - accelerating voltage
*	    theta  - electron scattering angle
*	    ratio  - ratio of Mott Cross-section to Rutherford CS returned
*
*        REMARKS
*	    real*16
*
*        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
*	    EULER  - euler transformation
*	    LGDRE  - legendre function
*	    D	   - coefficient precursors
*	    LNGAM  - natural logarithm of complex gamma functions
*
*        METHOD
*	    Direct evaluation of the series
*	    G1 and F1 convergence accelerated by use of
*	    the euler transformation
*	 see:
*	    McKinley & Feshbach
*	    Phys Rev 74,1759(1948)
*	 and;
*	    Feshbach
*	    Phys Rev 88,295(1952)
*	 and;
*	    Doggett & Spencer
*	    Phys Rev 103,1597(1956)
*	 and;
*	    O. S. Oen
*	    CROSS-SECTIONS FOR ATOMIC DISPLACEMENTS
*	    IN SOLIDS BY FAST ELECTRONS
*	    ORNL Report # TM-4897,(1973)
*
*     ..................................................................
*

	subroutine sigmr(Z,AV,theta,ratio)
   
	COMMON /CONST/elmass,amu,evpKeV,fsc,CLrad,hbar,c,pi,one
	COMMON /REPEAT/AVold

	parameter(maxord = 120)	! max number of coefficients to evalutate
				! near 125 the parameter 2**(i+1) in subroutine
				! euler will overflow.
*	variables used

	character*1 comma	! comma

	integer*2 k		! coef order
	integer*2 nterms	! num terms to use in sums

	real*16 Z		! atomic number of target
	real*16 AV		! incident beam K.E.
	real*16 theta		! scattering angle
	real*16 ratio		! ratio of Mott CS to Rutherford
	real*16 x		! cosine of the scattering angle 
	real*16 kr		! real translation of k
	real*16 P(0:maxord)	! values of legendre fcns
	real*16 zero		! constants
	real*16 alpha		! fsc * Z
	real*16 beta		! relativistic beta = v/c
	real*16 q		! McKinley Feshbach variable, alpha/beta
	real*16 kpsqr		! (momentum/hbar)**2
				! in E-12 cm -- e**2/mc**2
	real*16 y(maxord)	! input array - terms in series to sum
	real*16 sumy		! estimated sum
	real*16 zz		! auxiliary variable
	real*16 trad		! theta in radians

*	DECLARATIONS FOR COMMON BLOCK REPEAT
	real*16 AVold	! repeated energy ?

*	DECLARATIONS FOR COMMON BLOCK CONST - CONSTANTS
	real*16 elmass	! rest mass of electron (KeV)
	real*16 amu	! atomic mass unit (KeV)
	real*16 evpKeV	! ev/KeV conversion
	real*16 fsc	! fine structure constant
	real*16 CLrad	! classical radius of the electron
			! in E-12 cm -- e**2/mc**2
	real*16 hbar	! h bar quantum mechanical const.(eV s)
	real*16 c	! speed of light(cm/s)
	real*16 pi	! constant
	real*16 one	! constant

	complex*16 coschi	! cos(CHI) in small angle asymtotic reln
	complex*16 arga,argb,ga,gb	! utility variables
	complex*16 Im		! i
	complex*16 Aco,B,S,V	! partial evaluation of coef in series
	complex*16 F0,G0	! terms in Mott Cross-Section
	complex*16 F1(0:maxord)	! terms in Mott Cross-Section
	complex*16 G1(0:maxord)	! terms in Mott Cross-Section
	complex*16 Fcoef(0:maxord)	! F1 = f0 + f1*P1 + f2*P2
	complex*16 Gcoef(0:maxord)	! similarly
	complex*16 sumF1	! F1 sum
	complex*16 sumG1	! G1 sum
	complex*16 F,G		! Mott series sums
	complex*16 D		! series coefficient precursors
	complex*16 lngam	! ln gamma fcns

	external euler		! use Euler Transformation to sum series
	external lgdre		! subroutine to calc legendre fcns
	external D		! function to calc coefficients
	external lngam		! natural log of the complex gamma fcn

*---------------------------------

*	INITIALIZE VARIABLES

	comma = ','
	zero = 0.
	Im = dcmplx(zero,one)
	alpha = Z*fsc
	beta = sqrt(1.-1./(1.+AV/elmass)**2)
	q = alpha/beta
	kpsqr = ((AV+elmass)**2 - elmass**2)*evpKeV**2/((hbar*c)**2) ! (cm**-2)
	kpsqr = kpsqr*1E-24		! convert units to barns**-1


*	AVOID DIVERGENCE AT 180 DEGREES
	if(theta.eq.180.) theta = 178.

	x = qcosd(theta)
	call lgdre(maxord,x,P)		! generate legendre polynomials

*	GENERATE RATIO USING ASYMPTOTIC RELN
*	cited by Doggett & Spencer
*	if angle less than 0.5 degrees

	arga = 0.5 - Im*q
	argb = one + Im*q
	ga = lngam(arga)
	gb = lngam(argb)
	coschi = cdexp(ga+gb-conjg(ga)-conjg(gb))
	trad = theta*pi/180.
	if(theta.le.0.5) then
	  ratio = dreal(coschi)*pi*alpha*beta*trad/2. + 1.
*	  type*,'asymtotic reln used'
	  goto 250
	 else
	endif

	nterms = 75		! fixed number of terms to evaluate

*	EVALUATE SERIES

	zz = (qsind(theta/2.))**2
	F0 = Im*0.5*cdexp(Im*q*qlog(zz))
     @     *cdexp(lngam(one-Im*q))/cdexp(lngam(one+Im*q))

	zz = (qtand(theta/2.))**2
	G0 = -1.*Im*q*F0*(1.- x)/zz	! (1-cos)G0

	if(AV.eq.AVold) goto 150	! if accel voltage unchanged
	AVold = AV		! do not re-evaluate coeff.

*	Evaluate Fcoef, Gcoef
*	The series G1 has been multipied by (1-cos(theta))
*	and the result expanded in Legendre Polynomials.
*	The coefficients are expressed in terms of the coefficients
*	D sub k .  This results in faster convergence near 5 to 15 degrees

	k = 0		! zeroth term
	  Aco = 0.
	  B = D(k,Z,AV)
	  S = D(k+1,Z,AV)
	do 100 k=0,nterms
	  kr = qfloat(k)
	  V = D(k+2,Z,AV)
	  Fcoef(k) = B*kr+S*(kr+1.)
	  Gcoef(k) = (kr-1.)**2*kr*Aco/(2.*kr-1.)
	    Gcoef(k) = Gcoef(k) + kr**2*(kr-1.)*B/(2.*kr-1.)
	    Gcoef(k) = Gcoef(k) - (kr+1.)**2*(kr+2.)*S/(2.*kr+3.)
	    Gcoef(k) = Gcoef(k) - (kr+2.)**2*(kr+1.)*V/(2.*kr+3.)
	  Aco = B		! increment functions
	  B = S		! should be faster than recomputing
	  S = V		! functions each time
100	continue
150	continue
	do 200 k=0,nterms
	G1(k) = Gcoef(k)*P(k)
200	F1(k) = Fcoef(k)*P(k)

*	CALCULATE SUM F1 AND SUM G1

	do 550 i=1,maxord		! Im(sumF1)
550	  y(i) = dreal(F1(i-1))
	call euler(y,nterms,sumy)
	sumF1 = Im*0.5*dcmplx(sumy,zero)

	do 575 i=1,maxord		! Re(sumF1)
575	  y(i) = dimag(F1(i-1))
	call euler(y,nterms,sumy)
	sumF1 = sumF1 - 0.5*sumy

	do 650 i=1,maxord		! Im(sumG1)
650	y(i) = dreal(G1(i-1))
	call euler(y,nterms,sumy)
	sumG1 = Im*0.5*dcmplx(sumy,zero)

	do 675 i=1,maxord		! Re(sumG1)
675	  y(i) = dimag(G1(i-1))
	call euler(y,nterms,sumy)
	sumG1 = sumG1 - 0.5*sumy

*	SUM SERIES

	F = F0 + sumF1
	G = G0 + sumG1

*	Evaluate the ratio of the Mott cross-section
*	to the Rutherford cross-section

	ratio = q**2*conjg(F)*F*(1.-x) + conjg(G)*G/((1.+x)*(1.-beta**2))
	ratio = 2.*beta**4*ratio/(z**2*CLrad**2*kpsqr)

1000	continue

250	continue
	return
	END