clc; close all;

N_plot = 15; % should be chosen based on how fine the grid is
N_p_plot = (N_plot * N_plot + N_plot) / 2;
r_1 = linspace(-1,1,N_plot); % nodes on the face of I on axis x
s_1 = linspace(-1,1,N_plot); % nodes on the face of I on axis y
r_plot = r_1;
s_plot = -1 * ones(1,N_plot);
for i = 2:N_plot
    r_old = r_plot;
    s_old = s_plot;
    clear r_plot s_plot;
    r_temp = r_1(1:(N_plot - (i - 1)));
    s_temp = s_1(i) * ones(1,N_plot - (i - 1));
    r_plot = [r_old, r_temp];
    s_plot = [s_old, s_temp];
end

[a_plot,b_plot] = rstoab(r_plot,s_plot);
values_mode = zeros(N_p,N_p_plot);
temp = 1;
for i = 0:N
    for j = 0:(N - i)
        values_mode(temp,:) = ModalBaseValues2D(a_plot,b_plot,i,j);
        temp = temp + 1;
    end
end
values_node = inv(V') * values_mode;

mesh_size_string = sprintf('%0.2f',mesh_size_parametr);
mesh_size_string = strrep(mesh_size_string,'.',',');
interval = 5;
n = length(p_states);

for state = 1:n
    p = p_states{state};
    
    fig_p_CS = figure('Name','Pressure computed solution');
    
    for i = 1:K
        element = E(:,i); v_1_ID = element(1); v_2_ID = element(2); v_3_ID = element(3);
        x_plot = 0.5 * (-(r_plot + s_plot)' * P(1,v_1_ID) + (r_plot + 1)' * P(1,v_2_ID) + (s_plot + 1)' * P(1,v_3_ID));
        y_plot = 0.5 * (-(r_plot + s_plot)' * P(2,v_1_ID) + (r_plot + 1)' * P(2,v_2_ID) + (s_plot + 1)' * P(2,v_3_ID));
        tri = delaunay(x_plot,y_plot);
        
        coef = p(:,i)';
        values = coef * values_node;
        
        trisurf(tri, x_plot, y_plot, values); hold on;
    end
    
    shading interp;
    
    xlabel('x'); ylabel('y'); zlabel('p');
    
    colorbar('Location','EastOutside');
    view(0, 90);
    axis equal;
    xlim(x_axis_limits); ylim(y_axis_limits);

    % skalovani barev
    cmax = 0.5;
    cmin = -cmax;
    caxis([cmin cmax]);

    zlim([-1 1]);
    plot3([x_l x_u x_u x_l x_l],[y_l y_l y_u y_u y_l],[0.9 0.9 0.9 0.9 0.9],'LineWidth',1,'Color','black');
    
    filename_string = ['Bailly_Juve_monopole_subsonic_p_solution_for_h=',mesh_size_string,',_N=',num2str(N),',_T=',num2str(state * interval)];
    % saveas(fig_p_CS,filename_string,'pdf'); % neni ted potreba, pro nazornost png staci
    saveas(fig_p_CS,filename_string,'png');
    
    close all;
end