function ch07_principal_value
%==========================================================================
% ch07_principal_value.m
%
% Cauchy principal value integral using analytical subtraction.
%
% The program evaluates
%
%                 1
%                 /
%                |       cos(x)
%      PV        |     --------- dx
%                |      x - 0.2
%                /
%                -1
%
% The singularity is removed analytically:
%
%      f(x)/(x-x0)
%        = [f(x)-f(x0)]/(x-x0) + f(x0)/(x-x0).
%
% Hence
%
%      PV integral
%        = integral [f(x)-f(x0)]/(x-x0) dx
%          + f(x0) log[(b-x0)/(x0-a)].
%
% The regular part is evaluated with Gauss-Legendre quadrature on
% [a,x0] and [x0,b].
%
% Based on the Fortran implementation by Alexander Godunov.
% MATLAB version prepared for the companion website, 2026.
%
% A plot of the integrand (or integration domain) is shown before the
% numerical calculation to make the character of the problem visible.
%==========================================================================

a = -1.0;
b = 1.0;
x0 = 0.2;

reference = -0.5912784964342436;
nvalues = [4, 8, 16, 32];

% === Plot the integrands ================================================
% Plot the original singular integrand on the two sides of the pole.
gap = 0.01;
xleft_plot  = linspace(a,x0-gap,1000);
xright_plot = linspace(x0+gap,b,1000);

figure('Name','Cauchy principal-value integrand')
plot(xleft_plot, f(xleft_plot)./(xleft_plot-x0), 'LineWidth',1.5)
hold on
plot(xright_plot,f(xright_plot)./(xright_plot-x0),'LineWidth',1.5)
hold off
xlabel('x')
ylabel('f(x)/(x-x_0)')
title('Original integrand with a pole at x_0 = 0.2')
grid on

% Analytical subtraction removes the pole.  Plot the regularized
% integrand that is actually integrated numerically.
xplot = linspace(a,b,2001);
f0plot = f(x0);
gplot = (f(xplot)-f0plot)./(xplot-x0);

% At x = x0 the removable limit is f''s derivative, -sin(x0).
mask = abs(xplot-x0) < 1.0e-12;
gplot(mask) = -sin(x0);

figure('Name','Regularized principal-value integrand')
plot(xplot,gplot,'LineWidth',1.5)
xlabel('x')
ylabel('[f(x)-f(x_0)]/(x-x_0)')
title('Regularized integrand after analytical subtraction')
grid on
% === End plots ==========================================================

fprintf('Cauchy principal value integral\n');
fprintf('PV integral of cos(x)/(x-0.2) from -1 to 1\n');
fprintf('Reference value      = %15.7e\n\n', reference);

fprintf('        n       Principal value\n');
fprintf('--------------------------------\n');

for n = nvalues
    result = principal_value(@f, a, b, x0, n);
    fprintf('%9d%21.7e\n', n, result);
end

end


function y = f(x)
%--------------------------------------------------------------------------
% Numerator function in PV integral f(x)/(x-x0).
%--------------------------------------------------------------------------
y = cos(x);
end


function result = principal_value(fun, a, b, x0, n)
%--------------------------------------------------------------------------
% Cauchy principal value using analytical subtraction.
%--------------------------------------------------------------------------

if x0 <= a || x0 >= b
    error('principal_value: x0 must lie inside (a,b).');
end

if n < 1
    error('principal_value: n must be positive.');
end

f0 = fun(x0);

left  = gauss_legendre_regular(fun, a,  x0, x0, f0, n);
right = gauss_legendre_regular(fun, x0, b,  x0, f0, n);

singular_part = f0*log((b-x0)/(x0-a));

result = left + right + singular_part;

end


function integral = gauss_legendre_regular(fun, a, b, x0, f0, n)
%--------------------------------------------------------------------------
% Gauss-Legendre quadrature for the regularized function
%
%             [f(x)-f(x0)]/(x-x0).
%
% In this application x0 is an endpoint of the subinterval and is
% therefore not sampled by Gauss-Legendre quadrature.
%--------------------------------------------------------------------------

tol = 1.0e-14;
max_iter = 100;

midpoint = (a+b)/2.0;
halfwidth = (b-a)/2.0;
integral = 0.0;

m = floor((n+1)/2);

for i = 1:m

    z = cos(pi*(i-0.25)/(n+0.5));
    converged = false;

    for iter = 1:max_iter
        [pn,pnm1] = legendre_pair(n,z);

        if n == 1
            dp = 1.0;
        else
            dp = n*(z*pn-pnm1)/(z*z-1.0);
        end

        zold = z;
        z = zold - pn/dp;

        if abs(z-zold) <= tol
            converged = true;
            break;
        end
    end

    if ~converged
        error('gauss_legendre_regular: Newton iteration failed.');
    end

    [pn,pnm1] = legendre_pair(n,z);

    if n == 1
        dp = 1.0;
    else
        dp = n*(z*pn-pnm1)/(z*z-1.0);
    end

    weight = 2.0/((1.0-z*z)*dp*dp);

    xleft  = midpoint-halfwidth*z;
    xright = midpoint+halfwidth*z;

    gleft  = (fun(xleft)-f0)/(xleft-x0);
    gright = (fun(xright)-f0)/(xright-x0);

    if abs(z) <= tol
        integral = integral + ...
            weight*(fun(midpoint)-f0)/(midpoint-x0);
    else
        integral = integral + weight*(gleft+gright);
    end
end

integral = halfwidth*integral;

end


function [pn,pnm1] = legendre_pair(n,z)
if n == 1
    pn = z;
    pnm1 = 1.0;
    return;
end

p0 = 1.0;
p1 = z;

for j = 2:n
    p2 = ((2.0*j-1.0)*z*p1 - (j-1.0)*p0)/j;
    p0 = p1;
    p1 = p2;
end

pn = p1;
pnm1 = p0;
end