General calculation of the poloidal magnetic field radius of curvature (see section 2.3 of "GS2 analytic geometry specification") ! JRB
we actually calculate 1/Rpol as this is what we actually use.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real, | intent(in), | dimension(-ntgrid:) | :: | theta | ||
| real, | intent(in), | dimension(-ntgrid:) | :: | th_bish | ||
| real, | intent(in), | dimension(-ntgrid:) | :: | ltheta | ||
| real, | intent(out), | dimension(-ntgrid:) | :: | invRpol | ||
| integer, | intent(in) | :: | nth |
pure subroutine inv_Rpol(theta, th_bish, ltheta, invRpol, nth)
use constants, only: twopi, pi
implicit none
integer, intent (in) :: nth
real, dimension(-ntgrid:), intent (in) :: th_bish, theta, ltheta
real, dimension(-ntgrid:), intent (out) :: invRpol
real, dimension(-ntgrid:ntgrid) :: dthdl
integer :: i
real :: dth_bish, dth
! Calculate du/dtheta (where u is defined according to Miller PoP 5 973 (1998))
! Note that we must properly account for discontinuities in u. u is th_bish
do i = -nth + 1, nth - 1
dth_bish = th_bish(i + 1) - th_bish(i - 1) ; dth = theta(i + 1) - theta(i - 1)
if ( dth_bish > pi) then
invRpol(i) = (dth_bish - twopi) / dth
else if(dth_bish < -pi) then
invRpol(i) = (dth_bish + twopi) / dth
else
invRpol(i) = dth_bish / dth
end if
end do
invRpol(nth) = (th_bish(-nth+1) - th_bish(nth-1)) / (theta(-nth+1) - theta(nth-1) + twopi)
invRpol(-nth) = invRpol(nth)
! Calculate dtheta/dl (where l is the poloidal arc length)
call gradl(ltheta, theta, dthdl, -twopi, nth)
! Scale by dthdl to complete calculation of 1/Rpol
! (the inverse radius of curvature of the poloidal magnetic field)
invRpol(-nth:nth) = invRpol(-nth:nth) * dthdl(-nth:nth)
end subroutine inv_Rpol