function areaq=quadriarea3D(xcoord1,ycoord1,zcoord1);
%function to calculate the area of a quadrilateral in 3D (may be non
%planar), by decomposing it into 2 triangles for which the areas are added
%to give the wuadirlateral area.
if nargin<3
   zcoord1=zeros(4,1); 
end

vec1=[xcoord1(2)-xcoord1(1),ycoord1(2)-ycoord1(1),zcoord1(2)-zcoord1(1)];
vec2=[xcoord1(4)-xcoord1(1),ycoord1(4)-ycoord1(1),zcoord1(4)-zcoord1(1)];
Areatri1=norm(cross(vec1,vec2))/2;
vec3=[xcoord1(2)-xcoord1(3),ycoord1(2)-ycoord1(3),zcoord1(2)-zcoord1(3)];
vec4=[xcoord1(4)-xcoord1(3),ycoord1(4)-ycoord1(3),zcoord1(4)-zcoord1(3)];
Areatri2=norm(cross(vec3,vec4))/2;
areaq=Areatri1+Areatri2;