module globals
	implicit none
!’PˆÊ‚Ícm,g,MeV
	real(8) ,parameter ::x=1.d0
	real(8) ,parameter ::Z=1.d0,A=1.d0,rho=2.d0/18.d0
	real(8) ,parameter ::Lrad=716.4d0&
				*A/(Z*(Z+1.d0)*log(287.d0/sqrt(Z)))
	real(8) ,parameter ::E0=100.d0	!100MeV
	real(8) ,parameter ::Me=0.511d0	!0.511MeV
	real(8) ,parameter ::re=2.817d-13,Na=6.022d23,const=0.1535d0
	real(8) ,save ::beta
end module globals

module subprogs
	use globals
	implicit none
contains
	function I() result(y)
		real(8) y
		if(Z<13.d0) then
			y=Z*(12.d0+7.d0/Z)
		else
			y=Z*(9.76d0+58.8d0*Z**(-1.19d0))
		endif
	endfunction
	function F(epsilon) result(y)
		real(8) epsilon,y,tau,beta
		tau=epsilon/Me
		beta=sqrt(1.d0-1.d0/tau**2)
		y=1.d0-beta**2+( tau**2/8.d0-(2*tau+1.d0)*log(2.d0) )&
				/(tau+1.d0)**2
	endfunction
	function B(epsilon) result(y)
		real(8) epsilon,y,tau,beta
		tau=epsilon/Me
		beta=sqrt(1.d0-1.d0/tau**2)
		y=const*rho*Z/(A*beta**2)&
		*(log(tau**2*(tau+2.d0)/(2.d0*(I()/Me)**2))+F(epsilon))
	endfunction
	function G(epsilon) result(y)
		real(8) epsilon,y
		y=-epsilon/Lrad-B(epsilon)
	endfunction
end module subprogs

program e_loss
	use globals
	use subprogs
	implicit none
	integer ::j,N=10000
	real(8) E,dE,dx
	real(8) k1,k2,k3,k4
	dx=x/N
	E=E0
	open(10,file='e_loss.dat')

	do j=1,N
		k1=dx*G(E)
		k2=dx*G(E+k1/2.d0)
		k3=dx*G(E+k2/2.d0)
		k4=dx*G(E+k3)
		E=E+1.d0/6.d0*(k1+2*k2+2*k3+k4)
		if(mod(j,10)==0) then
			write(10,*) j,E
		endif
	enddo
		dE=E-E0
	write(*,*) "Lrad=",Lrad/rho,"E=",E,"dE=",dE
endprogram