clc;
clear;
close all;

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 = 1.75, 1.50, 1.25, 1.00, 0.75, 0.50
K_vektor = [372, 520, 690, 1024, 1934, 4272];

% 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);

n_h = length(sqrtK_vektor); % pocet ruznych jemnosti diskretizace skutecne uvazovanych
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);
% casovy krok 10x kratsi
% % N = 2 - casovy krok 10x kratsi
% max_chyby_matice(1,:) = [0.1307, 0.0794, 0.0603, 0.0377, 0.0188, 0.0094];
% avg_chyby_matice(1,:) = [0.0019, 0.0013, 6.8617e-04, 4.0122e-04, 1.1764e-04, 3.7801e-05];
% N = 2 - casovy krok nezkracovan, az od h = 0.75 2x kratsi
max_chyby_matice(1,:) = [1.3065e-01, 7.9388e-02, 6.0272e-02, 3.7670e-02, 1.8830e-02, 9.3700e-03];
avg_chyby_matice(1,:) = [1.8914e-03, 1.2951e-03, 6.8564e-04, 4.0110e-04, 1.1761e-04, 3.7797e-05];
% % N = 4 - casovy krok 10x kratsi
% max_chyby_matice(2,:) = [0.0096, 0.0049, 0.0024, 0.0021, 3.8458e-04, 1.0822e-04];
% avg_chyby_matice(2,:) = [7.9298e-05, 2.9079e-05, 1.3418e-05, 7.4753e-06, 9.7391e-07, 1.9352e-07];
% N = 4 - casovy krok nezkracovan, az od h = 0.75 2x kratsi
max_chyby_matice(2,:) = [9.6249e-03, 4.8879e-03, 2.3657e-03, 2.1191e-03, 3.8459e-04, 1.0822e-04];
avg_chyby_matice(2,:) = [7.9296e-05, 2.9074e-05, 1.3414e-05, 7.4749e-06, 9.7392e-07, 1.9352e-07];
% % N = 6 - casovy krok 10x kratsi
% max_chyby_matice(3,:) = [3.8946e-04, 1.5211e-04, 6.8580e-05, 6.3249e-05, 4.0270e-06, 4.6996e-07];
% avg_chyby_matice(3,:) = [2.4961e-06, 7.1816e-07, 3.4740e-07, 1.0642e-07, 7.8802e-09, 7.9630e-10];
% N = 6 - casovy krok nezkracovan, az od h = 0.75 2x kratsi
max_chyby_matice(3,:) = [3.8942e-04, 1.5206e-04, 6.8560e-05, 6.3243e-05, 4.0268e-06, 4.6989e-07];
avg_chyby_matice(3,:) = [2.4959e-06, 7.1817e-07, 3.4737e-07, 1.0645e-07, 7.8802e-09, 7.9620e-10];
% % N = 8 - casovy krok 10x kratsi
% max_chyby_matice(4,:) = [2.7288e-05, 2.6705e-06, 1.5066e-06, 1.1925e-06, 4.7741e-08, 1.4508e-09];
% avg_chyby_matice(4,:) = [9.5408e-08, 1.5839e-08, 7.3130e-09, 1.3191e-09, 5.2328e-11, 2.9710e-12];
% N = 8 - casovy krok nezkracovan, az od h = 0.75 2x kratsi
max_chyby_matice(4,:) = [2.7290e-05, 2.6667e-06, 1.5090e-06, 1.1933e-06, 4.7766e-08, 1.4578e-09];
avg_chyby_matice(4,:) = [9.5395e-08, 1.5838e-08, 7.3102e-09, 1.3176e-09, 5.2386e-11, 2.9779e-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);

% predalokovana pole pro urcovani radu pomoci podilu
rady_max_vektor_K = zeros(1,n_h - 1);
rady_avg_vektor_K = zeros(1,n_h - 1);

% 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('Advection2D_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);
    % urceni radu pomoci podilu
    for i = 1:(n_h - 1)
        rady_max_vektor_K(i) = log10(max_chyby_vektor(i) / max_chyby_vektor(i + 1)) / log10(sqrtK_vektor(i) / sqrtK_vektor(i + 1));
        rady_avg_vektor_K(i) = log10(avg_chyby_vektor(i) / avg_chyby_vektor(i + 1)) / log10(sqrtK_vektor(i) / sqrtK_vektor(i + 1));
    end
    rad_podily_max_K = mean(rady_max_vektor_K);
    rad_podily_avg_K = mean(rady_avg_vektor_K);
    if vystup_konzole == 1
        rad_podily_max_K
        rad_podily_avg_K
        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^-9, 10^0]);
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^-12, 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 = 'Advection2D_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('Advection2D_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);