module double
integer, parameter :: dp = selected_real_kind(16)
end module double
module constants
use double
real(kind=dp), parameter :: mass = 9.109E-31
real(kind=dp), parameter :: charge = 1.6E-19
real(kind=dp), parameter :: epsilon = 8.85E-12
real(kind=dp), parameter :: pi=3.14
real(kind=dp), parameter :: L = 2.0
real(kind=dp), parameter :: p_energy = 1000.0
real(kind=dp), parameter :: c = 3.0E+8
real(kind=dp), parameter :: a=4.5E-3
real(kind=dp), parameter :: b=7.62E-3
end module constants
module density_mod
interface density_int
function density(current, energy, switch)
use constants
use double
implicit none
integer, optional :: switch
real :: current, density, energy,temp ! c - mA, e - keV
end function density
end interface density_int
end module density_mod
program test
use double
use density_mod
implicit none
real(kind=dp) :: new
new = density(2490.0,9.0)
write(*,*) new
end program test
! density function
function density(current, energy, switch)
use double
use constants
implicit none
integer,optional :: switch
real :: current, energy, density, temp
temp = sqrt(mass/(2.0*charge))*1.0/(2.0*pi)*current*1.E-3*1.0/sqrt(energy*1.E+3)
if(present(switch)) then
select case(switch)
case(1)
density=temp*2.0/a**2
return
case(2)
density=temp*(b**2-a**2)
return
case(3)
density=temp*1.0/(2.0/9.0*a**2)*(1.0-exp(-9.0/2.0))
return
case default
density=temp*1.0/(2.0/9.0*a**2)*(1.0-exp(-9.0/2.0))
return
end select
else
density=temp*1.0/(2.0/9.0*a**2)*(1.0-exp(-9.0/2.0))
return
end if
end function density
