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 = 2.00, 1.75, 1.50, 1.25, 1.00, 0.75, 0.50
K_vektor = [224, 342, 458, 616, 944, 1728, 3718];

% 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);
% N = 2 - casovy krok nezkracovan
max_chyby_matice(1,:) = [1.9550e-01, 1.3403e-01, 1.0421e-01, 6.5356e-02, 3.6598e-02, 1.8360e-02, 7.2501e-03];
avg_chyby_matice(1,:) = [3.2007e-03, 1.8792e-03, 1.1354e-03, 7.9014e-04, 3.2492e-04, 1.1877e-04, 3.4265e-05];
% N = 4 - casovy krok 2x kratsi
max_chyby_matice(2,:) = [1.6805e-02, 1.0586e-02, 6.0507e-03, 2.9393e-03, 1.2057e-03, 3.8054e-04, 6.3874e-05];
avg_chyby_matice(2,:) = [1.4037e-04, 4.7092e-05, 1.9030e-05, 1.2033e-05, 3.7279e-06, 8.2749e-07, 1.2068e-07];
% N = 6 - casovy krok 2x kratsi
max_chyby_matice(3,:) = [1.4096e-03, 4.2847e-04, 2.1477e-04, 6.3698e-05, 2.0107e-05, 3.2897e-06, 2.6365e-07];
avg_chyby_matice(3,:) = [4.8226e-06, 1.3715e-06, 4.5310e-07, 2.0192e-07, 4.0190e-08, 5.0215e-09, 3.2645e-10];
% N = 8 - casovy krok 2x kratsi
max_chyby_matice(4,:) = [7.9159e-05, 1.4940e-05, 4.5516e-06, 8.2401e-07, 1.9200e-07, 2.3351e-08, 7.5329e-10];
avg_chyby_matice(4,:) = [1.8208e-07, 3.2528e-08, 8.0024e-09, 2.5414e-09, 3.7155e-10, 2.4089e-11, 8.4053e-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);

% 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_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);
    % 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^-10, 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^-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 = 'Advection2D_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('Advection2D_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);