clc;
clear;
close all;

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

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

vynechane_cleny = [1];

% interval, na kterem ulohu resime, a pocet elementu
interval = 20;
K_vektor = [10, 15, 20, 40, 60, 80];

% pro ulozeni do souboru
K_vektor_print = K_vektor;

n_K = length(K_vektor); zustavajici_cleny = setdiff(1:n_K,vynechane_cleny);
h_vektor = interval ./ K_vektor;

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

n_h = length(h_vektor); % pocet ruznych jemnosti diskretizace skutecne uvazovanych
log_h_vektor = log10(h_vektor);
log_K_vektor = log10(K_vektor);

% E_max a E_avg metody pro jednoliva N a K
max_chyby_matice = zeros(n_N,n_K);
avg_chyby_matice = zeros(n_N,n_K);
% % u N = 2, 4 nezkracovano; u N = 6, 8 casovy krok 10 kratsi
% % N = 2
% max_chyby_matice(1,:) = [0.14905, 0.057242, 0.030397, 0.0043693, 0.0012970, 5.7830e-04];
% avg_chyby_matice(1,:) = [0.032319, 0.010305, 0.0037791, 2.9964e-04, 7.6836e-05, 3.0999e-05];
% % N = 4
% max_chyby_matice(2,:) = [0.013484, 0.0030931, 8.5616e-04, 2.3486e-05, 3.6468e-06, 8.3908e-07];
% avg_chyby_matice(2,:) = [0.0015983, 1.7472e-04, 4.0923e-05, 1.0055e-06, 1.2789e-07, 3.0509e-08];
% % N = 6
% max_chyby_matice(3,:) = [3.6778e-04, 7.6329e-05, 1.1562e-05, 7.9363e-08, 5.7937e-09, 7.5751e-10];
% avg_chyby_matice(3,:) = [2.8858e-05, 2.5682e-06, 5.4227e-07, 3.1383e-09, 1.7525e-10, 2.3030e-11];
% % N = 8
% max_chyby_matice(4,:) = [9.3538e-06, 1.5501e-06, 1.2389e-07, 2.2523e-10, 6.7669e-12, 5.6632e-13];
% avg_chyby_matice(4,:) = [8.0318e-07, 3.9753e-08, 4.5657e-09, 6.7361e-12, 1.7691e-13, 1.4683e-14];
% casovy krok nezkracovan - u N = 4, 6 u poslednich dvou 4x kratsi a u N = 8 casovy krok 10x kratsi
% N = 2
max_chyby_matice(1,:) = [1.4905e-01, 5.7601e-02, 3.0397e-02, 4.3693e-03, 1.2970e-03, 5.7830e-04];
avg_chyby_matice(1,:) = [3.2319e-02, 1.0340e-02, 3.7791e-03, 2.9964e-04, 7.6836e-05, 3.0999e-05];
% N = 4
max_chyby_matice(2,:) = [1.3484e-02, 3.0931e-03, 8.5616e-04, 2.3486e-05, 3.5695e-06, 8.1582e-07];
avg_chyby_matice(2,:) = [1.5983e-03, 1.7472e-04, 4.0923e-05, 1.0053e-06, 1.2721e-07, 3.0088e-08];
% N = 6
max_chyby_matice(3,:) = [3.6778e-04, 7.6329e-05, 1.1562e-05, 7.9338e-08, 5.7739e-09, 7.5278e-10];
avg_chyby_matice(3,:) = [2.8858e-05, 2.5682e-06, 5.4227e-07, 3.1569e-09, 1.7544e-10, 2.3106e-11];
% N = 8
max_chyby_matice(4,:) = [9.3538e-06, 1.5501e-06, 1.2390e-07, 2.2523e-10, 6.8378e-12, 5.6866e-13];
avg_chyby_matice(4,:) = [8.0318e-07, 3.9753e-08, 4.5658e-09, 6.7377e-12, 1.7749e-13, 1.5215e-14];

% 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_h = zeros(1,n_h - 1);
rady_avg_vektor_h = zeros(1,n_h - 1);
rady_max_vektor_K = zeros(1,n_h - 1);
rady_avg_vektor_K = zeros(1,n_h - 1);

% handles na figures
fig_max_h = figure;
fig_avg_h = figure;
fig_max_K = figure;
fig_avg_K = 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_h = zeros(1,n_N);
N_avg_h = zeros(1,n_N);
N_max_K = zeros(1,n_N);
N_avg_K = zeros(1,n_N);

fileID = fopen('Advection1D_konvergence_primky.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_h_max = polyfit(log_h_vektor,log_max_chyby_vektor,1);
    primka_h_avg = polyfit(log_h_vektor,log_avg_chyby_vektor,1);
    primka_K_max = polyfit(log_K_vektor,log_max_chyby_vektor,1);
    primka_K_avg = polyfit(log_K_vektor,log_avg_chyby_vektor,1);
    % urceni radu pomoci podilu
    for i = 1:(n_h - 1)
        rady_max_vektor_h(i) = log10(max_chyby_vektor(i) / max_chyby_vektor(i + 1)) / log10(h_vektor(i) / h_vektor(i + 1));
        rady_avg_vektor_h(i) = log10(avg_chyby_vektor(i) / avg_chyby_vektor(i + 1)) / log10(h_vektor(i) / h_vektor(i + 1));
        rady_max_vektor_K(i) = log10(max_chyby_vektor(i) / max_chyby_vektor(i + 1)) / log10(K_vektor(i) / K_vektor(i + 1));
        rady_avg_vektor_K(i) = log10(avg_chyby_vektor(i) / avg_chyby_vektor(i + 1)) / log10(K_vektor(i) / K_vektor(i + 1));
    end
    rad_podily_max_h = mean(rady_max_vektor_h);
    rad_podily_avg_h = mean(rady_avg_vektor_h);
    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_h
        rad_podily_avg_h
        rad_podily_max_K
        rad_podily_avg_K
        primka_h_max
        primka_h_avg
        primka_K_max
        primka_K_avg
    end
    if zapis_soubor == 1
        fprintf(fileID, 'max: ');
        fprintf(fileID, '%f ', primka_h_max);
        fprintf(fileID,'\n');
        fprintf(fileID, 'avg: ');
        fprintf(fileID, '%f ', primka_h_avg);
        fprintf(fileID,'\n\n');
    end
    
    color_spec = color_spec_matice(i_N,:);
    set(0, 'CurrentFigure', fig_max_h);
    scatter(log_h_vektor, log_max_chyby_vektor, marker_size, 'o', 'MarkerEdgeColor', 'black', 'MarkerFaceColor', color_spec); hold on;
    N_max_h(i_N) = plot(log_h_vektor, polyval(primka_h_max,log_h_vektor), 'Color', color_spec,'LineWidth', 2); hold on;
    set(0, 'CurrentFigure', fig_avg_h);
    scatter(log_h_vektor, log_avg_chyby_vektor, marker_size, 'o', 'MarkerEdgeColor', 'black', 'MarkerFaceColor', color_spec); hold on;
    N_avg_h(i_N) = plot(log_h_vektor, polyval(primka_h_avg,log_h_vektor), 'Color', color_spec, 'LineWidth', 2); hold on;
    set(0, 'CurrentFigure', fig_max_K);
    scatter(K_vektor, max_chyby_vektor, marker_size, 'o', 'MarkerEdgeColor', 'black', 'MarkerFaceColor', color_spec); hold on;
    N_max_K(i_N) = plot(K_vektor, 10 .^ polyval(primka_K_max,log_K_vektor), 'Color', color_spec,'LineWidth', 2); hold on;
    set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
    set(0, 'CurrentFigure', fig_avg_K);
    scatter(K_vektor, avg_chyby_vektor, marker_size, 'o', 'MarkerEdgeColor', 'black', 'MarkerFaceColor', color_spec); hold on;
    N_avg_K(i_N) = plot(K_vektor, 10 .^ polyval(primka_K_avg,log_K_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_h);
xlabel('log_{10}(h)'); ylabel('log_{10}(E_{max})');
legend(N_max_h, 'N = 2', 'N = 4', 'N = 6', 'N = 8', 'Location', 'SouthEast');
set(0, 'CurrentFigure', fig_avg_h);
xlabel('log_{10}(h)'); ylabel('log_{10}(E_{avg})');
legend(N_avg_h, 'N = 2', 'N = 4', 'N = 6', 'N = 8', 'Location', 'SouthEast');
set(0, 'CurrentFigure', fig_max_K);
xlim([10^1, 10^2]); ylim([10^-13, 10^-1]);
xlabel('K'); ylabel('E_{max}');
legend(N_max_K, 'N = 2', 'N = 4', 'N = 6', 'N = 8', 'Location', 'SouthWest');
set(0, 'CurrentFigure', fig_avg_K);
xlim([10^1, 10^2]); ylim([10^-14, 10^-1]);
xlabel('K'); ylabel('E_{avg}');
% xlabel('K'); ylabel('');
legend(N_avg_K, 'N = 2', 'N = 4', 'N = 6', 'N = 8', 'Location', 'SouthWest');

if save_figures == 1
    % E_max against K
    filename_string = 'Advection1D_convergence_analysis_Emax_against_K';
    saveas(fig_max_K,filename_string,'pdf');
    saveas(fig_max_K,filename_string,'png');
    
    % E_avg against K
    filename_string = 'Advection1D_convergence_analysis_Eavg_against_K';
    saveas(fig_avg_K,filename_string,'pdf');
    saveas(fig_avg_K,filename_string,'png');
end


% zapsani do diplomky na 2 desetinna mista
fileID = fopen('Advection1D_errors.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);