subroutine findb0(me, stps, bdel, value, bequ, rr0)
class(shellig_type), intent(inout) :: me
real(wp), intent(in) :: stps
real(wp), intent(inout) :: bdel
real(wp), intent(out) :: bequ
logical, intent(out) :: value
real(wp), intent(out) :: rr0
real(wp) :: b, bdelta, bmin, bold, bq1, &
bq2, bq3, p(8, 4), r1, r2, r3, &
rold, step, step12, zz
integer :: i, irun, j, n
step = stps
irun = 0
rold = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings
main: do
irun = irun + 1
if (irun > 5) then
value = .false.
exit main
end if
! first three points
p(1, 2) = me%sp(1)
p(2, 2) = me%sp(2)
p(3, 2) = me%sp(3)
step = -sign(step, p(3, 2))
call me%stoer(p(1, 2), bq2, r2)
p(1, 3) = p(1, 2) + 0.5_wp * step * p(4, 2)
p(2, 3) = p(2, 2) + 0.5_wp * step * p(5, 2)
p(3, 3) = p(3, 2) + 0.5_wp * step
call me%stoer(p(1, 3), bq3, r3)
p(1, 1) = p(1, 2) - step * (2.0_wp * p(4, 2) - p(4, 3))
p(2, 1) = p(2, 2) - step * (2.0_wp * p(5, 2) - p(5, 3))
p(3, 1) = p(3, 2) - step
call me%stoer(p(1, 1), bq1, r1)
p(1, 3) = p(1, 2) + step * (20.0_wp * p(4, 3) - 3.*p(4, 2) + p(4, 1)) / 18.0_wp
p(2, 3) = p(2, 2) + step * (20.0_wp * p(5, 3) - 3.*p(5, 2) + p(5, 1)) / 18.0_wp
p(3, 3) = p(3, 2) + step
call me%stoer(p(1, 3), bq3, r3)
! invert sense if required
if (bq3 > bq1) then
step = -step
r3 = r1
bq3 = bq1
do i = 1, 5
zz = p(i, 1)
p(i, 1) = p(i, 3)
p(i, 3) = zz
end do
end if
! initialization
step12 = step / 12.0_wp
value = .true.
bmin = 1.0e4_wp
bold = 1.0e4_wp
! corrector (field line tracing)
n = 0
corrector: do
p(1, 3) = p(1, 2) + step12 * (5.0_wp * p(4, 3) + 8.0_wp * p(4, 2) - p(4, 1))
n = n + 1
p(2, 3) = p(2, 2) + step12 * (5.0_wp * p(5, 3) + 8.0_wp * p(5, 2) - p(5, 1))
! predictor (field line tracing)
p(1, 4) = p(1, 3) + step12 * (23.0_wp * p(4, 3) - 16.0_wp * p(4, 2) + 5.0_wp * p(4, 1))
p(2, 4) = p(2, 3) + step12 * (23.0_wp * p(5, 3) - 16.0_wp * p(5, 2) + 5.0_wp * p(5, 1))
p(3, 4) = p(3, 3) + step
call me%stoer(p(1, 4), bq3, r3)
do j = 1, 3
do i = 1, 8
p(i, j) = p(i, j + 1)
end do
end do
b = sqrt(bq3)
if (b < bmin) bmin = b
if (b > bold) exit corrector
bold = b
rold = 1.0_wp / r3
me%sp(1) = p(1, 4)
me%sp(2) = p(2, 4)
me%sp(3) = p(3, 4)
end do corrector
if (bold /= bmin) value = .false.
bdelta = (b - bold) / bold
if (bdelta <= bdel) exit main
step = step / 10.0_wp
end do main
rr0 = rold
bequ = bold
bdel = bdelta
end subroutine findb0