Title                  :TOTCS
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
Abstract:
	This program calculates displacement cross-sections that are useful
in estimating the amount of sputtering and displacement damage due to
collisions of relativistic electrons with nuclei in elemental samples. The
total displacement cross-section is obtained by evaluating the differential
Mott cross-section and numerically integrating over the range of allowed
angles.  An upper limit for cascade cross-sections is evaluated in version 2.0. 
For the medium and high voltage 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. This program allows
the analyst to estimate the potential adverse effects on the specimen prior
to study in the TEM.
Title                  :TOTCS
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:
	2.0

Purpose:
	Compute the total Mott cross-section for a range of accelerating
	voltages running up to 1500KV.	The calculation method of
	Doggett and Spencer is used and results agree with O. S. Oen's
	implementation of that calculation to within 0.5%. The program
	requires as input Z - the atomic number, A - the atomic mass,
	and Ed the assumed displacement energy. Preliminary results
	of this study are to be published in:

	Proceedings of the Workshop on Analytical Electron Microscopy
	San Francisco Press, Fall 1987


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

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

Input/Output:
	UNIT 1	input	(to be assigned)
	UNIT 2	output	(to be assigned)

Test Data: 
	This routine is designed to read input from a data file
	with the format (3F10.3).

Input file format:
	Z,A,Es
	(EOF)
	where;	Z = atomic number
		A = atomic mass(amu)
		Es = Displacement/sputtering energy(eV).  The energy required
	to produce the defect or sputter an atom from the surface.
	The atom is bound by an isotropic square well of depth Es.

example:
	13.000,26.980,3.420
	(EOF)

General Comments: 
	To run the program it is convienient to create a DCL command
	file which assigns fortran device 1 to the input file and
	device 2 to an output file.
  
	$ASSIGN filename.INP  FOR001
	$ASSIGN filename.OUT  FOR002
	$RUN TOTCS

	Since the program takes a substantial amount of CPU time the
	program is usually submitted to a batch queue.  This can be done
	in VAX/VMS by typing:

	submit/noprint/notify <fn>

	where <fn> is the above command file's name.


	VAX/VMS CPU time is approximately 30 seconds for each point calculated.
	The number of points varies since the threshold is a function of the
	atomic mass, and Z.  For the example here 1 hr 35 min 56.61 sec
	of CPU time was used.

	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.
		QATRGL a 16-point Gauss-Legendre numerical quadrature routine

Sample Output: 
	The following is output generated by the above sample input
	which corresponds to Al with an assumed Es=3.420 eV.

 Atomic number of target =       13.
   Atomic mass of target =    26.980
       Sputtering Energy =     3.420
  Threshold Energy (KeV) =    40.436
  
 Accelerating	Cross-Section	Cascade CS
 Voltage(KV)	(Barns/atom)	(Barns/atom)
  
     40.436	     0.000	     0.000
     41.000	    22.380	    22.380
     42.000	    59.264	    59.264
     43.000	    92.871	    92.871
     44.000	   123.530	   123.530
     45.000	   151.533	   151.533
     46.000	   177.140	   177.140
     47.000	   200.578	   200.578
     48.000	   222.052	   222.052
     49.000	   241.743	   241.743
     50.000	   259.813	   259.813
     55.000	   330.434	   330.434
     60.000	   377.058	   377.058
     65.000	   407.918	   407.918
     70.000	   428.164	   428.164
     75.000	   441.128	   441.128
     80.000	   449.016	   449.143
     85.000	   453.328	   454.798
     90.000	   455.099	   458.956
     95.000	   455.058	   461.966
    100.000	   453.724	   464.089
    110.000	   448.570	   466.425
    120.000	   441.574	   467.070
    130.000	   433.788	   466.693
    140.000	   425.790	   465.713
    150.000	   417.901	   464.395
    160.000	   410.290	   462.909
    170.000	   403.045	   461.366
    180.000	   396.203	   459.837
    190.000	   389.773	   458.368
    200.000	   383.746	   456.988
    210.000	   378.107	   455.714
    220.000	   372.833	   454.556
    230.000	   367.901	   453.517
    240.000	   363.288	   452.599
    250.000	   358.970	   451.799
    260.000	   354.926	   451.114
    270.000	   351.133	   450.541
    280.000	   347.573	   450.074
    290.000	   344.228	   449.708
    300.000	   341.082	   449.438
    310.000	   338.120	   449.259
    320.000	   335.327	   449.165
    330.000	   332.692	   449.152
    340.000	   330.202	   449.215
    350.000	   327.847	   449.350
    360.000	   325.618	   449.551
    370.000	   323.505	   449.815
    380.000	   321.501	   450.138
    390.000	   319.598	   450.517
    400.000	   317.790	   450.947
    450.000	   309.961	   453.772
    500.000	   303.738	   457.487
    550.000	   298.703	   461.847
    600.000	   294.563	   466.672
    650.000	   291.114	   471.828
    700.000	   288.205	   477.219
    750.000	   285.726	   482.769
    800.000	   283.592	   488.422
    850.000	   281.742	   494.137
    900.000	   280.124	   499.880
    950.000	   278.700	   505.626
   1000.000	   277.438	   511.356
   1050.000	   276.314	   517.054
   1100.000	   275.307	   522.709
   1150.000	   274.401	   528.313
   1200.000	   273.580	   533.857
   1250.000	   272.836	   539.338
   1300.000	   272.156	   544.751
   1350.000	   271.534	   550.095
   1400.000	   270.963	   555.366
   1450.000	   270.437	   560.565
   1500.000	   269.951	   565.691
Title                  :TOTCS
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 TOTCS
*
*	Version		Date		Initials	Reason
*-----------------------------------------------------------------------
*	 1.0		6/25/87		CRB		Creation
*	 1.1		7/21/87		CRB		Modify AV selected
*	 2.0		8/31/87		CRB		Add Kinchin and Pease
*							cascade
*
*        Purpose
*           Compute the total mott cross-section
*	    for a range of accelerating voltages
*	    running up to AVmax set here to 1500KV
*
*        Subroutines and function subprograms required
*	    are included in this listing
*
*	 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 totcs

	COMMON Z,AV,A,Ed
	COMMON /CONST/elmass,amu,evpKeV,fsc,CLrad,hbar,c,pi,one
	COMMON /REPEAT/AVold

	integer*2 ke	! utility counter

	real*8 fun	! function to integrate
	real*8 fcn	! another function to integrate
	real*8 twoEd	! integration limit
	real*8 minusz
	real*8 xl 	! lower limit
	real*8 xu	! upper limit
	real*8 y	! integral approximation
	real*8 halfxl
	real*8 zero
	real*8 xulog
	real*8 intgl	! integral approximation
	real*8 intgl2	! integral approximation
	real*8 rth	! Rutherford C-S
	real*8 CS	! cross-section in barns (E-24 cm**2)
	real*8 CASCS	! apparent cross-section in barns (E-24 cm**2)
			! after cascade
	real*16 Z	! atomic number of target
	real*16 AV	! accelerating voltage
	real*16 AVmax	! max accelerating voltage
	real*16 Th	! threshold incident KE that delives Ed
	real*16 AVbase	! accelerating voltage base
	real*16 AVb2	! accelerating voltage base
	real*16 AVb3	! accelerating voltage base
	real*16 AVb4	! accelerating voltage base
	real*16 AVb5	! accelerating voltage base
	real*16 Ed	! Displacement or sputtering energy
	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

*	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

	external qatrgl		! integration routine
	external fun		! integrand
	external fcn		! integrand

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

	read(1,20) Z,A,Ed
20	format(3F10.3)

*---------OUTPUT TO DEVICE 2 ------------------------
	write(2,40) Z
40	format(' ','Atomic number of target =',F10.0)
	write(2,50) A
50	format(' ','  Atomic mass of target =',F10.3)
	write(2,60) Ed
60	format(' ','      Sputtering Energy =',F10.3)

*-------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

	alpha = Z * fsc

	Th = -1.*elmass + qsqrt(elmass**2 + A*amu*Ed/(2.*evpKeV))
*	type*,'th =',Th
	write(2,70) Th
70	format(' ',' Threshold Energy (KeV) =',F10.3)
	write(2,80)
80	format(' ',' ')
	write(2,85)
85	format(' ','Accelerating	Cross-Section	Cascade CS')
	write(2,90)
90	format(' ','Voltage(KV)	(Barns/atom)	(Barns/atom)')
	write(2,80)

	AVbase = qfloat(int(Th)+1)
	AVb2 = 10.*qfloat(int(AVbase/10.)+1)
	AVb3 = 20.*qfloat(int(AVb2/20.)+3)
	AVb4 = 50.*qfloat(int(AVb3/50.)+4)
	AVb5 = 100.*qfloat(int(AVb4/100.)+1)

	ke = 0
	AVmax = 1500.0

	do 2000 while (AV.lt.AVmax)	! carry calculation to 1500.0 kV
	if(ke.eq.0) then
	  AV = Th			! energy in KeV
	  ke = ke + 1
	else if(ke.eq.1)then
	  AV = AVbase
	  ke = ke + 1
	  if (AV.gt.AVb5) then
	   AV = 50.*qfloat(int(AV/50.)+1)
	   else
	   endif
	else if(AV.ge.AVb5.and.AV.le.AVmax) then
	  AV = AV + 50.0
	else if(AV.ge.AVb3) then
	  AV = AV + 10.0
	else if(AV.ge.AVb2.and.AV.lt.AVb3) then
	  AV = AV + 5.0
	else if(AV.ge.AVbase.and.AV.lt.AVb2) then
	  AV = AV + 1.0 
	endif

	beta = sqrt(1. - 1./(1.+ AV/elmass)**2)
	Em = (2./amu)*evpKeV * AV * (AV + 2.*elmass)/A	! threshold

	 xl = Em/Ed	! Integration limits
	 twoEd = 2.*Ed
	 halfxl = 0.5*xl
	 xu = 1.0
	 zero = 0.0
	 xulog = dlog(halfxl)
	 Rth = pi*(Z*Clrad)**2*(1.-beta**2)/beta**4

	IF(Em.le.Ed) then
	 CASCS = 0.0		! cross-section is zero if max energy
	 cs = 0.0		! transfer <= the sputtering energy
	else
	 call qatrgl(xu,xl,fun,intgl)
	 cs = intgl * Rth
	endif

	if(Em.gt.Ed.and.Em.le.twoEd) then
	 cascs = cs			! no cascade

	else if(Em.gt.twoEd) then
	 call qatrgl(halfxl,xl,fun,intgl)
	 call qatrgl(zero,xulog,fcn,intgl2)
	 intgl = intgl + 0.5*xl*intgl2
	 cascs = intgl*Rth
	endif


	write(2,95) AV , cs, CASCS
95	format(' ',F10.3,2('	',F10.3))
2000	continue

3000	continue
	call exit
	end

*
*     ..................................................................
*
*        FUNCTION FCN
*
*        PURPOSE
*           Compute one of the integrands required by
*           integration routine qatrgl
*
*        USAGE
*           X = FCN(W)
*
*        DESCRIPTION OF PARAMETERS
*           W      - Arguement of the integrand
*
*        REMARKS
*	    returns real*8 value
*
*        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
*           SIGMR
*
*        METHOD
*	    Implements two simple variable changes as
*	    suggested by O. S. OEN
*	    ORNL REPORT # TM-4897,(1973)
*	    and then calls subroutine sigmr
*     ..................................................................
*
	function fcn(W)

	COMMON Z,AV,A,Ed
	COMMON /CONST/elmass,amu,evpKeV,fsc,CLrad,hbar,c,pi,one
	COMMON /REPEAT/AVold

	real*8 fcn
	real*8 x
	real*8 W

	real*16 Z		! atomic number of target
	real*16 AV		! incident beam K.E.
	real*16 A		! atomic mass
	real*16 Ed		! sputtering energy
	real*16 theta		! scattering angle
	real*16 ratio		! ratio of Mott CS to Rutherford

*	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

	external sigmr	! calc differential mott CS

*	Initialize values

	x = 1.0/dexp(W)
	theta = qextd(2.*dasind(dsqrt(x)))	! chg variables again

	call sigmr(Z,AV,theta,ratio)

	fcn = ratio

	return
	end

*
*     ..................................................................
*
*        FUNCTION FUN
*
*        PURPOSE
*           Compute the integrand required by
*           integration routine qatrgl
*
*        USAGE
*           X = FUN(Y)
*
*        DESCRIPTION OF PARAMETERS
*           Y      - Arguement of the integrand
*
*        REMARKS
*	    returns real*8 value
*
*        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
*           SIGMR
*
*        METHOD
*	    Implements two simple variable changes as
*	    suggested by O. S. OEN
*	    ORNL REPORT # TM-4897,(1973)
*	    and then calls subroutine sigmr
*     ..................................................................
*
	function fun(y)

	COMMON Z,AV,A,Ed
	COMMON /CONST/elmass,amu,evpKeV,fsc,CLrad,hbar,c,pi,one
	COMMON /REPEAT/AVold

	real*8 fun
	real*8 x
	real*8 y

	real*16 Z		! atomic number of target
	real*16 AV		! incident beam K.E.
	real*16 A		! atomic mass
	real*16 Ed		! sputtering energy
	real*16 theta		! scattering angle
	real*16 ratio		! ratio of Mott CS to Rutherford

*	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

	external sigmr	! calc differential mott CS

*	Initialize values

	x = 1./y				! chg variables
	theta = qextd(2.*dasind(dsqrt(x)))	! chg variables again

	call sigmr(Z,AV,theta,ratio)

	fun = ratio

	return
	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 = 179.99

	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
*
*     ..................................................................
*
*        SUBROUTINE QATRGL
*
*        PURPOSE
*           Compute an approximation for integral(fct(x), summed
*           over x from xl to xu).
*
*        USAGE
*           CALL QATR (XL,XU,FCN,INT)
*           parameter fcn requires an external statement.
*
*        DESCRIPTION OF PARAMETERS
*           XL     - The lower bound of the interval.
*           XU     - The upper bound of the interval.
*           FCN    - The name of the external function subprogram used.
*           INT    - The resulting approximation for the integral value.
*
*        REMARKS
*
*        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
*           the external function subprogram fcn(x) must be coded by
*           the user.
*
*        METHOD
*           Uses gauss legendre quadrature
*           SEE:
*           Abromowitz and Stegun
*           HANDBOOK OF MATHEMATICAL FUNCTIONS
*           25.4.30, AND PG 916
*           Dover
*           and,
*           NUMERICAL ANALYSIS
*           Francis Scheid
*           SCHAUM'S OUTLINE SERIES
*           MCGRAW HILL CO
*     ..................................................................
*
      SUBROUTINE QATRGL(XL,XU,FCN,INT)
*
	COMMON /REPEAT/AVold

	real*8 xl	! lower limit
	real*8 xu	! upper limit
	real*8 int	! integral value

	real*8 x(8)	! evaluation points
	real*8 w(8)	! weights
	real*8 y	! utility variable
*
*	DECLARATIONS FOR COMMON BLOCK REPEAT
	real*16 AVold	! repeated energy ?

	data x/0.095012509837637440185,0.281603550779258913230,
     @  0.458016777657227386342,0.617876244402643748447,
     @  0.755404408355003033895,0.865631202387831743880,
     @  0.944575023073232576078,0.989400934991649932596/

	data w/0.189450610455068496285,0.182603415044923588867,
     @  0.169156519395002538189,0.149595988816576732081,
     @  0.124628971255533872052,0.095158511682492784810,
     @  0.062253523938647892863,0.027152459411754094852/
*
	int = 0.0	! zero integral
	do 100 i=1,8	! 16 point integration
	do 100 k=1,2
	y = (xu-xl)*x(i)*(-1.)**k/2. + (xu + xl)/2.
100	int = int+fcn(y)*w(i)
	int = (xu - xl)*int/2.0
*
      RETURN
      END
*
*	The following is a two point integration
*	which should provide only modest results,
*	but can be patched in when debugging to
*	make the program run more quickly.
*	The patch is effected by changing the
*	name of the subroutine in the call statement
*	from QATRGL to xQATRGL.
*	The two point integration is identical
*	except for the printing of a warning.
*
      SUBROUTINE xQATRGL(XL,XU,FCN,INT)
*
	COMMON /REPEAT/AVold

	real*8 xl	! lower limit
	real*8 xu	! upper limit
	real*8 int	! integral value
	real*8 x(1)	! evaluation points
	real*8 w(1)	! weights
	real*8 y	! utility variable
*
*	DECLARATIONS FOR COMMON BLOCK REPEAT
	real*16 AVold	! repeated energy ?

	data x/0.577350269189626/

	data w/1.0/
*
	write(2,200)
200	format(' ','WARNING - two point integration in use')
*
	int = 0.0	! zero integral
	do 100 k=1,2
	y = (xu-xl)*x(i)*(-1.)**k/2. + (xu + xl)/2.
100	int = int+fcn(y)*w(i)
	int = (xu - xl)*int/2.0
*
      RETURN
      END