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~ndmc ndmc proc~ndmc->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~ndmc program~fidasim->proc~fida_weights_mc program~fidasim->proc~fida_f program~fidasim->proc~neutron_mc proc~pfida_mc pfida_mc program~fidasim->proc~pfida_mc proc~pfida_f pfida_f program~fidasim->proc~pfida_f proc~neutron_f neutron_f program~fidasim->proc~neutron_f proc~fida_mc fida_mc program~fidasim->proc~fida_mc proc~pfida_mc->proc~gyro_correction proc~pfida_f->proc~gyro_correction 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, vpar, term1, term2

    if(inputs%flr.ge.1) then
        uvw = fields%uvw
        R = sqrt(uvw(1)**2 + uvw(2)**2)
        phi = atan2(uvw(2),uvw(1))
        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, in beam coordinates

        if(inputs%flr.ge.2) then
            !! convert the r_gyro vector to machine coordiantes
            if(fields%coords.eq.0) then
                rg_uvw=  matmul(beam_grid%basis, r_gyro)
            endif
            if(fields%coords.eq.1) then
                rg_uvw = r_gyro
            endif

            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) = (1./R*fields%dbz_dphi-fields%dbt_dz)/fields%b_abs
            cuvrxb(2) = (fields%dbr_dz - fields%dbz_dr)/fields%b_abs
            cuvrxb(3) = (1.0/R*fields%bt + fields%dbt_dr - 1.0/R*fields%dbr_dphi)/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) = 1.0/R*(fields%br*fields%dbr_dphi + fields%bt * fields%dbt_dphi + fields%bz*fields%dbz_dphi)/&
                        fields%b_abs
            grad_B(3) = (fields%br*fields%dbr_dz + fields%bt * fields%dbt_dz + fields%bz*fields%dbz_dz)/&
                        fields%b_abs
            !convert rg_uvw vector to cylindrical coordiantes
            rg_rtz(1) = rg_uvw(1)*cos(phi) + rg_uvw(2)*sin(phi)
            rg_rtz(2) = -rg_uvw(1)*sin(phi) + rg_uvw(2)*cos(phi)
            rg_rtz(3) = rg_uvw(3)
            term2 = -1.0 / (2.0 * fields%b_abs)*dot_product(rg_rtz,grad_B)
        else
            term1 = 0.0
            term2 = 0.0
        endif

        r_gyro = r_gyro * (1.0 - term1 - term2)
        if ((1.0 - term1 - term2 .le. 0.0) .or. (1.0 - term1 - term2 .ge. 2.0) ) then
            write(*,*) 'GYRO_STEP: Gyro correction results in negative distances or too large shift: ', &
                          1.0-term1-term2
            stop
        endif
    else
        r_gyro = 0.d0
    endif

end subroutine gyro_step