function [excess]=excess_of_quad(V1,V2,V3,V4)
% excess=excess_of_quad(V1,V2,V3,V4);
%
% where Vn are vectors

plane1=plane_normal(V1,V2);
plane2=plane_normal(V2,V3);
plane3=plane_normal(V3,V4);
plane4=plane_normal(V4,V1);

angle12=pi-angle_between_vectors(plane2,plane1);
angle23=pi-angle_between_vectors(plane3,plane2);
angle34=pi-angle_between_vectors(plane4,plane3);
angle41=pi-angle_between_vectors(plane1,plane4);

excess=angle12+angle23+angle34+angle41-2*pi;

%%%

function [plane]=plane_normal(P1,P2);

if ndims(P1)==4
 plane(:,:,:,1)=P1(:,:,:,2).*P2(:,:,:,3)-P1(:,:,:,3).*P2(:,:,:,2);
 plane(:,:,:,2)=P1(:,:,:,3).*P2(:,:,:,1)-P1(:,:,:,1).*P2(:,:,:,3);
 plane(:,:,:,3)=P1(:,:,:,1).*P2(:,:,:,2)-P1(:,:,:,2).*P2(:,:,:,1);
 mag=sqrt(plane(:,:,:,1).^2 + plane(:,:,:,2).^2 + plane(:,:,:,3).^2);
 plane(:,:,:,1)=plane(:,:,:,1)./mag;
 plane(:,:,:,2)=plane(:,:,:,2)./mag;
 plane(:,:,:,3)=plane(:,:,:,3)./mag;
elseif ndims(P1)==3
 plane(:,:,1)=P1(:,:,2).*P2(:,:,3)-P1(:,:,3).*P2(:,:,2);
 plane(:,:,2)=P1(:,:,3).*P2(:,:,1)-P1(:,:,1).*P2(:,:,3);
 plane(:,:,3)=P1(:,:,1).*P2(:,:,2)-P1(:,:,2).*P2(:,:,1);
 mag=sqrt(plane(:,:,1).^2 + plane(:,:,2).^2 + plane(:,:,3).^2);
 plane(:,:,1)=plane(:,:,1)./mag;
 plane(:,:,2)=plane(:,:,2)./mag;
 plane(:,:,3)=plane(:,:,3)./mag;
else
 error('Wrong ndims')
end