exponential_torus Subroutine

private subroutine exponential_torus(params, spp)

@brief Subroutine that generates a exponentially decaying radial distribution of particles in a circular cross-section torus of major and minor radi $R_0$ and $r_0$, respectively. @details We generate this exponentially decaying radial distribution $f(r)$ following the same approach as in \ref korc_spatial_distribution.disk, but this time, the radial distribution is given by:

$$f(r) = \left\{ \begin{array}{ll} \frac{k_0^2}{4\pi^2 R_0}\frac{\exp{\left( -k_0 r \right)}}{1 - \exp{\left( -k_0r_0\right)}\left[ 1 + k_0 r_0 \right]} &\ r<r_0 \\ 0 & r>r_0 \end{array} \right.$$

The radial position of the particles $r$ is obtained using the Inverse Trasnform Sampling method, finding $r$ numerically through the Newton-Raphson method. First, we calculate the particles' radial distribution in a disk centered at $(R,Z) = (0,0)$. Then, we transfor to a new set of coordinates where the disk is centered at $(R,Z) = (R_0,Z_0)$. Finally, we generate the toroidal distribution by givin each particle a toroidal angle $\zeta$ which follows a uniform distribution in the interval $[0,2\pi]$.

@param[in] params Core KORC simulation parameters. @param[in,out] spp An instance of the derived type SPECIES containing all the parameters and simulation variables of the different species in the simulation. @param fl Variable used in the Newton-Raphson method for finding the radial position of each particle. @param fr Variable used in the Newton-Raphson method for finding the radial position of each particle. @param fm Variable used in the Newton-Raphson method for finding the radial position of each particle. @param rl Variable used in the Newton-Raphson method for finding the radial position of each particle. @param rr Variable used in the Newton-Raphson method for finding the radial position of each particle. @param rm Variable used in the Newton-Raphson method for finding the radial position of each particle. @param relerr Tolerance used to determine when to stop iterating the Newton-Raphson method for finding $r$. @param r Radial position of the particles $r$. @param theta Uniform deviates in the range $[0,2\pi]$ representing the uniform poloidal angle $\theta$ distribution of the particles. @param zeta Uniform deviates in the range $[0,2\pi]$ representing the uniform toroidal angle $\zeta$ distribution of the particles. @param pp Particle iterator.

Arguments

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