FIXME : Add documentation
If vcrit is specified, there is no way of knowing its radial derivative, so set dvc/drho = 0
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(in) | :: | is | |||
real, | intent(in) | :: | v | |||
real, | intent(inout) | :: | f0 | |||
real, | intent(inout) | :: | df0dE_out | |||
real, | intent(inout) | :: | f0prim |
subroutine eval_f0_sdanalytic(is,v,f0,df0dE_out,f0prim)
use constants, only: pi
use mp, only: mp_abort
implicit none
integer,intent(in):: is
real,intent(in):: v
real,intent(inout):: f0,df0dE_out, f0prim
real:: A, df0dv,vcprim, vcrit
if (spec(is)%vcrit .LE. 0.0) then
call calculate_slowingdown_parameters(is,vcrit,vcprim)
if (abs(spec(is)%vcprim + 999.0) .LT. 1.0e-3) vcprim = spec(is)%vcprim
else
!> If vcrit is specified, there is no way of knowing its radial derivative, so set dvc/drho = 0
vcrit = spec(is)%vcrit
if (abs(spec(is)%vcprim + 999.0) .LT. 1.0e-3) vcprim = 0.0
end if
A = (3.0/(4.0*pi))/log( (vcrit**3 + 1.0)/(vcrit**3))
f0 = A/(vcrit**3 + v**3)
df0dv = -A*3.0*v**2/(vcrit**3 + v**3)**2
df0dE_out = (0.5/v)*df0dv / ( f0 * spec(is)%temp )
f0prim = -spec(is)%fprim - (-(vcrit**3/(vcrit**3 + v**3)) + (1.0/( log(1.0 + vcrit**(-3)) * (1.0 + vcrit**(3)))))*vcprim
if (v .GT. 1.0) then
f0 = 0.0
f0prim = 0.0
df0dE_out = 1.0
end if
! write(*,*) v, vcrit,vcprim,f0prim
end subroutine eval_f0_sdanalytic