% ============================================================= % 
% 
%  This script compares different PSWF-evaluation methods 
% 
%   {Use Fredholm determinants to evaluate something ... }
% 
% ============================================================= % 

clear; close all; clc; 

W1           = 4; 
Wc          = 0; % Up-conversion frequency parameter 
b1           = 10; 
a1           = -10; 
h           = 0.02; 
saveFigs    = true; 

t_arr = ((3*a1)+h/2:h:(3*b1)-h/2); 
Nmax = 8; 

% % Method 1 
% [lambdas1, phi_q1] = getShifted1DPSWFs(t_arr, a1, b1, Wc, W1, h, ceil(3*Nmax+1), Nmax); 
% 
% % Method 2
% [lambdas2, phi_q2] = get1DPSWFs_numerical(t_arr, a1, b1, Wc, W1, h, ceil(3*Nmax+1), 12); 

% Multi-interval tests 
intervals = {[-3, -2]/10, [1, 2]/10, [5, 6.5]/10}; 
[lambdas3, phi_q3] = get1DPSWFs_multipleFreqInterval(t_arr, a1, b1, intervals, h, 12); 


%% Visualization 
close all; clc; 

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

fig1 = figure(1); 
subplot(4, 2, 1); 
showPSWFTimeDomainConcentration(t_arr, phi_q3(1, :), a1, b1); 

subplot(4, 2, 2); 
showPSWFSpectralConcentration(phi_q3(1,:), intervals, 16384, h);  

subplot(4, 2, 3); 
showPSWFTimeDomainConcentration(t_arr, phi_q3(2, :), a1, b1); 

subplot(4, 2, 4); 
showPSWFSpectralConcentration(phi_q3(2,:), intervals, 16384, h); 

subplot(4, 2, 5); 
showPSWFTimeDomainConcentration(t_arr, phi_q3(3, :), a1, b1); 

subplot(4, 2, 6); 
showPSWFSpectralConcentration(phi_q3(3,:), intervals, 16384, h); 

subplot(4, 2, 7); 
showPSWFTimeDomainConcentration(t_arr, phi_q3(4, :), a1, b1); 

subplot(4, 2, 8); 
showPSWFSpectralConcentration(phi_q3(4,:), intervals, 16384, h); 

if saveFigs
    exportgraphics(fig1, './results/PSWFGraphics/PSWF_space_wavenumber.pdf', 'ContentType', 'vector'); 
    fprintf('Files saved.\n'); 
end



% subplot(3, 2, 5); 
% stem((0:length(lambdas3)-1).', lambdas3); 

% plot((0:length(lambdas1)-1).', lambdas1); hold on; 
% plot((0:length(lambdas2)-1).', lambdas2); 
% grid on; 
% 
% xlabel('Eigenvalue index $n$', 'Interpreter', 'latex'); 
% ylabel('PSWF Eigenvalue', 'Interpreter', 'latex');  
% 
% % yyaxis right; 
% % plot((0:length(lambdas2)-1).', abs(lambdas2-lambdas1(1:length(lambdas2)))); 

% We can plot the angular power distribution of the generated PSWFs 
% Fourier transform convention: *exp(-i*2pi*f*t)


fprintf('End of this script.\n'); 


function showPSWFTimeDomainConcentration(t_arr, x_sig, a, b)
    h1 = plot(t_arr, real(x_sig), 'Color', 'r'); hold on; 
    h2 = plot(t_arr, imag(x_sig), 'Color', 'b');  
    h3 = plot(t_arr, abs(x_sig), 'Color', 'k', 'LineStyle', ':'); 
    
    xlabel('Space $x/\lambda$', 'Interpreter', 'latex'); 
    ylabel('PSWF value'); 

    area_arr = zeros(1, length(t_arr)); 
    area_arr(and(t_arr <= b, t_arr >= a)) = 0.5; 
    area(t_arr, area_arr, "FaceAlpha", 0.1, "FaceColor", "b"); 
    area(t_arr, -area_arr, "FaceAlpha", 0.1, "FaceColor", "b"); 

    legend([h1, h2, h3], {"Real", "Imag", "Envelope"}); 
end

function showPSWFSpectralConcentration(x_sig, intervals, Nfft, h)
    Nx = length(x_sig); 
    
    x_sig = [x_sig, zeros(1, Nfft-Nx)];     % zero-padding to increase spectral resolution 
    blk_window = blackman(Nfft, "periodic").'; 
    X_spec = fftshift(fft((x_sig.*blk_window)*h)); 
    
    freq_arr = (-Nfft/2:1:(Nfft/2-1))/Nfft/h;   
    
    % Show the spectrum in units of 1/(b-a). Spectral analysis.  
    
    area_arr = zeros(1, Nfft); 
    for idx = 1:length(intervals)
        interval = intervals{idx}; 
        area_arr(and(freq_arr <= interval(2), freq_arr >= interval(1))) = 1; 
    end
    area_arr = area_arr * max(abs(X_spec)); 
    area(freq_arr, area_arr, 'FaceAlpha', 0.1, 'FaceColor', 'red', 'LineStyle',':'); hold on; 
    plot(freq_arr, abs(X_spec), 'Color', 'r'); hold on; 
    xlim([-1.2, 1.2]); 
    xlabel('Wavenumber $k/k_c$', 'Interpreter', 'latex'); 
    ylabel('PSWF spectrum'); 
end

