zbrent Function

  real function zbrent(self, x1, x2, rootval, thetroot, tol)
    use mp, only: mp_abort
    implicit none
    class(abstract_eqfile_cart_geo_type), intent(in) :: self
    real, intent (in) :: x1, x2, rootval, thetroot, tol
    real, parameter :: eps = 3.e-8, eps1 = 2.e-8
    integer, parameter :: itmax = 100
    real :: a, b, c, fa, fb, fc, d, e, tol1, xm, q, p, r, s
    integer :: iter
    zbrent = 0
    a = x1 ; b = x2
    fa = self%psi(a, thetroot) - rootval ; fb = self%psi(b, thetroot) - rootval
    if (fb * fa > 0.) then !If we're not bracketing the root
       ! check to see if fa or fb is close to zero already:
       if (abs(fa) < eps1) then
          zbrent = a ; return
       else if (abs(fb) < eps1) then
          zbrent = b ; return
       else
          write(*,*) a, b, fa, fb
          call mp_abort('root must be bracketed for zbrent.', .true.)
       end if
    end if

    c = 0. ; fc = fb
    do iter = 1, itmax
       if (fb * fc > 0.) then ! If root not bracketed by b and c swap a and c
          c = a ; fc = fa
          d = b - a ; e = d
       end if
       if (abs(fc) < abs(fb)) then !If c closer than b cyclic shuffle a<-b<-c<-a
          a = b ; b = c ; c = a
          fa = fb ; fb = fc ; fc = fa
       end if
       tol1 = 2. * eps * abs(b) + 0.5 * tol
       xm = 0.5 * (c - b) ! Bisect c and b
       if (abs(xm) <= tol1 .or. is_zero(fb)) then
          zbrent = b ; return
       end if
       if (abs(e) >= tol1 .and. abs(fa) > abs(fb)) then !e might not be set by this point?
          s = fb / fa
          if (exactly_equal(a, c)) then
             p = 2. * xm * s
             q = 1. - s
          else
             q = fa / fc
             r = fb / fc
             p = s * (2. * xm * q * (q - r) - (b - a) * (r - 1.))
             q = (q - 1.) * (r - 1.) * (s - 1.)
          end if
          if (p > 0.) q = -q
          p = abs(p)
          if (2. * p  <  min(3. * xm * q - abs(tol1 * q), abs(e * q))) then
             e = d
             d = p / q
          else
             d = xm
             e = d
          end if
       else
          d = xm
          e = d
       end if
       a = b ; fa = fb
       if (abs(d)  >  tol1) then
          b = b + d
       else
          b = b + sign(tol1, xm)
       end if
       fb = self%psi(b, thetroot) - rootval
    end do
    call mp_abort('zbrent exceeding maximum iterations.',.true.)
  end function zbrent