function out=rjsurfacenormal(X,Y,Z,varargin)
%RJSURFACENORMAL
% Calculates the vertex normals (non-normalized) over the whole surface.
% If the dimension of X is (m x n) -- the result is returned in
% the form of a (m x n x 3) matrix
%
% The algorithm is not the best imaginable - it is not quite equivalent
% to the one in rjfindnormal - but that one is not optimal either.
% One would like to have it such that when working with a sphere - but
% with different steps in phi and theta - one would get the exact
% normal. This only gets the exact normal if the steps in phi and theta
% are uniform - but it is fast:).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Extending the matrixes linearly, to cope with edges
X=rjextendmat(X);
Y=rjextendmat(Y);
Z=rjextendmat(Z);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pos(:,:,1)=X;
pos(:,:,2)=Y;
pos(:,:,3)=Z;
s1=size(X,1);
s2=size(X,2);
temp=[pos zeros(s1,2,3)] - [zeros(s1,2,3) pos];
diff_column=temp(2:size(temp,1)-1,3:(size(temp,2)-2),:);
temp=[pos;zeros(2,s2,3)] - [zeros(2,s2,3);pos];
diff_row=temp(3:(size(temp,1)-2),2:(size(temp,2)-1),:);
out=cross(diff_column,diff_row,3);