subroutine exponential_torus(params,spp)
TYPE(KORC_PARAMS), INTENT(IN) :: params
TYPE(SPECIES), INTENT(INOUT) :: spp
REAL(rp) :: fl
REAL(rp) :: fr
REAL(rp) :: fm
REAL(rp) :: rl
REAL(rp) :: rr
REAL(rp) :: rm
REAL(rp) :: relerr
REAL(rp), DIMENSION(:), ALLOCATABLE :: r
REAL(rp), DIMENSION(:), ALLOCATABLE :: theta
REAL(rp), DIMENSION(:), ALLOCATABLE :: zeta
INTEGER :: pp
ALLOCATE( theta(spp%ppp) )
ALLOCATE( zeta(spp%ppp) )
ALLOCATE( r(spp%ppp) )
! Initial condition of uniformly distributed particles on a
! disk in the xz-plane
! A unique velocity direction
call init_u_random(10986546_8)
call init_random_seed()
call RANDOM_NUMBER(theta)
theta = 2.0_rp*C_PI*theta
call init_random_seed()
call RANDOM_NUMBER(zeta)
zeta = 2.0_rp*C_PI*zeta
! Uniform distribution on a disk at a fixed azimuthal theta
call init_random_seed()
call RANDOM_NUMBER(r)
! Newton-Raphson applied here for finding the radial distribution
do pp=1_idef,spp%ppp
rl = 0.0_rp
rr = spp%r_outter
fl = fzero(rl,spp%r_outter,spp%falloff_rate,r(pp))
fr = fzero(rr,spp%r_outter,spp%falloff_rate,r(pp))
if (fl.GT.korc_zero) then
relerr = 100*ABS(fl - fr)/fl
else
relerr = 100*ABS(fl - fr)/fr
end if
do while(relerr.GT.1.0_rp)
rm = 0.5_rp*(rr - rl) + rl
fm = fzero(rm,spp%r_outter,spp%falloff_rate,r(pp))
if (SIGN(1.0_rp,fm).EQ.SIGN(1.0_rp,fr)) then
rr = rm
else
rl = rm
end if
fl = fzero(rl,spp%r_outter,spp%falloff_rate,r(pp))
fr = fzero(rr,spp%r_outter,spp%falloff_rate,r(pp))
if (fl.GT.korc_zero) then
relerr = 100*ABS(fl - fr)/fl
else
relerr = 100*ABS(fl - fr)/fr
end if
end do
r(pp) = rm
end do
spp%vars%X(:,1) = ( spp%Ro + r*COS(theta) )*SIN(zeta)
spp%vars%X(:,2) = ( spp%Ro + r*COS(theta) )*COS(zeta)
spp%vars%X(:,3) = spp%Zo + r*SIN(theta)
DEALLOCATE(theta)
DEALLOCATE(zeta)
DEALLOCATE(r)
end subroutine exponential_torus