clc; close all;

save_figures = 1;

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

t = 40;

y = 0;
x_step = 0.2; x_vect = x_l:x_step:x_u;
n = length(x_vect);
values_rho = zeros(1,n);
values_p = 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 = rho(:,j)';
    value = coef * values_node;
    values_rho(i) = value;
    coef = p(:,j)';
    value = coef * values_node;
    values_p(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_rho_exact = zeros(1,n);
values_p_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_rho_exact(i) = rho_point;
    values_p_exact(i) = p_point;
end
toc;

% vykresleni numerickeho reseni bodove
fig_rho = figure;
plot(x_vect_exact,values_rho_exact,'Color','red','LineWidth',2); hold on;
plot(x_vect,values_rho,'*');
xlim([x_l, x_u]);
ylim([-0.06, 0.14]);
legend('analytical','numerical','Location','NorthEast');
xlabel('x'); ylabel('\rho''');

fig_p = figure;
plot(x_vect_exact,values_p_exact,'Color','red','LineWidth',2); hold on;
plot(x_vect,values_p,'*');
xlim([x_l, x_u]);
ylim([-0.06, 0.13]);
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 = ['CAA_Workshop1_Category3_Problem1_rho_solution_cut_y=0_for_h=',mesh_size_string,',_N=',num2str(N),',_t=',num2str(t)];
    saveas(fig_rho,filename_string,'pdf');
    saveas(fig_rho,filename_string,'png');
    
    filename_string = ['CAA_Workshop1_Category3_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