subroutine zero_chandrupatla ( f, x1, x2, xm, fm, calls )

c*********************************************************************72
c
cc zero_chandrupatla() seeks a zero of a function using Chandrupatla's algorithm.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    19 March 2024
c
c  Author:
c
c    Original QBASIC version by Tirupathi Chandrupatla.
c    This version by John Burkardt.
c
c  Reference:
c
c    Tirupathi Chandrupatla,
c    A new hybrid quadratic/bisection algorithm for finding the zero of a
c    nonlinear function without using derivatives,
c    Advances in Engineering Software,
c    Volume 28, Number 3, pages 145-149, 1997.
c
c  Input:
c
c    function f ( x ): the name of the user-supplied function.
c
c    double precision a, b: the endpoints of the change of sign interval.
c
c  Output:
c
c    double precision z, fz: the estimated root and its function value.
c
c    integer calls: the number of function calls.
c
      implicit none

      double precision a
      double precision al
      double precision b
      double precision c
      integer calls
      double precision d
      double precision delta
      double precision eps
      double precision, external :: f 
      double precision f0
      double precision f1
      double precision f2
      double precision f3
      double precision fh
      double precision fl
      double precision fm
      double precision ph
      double precision t
      double precision tl
      double precision tol
      double precision x0
      double precision x1
      double precision x2
      double precision x3
      double precision xi
      double precision xm

      eps = 1.0D-10
      delta = 0.00001

      f1 = f ( x1 )
      f2 = f ( x2 )
      calls = 2

      t = 0.5

      do while ( .true. )

        x0 = x1 + t * ( x2 - x1 )
        f0 = f ( x0 )
        calls = calls + 1
c
c  Arrange 2-1-3: 2-1 Interval, 1 Middle, 3 Discarded point.
c
        if ( ( 0.0 < f0 ) .eqv. ( 0.0 < f1 ) ) then
          x3 = x1
          f3 = f1
          x1 = x0
          f1 = f0
        else
          x3 = x2
          f3 = f2
          x2 = x1
          f2 = f1
          x1 = x0
          f1 = f0
        end if
c
c  Identify the one that approximates zero.
c
        if ( abs ( f2 ) < abs ( f1 ) ) then
          xm = x2
          fm = f2
        else
          xm = x1
          fm = f1
        end if

        tol = 2.0 * eps * abs ( xm ) + 0.5 * delta
        tl = tol / abs ( x2 - x1 )

        if ( 0.5 < tl .or. fm == 0.0 ) then
          exit
        end if
c
c  If inverse quadratic interpolation holds, use it.
c
        xi = ( x1 - x2 ) / ( x3 - x2 )
        ph = ( f1 - f2 ) / ( f3 - f2 )
        fl = 1.0 - sqrt ( 1.0 - xi )
        fh = sqrt ( xi )

        if ( fl < ph .and. ph < fh ) then
          al = ( x3 - x1 ) / ( x2 - x1 )
          a = f1 / ( f2 - f1 )
          b = f3 / ( f2 - f3 )
          c = f1 / ( f3 - f1 )
          d = f2 / ( f3 - f2 )
          t = a * b + c * d * al
        else
          t = 0.5
        end if
c
c  Adjust T away from the interval boundary.
c
        t = max ( t, tl )
        t = min ( t, 1.0 - tl )

      end do

      return
      end