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 = [6]; % pro N = 8 a h = 1.00 chyba uplne mimo

% pocet elementu
% pocitano pro h = 3.50, 3.00, 2.50, 2.00, 1.50, 1.00, 0.75
K_vektor = [368, 518, 678, 1030, 1936, 4272, 7688];

% 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 10x kratsi
% max_chyby_matice(1,:) = [2.3091e-02, 2.9135e-02, 2.3790e-02, 1.3654e-02, 8.3578e-03, 2.2009e-03]; % 9.8200e-04
% avg_chyby_matice(1,:) = [4.2612e-03, 1.3892e-03, 1.0141e-03, 6.1554e-04, 2.0102e-04, 5.7369e-05]; % 2.2861e-05
% N = 2 - casovy krok nezkracovan
max_chyby_matice(1,:) = [2.3086e-02, 2.9158e-02, 2.3781e-02, 1.3678e-02, 8.3586e-03, 2.2009e-03, 9.8162e-04];
avg_chyby_matice(1,:) = [4.2538e-03, 1.3878e-03, 1.0123e-03, 6.1484e-04, 2.0088e-04, 5.7355e-05, 2.2860e-05];
% N = 4 - casovy krok 10x kratsi
% max_chyby_matice(2,:) = [3.7568e-03, 3.5955e-03, 2.3409e-03, 6.0863e-04, 2.3495e-04, 2.2047e-05];
% avg_chyby_matice(2,:) = [1.2369e-04, 3.9545e-05, 2.3110e-05, 7.4786e-06, 1.4877e-06, 1.9464e-07];
% N = 4 - casovy krok nezkracovan, az od h = 1.00 2x kratsi
max_chyby_matice(2,:) = [3.7553e-03, 3.5947e-03, 2.3405e-03, 6.0869e-04, 2.3492e-04, 2.2047e-05, 6.4371e-06];
avg_chyby_matice(2,:) = [1.2356e-04, 3.9533e-05, 2.3103e-05, 7.4770e-06, 1.4876e-06, 1.9464e-07, 4.5008e-08];
% N = 6 - casovy krok 10x kratsi
% max_chyby_matice(3,:) = [2.2600e-04, 2.2115e-04, 8.0096e-05, 1.3562e-05, 3.6162e-06, 1.6458e-07];
% avg_chyby_matice(3,:) = [3.3211e-06, 9.8937e-07, 4.4173e-07, 8.9728e-08, 1.0835e-08, 6.1203e-10];
% N = 6 - casovy krok 2x kratsi, od h = 1.50 4x kratsi
max_chyby_matice(3,:) = [2.2600e-04, 2.2115e-04, 8.0097e-05, 1.3562e-05, 3.6162e-06, 1.6458e-07, 2.6427e-08];
avg_chyby_matice(3,:) = [3.3210e-06, 9.8936e-07, 4.4172e-07, 8.9736e-08, 1.0835e-08, 6.1203e-10, 7.6270e-11];
% N = 8 - casovy krok 10x kratsi
% max_chyby_matice(4,:) = [1.1162e-05, 8.7087e-06, 2.1544e-06, 3.3304e-07, 1.6272e-07, 6.5546e-07];
% avg_chyby_matice(4,:) = [1.3071e-07, 2.2014e-08, 7.7000e-09, 9.9954e-10, 7.7271e-11, 2.8551e-11];
% N = 8 - casovy krok 4x kratsi (prepocitano jen pro h = 3.50, 1.00 a 0.75)
% pro h = 1.00 odlehla hodnota, kazi jinak optimalni rad konvergence
max_chyby_matice(4,:) = [1.1162e-05, 8.7087e-06, 2.1544e-06, 3.3304e-07, 1.6272e-07, 6.5546e-07, 1.1496e-08];
avg_chyby_matice(4,:) = [1.3071e-07, 2.2014e-08, 7.7000e-09, 9.9954e-10, 7.7271e-11, 2.8551e-11, 4.6948e-13];
% N = 8 - casovy krok 4x kratsi - spocteno pro h = 3.50 az 1.50
% max_chyby_matice(4,:) = [1.1162e-05, 8.7087e-06, 2.1544e-06, 3.3305e-07, 1.6271e-07, 6.5546e-07, 1.1496e-08];
% avg_chyby_matice(4,:) = [1.3071e-07, 2.2014e-08, 7.7674e-09, 1.0685e-09, 1.3398e-10, 2.8551e-11, 4.6948e-13];

% 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_Delaunay.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]);
legend(N_avg_sqrtK, 'N = 2', 'N = 4', 'N = 6', 'N = 8', 'Location', 'SouthWest');
% legend(N_avg_sqrtK, 'N = 2', 'N = 4', 'N = 6', 'Location', 'SouthWest'); % bez N 8 
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]);
legend(N_avg_sqrtK, 'N = 2', 'N = 4', 'N = 6', 'N = 8', 'Location', 'SouthWest');
% xlim([10^1, 10^2]); ylim([10^-11, 10^-2]); % bez N 8 
% legend(N_avg_sqrtK, 'N = 2', 'N = 4', 'N = 6', 'Location', 'SouthWest'); % bez N 8 
xlabel('sqrt(K)'); ylabel('E_{avg}');


if save_figures == 1
    % E_avg against sqrtK
    filename_string = 'Acoustic_pulse_propagation_convergence_analysis_Eavg_against_sqrtK_Delaunay';
    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_Delaunay.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);