module globals
	implicit none
!単位はcm,g,MeV
	real(8) ,parameter ::x=8.d-1,rho0=3.3d0
	real(8) ,parameter ::Z1=20.d0,A1=40.d0,rho1=rho0*40.d0/56.d0
	real(8) ,parameter ::Z2=8.d0,A2=16.d0,rho2=rho0*16.d0/56.d0
	real(8) ,parameter ::w1=1.d0*A1/(A1+A2),w2=1.d0*A2/(A1+A2)
			!wi=ai*Ai/Atotal(aiは分子内のiの個数)
	real(8) ,parameter ::Lrad1=716.4d0&
				*A1/(Z1*(Z1+1.d0)*log(287.d0/sqrt(Z1)))
	real(8) ,parameter ::Lrad2=716.4d0&
				*A2/(Z2*(Z2+1.d0)*log(287.d0/sqrt(Z2)))
	real(8) ,parameter ::Lrad=1.d0/( w1/Lrad1+w2/Lrad2 )
	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(Z) result(y)
		real(8) Z,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(Z,A,rho,epsilon) result(y)
		real(8) Z,A,rho,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(Z)/Me)**2))+F(epsilon))
	endfunction
	function G(epsilon) result(y)
		real(8) epsilon,y
		y=-epsilon/Lrad-w1*rho0/rho1*B(Z1,A1,rho1,epsilon)&
			-w2*rho0/rho2*B(Z2,A2,rho2,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_loss2.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(*,*) Lrad1/rho1,Lrad2/rho2,"Lrad=",Lrad/rho0,"E=",E,"dE=",dE
endprogram