gyro_step Subroutine

  1. fidasim.f90
  2. libfida
  3. gyro_step

public subroutine gyro_step(vi, fields, r_gyro)

Calculates gyro-step

References

  • Belova, E. V., N. N. Gorelenkov, and C. Z. Cheng. "Self-consistent equilibrium model of low aspect-ratio toroidal plasma with energetic beam ions." Physics of Plasmas (1994-present) 10.8 (2003): 3240-3251. Appendix A: Last equation

Arguments

Type IntentOptional AttributesName
real(kind=Float64), intent(in), dimension(3):: vi

Ion velocity
type(LocalEMFields), intent(in) :: fields

Electro-magnetic fields
real(kind=Float64), intent(out), dimension(3):: r_gyro

Gyro-step Gyro-radius vector from particle position to guiding center

Calls

proc~~gyro_step~~CallsGraph proc~gyro_step gyro_step proc~cross_product cross_product proc~gyro_step->proc~cross_product

Called by

proc~~gyro_step~~CalledByGraph proc~gyro_step gyro_step proc~gyro_correction gyro_correction proc~gyro_correction->proc~gyro_step proc~npa_weights npa_weights proc~npa_weights->proc~gyro_step proc~gyro_radius gyro_radius proc~gyro_radius->proc~gyro_step proc~fida_weights_mc fida_weights_mc proc~fida_weights_mc->proc~gyro_correction proc~fida_f fida_f proc~fida_f->proc~gyro_correction proc~neutron_mc neutron_mc proc~neutron_mc->proc~gyro_correction program~fidasim fidasim program~fidasim->proc~npa_weights program~fidasim->proc~fida_weights_mc program~fidasim->proc~fida_f program~fidasim->proc~neutron_mc proc~neutron_f neutron_f program~fidasim->proc~neutron_f proc~fida_mc fida_mc program~fidasim->proc~fida_mc proc~neutron_f->proc~gyro_correction proc~fida_mc->proc~gyro_correction

Contents

Source Code

gyro_step

Source Code

subroutine gyro_step(vi, fields, r_gyro)
    !+ Calculates gyro-step
    !+
    !+###References
    !+* Belova, E. V., N. N. Gorelenkov, and C. Z. Cheng. "Self-consistent equilibrium model of low aspect-ratio toroidal plasma with energetic beam ions." Physics of Plasmas (1994-present) 10.8 (2003): 3240-3251. Appendix A: Last equation
    real(Float64), dimension(3), intent(in)  :: vi
        !+ Ion velocity
    type(LocalEMFields), intent(in)          :: fields
        !+ Electro-magnetic fields
    real(Float64), dimension(3), intent(out) :: r_gyro
        !+ Gyro-step
        !+ Gyro-radius vector from particle position to guiding center

    real(Float64), dimension(3) :: vxB, rg_uvw, uvw, cuvrxb, b_rtz, grad_B, rg_rtz
    real(Float64) :: one_over_omega, phi, R, rg_r, vpar, term1, term2

    if(inputs%no_flr.eq.0) then
        one_over_omega=inputs%ab*mass_u/(fields%b_abs*e0)
        vxB = cross_product(vi,fields%b_norm)
        vpar =  dot_product(vi,fields%b_norm)
        r_gyro = vxB*one_over_omega !points towards gyrocenter

        !! Second order correction approximation derived from
        !! Belova, E. V., N. N. Gorelenkov, and C. Z. Cheng.
        !! "Self-consistent equilibrium model of low aspect-ratio
        !! toroidal plasma with energetic beam ions."
        !! Physics of Plasmas (1994-present) 10.8 (2003): 3240-3251.
        !! Appendix A: Last equation
        uvw = fields%uvw
        R = sqrt(uvw(1)**2 + uvw(2)**2)
        phi = atan2(uvw(2),uvw(1))
        if(fields%machine_coords) then
            rg_uvw = r_gyro
        else
            rg_uvw = matmul(beam_grid%basis,r_gyro)
        endif
        rg_r = rg_uvw(1)*cos(phi) + rg_uvw(2)*sin(phi)
        b_rtz(1) = fields%br/fields%b_abs
        b_rtz(2) = fields%bt/fields%b_abs
        b_rtz(3) = fields%bz/fields%b_abs
        cuvrxb(1) = -fields%dbt_dz/fields%b_abs
        cuvrxb(2) = (fields%dbr_dz - fields%dbz_dr)/fields%b_abs
        cuvrxb(3) = fields%dbt_dr/fields%b_abs
        term1 = vpar*one_over_omega*dot_product(b_rtz,cuvrxb)
        grad_B(1) = (fields%br*fields%dbr_dr + fields%bt * fields%dbt_dr + fields%bz*fields%dbz_dr)/&
                    fields%b_abs
        grad_B(2) = 0.0
        grad_B(3) = (fields%br*fields%dbr_dz + fields%bt * fields%dbt_dz + fields%bz*fields%dbz_dz)/&
                    fields%b_abs
        rg_rtz(1) = rg_uvw(1)*cos(phi) + rg_uvw(2)*sin(phi)
        rg_rtz(2) = 0.0
        rg_rtz(3) = rg_uvw(3)
        term2 = -1.0 / (2.0 * fields%b_abs)*dot_product(rg_rtz,grad_B)
        r_gyro = r_gyro * (1.0 - term1 - term2)
        if (1.0 - term1 - term2 .le. 0.0) then
            write(*,*) 'GYRO_STEP: Gyro correction results in negative distances: ', &
                          1.0-term1-term2
            stop
        endif
    else
        r_gyro = 0.d0
    endif

end subroutine gyro_step