clc;
clear;
close all;

% provadi se analyza pro chybu ve slozce p (tlak)

vystup_konzole = 1;
zapis_soubor = 0;
save_figures = 0;

% rady aproximace
N_vektor = [2, 4, 6, 8]; n_N = length(N_vektor);

vynechane_cleny = [];

% pocet elementu
% pocitano pro h = 4.00, 3.50, 3.00, 2.50, 2.00, 1.50, 1.00, 0.75
K_vektor = [246, 344, 460, 618, 942, 1734, 3726, 6750];

% pro ulozeni do souboru
K_vektor_print = K_vektor;

n_K = length(K_vektor); zustavajici_cleny = setdiff(1:n_K,vynechane_cleny);
sqrtK_vektor = sqrt(K_vektor);

% pripadne vynechani nekterych clenu
K_vektor = K_vektor(zustavajici_cleny);
sqrtK_vektor = sqrtK_vektor(zustavajici_cleny);

log_K_vektor = log10(K_vektor);
log_sqrtK_vektor = log10(sqrtK_vektor);

% E_max a E_avg metody pro jednotliva N a K
max_chyby_matice = zeros(n_N,n_K);
avg_chyby_matice = zeros(n_N,n_K);
% N = 2 - casovy krok nezkracovan
max_chyby_matice(1,:) = [3.1359e-02, 2.6246e-02, 1.3770e-02, 1.2507e-02, 7.5214e-03, 4.0132e-03, 1.3635e-03, 5.7294e-04];
avg_chyby_matice(1,:) = [4.3479e-03, 2.4920e-03, 2.0242e-03, 9.8412e-04, 5.1459e-04, 2.0169e-04, 6.0197e-05, 2.4142e-05];
% N = 4 - casovy krok 2x kratsi
max_chyby_matice(2,:) = [7.2047e-03, 3.4516e-03, 1.5245e-03, 8.7281e-04, 2.9358e-04, 6.8807e-05, 9.5886e-06, 2.1662e-06];
avg_chyby_matice(2,:) = [1.8752e-04, 7.2488e-05, 4.0577e-05, 1.7293e-05, 5.8318e-06, 1.3307e-06, 1.9756e-07, 4.4360e-08];
% N = 6 - casovy krok 2x kratsi
max_chyby_matice(3,:) = [5.0127e-04, 1.6431e-04, 4.5018e-05, 2.1232e-05, 4.5973e-06, 5.8943e-07, 1.6345e-07, 8.7338e-08];
avg_chyby_matice(3,:) = [6.7099e-06, 1.9886e-06, 7.5880e-07, 2.7365e-07, 6.1280e-08, 7.7720e-09, 5.1835e-10, 6.5117e-11];
% N = 8 - casovy krok 4x kratsi
max_chyby_matice(4,:) = [2.3659e-05, 6.5315e-06, 1.4698e-06, 3.6804e-07, 1.6456e-07, 1.1502e-08, 1.6277e-07, 6.5551e-07];
avg_chyby_matice(4,:) = [2.2836e-07, 4.7629e-08, 1.2306e-08, 3.5377e-09, 5.0908e-10, 3.5093e-11, 1.3209e-11, 7.0027e-12];

% pro ulozeni do souboru
max_chyby_matice_print = max_chyby_matice;
avg_chyby_matice_print = avg_chyby_matice;

% pripadne vynechani nekterych clenu
max_chyby_matice = max_chyby_matice(:,zustavajici_cleny);
avg_chyby_matice = avg_chyby_matice(:,zustavajici_cleny);

% handles na figures
fig_max_sqrtK = figure;
fig_avg_sqrtK = figure;
% barvy pro jednotliva N
color_spec_matice = zeros(n_N,3);
color_spec_matice(1,:) = [0.8, 0, 0];
color_spec_matice(2,:) = [0, 0, 0.8];
color_spec_matice(3,:) = [0, 0.6, 0];
color_spec_matice(4,:) = [0.7, 0, 0.7];
% marker size pro funkci scatter
marker_size = 40;
% pole pro ulozeni odkazu na plots
N_max_sqrtK = zeros(1,n_N);
N_avg_sqrtK = zeros(1,n_N);

fileID = fopen('Acoustic_pulse_propagation_konvergence_primky_Frontal.txt','w');
for i_N = 1:n_N
    N = N_vektor(i_N);
    if vystup_konzole == 1
        fprintf('N = %d\n', N);
    end
    if zapis_soubor == 1
        fprintf(fileID, 'N = %d\n', N);
    end
    
    max_chyby_vektor = max_chyby_matice(i_N,:);
    avg_chyby_vektor = avg_chyby_matice(i_N,:);
    
    log_max_chyby_vektor = log10(max_chyby_vektor);
    log_avg_chyby_vektor = log10(avg_chyby_vektor);
    
    % urceni radu prolozenim primky
    primka_K_max = polyfit(log_sqrtK_vektor,log_max_chyby_vektor,1);
    primka_K_avg = polyfit(log_sqrtK_vektor,log_avg_chyby_vektor,1);
    
    if vystup_konzole == 1
        primka_K_max
        primka_K_avg
    end
    if zapis_soubor == 1
        fprintf(fileID, 'max: ');
        fprintf(fileID, '%f ', primka_K_max);
        fprintf(fileID,'\n');
        fprintf(fileID, 'avg: ');
        fprintf(fileID, '%f ', primka_K_avg);
        fprintf(fileID,'\n\n');
    end
    
    color_spec = color_spec_matice(i_N,:);
    set(0, 'CurrentFigure', fig_max_sqrtK);
    scatter(sqrtK_vektor, max_chyby_vektor, marker_size, 'o', 'MarkerEdgeColor', 'black', 'MarkerFaceColor', color_spec); hold on;
    N_max_sqrtK(i_N) = plot(sqrtK_vektor, 10 .^ polyval(primka_K_max,log_sqrtK_vektor), 'Color', color_spec,'LineWidth', 2); hold on;
    set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
    set(0, 'CurrentFigure', fig_avg_sqrtK);
    scatter(sqrtK_vektor, avg_chyby_vektor, marker_size, 'o', 'MarkerEdgeColor', 'black', 'MarkerFaceColor', color_spec); hold on;
    N_avg_sqrtK(i_N) = plot(sqrtK_vektor, 10 .^ polyval(primka_K_avg,log_sqrtK_vektor), 'Color', color_spec, 'LineWidth', 2); hold on;
    set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
end
fclose(fileID);

% adjustments to the figures
set(0, 'CurrentFigure', fig_max_sqrtK);
xlim([10^1, 10^2]); ylim([10^-8, 10^-1]);
xlabel('sqrt(K)'); ylabel('E_{max}');
% legend(N_max_sqrtK, 'N = 2', 'N = 4', 'N = 6', 'N = 8', 'Location', 'SouthWest');
set(0, 'CurrentFigure', fig_avg_sqrtK);
xlim([10^1, 10^2]); ylim([10^-13, 10^-2]);
xlabel('sqrt(K)'); ylabel('E_{avg}');
% legend(N_avg_sqrtK, 'N = 2', 'N = 4', 'N = 6', 'N = 8', 'Location', 'SouthWest'); 


if save_figures == 1
    % E_avg against sqrtK
    filename_string = 'Acoustic_pulse_propagation_convergence_analysis_Eavg_against_sqrtK_Frontal';
    saveas(fig_avg_sqrtK,filename_string,'pdf');
    saveas(fig_avg_sqrtK,filename_string,'png');
end


% zapsani do diplomky na 2 desetinna mista
fileID = fopen('Acoustic_pulse_propagation_errors_Frontal.txt','w');
fprintf(fileID,'K:\n');
fprintf(fileID,'%d',K_vektor_print(1));
for i = 2:n_K
    fprintf(fileID,' %d',K_vektor_print(i));
end
fprintf(fileID,'\n\n');
% 2 des. mista
fprintf(fileID,'E_avg:\n');
for i = 1:n_N
    fprintf(fileID,'%0.2e',avg_chyby_matice_print(i,1));
    for j = 2:n_K
        fprintf(fileID,' %0.2e',avg_chyby_matice_print(i,j));
    end
    fprintf(fileID,'\n');
end
fprintf(fileID,'\n');
% 3 des. mista
fprintf(fileID,'E_avg:\n');
for i = 1:n_N
    fprintf(fileID,'%0.3e',avg_chyby_matice_print(i,1));
    for j = 2:n_K
        fprintf(fileID,' %0.3e',avg_chyby_matice_print(i,j));
    end
    fprintf(fileID,'\n');
end
fprintf(fileID,'\n');
% 4 des. mista
fprintf(fileID,'E_avg:\n');
for i = 1:n_N
    fprintf(fileID,'%0.4e',avg_chyby_matice_print(i,1));
    for j = 2:n_K
        fprintf(fileID,' %0.4e',avg_chyby_matice_print(i,j));
    end
    fprintf(fileID,'\n');
end
fclose(fileID);