Returns the distance to the magnetic axis as a function of rp (the normalised poloidal flux variable) and theta.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
class(abstract_eqfile_cart_geo_type), | intent(in) | :: | self | |||
real, | intent(in) | :: | r | |||
real, | intent(in) | :: | theta |
real function rfun(self, r, theta)
implicit none
class(abstract_eqfile_cart_geo_type), intent(in) :: self
real, intent(in) :: r, theta
integer, parameter :: jmax = 5, kmax = 5
real, parameter :: xerrsec = 1.e-9
integer :: i, j, k, imax
real :: fa, fb, bmult, thetroot, a, b, local_broot
local_broot = self%bound(theta) ; thetroot = theta
a = 0. ; b = local_broot ; bmult = 0.01
fa = self%psi(a, thetroot) - r ; fb = self%psi(b, thetroot) - r
! need to bracket root, if we don't then try to search for a different b which does bracket
if (fa * fb > 0.) then
outer: do k = 1, kmax !Each k increases the distance we scan by (bmult)
bmult = bmult * 2. ; imax = 10
do j = 1, jmax !Each j increases the number of steps (imax)
imax = imax * 2
do i = 1, imax !Scan from local_broot to local_broot*(1 + bmult) in imax steps
b = local_broot * (1 + bmult * float(i) / float(imax))
fb = self%psi(b, thetroot) - r
if (fa * fb <= 0.) exit outer !If a ; b brackets root skip ahead
end do
do i = 1, imax !Scan from local_broot to local_broot/(1 + bmult) in imax steps
b = local_broot / (1 + bmult * float(i) / float(imax))
fb = self%psi(b, thetroot) - r
if (fa * fb <= 0.) exit outer !If a ; b brackets root skip ahead
end do
end do
end do outer !No promise that we have found a bracket here, but zbrent will report this.
end if
rfun = self%zbrent(a, b, r, thetroot, xerrsec)
end function rfun