function [ xm, fm, calls ] = zero_chandrupatla ( f, x1, x2 )

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

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

  t = 0.5;

  while ( true )

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

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

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

    if ( fl < ph && ph < fh )
      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
%
%  Adjust T away from the interval boundary.
%
    t = max ( t, tl );
    t = min ( t, 1.0 - tl );

  end

  return
end

