Root_finding

% Experimental script to find the roots of the Laplacian spectrum function
% using different methods; for testing purposes only, see the eigfinder
% function for the current working version.

clear
clc
%close all
set(0,'defaultlegendFontSize',14)
set(0,'defaulttextInterpreter','latex') %latex axis labels
set(0,'defaultlegendInterpreter','latex') %latex axis labels
set(0,'defaultAxesFontSize',14)
set(0,'defaultAxesTickLabelInterpreter','latex')

%% Problem conditions %%
global m kappa D l Kp Km
m = 3;
kappa = ones(1,m)*0.1;
D = ones(1,m);
l = ones(1,m);
Kp = 0;
Km = 0;

%% Calculate function %%
N = 100;
maxeig=100;
eta = linspace(10^-3,maxeig,N);
F = zeros(1,N); % Using my function
F_Moutal = zeros(1,N); % Using Moutal's function
for i = 1:N
    F(i) = function_F(eta(i));
    F_Moutal(i) = general_condition(sqrt(eta(i)),sqrt(D),kappa,l,Km,Kp);
end

figure
plot(eta,F,'LineWidth',1.5)
hold on
plot(eta,F_Moutal,'LineWidth',1.5,'LineStyle','--')
plot([0 maxeig],[0 0],'k','LineWidth',1,'LineStyle',':')
xlabel('Trial value $\eta$')
ylabel('Function evaluation $F(\eta)$')

%% Find preliminary zeros %%
s = spline(eta,F);
z = fnzeros(s,[0 maxeig]); % Finds the eigenvalues within range
spectrum = [0 z(1,:)];

% Moutal's result
[spectrum_M,~]=function_modes(D, kappa, l, Km, Kp, 10, maxeig);

scatter(spectrum,zeros(1,numel(spectrum)),'LineWidth',1.5)
hold on
scatter(spectrum_M,zeros(1,numel(spectrum_M)),'LineWidth',1.5)
legend('Function','Zero','In-house','Moutal')

%% check if zeros are actually zeros %%
Z_Moutal = zeros(1,numel(spectrum_M));
for i = 1:numel(spectrum_M)
    Z_Moutal(i) = general_condition(sqrt(spectrum_M(i)),sqrt(D),kappa,l,Km,Kp);
    %Z_Moutal(i) = function_F(spectrum_M(i));
end
for i = numel(spectrum)
    Z(i) = function_F(spectrum(i));
end
figure
scatter(spectrum_M,spectrum_M)
hold on
scatter(spectrum,spectrum)

%% Get zeros fancy %%
clearvars F
maxeig=1000;
N = maxeig*10;
tic
eta = linspace(10^-3,maxeig,N);
for i = 1:N
    F(i) = function_F(eta(i));
end
threshold = 2;
dfine = 10^-2;
% Now we will try to identify the intervals close to or already including a
% zero. Note that using the logical arrays, we will identify the intervals
% to be explored further with ZEROS and the areas to be scrapped with ONES.
range = ((abs(F)<threshold) == 0);
signchange = (zeros(1,N)==0);
signchange(1:N-1) = (diff(sign(F)) == 0);
signchange(2:N) = signchange(2:N).*(diff(sign(F)) == 0);
% Note: two previous lines insure that if a sign change was missed, a new
% INTERVAL (i.e. two consecutive logical zeros) are added and not a single
% point (remember diff uses two consecutive points).
range = range.*signchange; % will add an interval if a signchange occurs outside of bounds
% Add the first zero
IDbounds = [1 find(diff(range))];
if (-1)^numel(IDbounds) ~= 1 % Essentially closes last bound if left open
    IDbounds = [IDbounds N];
end

plot(eta,F,'LineWidth',1,'Color','k')
hold on
I = 2;
roots = 0;
for i = 1:2:numel(IDbounds)
    lowb = eta(IDbounds(i));
    upb = eta(IDbounds(i+1));
    etafine = lowb:dfine:upb;
    Ffine = zeros(1,length(etafine));
    for j = 1:numel(etafine)
        Ffine(j) = function_F(etafine(j));
    end
    s = spline(etafine,Ffine);
    z = fnzeros(s,[min(etafine) max(etafine)]);
    roots = [roots z(1,:)];
    plot(etafine,Ffine,'LineWidth',4);
end
toc

tic
roots2 = eigfinder(maxeig,2);
toc

tic
roots_M=function_eigenvalues(D, kappa, l, Km, Kp, maxeig);
toc

plot([0 maxeig],[0 0],'k','LineWidth',1,'LineStyle',':')
scatter(roots2,zeros(1,numel(roots2)),35,'MarkerEdgeColor','black','LineWidth',1.5)
xlabel('$\eta$')
ylabel('$F(\eta)$')

disp(sum(abs(roots2-roots_M)))