function [n_x,n_y,n_z,surfJ] = Normals3D_my()

Globals3D;

n_x = zeros(N_faces * N_fp,K); n_y = zeros(N_faces * N_fp,K); n_z = zeros(N_faces * N_fp,K);
faceIDs = [(1:N_fp)', ((N_fp + 1):(2 * N_fp))', ((2 * N_fp + 1):(3 * N_fp))', ((3 * N_fp + 1):(4 * N_fp))'];
face_def = [1, 2, 3; 1, 2, 4; 2, 3, 4; 1, 3, 4];

for i = 1:N_faces
    face = face_def(i,:);
    
    % smerove a normalove vektory
    t_matrix1 = [P(1,E(face(1),:)) - P(1,E(face(2),:)); P(2,E(face(1),:)) - P(2,E(face(2),:)); P(3,E(face(1),:)) - P(3,E(face(2),:))];
    t_matrix2 = [P(1,E(face(1),:)) - P(1,E(face(3),:)); P(2,E(face(1),:)) - P(2,E(face(3),:)); P(3,E(face(1),:)) - P(3,E(face(3),:))];
    n_matrix = cross(t_matrix1,t_matrix2,1);
    
    % otoceni vnitrnich normal
    const_vect = -(n_matrix(1,:) .* P(1,E(face(1),:)) + n_matrix(2,:) .* P(2,E(face(1),:)) + n_matrix(3,:) .* P(3,E(face(1),:)));
    pointID = setdiff(1:4,face);
    points = [P(1,E(pointID,:)); P(2,E(pointID,:)); P(3,E(pointID,:))];
    signs = n_matrix(1,:) .* points(1,:) + n_matrix(2,:) .* points(2,:) + n_matrix(3,:) .* points(3,:) + const_vect;
    IDs = (signs > 0);
    n_matrix(:,IDs) = -n_matrix(:,IDs);
    
    % zapsani do n_x a n_y
    n_x(faceIDs(:,i),:) = repmat(n_matrix(1,:),N_fp,1);
    n_y(faceIDs(:,i),:) = repmat(n_matrix(2,:),N_fp,1);
    n_z(faceIDs(:,i),:) = repmat(n_matrix(3,:),N_fp,1);
end

surfJ = sqrt(n_x .* n_x + n_y .* n_y + n_z .* n_z);
n_x = n_x ./ surfJ; n_y = n_y ./ surfJ; n_z = n_z ./ surfJ;
surfJ = surfJ / 2; % protoze puvodne byly smerove vektory jen poloviny stran
surfJ = surfJ .* J(Fmask(:),:); % dle vzorce na str. 426 (Hesthaven dG)

end