clc; close all;

save_figures = 1;

% srovnani presneho a numerickeho reseni v case t = 10 v rezech

t = 10;

y = 0;
x_step = 0.2; x_vect = x_l:x_step:x_u;
x = 0;
y_step = 0.2; y_vect = y_l:y_step:y_u;
n = length(x_vect);
values_u = zeros(1,n);
values_v = zeros(1,n);
values_p = zeros(1,n);

tic;
for i = 1:n
    x_point = 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_point - 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_point - vertices(1,3)) + (vertices(1,3) - vertices(1,2)) * (y - vertices(2,3))) / det_T;
    lambda3 = ((vertices(2,3) - vertices(2,1)) * (x_point - 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 = u(:,j)';
    values_u(i) = coef * values_node;
    coef = p(:,j)';
    values_p(i) = coef * values_node;
end
for i = 1:n
    y_point = y_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_point - 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_point - vertices(2,3))) / det_T;
    lambda3 = ((vertices(2,3) - vertices(2,1)) * (x - vertices(1,3)) + (vertices(1,1) - vertices(1,3)) * (y_point - 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 = v(:,j)';
    values_v(i) = coef * values_node;
end
toc;


x_step_exact = 0.01; x_vect_exact = x_l:x_step_exact:x_u;
y_step_exact = 0.01; y_vect_exact = y_l:y_step_exact:y_u;
n = length(x_vect_exact);
values_u_exact = zeros(1,n);
values_v_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_u_exact(i) = u_point;
    values_p_exact(i) = p_point;
    [rho_point,u_point,v_point,p_point] = ExactSolution(x,y_vect_exact(i),t);
    values_v_exact(i) = v_point;
end
toc;

% vykresleni numerickeho reseni bodove
fig_u = figure;
plot(x_vect_exact,values_u_exact,'Color','red','LineWidth',2); hold on;
plot(x_vect,values_u,'*');
legend('analytical','numerical','Location','NorthEast');
ylim([-0.2,0.25]);
xlabel('x'); ylabel('u''');

fig_v = figure;
plot(y_vect_exact,values_v_exact,'Color','red','LineWidth',2); hold on;
plot(y_vect,values_v,'*');
legend('analytical','numerical','Location','NorthEast');
ylim([-0.2,0.25]);
xlabel('y'); ylabel('v''');

fig_p = figure;
plot(x_vect_exact,values_p_exact,'Color','red','LineWidth',2); hold on;
plot(x_vect,values_p,'*');
legend('analytical','numerical','Location','NorthEast');
ylim([-0.1,0.20]);
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 = ['Acoustic_pulse_propagation_u_solution_cut_y=0_for_h=',mesh_size_string,',_N=',num2str(N),',_t=',num2str(t)];
    saveas(fig_u,filename_string,'pdf');
    saveas(fig_u,filename_string,'png');
    
    filename_string = ['Acoustic_pulse_propagation_v_solution_cut_y=0_for_h=',mesh_size_string,',_N=',num2str(N),',_t=',num2str(t)];
    saveas(fig_v,filename_string,'pdf');
    saveas(fig_v,filename_string,'png');
    
    filename_string = ['Acoustic_pulse_propagation_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