function [sx,sy,s]=slope(md)
%SLOPE - compute the surface slope
%
%   Usage:
%      [sx,sy,s]=slope(md)

%load some variables (it is much faster if the variab;es are loaded from md once for all) 
if strcmpi(md.type,'2d'),
	numberofelements=md.numberofelements;
	numberofgrids=md.numberofgrids;
	index=md.elements;
	x=md.x; y=md.y; z=md.z;
else
	numberofelements=md.numberofelements2d;
	numberofgrids=md.numberofgrids2d;
	index=md.elements2d;
	x=md.x2d; y=md.y2d; z=md.z2d;
end

%initialization
alpha=zeros(numberofelements,3);
beta=zeros(numberofelements,3);

%build some usefull variables
summation=1/3*ones(3,1);
x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));

%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
invdet=1./(x1.*(y2-y3)-x2.*(y1-y3)+x3.*(y1-y2));
alpha=[invdet.*(y2-y3) invdet.*(y3-y1) invdet.*(y1-y2)];
beta=[invdet.*(x3-x2) invdet.*(x1-x3) invdet.*(x2-x1)];

summation=[1;1;1];
sx=(md.surface(index).*alpha)*summation;
sy=(md.surface(index).*beta)*summation;
s=sqrt(sx.^2+sy.^2);

if strcmpi(md.type,'3d'),
	sx=project3d(md,sx,'element');
	sy=project3d(md,sy,'element');
	s=sqrt(sx.^2+sy.^2);
end
