subroutine gaussian_torus(params,spp)
TYPE(KORC_PARAMS), INTENT(IN) :: params
TYPE(SPECIES), INTENT(INOUT) :: spp
REAL(rp), DIMENSION(:), ALLOCATABLE :: theta
REAL(rp), DIMENSION(:), ALLOCATABLE :: zeta
REAL(rp), DIMENSION(:), ALLOCATABLE :: r ! temporary vars
REAL(rp) :: sigma
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)
sigma = 1.0_rp/SQRT(2.0_rp*(spp%falloff_rate/params%cpp%length))
sigma = sigma/params%cpp%length
r = sigma*SQRT(-2.0_rp*LOG(1.0_rp - (1.0_rp - &
EXP(-0.5_rp*spp%r_outter**2/sigma**2))*r))
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 gaussian_torus