%=======================================================================
% Driver script: 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).
%=======================================================================

clear;
clc;

number_of_data = 31;
x_test = [0.00, 0.20, 0.70, 1.25, 1.95, 2.75, 3.00];

x_data = zeros(1,number_of_data);
y_data = zeros(1,number_of_data);

%-----------------------------------------------------------------------
% 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.
%-----------------------------------------------------------------------
for i = 1:number_of_data
    fraction = (i-1)/(number_of_data-1);
    x_data(i) = 3.0*fraction^1.15;
    y_data(i) = exp(x_data(i));
end

%-----------------------------------------------------------------------
% Step 2: Calculate the derivatives at the requested points.
%-----------------------------------------------------------------------
fprintf('Numerical derivatives from tabulated data\n');
fprintf('Local four-point Lagrange interpolation\n\n');
fprintf(['        x      exact f''          numerical f''', ...
         '      exact f''''         numerical f''''\n']);

for i = 1:length(x_test)
    x_eval = x_test(i);

    [first_derivative,status_first] = ...
        derivative_from_data(x_eval,x_data,y_data,1);

    [second_derivative,status_second] = ...
        derivative_from_data(x_eval,x_data,y_data,2);

    if status_first ~= 0 || status_second ~= 0
        fprintf(['Derivative could not be evaluated at x = %.3f; ', ...
                 'status = %d %d\n'], ...
                x_eval,status_first,status_second);
        continue
    end

    exact = exp(x_eval);
    fprintf('%9.3f%18.8e%18.8e%18.8e%18.8e\n', ...
            x_eval,exact,first_derivative,exact,second_derivative);
end


function [derivative,status] = ...
    derivative_from_data(x_eval,x_data,y_data,derivative_order)
%=======================================================================
% 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
%=======================================================================

stencil_size = 4;
derivative = 0.0;
status = 0;
number_of_points = length(x_data);

% Check the input before selecting the interpolation stencil.
if length(y_data) ~= number_of_points
    status = 1;
    return
end

if number_of_points < stencil_size
    status = 2;
    return
end

if derivative_order < 1 || derivative_order > 2
    status = 3;
    return
end

if any(diff(x_data) <= 0)
    status = 4;
    return
end

if x_eval < x_data(1) || x_eval > x_data(end)
    status = 5;
    return
end

% 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,stencil_size);
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


function first = select_stencil(x_eval,x_data,stencil_size)
%=======================================================================
% 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.
%=======================================================================

number_of_points = length(x_data);
left = 1;
right = number_of_points;

while right > left+1
    middle = floor((left+right)/2);

    if x_eval < x_data(middle)
        right = middle;
    else
        left = middle;
    end
end

first = left-floor((stencil_size-2)/2);
first = max(1,min(first,number_of_points-stencil_size+1));
end


function value = ...
    lagrange_derivative(x_eval,x,y,derivative_order)
%=======================================================================
% 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.
%=======================================================================

number_of_points = length(x);
value = 0.0;

for j = 1:number_of_points
    % Denominator of the j-th Lagrange basis polynomial.
    denominator = 1.0;
    for k = 1:number_of_points
        if k ~= j
            denominator = denominator*(x(j)-x(k));
        end
    end

    if derivative_order == 1
        % In each term of L'_j(x), one factor is differentiated.
        for k = 1:number_of_points
            if k == j
                continue
            end

            term = 1.0;
            for i = 1:number_of_points
                if i ~= j && i ~= k
                    term = term*(x_eval-x(i));
                end
            end
            value = value+y(j)*term/denominator;
        end

    else
        % In each term of L''_j(x), two factors are differentiated.
        for k = 1:number_of_points
            if k == j
                continue
            end

            for ell = 1:number_of_points
                if ell == j || ell == k
                    continue
                end

                term = 1.0;
                for i = 1:number_of_points
                    if i ~= j && i ~= k && i ~= ell
                        term = term*(x_eval-x(i));
                    end
                end
                value = value+y(j)*term/denominator;
            end
        end
    end
end
end