!=======================================================================
! Driver program: Calculates first and second derivatives from
! tabulated data using local four-point Lagrange interpolation.
!
! The data points may be unequally spaced. Near either endpoint, the
! four-point stencil is shifted automatically so that it remains inside
! the tabulated interval.
!
! Based on deriv_Lg.f90 by Alexander Godunov (October 2009).
! Revised for the companion website (September 2026).
!=======================================================================
program ch06_diff_data
  use iso_fortran_env, only : real64
  implicit none

  integer, parameter :: number_of_data  = 31
  integer, parameter :: number_of_tests = 7
  integer, parameter :: stencil_size    = 4

  real(real64), parameter :: x_test(number_of_tests) = [ &
       0.00_real64, 0.20_real64, 0.70_real64, 1.25_real64, &
       1.95_real64, 2.75_real64, 3.00_real64 ]

  real(real64) :: x_data(number_of_data)
  real(real64) :: y_data(number_of_data)
  real(real64) :: first_derivative, second_derivative
  real(real64) :: exact, fraction
  integer :: i, status_first, status_second

  !---------------------------------------------------------------------
  ! Step 1: Generate sample tabulated data.
  !
  ! The example uses f(x)=exp(x) on a mildly nonuniform grid. Both exact
  ! derivatives equal exp(x), so the numerical results are easy to check.
  !
  ! For measured or previously calculated data, replace this loop with
  ! assignments to x_data and y_data or with suitable file-input code.
  !---------------------------------------------------------------------
  do i = 1, number_of_data
    fraction = real(i - 1, real64)/real(number_of_data - 1, real64)
    x_data(i) = 3.0_real64*fraction**1.15_real64
    y_data(i) = exp(x_data(i))
  end do

  !---------------------------------------------------------------------
  ! Step 2: Calculate the derivatives at the requested points.
  !---------------------------------------------------------------------
  write (*,'(a)') 'Numerical derivatives from tabulated data'
  write (*,'(a)') 'Local four-point Lagrange interpolation'
  write (*,'(a)') ''
  write (*,'(a)') &
       "       x       exact f'          numerical f'      exact f''         numerical f''"

  do i = 1, number_of_tests
    call derivative_from_data(x_test(i), x_data, y_data, 1, &
                              first_derivative, status_first)

    call derivative_from_data(x_test(i), x_data, y_data, 2, &
                              second_derivative, status_second)

    if (status_first /= 0 .or. status_second /= 0) then
      write (*,'(a,f8.3,a,2(i0,1x))') &
           'Derivative could not be evaluated at x = ', x_test(i), &
           '; status = ', status_first, status_second
      cycle
    end if

    exact = exp(x_test(i))
    write (*,'(f9.3,4es18.8)') x_test(i), exact, first_derivative, &
                               exact, second_derivative
  end do

contains

  subroutine derivative_from_data(x_eval, x_data, y_data, derivative_order, &
                                  derivative, status)
    !-------------------------------------------------------------------
    ! Estimate a first or second derivative from tabulated data.
    !
    ! A cubic Lagrange polynomial is constructed from four consecutive
    ! data points near x_eval. Differentiating this polynomial gives the
    ! requested derivative.
    !
    ! Input
    !   x_eval            point where the derivative is required
    !   x_data(:)         strictly increasing data abscissas
    !   y_data(:)         tabulated function values
    !   derivative_order  1 for the first derivative, 2 for the second
    !
    ! Output
    !   derivative        numerical value of the requested derivative
    !   status            0: success
    !                     1: x_data and y_data have different sizes
    !                     2: fewer than four data points were supplied
    !                     3: derivative_order is not 1 or 2
    !                     4: x_data is not strictly increasing
    !                     5: x_eval is outside the tabulated interval
    !-------------------------------------------------------------------
    real(real64), intent(in)  :: x_eval
    real(real64), intent(in)  :: x_data(:), y_data(:)
    integer,      intent(in)  :: derivative_order
    real(real64), intent(out) :: derivative
    integer,      intent(out) :: status

    real(real64) :: x_local(stencil_size)
    real(real64) :: y_local(stencil_size)
    integer :: first, i, n

    derivative = 0.0_real64
    status = 0
    n = size(x_data)

    ! Check the input before selecting the interpolation stencil.
    if (size(y_data) /= n) then
      status = 1
      return
    end if

    if (n < stencil_size) then
      status = 2
      return
    end if

    if (derivative_order < 1 .or. derivative_order > 2) then
      status = 3
      return
    end if

    do i = 1, n - 1
      if (x_data(i + 1) <= x_data(i)) then
        status = 4
        return
      end if
    end do

    if (x_eval < x_data(1) .or. x_eval > x_data(n)) then
      status = 5
      return
    end if

    ! Select four nearby points. At a boundary, the stencil is shifted
    ! automatically so that all four points remain in the data range.
    first = select_stencil(x_eval, x_data)
    x_local = x_data(first:first + stencil_size - 1)
    y_local = y_data(first:first + stencil_size - 1)

    derivative = lagrange_derivative(x_eval, x_local, y_local, &
                                     derivative_order)
  end subroutine derivative_from_data


  integer function select_stencil(x_eval, x_data) result(first)
    !-------------------------------------------------------------------
    ! Locate the first point of the local four-point stencil.
    !
    ! The binary search finds left such that
    !
    !       x_data(left) <= x_eval < x_data(left+1).
    !
    ! The stencil is centered on this interval as closely as possible.
    !-------------------------------------------------------------------
    real(real64), intent(in) :: x_eval
    real(real64), intent(in) :: x_data(:)

    integer :: left, right, middle, n

    n = size(x_data)
    left = 1
    right = n

    do while (right > left + 1)
      middle = (left + right)/2
      if (x_eval < x_data(middle)) then
        right = middle
      else
        left = middle
      end if
    end do

    first = left - (stencil_size - 2)/2
    first = max(1, min(first, n - stencil_size + 1))
  end function select_stencil


  real(real64) function lagrange_derivative(x_eval, x, y, &
                                             derivative_order) result(value)
    !-------------------------------------------------------------------
    ! Differentiate the local Lagrange interpolating polynomial.
    !
    ! For the local data (x_j,y_j),
    !
    !       P(x) = sum_j y_j L_j(x),
    !
    ! where
    !
    !       L_j(x) = product_{k /= j} (x-x_k)/(x_j-x_k).
    !
    ! The loops evaluate the first or second derivative of each basis
    ! polynomial L_j directly. No constant grid spacing is assumed.
    !-------------------------------------------------------------------
    real(real64), intent(in) :: x_eval
    real(real64), intent(in) :: x(:), y(:)
    integer,      intent(in) :: derivative_order

    real(real64) :: denominator, term
    integer :: i, j, k, l, n

    n = size(x)
    value = 0.0_real64

    do j = 1, n
      ! Denominator of the j-th Lagrange basis polynomial.
      denominator = 1.0_real64
      do k = 1, n
        if (k /= j) denominator = denominator*(x(j) - x(k))
      end do

      if (derivative_order == 1) then
        ! In each term of L'_j(x), one factor has been differentiated.
        do k = 1, n
          if (k == j) cycle
          term = 1.0_real64
          do i = 1, n
            if (i /= j .and. i /= k) term = term*(x_eval - x(i))
          end do
          value = value + y(j)*term/denominator
        end do

      else
        ! In each term of L''_j(x), two factors have been differentiated.
        do k = 1, n
          if (k == j) cycle
          do l = 1, n
            if (l == j .or. l == k) cycle
            term = 1.0_real64
            do i = 1, n
              if (i /= j .and. i /= k .and. i /= l) then
                term = term*(x_eval - x(i))
              end if
            end do
            value = value + y(j)*term/denominator
          end do
        end do
      end if
    end do
  end function lagrange_derivative

end program ch06_diff_data