cart_to_tor_check_if_confined Subroutine

public subroutine cart_to_tor_check_if_confined(X, F, Xtor, flag)

The coordinate transformation is given by:

where is the radial position of the magnetic axis.

Arguments

Type IntentOptional AttributesName
real(kind=rp), intent(in), DIMENSION(:,:), ALLOCATABLE:: X

Particles' position in Cartesian coordinates. X(1,:) = , X(2,:) = , X(3,:) =

type(FIELDS), intent(in) :: F

An instance of the KORC derived type FIELDS.

real(kind=rp), intent(inout), DIMENSION(:,:), ALLOCATABLE:: Xtor

Particles' position in cylindrical coordinates. Xtor(1,:) = , Xtor(2,:) = , Xtor(3,:) =

integer(kind=is), intent(inout), DIMENSION(:), ALLOCATABLE:: flag

Called by

proc~~cart_to_tor_check_if_confined~~CalledByGraph proc~cart_to_tor_check_if_confined cart_to_tor_check_if_confined proc~get_analytical_fields get_analytical_fields proc~get_analytical_fields->proc~cart_to_tor_check_if_confined proc~get_profiles get_profiles proc~get_profiles->proc~cart_to_tor_check_if_confined proc~avalanche_4d Avalanche_4D proc~avalanche_4d->proc~get_profiles proc~get_fields get_fields proc~get_fields->proc~get_analytical_fields proc~mh_psi MH_psi proc~mh_psi->proc~get_fields proc~get_avalanche_4d get_Avalanche_4D proc~get_avalanche_4d->proc~avalanche_4d proc~sample_hollmann_distribution_3d_psi sample_Hollmann_distribution_3D_psi proc~sample_hollmann_distribution_3d_psi->proc~get_fields proc~unitvectors unitVectors proc~unitvectors->proc~get_fields proc~get_hollmann_distribution_3d_psi get_Hollmann_distribution_3D_psi proc~get_hollmann_distribution_3d_psi->proc~sample_hollmann_distribution_3d_psi proc~gyro_distribution gyro_distribution proc~gyro_distribution->proc~unitvectors

Contents


Source Code

  subroutine cart_to_tor_check_if_confined(X,F,Xtor,flag)
    !! @note Subroutine that converts the position of simulated particles
    !! from Cartesian \((x,y,z)\) to toroidal \((r,\theta, \zeta)\) coordinates.
    !! In addition to performing the coordinate transformation, this
    !! subroutine checks whether a given particle is within the plasma or not.
    !! A particle is not longer considered to be within the plasma if its
    !! minor radius \(r > r_{edge}\), where \(r_{edge}\) is the radial
    !! distance to the plasma edge as measured from the magnetic axis. For
    !! more details see the analytical model of the magnetic field in
    !! [[korc_types]] and [[korc_fields]].
    !!
    !! The coordinate transformation is given by:
    !!
    !! $$r = \sqrt{ \left[\sqrt{x^2 + y^2}-R_0\right]^2 + z^2 },$$
    !! $$\theta = \arctan{\left( \frac{z}{\sqrt{x^2 + y^2}-Ro} \right)}.$$
    !! $$\zeta = \arctan{\left( \frac{x}{y} \right)},$$
    !!
    !! where \(R_0\) is the radial position of the magnetic axis.
    REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN)      :: X
    !! Particles' position in Cartesian coordinates. X(1,:) = \(x\), X(2,:)
    !! = \(y\), X(3,:) = \(z\)
    TYPE(FIELDS), INTENT(IN)                               :: F
    !! An instance of the KORC derived type FIELDS.
    REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT)   :: Xtor
    !! Particles' position in cylindrical coordinates. Xtor(1,:) =
    !! \(r\), Xtor(2,:) = \(\theta\), Xtor(3,:) = \(\zeta\)
    INTEGER(is), DIMENSION(:), ALLOCATABLE, INTENT(INOUT)  :: flag
    REAL(rp)                                               :: a
    !! Distance to plasma edge as measured from the magnetic axis.
    REAL(rp)                                               :: Ro
    !! Radial position of the magnetic axis.
    INTEGER                                                :: pp
    !! Iterator.
    INTEGER                                                :: ss
    !! Iterator.

    if (size(X,1).eq.1) then
       ss = size(X,1)
    else
       if (X(2,1).eq.0) then
          ss=1_idef
       else
          ss = size(X,1)
       end if
    endif
    
    a = F%AB%a
    Ro = F%AB%Ro

!    write(output_unit_write,'("X c2tor: ",E17.10)') X(1,:)
    
    !$OMP PARALLEL DO FIRSTPRIVATE(ss,a,Ro) PRIVATE(pp) SHARED(X,Xtor,flag)
    do pp=1_idef,ss
       if ( flag(pp) .EQ. 1_is ) then
          Xtor(pp,1) = SQRT( (SQRT(X(pp,1)**2 + X(pp,2)**2) - Ro)**2 + &
               X(pp,3)**2 )
          Xtor(pp,2) = ATAN2(X(pp,3), SQRT(X(pp,1)**2 + X(pp,2)**2) - Ro)
          Xtor(pp,2) = MODULO(Xtor(pp,2),2.0_rp*C_PI)
          Xtor(pp,3) = ATAN2(X(pp,1),X(pp,2))
          Xtor(pp,3) = MODULO(Xtor(pp,3),2.0_rp*C_PI)

!          write(output_unit_write,'("r: ",E17.10)') Xtor(1,1)
!          write(output_unit_write,'("a: ",E17.10)') a
!          write(output_unit_write,'("Ro: ",E17.10)') Ro
          
          if (Xtor(pp,1) .GT. F%AB%a) then
             flag(pp) = 0_is

!             stop 'error in dist init'
          end if
       end if
    end do
    !$OMP END PARALLEL DO
  end subroutine cart_to_tor_check_if_confined