% ============================= %
% 
%   This script evaluates the 1D-PSWFs within arbitrary point in [-1,1].  
% 
%
% ============================= %
clear; close all; clc; 

set(0,'DefaultLineMarkerSize',  6);
set(0,'DefaultTextFontSize',    14);
set(0,'DefaultAxesFontSize',    12);
set(0,'DefaultLineLineWidth',   1.4);
set(0,'defaultfigurecolor',     'w');



%% Directly test the 1DPSWF function. 

L           = 10; 
b           = 0.3; 
a           = -0.3; 
c           = pi*(b-a)*L/2; 
saveFigs    = false; 

NyqNumberEst = ceil(2*c/pi);

nmax    = 40;               % highest Legendre order. 
h       = 0.01; 


[mu, NyquistNumber] = get1DPSWFEigenvalus(c, nmax, h); 




%% Self-tests. 

Pbar = getNormalizedLegendrePolynomials(nmax); 
x = -1+h/2:h:1-h/2;

% Construct A matrix. 
A = zeros(nmax+1); 
for ii = 1:nmax+1
    k = ii-1; 
    A(ii, ii) = k*(k+1)+(2*k*(k+1)-1)/((2*k+3)*(2*k-1))*c^2;
end

for ii = 1:nmax-1
    k = ii-1;
    A(ii, ii+2) = (k+2)*(k+1)/((2*k+3)*sqrt((2*k+1)*(2*k+5)))*c^2; 
end

for ii = 1:nmax-1
    k = ii-1; 
    A(ii+2, ii) = (k+2)*(k+1)/((2*k+3)*sqrt((2*k+1)*(2*k+5)))*c^2; 
end

Aeven = A(1:2:nmax+1, 1:2:nmax+1); 
Aodd = A(2:2:nmax+1, 2:2:nmax+1); 

% Assume the diagonals of Ds are in increasing order. 
[Veven, ~] = eig(Aeven); 
[Vodd, ~] = eig(Aodd); 



% Compute the first 5 1D-PSWFs on the interval [-2,2] 
x_ext = [-2+h/2:h:-1, x, 1:h:2-h/2]; 

Ny = 5;
C = linspecer(Ny); 
Ys = zeros(Ny, length(x_ext)); 
legends = {}; 

for idx = 1:Ny
    Ys(idx, :) = evalPSWF(Pbar, Veven, Vodd, idx-1, x_ext); 
    legends{idx} = sprintf('$\\psi_{%d}(x)$', idx-1); 
    % pswfEnergy = numericalInnerProd(Ys(idx, :), Ys(idx, :), h); 
end


%% ======== Plot PSWF eigenvalues and values within interval [a,b] ========
figure(1); 
Ndisp = 20; 
stem((1:Ndisp).', mu(1:Ndisp), 'Color', 'b'); grid on; 
xlabel(' Eigenvalue Index n', 'FontSize', 14); 
ylabel('$\lambda_n$', 'Interpreter', 'latex', 'FontSize', 14); 
ylim([0, 1.2]); 

if saveFigs
    exportgraphics(gcf, 'results/PSWFGraphics/Eigenvalues.pdf', 'ContentType','vector'); 
end

figure(2); 
plot((1:Ndisp).', mu(1:Ndisp), 'Marker','o', 'Color', 'b'); grid on; 
xlabel('n'); 
ylabel('$\lambda_n$', 'Interpreter','latex'); 
set(gca, 'yscale', 'log'); 

if saveFigs
    exportgraphics(gcf, 'results/PSWFGraphics/EigenvaluesLogScale.pdf', 'ContentType','vector'); 
end

figure(3); 
for idx = 1:Ny
    plot(x_ext, Ys(idx, :), 'Color', C(idx, :)); hold on; grid on; box on; 
end

legend(legends, 'Location','best', 'Interpreter','latex'); 
xlabel('Variable x'); 
ylabel('PSWF value'); 

if saveFigs
    exportgraphics(gcf, 'results/PSWFGraphics/PSWFs1D.pdf', 'ContentType','vector'); 
end




%% Utils. 



function Pbar = getNormalizedLegendrePolynomials(n)

    Pbar = zeros(n+1); 

    Pbar(1,1) = 1;  % P0(x) = 1
    Pbar(2,2) = 1;  % P1(x) = x

    for ord = 2:n
        Pbar(ord+1, :) = (2*ord-1)/(ord)*[0, Pbar(ord, 1:n)] - (ord-1)/ord*Pbar(ord-1,:); 
    end

    for idx = 1:n+1
        Pbar(idx, :) = Pbar(idx, :)*sqrt(1/2+idx-1); 
    end

end

function y = evalPoly(p, x)
    ord = length(p)-1;
    y = zeros(size(x));
    t = ones(size(x));
    for k = 0:ord
        y = y + p(k+1)*t;
        t = t .* x; 
    end
end

function y = evalPSWF(Pbar, Veven, Vodd, ordPSWF, x)
    [nmax, ~] = size(Pbar); 
    nmax = nmax - 1;  

    y = zeros(size(x)); 
    
    if mod(ordPSWF, 2) == 0
        v = Veven(:, ordPSWF/2+1); 
        Ps = Pbar(1:2:end, :); 
    else
        v = Vodd(:, (ordPSWF+1)/2); 
        Ps = Pbar(2:2:end, :); 
    end

    for idx = 1:length(v)
        y = y + v(idx)*evalPoly(Ps(idx, :), x); 
    end

end

function y = evalDerPSWF(Pbar, Veven, Vodd, ordPSWF, x)
    [nmax, ~] = size(Pbar); 
    nmax = nmax - 1;  

    y = zeros(size(x)); 
    
    if mod(ordPSWF, 2) == 0
        v = Veven(:, ordPSWF/2+1); 
        Ps = Pbar(1:2:end, :); 
    else
        v = Vodd(:, (ordPSWF+1)/2); 
        Ps = Pbar(2:2:end, :); 
    end

    dPs = Ps(:,2:end).*(1:nmax); 

    for idx = 1:length(v)
        y = y + v(idx)*evalPoly(dPs(idx, :), x); 
    end

end


function drawPoly(p, x)
    y = evalPoly(p, x); 
    plot(x, y);
end

function I = numericalInnerProd(x, y, h)
    I = sum(x.*y)*h;
end

function y = numericalBandlimitedOperatorApply(phi, c, x, h)
    L = length(x); 
    y = zeros(1, L);

    for idx = 1:L
        y(idx) = numericalInnerProd(exp(1i*c*x(idx)*x), phi, h); 
    end

end




