clc; close all;

save_figures = 0;

% obrazek z Bailly, Juve - srovnani presneho a numerickeho reseni v case t = 270 v rezu y = 0
interval = 5;
state = 6; % state = 54 odpovida casu t = 270, protoze interval je 5
p = p_states{state};
t = state * interval; % tj t = 270
K = max(elements_Euler);

y = 0;
x_step = 0.5; x_vect = x_l:x_step:x_u; % pro sve analytical solution
% x_step = 0.1; x_vect = x_l:x_step:x_u; % pro Nogueira analytical solution
n = length(x_vect);
values = zeros(1,n);

tic;
for i = 1:n
    x = x_vect(i);
    for j = 1:K
        vertices = P(:,E(:,j));
        A = [vertices(1,2) - vertices(1,1), vertices(1,3) - vertices(1,1); vertices(2,2) - vertices(2,1), vertices(2,3) - vertices(2,1)];
        b = [x - vertices(1,1); y - vertices(2,1)];
        koef_temp = A\b;
        koef = [1 - sum(koef_temp); koef_temp];
        
        if (min(koef) >= 0)
            break;
        end
    end
    
    T_matrix = [vertices(1,1) - vertices(1,3), vertices(1,2) - vertices(1,3); vertices(2,1) - vertices(2,3), vertices(2,2) - vertices(2,3)];
    det_T = det(T_matrix);
    lambda2 = ((vertices(2,2) - vertices(2,3)) * (x - vertices(1,3)) + (vertices(1,3) - vertices(1,2)) * (y - vertices(2,3))) / det_T;
    lambda3 = ((vertices(2,3) - vertices(2,1)) * (x - vertices(1,3)) + (vertices(1,1) - vertices(1,3)) * (y - vertices(2,3))) / det_T;
    lambda1 = 1 - lambda2 - lambda3;
    
    r = -lambda2 + lambda3 - lambda1;
    s = -lambda2 - lambda3 + lambda1;
    
    [a,b] = rstoab(r,s);
    values_mode = zeros(N_p,1);
    temp = 1;
    for k = 0:N
        for l = 0:(N - k)
            values_mode(temp,:) = ModalBaseValues2D(a,b,k,l);
            temp = temp + 1;
        end
    end
    values_node = inv(V') * values_mode;
    
    coef = p(:,j)';
    value = coef * values_node;
    values(i) = value;
end
toc;


x_step_exact = 0.05; x_vect_exact = x_l:x_step_exact:x_u;
n = length(x_vect_exact);
values_exact = zeros(1,n);

tic;
for i = 1:n
    p_point = ExactSolution(x_vect_exact(i),y,t);
    values_exact(i) = p_point;
end
toc;


% vykresleni numerickeho reseni bodove
fig_p = figure;
% Nogueira
% load('Nogueira_analytical_solution_subsonic.mat');
% plot(Nogueira_XData,Nogueira_YData,'*'); hold on;
% plot(x_vect,values,'LineWidth',2,'Color','red');
% xlim([x_l, x_u]);
% ylim([-0.4, 0.3]);
% legend('analytical','numerical','Location','NorthEast');
% xlabel('x'); ylabel('p');
% moje vlastni
plot(x_vect_exact,values_exact,'Color','red','LineWidth',2); hold on;
plot(x_vect,values,'*');
xlim([x_l, x_u]);
ylim([-0.4, 0.3]); % nastavenio pro state = 54
legend('analytical','numerical','Location','NorthEast');
xlabel('x'); ylabel('p''');

if save_figures == 1
    mesh_size_string = sprintf('%0.2f',mesh_size_parametr);
    mesh_size_string = strrep(mesh_size_string,'.',',');
    
    filename_string = ['Bailly_Juve_monopole_subsonic_p_solution_cut_y=0_for_h=',mesh_size_string,',_N=',num2str(N),',_t=',num2str(t)];
    saveas(fig_p,filename_string,'pdf');
    saveas(fig_p,filename_string,'png');
end