When calculating gyrofrequencies always the electron charge is used, not considering the charge of the ion. This can be an issue if we want to use FIDASIM in the future for helium (not likely because CX cross-section would be pretty small in the Fusion energy range of alpha particle, but maybe needed to analyse some data from AUG, that can use He beams)
Affected lines
omega= (fields%b_abs*e0)/(Ai*mass_u), in the routine given below
subroutine gyro_surface(fields, energy, pitch, Ai, gs)
!+ Calculates the surface of all possible trajectories
type(LocalEMFields), intent(in) :: fields
!+ Electromagnetic fields at guiding center
real(Float64), intent(in) :: energy
!+ Energy of particle
real(Float64), intent(in) :: pitch
!+ Particle pitch w.r.t the magnetic field
real(Float64), intent(in) :: Ai
!+ Ion atomic mass
type(GyroSurface), intent(out) :: gs
!+ Gyro-surface
integer :: i
real(Float64) :: alpha, vabs, omega
real(Float64), dimension(3,3) :: s
vabs = sqrt(energy/(v2_to_E_per_amu*Ai))
omega= (fields%b_abs*e0)/(Ai*mass_u)
alpha = vabs/omega
gs%omega = omega
gs%v = vabs
gs%axes(1) = alpha*sqrt(1-pitch**2)
gs%axes(2) = alpha*sqrt(1-pitch**2)
gs%axes(3) = pitch*alpha
s = 0.d0
s(1,1) = gs%axes(1)**(-2)
s(2,2) = gs%axes(2)**(-2)
s(3,3) = -gs%axes(3)**(-2)
gs%center = fields%pos
gs%basis(:,1) = fields%a_norm
gs%basis(:,2) = fields%c_norm
gs%basis(:,3) = fields%b_norm
gs%A = matmul(gs%basis,matmul(s,transpose(gs%basis)))
end subroutine gyro_surface
- `one_over_omega=Ai*mass_u/(fields%b_abs*e0)` in the routine below:
```fortran
subroutine gyro_step(vi, fields, Ai, 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), intent(in) :: Ai
!+ Atomic mass of ion
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=Ai*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
Proposed Solution
Include a namelist parameter, qFI, default to 1.0d0 to handle this
Issue
When calculating gyrofrequencies always the electron charge is used, not considering the charge of the ion. This can be an issue if we want to use FIDASIM in the future for helium (not likely because CX cross-section would be pretty small in the Fusion energy range of alpha particle, but maybe needed to analyse some data from AUG, that can use He beams)
Affected lines
omega= (fields%b_abs*e0)/(Ai*mass_u), in the routine given belowend subroutine gyro_surface
Proposed Solution