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);

















