clc; close all;

save_figures = 1;

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

t = 5;

mesh_size_string = sprintf('%0.3f',mesh_size_parametr);
mesh_size_string = strrep(mesh_size_string,'.',',');


% rovina z = 0
z = 0;
x_step = 0.05; x_vect = x_l:x_step:x_u;
y_step = 0.05; y_vect = y_l:y_step:y_u;

[X_z,Y_z] = meshgrid(x_vect,y_vect);
tri_z = delaunay(X_z,Y_z);

n = length(x_vect);
values_p_z = zeros(n,n);

tic;
for i_x = 1:n % jeden radek obrazku (pro 400 bodu) trva zpracovat cca 10 sekund
    x_point = x_vect(i_x);
    for i_y = 1:n
        y_point = y_vect(i_y);
%         for j = 1:K % toto lze napat efektivneji pomoci tsearchn - hnedani bodu v simplexu
%             vertices = P(:,E(:,j));
%             A = [vertices(1,2) - vertices(1,1), vertices(1,3) - vertices(1,1), vertices(1,4) - vertices(1,1); vertices(2,2) - vertices(2,1), vertices(2,3) - vertices(2,1), vertices(2,4) - vertices(2,1); vertices(3,2) - vertices(3,1), vertices(3,3) - vertices(3,1), vertices(3,4) - vertices(3,1)];
%             b = [x_point - vertices(1,1); y_point - vertices(2,1); z - vertices(3,1)];
%             koef_temp = A\b;
%             koef = [1 - sum(koef_temp); koef_temp];
%             
%             if (min(koef) >= 0)
%                 break;
%             end
%         end
        
        % pri vypoctu lambda je evidentne chyba
%         T_matrix = [vertices(1,1) - vertices(1,4), vertices(1,2) - vertices(1,4), vertices(1,3) - vertices(1,4); vertices(2,1) - vertices(2,4), vertices(2,2) - vertices(2,4), vertices(2,3) - vertices(2,4); vertices(3,1) - vertices(3,4), vertices(3,2) - vertices(3,4), vertices(3,3) - vertices(3,4)];
%         b_T = [x_point - vertices(1,4); y_point - vertices(2,4); z - vertices(3,4)];
%         lambdas = T\b_T;
%         lambda3 = lambdas(1); lambda4 = lambdas(2); lambda2 = lambdas(3);
%         lambda1 = 1 - lambda2 - lambda3 - lambda4;

        [j,lambdas] = tsearchn(P',E',[x_point,y_point,z]);
        lambda3 = lambdas(1); lambda4 = lambdas(2); lambda2 = lambdas(3); lambda1 = lambdas(4);

        r = -lambda3 + lambda4 - lambda2 - lambda1;
        s = -lambda3 - lambda4 + lambda2 - lambda1;
        t = -lambda3 - lambda4 - lambda2 + lambda1;
        
        [a,b,c] = rsttoabc(r,s,t);
        values_mode = zeros(N_p,1);
        temp = 1;
        for k = 0:N
            for l = 0:(N - k)
                for m = 0:(N - k - l)
                    values_mode(temp,:) = ModalBaseValues3D(a,b,c,k,l,m);
                    temp = temp + 1;
                end
            end
        end
        values_node = inv(V') * values_mode;
        
        coef = p(:,j)';
        values_p_z(i_y,i_x) = coef * values_node;
    end
    fprintf(2,'layer %d of %d processed\n',i_x,n);
end
toc;

figure;
surf(x_vect, y_vect, values_p_z, 'EdgeColor', 'none');
xlabel('x'); ylabel('y'); zlabel('p''');
colorbar('Location','EastOutside');
view(0, 90);
grid on;
axis equal;
xlim(x_axis_limits); ylim(y_axis_limits);

filename = sprintf('solution_p_cut_z=0_for_N=%d_and_h=%s',N,mesh_size_string);
save(filename,'values_p_z');

if save_figures == 1
    mesh_size_string = sprintf('%0.3f',mesh_size_parametr);
    mesh_size_string = strrep(mesh_size_string,'.',',');
    filename = sprintf('rez_z=0_pro_N=%d_and_h=%s',N,mesh_size_string);
    saveas(gcf,filename,'png');
    saveas(gcf,filename,'pdf');
end


% rovina x = 0
x = 0;
y_step = 0.05; y_vect = y_l:y_step:y_u;
z_step = 0.05; z_vect = z_l:z_step:z_u;

[Y_x,Z_x] = meshgrid(y_vect,z_vect);
tri_x = delaunay(Y_x,Z_x);

n = length(y_vect);
values_p_x = zeros(n,n);

tic;
for i_y = 1:n % jeden radek obrazku (pro 400 bodu) trva zpracovat cca 10 sekund
    y_point = y_vect(i_y);
    for i_z = 1:n
        z_point = z_vect(i_z);
%         for j = 1:K % toto lze napat efektivneji pomoci tsearchn - hnedani bodu v simplexu
%             vertices = P(:,E(:,j));
%             A = [vertices(1,2) - vertices(1,1), vertices(1,3) - vertices(1,1), vertices(1,4) - vertices(1,1); vertices(2,2) - vertices(2,1), vertices(2,3) - vertices(2,1), vertices(2,4) - vertices(2,1); vertices(3,2) - vertices(3,1), vertices(3,3) - vertices(3,1), vertices(3,4) - vertices(3,1)];
%             b = [x_point - vertices(1,1); y_point - vertices(2,1); z - vertices(3,1)];
%             koef_temp = A\b;
%             koef = [1 - sum(koef_temp); koef_temp];
%             
%             if (min(koef) >= 0)
%                 break;
%             end
%         end
        
        % pri vypoctu lambda je evidentne chyba
%         T_matrix = [vertices(1,1) - vertices(1,4), vertices(1,2) - vertices(1,4), vertices(1,3) - vertices(1,4); vertices(2,1) - vertices(2,4), vertices(2,2) - vertices(2,4), vertices(2,3) - vertices(2,4); vertices(3,1) - vertices(3,4), vertices(3,2) - vertices(3,4), vertices(3,3) - vertices(3,4)];
%         b_T = [x_point - vertices(1,4); y_point - vertices(2,4); z - vertices(3,4)];
%         lambdas = T\b_T;
%         lambda3 = lambdas(1); lambda4 = lambdas(2); lambda2 = lambdas(3);
%         lambda1 = 1 - lambda2 - lambda3 - lambda4;

        [j,lambdas] = tsearchn(P',E',[x,y_point,z_point]);
        lambda3 = lambdas(1); lambda4 = lambdas(2); lambda2 = lambdas(3); lambda1 = lambdas(4);

        r = -lambda3 + lambda4 - lambda2 - lambda1;
        s = -lambda3 - lambda4 + lambda2 - lambda1;
        t = -lambda3 - lambda4 - lambda2 + lambda1;
        
        [a,b,c] = rsttoabc(r,s,t);
        values_mode = zeros(N_p,1);
        temp = 1;
        for k = 0:N
            for l = 0:(N - k)
                for m = 0:(N - k - l)
                    values_mode(temp,:) = ModalBaseValues3D(a,b,c,k,l,m);
                    temp = temp + 1;
                end
            end
        end
        values_node = inv(V') * values_mode;
        
        coef = p(:,j)';
        values_p_x(i_z,i_y) = coef * values_node;
    end
    fprintf(2,'layer %d of %d processed\n',i_x,n);
end
toc;

figure;
surf(y_vect, z_vect, values_p_x, 'EdgeColor', 'none');
xlabel('y'); ylabel('z'); zlabel('p''');
colorbar('Location','EastOutside');
view(0, 90);
grid on;
axis equal;
xlim(y_axis_limits); ylim(z_axis_limits);

filename = sprintf('solution_p_cut_x=0_for_N=%d_and_h=%s',N,mesh_size_string);
save(filename,'values_p_x');

if save_figures == 1
    mesh_size_string = sprintf('%0.3f',mesh_size_parametr);
    mesh_size_string = strrep(mesh_size_string,'.',',');
    filename = sprintf('rez_x=0_pro_N=%d_and_h=%s',N,mesh_size_string);
    saveas(gcf,filename,'png');
    saveas(gcf,filename,'pdf');
end


% rovina y = 0
y = 0;
x_step = 0.05; x_vect = x_l:x_step:x_u;
z_step = 0.05; z_vect = z_l:z_step:z_u;

[X_y,Z_y] = meshgrid(x_vect,z_vect);
tri_y = delaunay(X_y,Z_y);

n = length(x_vect);
values_p_y = zeros(n,n);

tic;
for i_x = 1:n % jeden radek obrazku (pro 400 bodu) trva zpracovat cca 10 sekund
    x_point = x_vect(i_x);
    for i_z = 1:n
        z_point = z_vect(i_z);
%         for j = 1:K % toto lze napat efektivneji pomoci tsearchn - hnedani bodu v simplexu
%             vertices = P(:,E(:,j));
%             A = [vertices(1,2) - vertices(1,1), vertices(1,3) - vertices(1,1), vertices(1,4) - vertices(1,1); vertices(2,2) - vertices(2,1), vertices(2,3) - vertices(2,1), vertices(2,4) - vertices(2,1); vertices(3,2) - vertices(3,1), vertices(3,3) - vertices(3,1), vertices(3,4) - vertices(3,1)];
%             b = [x_point - vertices(1,1); y_point - vertices(2,1); z - vertices(3,1)];
%             koef_temp = A\b;
%             koef = [1 - sum(koef_temp); koef_temp];
%             
%             if (min(koef) >= 0)
%                 break;
%             end
%         end
        
        % pri vypoctu lambda je evidentne chyba
%         T_matrix = [vertices(1,1) - vertices(1,4), vertices(1,2) - vertices(1,4), vertices(1,3) - vertices(1,4); vertices(2,1) - vertices(2,4), vertices(2,2) - vertices(2,4), vertices(2,3) - vertices(2,4); vertices(3,1) - vertices(3,4), vertices(3,2) - vertices(3,4), vertices(3,3) - vertices(3,4)];
%         b_T = [x_point - vertices(1,4); y_point - vertices(2,4); z - vertices(3,4)];
%         lambdas = T\b_T;
%         lambda3 = lambdas(1); lambda4 = lambdas(2); lambda2 = lambdas(3);
%         lambda1 = 1 - lambda2 - lambda3 - lambda4;

        [j,lambdas] = tsearchn(P',E',[x_point,y,z_point]);
        lambda3 = lambdas(1); lambda4 = lambdas(2); lambda2 = lambdas(3); lambda1 = lambdas(4);

        r = -lambda3 + lambda4 - lambda2 - lambda1;
        s = -lambda3 - lambda4 + lambda2 - lambda1;
        t = -lambda3 - lambda4 - lambda2 + lambda1;
        
        [a,b,c] = rsttoabc(r,s,t);
        values_mode = zeros(N_p,1);
        temp = 1;
        for k = 0:N
            for l = 0:(N - k)
                for m = 0:(N - k - l)
                    values_mode(temp,:) = ModalBaseValues3D(a,b,c,k,l,m);
                    temp = temp + 1;
                end
            end
        end
        values_node = inv(V') * values_mode;
        
        coef = p(:,j)';
        values_p_y(i_y,i_x) = coef * values_node;
    end
    fprintf(2,'layer %d of %d processed\n',i_x,n);
end
toc;

figure;
surf(x_vect, z_vect, values_p_y, 'EdgeColor', 'none');
xlabel('x'); ylabel('z'); zlabel('p''');
colorbar('Location','EastOutside');
view(0, 90);
grid on;
axis equal;
xlim(x_axis_limits); ylim(z_axis_limits);

filename = sprintf('solution_p_cut_y=0_for_N=%d_and_h=%s',N,mesh_size_string);
save(filename,'values_p_y');

if save_figures == 1
    mesh_size_string = sprintf('%0.3f',mesh_size_parametr);
    mesh_size_string = strrep(mesh_size_string,'.',',');
    filename = sprintf('rez_y=0_pro_N=%d_and_h=%s',N,mesh_size_string);
    saveas(gcf,filename,'png');
    saveas(gcf,filename,'pdf');
end



% obrazek ve 3D
view3D = figure;
trisurf(tri_z,X_z,Y_z,zeros(size(X_z)),values_p_z); hold on;
trisurf(tri_x,zeros(size(Y_x)),Y_x,Z_x,values_p_x); hold on;
trisurf(tri_y,X_y,zeros(size(X_z)),Z_y,values_p_y); hold on;
shading interp;
axis equal;
xlabel('x'); ylabel('y'); zlabel('z');
xlim(x_axis_limits); ylim(y_axis_limits); zlim(z_axis_limits);