subroutine golden_section ( f, n, x_tol, a, b, it )

!*****************************************************************************80
!
!! golden_section() seeks a minimizer of a unimodal function in [a,b].
!
!  Discussion:
!
!    A unimodal function f(x) in [a,b] has a minimizer a <= c <= b such that
!    f(x) strictly decreases between a and c, and strictly increases between
!    c and b.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    23 February 2026
!
!  Author:
!
!    Original MATLAB code by J Nathan Kutz.
!    This version by John Burkardt.
!
!  Reference:
!
!    J Nathan Kutz,
!    Data Driven Modeling and Scientific Computation,
!    Oxford University Press, 2013,
!    ISBN: 978-0-19-966034-6.
!
!  Input:
!
!    f: the function to be minimized.
!
!    integer n: the maximum number of iterations allowed.
!    n must be at least 1.
!
!    real x_tol: a tolerance for the width of the x interval.
!
!    real a, b: the left and right endpoints.
!
!  Output:
!
!    real a, b: the current interval containing the minimizer.
!
!    integer it: the number of steps taken.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) a
  real ( kind = rk8 ) b
  real ( kind = rk8 ), external :: f
  real ( kind = rk8 ) f1
  real ( kind = rk8 ) f2
  real ( kind = rk8 ) g
  integer it
  integer n
  real ( kind = rk8 ) x1
  real ( kind = rk8 ) x2
  real ( kind = rk8 ) x_tol

  do it = 0, n

    if ( it == 0 ) then

      g = ( - 1.0D+00 + sqrt ( 5.0D+00 ) ) / 2.0D+00
      x1 =             g   * a + ( 1.0D+00 - g ) * b
      x2 = ( 1.0D+00 - g ) * a +             g   * b
      f1 = f ( x1 )
      f2 = f ( x2 )

    else
 
      if ( f1 < f2 ) then

        b = x2
        x2 = x1
        f2 = f1
        x1 = g * a + ( 1.0D+00 - g ) * b
        f1 = f ( x1 )

      else

        a = x1
        x1 = x2
        f1 = f2
        x2 = ( 1.0D+00 - g ) * a + g * b
        f2 = f ( x2 )

      end if

    end if

    if ( abs ( b - a ) <= x_tol ) then
      return
    end if

  end do

  return
end