clc; close all;

save_figures = 1;

% srovnani presneho a numerickeho reseni v case t = 10 v rezu y = 0

t = 10;

y = 0;
x_step = 0.1; x_vect = x_l:x_step:x_u;
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.01; 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
    [rho_point,u_point,v_point,p_point] = ExactSolution(x_vect_exact(i),y,t);
    values_exact(i) = p_point; 
end
toc;

% kontrola, ze bod je skutecne uvnitr trojuhelniku a ze r,s bylo napocitano dobre
% figure;
% trimesh(E(:,j)',P(1,:),P(2,:)); hold on;
% x = 0.5 * (-(r + s) * vertices(1,1) + (r + 1) * vertices(1,2) + (s + 1) * vertices(1,3));
% y = 0.5 * (-(r + s) * vertices(2,1) + (r + 1) * vertices(2,2) + (s + 1) * vertices(2,3));
% plot(x,y,'*');

% vykresleni numerickeho reseni bodove
fig_p = figure;
plot(x_vect_exact,values_exact,'Color','red','LineWidth',2); hold on;
plot(x_vect,values,'*');
legend('analytical','numerical','Location','NorthEast');

if save_figures == 1
    mesh_size_string = sprintf('%0.2f',mesh_size_parametr);
    mesh_size_string = strrep(mesh_size_string,'.',',');
    
    filename_string = ['Bajko_Problem1_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