subroutine gaussian_elliptic_torus(params,spp)
TYPE(KORC_PARAMS), INTENT(IN) :: params
TYPE(SPECIES), INTENT(INOUT) :: spp
REAL(rp), DIMENSION(:), ALLOCATABLE :: rotation_angle
REAL(rp), DIMENSION(:), ALLOCATABLE :: r
REAL(rp), DIMENSION(:), ALLOCATABLE :: theta
REAL(rp), DIMENSION(:), ALLOCATABLE :: zeta
REAL(rp), DIMENSION(:), ALLOCATABLE :: X
REAL(rp), DIMENSION(:), ALLOCATABLE :: Y
REAL(rp), DIMENSION(:), ALLOCATABLE :: X1
REAL(rp), DIMENSION(:), ALLOCATABLE :: Y1
REAL(rp) :: sigma
INTEGER :: pp
ALLOCATE(X1(spp%ppp))
ALLOCATE(Y1(spp%ppp))
ALLOCATE(X(spp%ppp))
ALLOCATE(Y(spp%ppp))
ALLOCATE( rotation_angle(spp%ppp) )
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)
Y = r*SIN(theta)
X = r*COS(theta) + spp%shear_factor*Y
rotation_angle = 0.5_rp*C_PI - ATAN(1.0_rp,1.0_rp + spp%shear_factor);
X1 = X*COS(rotation_angle) - Y*SIN(rotation_angle) + spp%Ro
Y1 = X*SIN(rotation_angle) + Y*COS(rotation_angle) + spp%Zo
spp%vars%X(:,1) = X1*SIN(zeta)
spp%vars%X(:,2) = X1*COS(zeta)
spp%vars%X(:,3) = Y1
DEALLOCATE(X1)
DEALLOCATE(Y1)
DEALLOCATE(X)
DEALLOCATE(Y)
DEALLOCATE(rotation_angle)
DEALLOCATE(theta)
DEALLOCATE(zeta)
DEALLOCATE(r)
end subroutine gaussian_elliptic_torus