function [normVectors,cosProjs]=plotAnglePlane(angles)
if size(angles,1)<size(angles,2)
    angles=angles';
end

% initialize normal vectors and rotation matrices array
% normVECs = zeros(size(angles,1),3);
rotMATs = zeros(3,3,size(angles,1));
cosProjs = zeros(size(angles,1),1);

phis = angles(:,1);
thetas = angles(:,2);
psis = angles(:,3);


figure; hold on; 
zlim([-1 1]);
xlim([-1 1]);

ylim([-1 1]);

axis equal;
grid on
view([1 1 1]);

normVectors=zeros(3,size(angles,1));

% calculate rotation matrices and normal vectors from the rotation matrices
for i=1:size(angles,1)
    vector1 = [0 0 1];
    rotmat1 = MatrixQuaternionRot(vector1,phis(i));
    
    vector2 = [0 1 0];
    rotmat2 = MatrixQuaternionRot(vector2,thetas(i));

    vector3 = [1 0 0];
    rotmat3 = MatrixQuaternionRot(vector3,psis(i));

    rotMATs(:,:,i) =  (rotmat1*rotmat2*rotmat3);

    init_normvec = [0 0 1];
    normVec = squeeze(rotMATs(:,:,i))*init_normvec';
    A=normVec(1);
    B=normVec(2);
    C=normVec(3);
    
    normVectors(:,i)=normVec;
    cosProjs(i)=acosd(C);
       
    [x y] = ndgrid(-1:0.1:1); % Generate x and y data
    z = -1/C*(A*x + B*y); % Solve for z vertices data

    surf(y,x,z) %Plot the surface
    vectarrow([0 0 0],normVec([2 1 3]));
%     pause

    
%     x = [1 -1 -1 1]; % Generate data for x vertices
%     y = [1 1 -1 -1]; % Generate data for y vertices
%     patch(x, y, z);

% pause
end
hold off

fprintf('The corresponding angles are:\n');
cosProjs