program energ_loss_double
	
* This program calcultaes the energy loss
* by electron after travelling through a
* specific material.
*
*
*	dEdx = dEdx_coll + dEdx_rad
*
*
* Draft: May 28, 1996
* Author: Paul Gueye
*
*
* Update: February 14, 1997
* The formula above is from:
* \item W. Leo, {\it Techniques for Nuclear and Particle Physics
* 	Experiments}, Springer Verlag, (1987).
* \item R. M. Sternheimer et al., {\it Atomic Data and Nuclear
* 	Data Tables}, {\bf 30-2}, p261-271, (1984)
* \item Review of Particle Physics, Part I, {\it Phys. Rev. D},
* 	{\bf vol54}, (July 1996)
*
* Comparison with Monte Carlo energy loss formula:
*	dE_e (MeV) = 0.1536(Z/A)t[19.26 + ln(t/rho)]
*	dE_p (MeV) = 0.3072(Z/A)(t/beta^2)
*		[(1/2)ln{(2m_e/16Z^0.9)(beta^2/(1-beta^2)} - beta^2]
*


	implicit double precision (a-h,o-z)
	
	real*8	A_med,Z_med,rho_med
	real*8	I_med,C_med,ap_med,m_med
	real*8	X_med,X0_med,X1_med,delta_med
	real*8	N_avog,pi,N_atoms,a
	real*8	r_e,m_e,p_e,T_e,T_emin,T_emax,beta_e,gamma_e
	real*8	K_coll,tau,F_tau1,F_tau2,F_tau3,F_tau
	real*8	dEdx1_coll,dEdx2_coll,dEdx3_coll,dEdx_coll
	real*8	dEdx1_rad,dEdx2_rad,dEdx3_rad,dEdx4_rad,dEdx_rad
	real*8	dEdx,alpha,K_rad,top_rad,phi_rad
	real*8	deltae,dT_e,thick
	real*8	neta,C1_shell,C2_shell,C3_shell,C4_shell
	real*8	C5_shell,C6_shell,C7_shell,C8_shell,C_shell
	real*8	F1_Z,F2_Z,F3_Z,F_Z
	real*8	dE
	integer	medium,restart
	
	T_e = 0.
		
* ==> Open files
	open(unit=1,file='energy_loss.dat')

* ==> Give parameters for electron and medium
1	write(6,*)
	write(6,*)
	write(6,*)'Thickness of material (in cm)'
	read(5,*)thick
	write(6,*)'Energy range and step [T_emin,T_emax,Deltae] (in MeV)'
	read(5,*)T_emin,T_emax,deltae

	if (T_e.lt.0) then
	write(6,*)'Negative kinetic energy impossible !'
	goto 1
	endif

2	write(6,*)
	write(6,*)
	write(6,*)'Choose Medium: Al(1), Air_gas(2), Cu(3), Scint.Plas(4)'
	write(6,*)'Pb(5), H2O(6), W(7)'
	read(5,*)medium
	
	if (medium.eq.1) goto 10
	if (medium.eq.2) goto 20
	if (medium.eq.3) goto 30
	if (medium.eq.4) goto 40
	if (medium.eq.5) goto 50
	if (medium.eq.6) goto 60
	if (medium.eq.7) goto 70
	write(6,*)'Medium does not exists in nature !!'
	goto 2

10	Z_med	= 13.d0
	A_med	= 26.98d0
	rho_med	= 2.7d0
	I_med	= 166.d0
	C_med	= -4.24d0
	ap_med	= 0.0802d0
	m_med	= 3.63d0
	X0_med	= 0.1708d0
	X1_med	= 3.01d0
	goto 100
20	Z_med	= 7.d0
	A_med	= 14.d0
	rho_med	= 1.205d0
	I_med	= 85.7d0
	C_med	= -10.6d0
	ap_med	= 0.1091d0
	m_med	= 3.4d0
	X0_med	= 1.742d0
	X1_med	= 4.28d0
	goto 100
30	Z_med	= 29.d0
	A_med	= 63.55d0
	rho_med	= 8.96d0
	I_med	= 322.d0
	C_med	= -4.42d0
	ap_med	= 0.1434d0
	m_med	= 2.9d0
	X0_med	= -0.0254d0
	X1_med	= 3.28d0
	goto 100
40	Z_med	= 5.577d0
	A_med	= 11.1081d0
	rho_med	= 1.049d0
	I_med	= 64.7d0
	C_med	= -3.20d0
	ap_med	= 0.161d0
	m_med	= 3.24d0
	X0_med	= 0.1464d0
	X1_med	= 2.49d0
	goto 100
50	Z_med	= 82.d0
	A_med	= 207.19d0
	rho_med	= 11.35d0
	I_med	= 823.d0
	C_med	= -6.20d0
	ap_med	= 0.0936d0
	m_med	= 3.16d0
	X0_med	= 0.3776d0
	X1_med	= 3.81d0
	goto 100
60	Z_med	= (2*1.+8)*1.d0
	A_med	= (2*1.+16)*1.d0
	rho_med	= 1.0d0
	I_med	= 75.d0
	C_med	= -3.5017d0
	ap_med	= 0.09116d0
	m_med	= 3.4773d0
	X0_med	= 0.24d0
	X1_med	= 2.8004d0
	goto 100
70	Z_med	= 74*1.d0
	A_med	= 183.851*1.d0
	rho_med	= 19.3d0
	I_med	= 727.d0
	C_med	= -5.4059d0
	ap_med	= 0.15509d0
	m_med	= 2.8447d0
	X0_med	= 0.2167d0
	X1_med	= 3.4960d0
	goto 100
	
	
* ==> Constants
100	pi	= 3.14592654d0
	N_avog	= 6.0221367d23
	r_e	= 2.81794092d-15
	m_e	= 0.510999d0
	alpha	= 1.d0/137.d0
	K_coll	= 2.d0 * pi * N_avog * r_e**2 * m_e * 1.d4
	T_e	= T_emin
	K_rad	= 4.d0 * Z_med**2 * r_e**2 * alpha
	top_rad	= 137.d0 * m_e * (Z_med)**(1.d0/3.d0)
	N_atoms	= rho_med * N_avog / A_med
	a	= Z_med / 137.d0

	write(1,101)
	write(1,102)

	write(1,*)' '
101	format(5x,'Energy',10x,'dE/dx_coll',5x,'dE/dx_rad',
     >		7x,'dE/dx',8x,'Einc-dE',8x,'dE')
102	format(5x,'[MeV]',10x,'[MeV/(g/cm2)]',1x,'[MeV/(g/cm2)]',
     >		1x,'[MeV/(g/cm2)]',6x,'[MeV]',8x,'[MeV]')

	write(6,*)
	write(6,*)
	
* ==> dEdx_coll calculation
* ==> p_e,beta_e
103	p_e	= dsqrt(T_e * (T_e + (2.d0 * m_e)))
	beta_e	= p_e/dsqrt(m_e**2 + p_e**2)
	gamma_e	= 1.d0 /dsqrt(1.d0 - beta_e**2)

* ==> X calcuation
	if (T_e.eq.0.d0)then
	X_med	= 0.d0
	else
	X_med	= dlog(beta_e * gamma_e )/dlog(10.d0)
	endif
	
	if (X_med.lt.X0_med) delta_med = 0.d0
	if (X_med.gt.X0_med .and. X_med.lt.X1_med)
     >		delta_med = (4.6052d0 * X_med) + C_med +
     >		(ap_med * (X1_med - X_med)**m_med)
	if (X_med.gt.X1_med) delta_med = (4.6052d0 * X_med) + C_med

* ==> Tau
	tau	= T_e/m_e

* ==> F_tau
	F_tau1	= 1.d0 - beta_e**2
	F_tau2	= (tau**2)/8.d0 - ((2.d0 * tau + 1.d0 ) * dlog(2.d0))
	F_tau3	= (tau + 1.d0)**2
	F_tau	= F_tau1 + (F_tau2/F_tau3)

* ==> neta
	neta = beta_e * gamma_e
	
	if (neta.ge.0.1d0) then
* ==> Shell correction
	C1_shell = 0.422377d0 * neta**(-2)
	C2_shell = 0.0304043d0 * neta**(-4)
	C3_shell = 0.00038106d0 * neta**(-6)
	C4_shell = 3.85019 d0* neta**(-2)
	C5_shell = 0.1667989d0 * neta**(-4)
	C6_shell = 0.00157955d0 * neta**(-6)
	C7_shell = (C1_shell + C2_shell - C3_shell)*1.d-6*(I_med**2)
	C8_shell = (C4_shell - C5_shell + C6_shell)*1.d-9*(I_med**3)
	C_shell = C7_shell + C8_shell
	
	goto 104
	
	else
	write(6,*)'neta < 0.1 for T_e =',T_e,' MeV'
	goto 111
	endif

* ==> dEdx_coll
104	dEdx1_coll = K_coll * (Z_med/A_med) * (1.d0/beta_e**2)
	dEdx2_coll = dlog((tau**2 * (tau + 2.))
     >		/(2.d0 * ((I_med * 1d-6)/m_e)**2))
	dEdx3_coll = F_tau - delta_med - (2.d0 * (C_shell/Z_med))
	dEdx_coll  = dEdx1_coll * (dEdx2_coll + dEdx3_coll)

* ==> F_Z calculation
	F1_Z = a**2
	F2_Z = (1.d0 + a**2)**(-1)
	F3_Z = 0.20206d0 - 0.0369d0*a**2 + 0.0083d0*a**4 - 0.002d0*a**6
	F_Z  = F1_Z * (F2_Z + F3_Z)

* ==> dEdx_rad calculation
	dEdx1_rad = dlog((2.d0 * T_e)/m_e)
	dEdx2_rad = -1.d0/3.d0 - F_Z
	dEdx3_rad = dlog(183.d0 * (Z_med)**(-1.d0/3.d0))
	dEdx4_rad = 1.d0/18.d0 - F_Z 

* ==> Flag
	if (T_e.ge.m_e .and. T_e.le.top_rad)
     >		phi_rad	= K_rad * (dEdx1_rad + dEdx2_rad)
	if (T_e.ge.top_rad)
     >		phi_rad	= K_rad * (dEdx3_rad + dEdx4_rad)
	
* ==> dEdx_rad
	dEdx_rad = N_atoms * T_e * phi_rad
	
* ==> dEdx
	dEdx	= dEdx_coll + dEdx_rad
	
* ==> dT_e = Einc-dE
	dE		= dEdx * rho_med * thick
	dT_e	= T_e - dE

* ==> Write output	
	write(1,110)T_e,dEdx_coll,dEdx_rad,dEdx,dT_e,dE
110	format(4(5x,e10.5),2(5x,e10.4))

	
111	T_e	= T_e + deltae
	if (T_e.gt.T_emax) goto 900
	
	goto 103
	
* ==> Restart ?
900	close(unit=1)

	write(6,*)
	write(6,*)
	write(6,*)'Do another calculation: yes (1) or no (2) ?'
	read(5,*)restart
	if (restart.eq.1) then
	goto 1
	elseif (restart.eq.2) then
	goto 1000
	else
	write(6,*)'Answer the question !'
	goto 900
	endif
	
	
1000	end