analytical_fields_p Subroutine

  subroutine analytical_fields_p(params,pchunk,F,X_X,X_Y,X_Z, &
       B_X,B_Y,B_Z,E_X,E_Y,E_Z,flag_cache)
    TYPE(KORC_PARAMS), INTENT(IN)                              :: params
    !! Core KORC simulation parameters.
    TYPE(FIELDS), INTENT(IN)                                   :: F
    INTEGER, INTENT(IN)  :: pchunk
    REAL(rp)      :: R0,B0,lam,q0,E0,ar
    REAL(rp),  INTENT(IN),DIMENSION(pchunk)      :: X_X,X_Y,X_Z
    REAL(rp),  INTENT(OUT),DIMENSION(pchunk)     :: B_X,B_Y,B_Z
    REAL(rp),  INTENT(OUT),DIMENSION(pchunk)     :: E_X,E_Y,E_Z
    INTEGER(is),  INTENT(INOUT),DIMENSION(pchunk)     :: flag_cache
    REAL(rp),DIMENSION(pchunk)     :: T_R,T_T,T_Z
    REAL(rp),DIMENSION(pchunk)                               :: Ezeta
    !! Toroidal electric field \(E_\zeta\).
    REAL(rp),DIMENSION(pchunk)                               :: Er
    !! changed YG Radial electric field \(E_\r\).
    REAL(rp),DIMENSION(pchunk)                               :: Bzeta
    !! Toroidal magnetic field \(B_\zeta\).
    REAL(rp),DIMENSION(pchunk)                              :: Bp
    !! Poloidal magnetic field \(B_\theta(r)\).
    REAL(rp),DIMENSION(pchunk)                               :: eta
    !! Aspect ratio \(\eta\).
    REAL(rp),DIMENSION(pchunk)                                :: q
    !! Safety profile \(q(r)\).
    REAL(rp),DIMENSION(pchunk)                             :: cT,sT,cZ,sZ
    INTEGER                                      :: cc
    !! Particle chunk iterator.
    REAL(rp) :: Er0,rrmn,sigmaamn
    
    B0=F%Bo
    E0=F%Eo
    lam=F%AB%lambda
    R0=F%AB%Ro
    q0=F%AB%qo
    ar=F%AB%a

    Er0=F%AB%Ero
    rrmn=F%AB%rmn
    sigmaamn=F%AB%sigmamn
    
    call cart_to_tor_check_if_confined_p(pchunk,ar,R0,X_X,X_Y,X_Z, &
         T_R,T_T,T_Z,flag_cache)

    !$OMP SIMD
    !    !$OMP& aligned(cT,sT,cZ,sZ,eta,q,Bp,Bzeta,B_X,B_Y,B_Z, &
    !    !$OMP& Ezeta,E_X,E_Y,E_Z,T_T,T_Z,T_R)
    do cc=1_idef,pchunk
       cT(cc)=cos(T_T(cc))
       sT(cc)=sin(T_T(cc))
       cZ(cc)=cos(T_Z(cc))
       sZ(cc)=sin(T_Z(cc))

       eta(cc) = T_R(cc)/R0
       q(cc) = q0*(1.0_rp + (T_R(cc)*T_R(cc)/(lam*lam)))
       Bp(cc) = -eta(cc)*B0/(q(cc)*(1.0_rp + eta(cc)*cT(cc)))
       !changed kappa  YG
       Bzeta(cc) = B0/( 1.0_rp + eta(cc)*cT(cc))


       B_X(cc) =  Bzeta(cc)*cZ(cc) - Bp(cc)*sT(cc)*sZ(cc)
       B_Y(cc) = -Bzeta(cc)*sZ(cc) - Bp(cc)*sT(cc)*cZ(cc)
       B_Z(cc) = Bp(cc)*cT(cc)

       !write(6,*) 'Ero ',Ero,'Er0 ',Er0
       !write(6,*) 'rmn ',rmn,'rrmn ',rrmn
       !write(6,*) 'sigmamn ',sigmamn,'sigmaamn ',sigmaamn
       !write(6,*) 'T_R ',T_R(cc)*params%cpp%length
       
       Ezeta(cc) = -E0/( 1.0_rp + eta(cc)*cT(cc))
       Er(cc) =Er0*(1/cosh((T_R(cc)-rrmn)/sigmaamn))

       E_X(cc) = Ezeta(cc)*cZ(cc)+Er(cc)*cT(cc)*sZ(cc)
       E_Y(cc) = -Ezeta(cc)*sZ(cc)+Er(cc)*cT(cc)*cZ(cc)
       E_Z(cc) = Er(cc)*sT(cc)

       !write(6,*) 'Er ',Er(cc)
    end do
    !$OMP END SIMD

  end subroutine analytical_fields_p