function [P0] = project2Plane(P,N,Q,m)

%  Let P be the m x 3 array of the 3D points to be projected, let Q be the
%  1 x 3 vector of the given point on the plane, let N be the 1 x 3 vector 
%  of the normal direction to the plane, and let P0 be the m x 3 array of 
%  points orthogonally projected from P onto the plane. Then do this:
% 
N = N/norm(N); % <-- do this if N is not normalized
% V = P-Q;
% D = dot(V,D);
% P0 = P-D*N


N2 = N.'*N;
P0 = P*(eye(3)-N2)+repmat(Q*N2,m,1);
%P0 = P0/norm(P0);

end