clear all; close all;

f = 0.99; % Portion of q_max that is reached when Delta_t is evaluated

zeta = 1;
r_p_to_lambda = [1, 2, 5, 10, 20, 50];
L_p_to_r_p = [5, 10, 20, 50];

result = zeros(length(r_p_to_lambda)*length(L_p_to_r_p), 3); % Empty matrix for the reults. First column: L_p_to_r_p, second column: r_p_to_lambda, third column: Delta_t

for i=1:length(r_p_to_lambda)

    r = r_p_to_lambda(1, i);

    for j=1:length(L_p_to_r_p)

        L = L_p_to_r_p(j);

        %% load simulation data
        data_conv = readmatrix("../t_and_q/t_q_zeta" + zeta + "_r" + r + "_L" + L + "_MPNPNS.csv");
        data_ref = readmatrix("../t_and_q/t_q_zeta" + zeta + "_r" + r + "_L" + L + "_MPNP.csv");

        %% find unique rows
        [rows, columns] = size(data_conv);

        [q_data_conv,ia_conv,ic_conv] = unique(data_conv(:, columns), 'stable');
        t_data_conv = data_conv(ia_conv, columns-1);

        [q_data_ref,ia_ref,ic_ref] = unique(data_ref(:, columns), 'stable');
        t_data_ref = data_ref(ia_ref, columns-1);

        %% interpolate data
        t_max = t_data_conv(length(t_data_conv));
        t = linspace(0, t_max, 5000000);

        q_conv = interp1(t_data_conv, q_data_conv, t, 'spline');
        q_ref = interp1(t_data_ref, q_data_ref, t, 'spline');

        %% calculate charging times
        [val1, idx1] = min(abs(q_conv-f*max(q_conv)));
        t_f_conv = t(idx1);

        [val2, idx2] = min(abs(q_ref-f*max(q_ref)));
        t_f_ref = t(idx2);

        %% calculate relative deviations in charging times
        Delta_t = (t_f_ref - t_f_conv)/t_f_conv;

        result(i+(j-1)*length(r_p_to_lambda), 1) = L;
        result(i+(j-1)*length(r_p_to_lambda), 2) = r;
        result(i+(j-1)*length(r_p_to_lambda), 3) = Delta_t;

    end;
end;

%% export results to a .txt-file adding a header for plotting the heatmap
header = {'L_to_r', 'r_to_lambda', 'Delta_t'};
exportWithHeader = [header; num2cell(result)];
writecell(exportWithHeader, "data_Deltat" + f*100 + "_zeta" + zeta + ".txt", 'Delimiter', 'tab');
