torus Subroutine

subroutine torus(params,spp)
  !! @note Subrotuine for generating a uniform torus/torus 
  !! shell as the initial spatial condition of a given species 
  !! of particles in the simulation.@endnote
  !! This distribution is generated using the Inverse Transform 
  !! Sampling method. This distribution follows the same radial 
  !! distribution of a uniform disk/ring distribution, see the 
  !! documentation of the [[disk]] subroutine.
  TYPE(KORC_PARAMS), INTENT(IN) 	:: params
    !! Core KORC simulation parameters.
  TYPE(SPECIES), INTENT(INOUT) 		:: spp
    !! An instance of the derived type SPECIES 
    !! containing all the parameters and simulation variables of the 
    !! different species in the simulation.
  REAL(rp), DIMENSION(:), ALLOCATABLE 	:: r
    !! Radial position of the particles \(r\).
  REAL(rp), DIMENSION(:), ALLOCATABLE 	:: theta
  !! Uniform deviates in the range \([0,2\pi]\) 
  !! representing the uniform poloidal angle \(\theta\)
  !! distribution of the particles.
  REAL(rp), DIMENSION(:), ALLOCATABLE 	:: zeta
  !! Uniform deviates in the range \([0,2\pi]\) representing 
  !! the uniform toroidal angle \(\zeta\) distribution of the particles.
  INTEGER,DIMENSION(33) :: seed=(/1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/)

  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)

  if (.not.params%SameRandSeed) then
     call init_random_seed()
  else
     call random_seed(put=seed)
  end if
  call RANDOM_NUMBER(theta)
  theta = 2.0_rp*C_PI*theta

  if (.not.params%SameRandSeed) then
     call init_random_seed()
  else
     call random_seed(put=seed)
  end if
  call RANDOM_NUMBER(zeta)
  zeta = 2.0_rp*C_PI*zeta

  ! Uniform distribution on a disk at a fixed azimuthal theta
  if (.not.params%SameRandSeed) then
     call init_random_seed()
  else
     call random_seed(put=seed)
  end if
  call RANDOM_NUMBER(r)

  r = SQRT((spp%r_outter**2 - spp%r_inner**2)*r + spp%r_inner**2)
  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 torus