matlab_class/gcmfaces_misc/gcmfaces_polygonangle.m

function [angsum]=gcmfaces_polygonangle(px,py,x,y);
%[angsum]=gcmfaces_polygonangle(px,py,x,y);
%object:    compute sum of interior angles for polygons (when input
%           is px,py) or points vs polygons (when input is px,py,x,y)
%inputs:    px,py are MxN matrices where each line specifies one polygon
%(optional) x,y are position vectors
%outputs:   ang are the corresponding sums of interior angles

M=size(px,1); N=size(px,2);
doPointsInPolygon=0; P=1;
if nargin>2;
    doPointsInPolygon=1;
    sizxy=size(x);
    x=reshape(x,[1 prod(sizxy)]);
    y=reshape(y,[1 prod(sizxy)]);
    P=length(x);
end;

angsum=zeros(M,P);
for ii=0:N-1;
    ppx=circshift(px,[0 -ii]);
    ppy=circshift(py,[0 -ii]);
    
    if ~doPointsInPolygon;
        %compute sum of corner angles
    v1x=ppx(:,2)-ppx(:,1); 
    v1y=ppy(:,2)-ppy(:,1);
    v2x=ppx(:,4)-ppx(:,1); 
    v2y=ppy(:,4)-ppy(:,1);
    else;
        %compute sum of sector angles        
        v1x=ppx(:,1)*ones(1,P)-ones(M,1)*x; 
        v1y=ppy(:,1)*ones(1,P)-ones(M,1)*y;
        v2x=ppx(:,2)*ones(1,P)-ones(M,1)*x; 
        v2y=ppy(:,2)*ones(1,P)-ones(M,1)*y;
    end;
    g_acos=acos( ( v1x.*v2x+v1y.*v2y )./sqrt( v1x.*v1x+v1y.*v1y )./sqrt( v2x.*v2x+v2y.*v2y ) );
    g_sin= ( v1x.*v2y-v1y.*v2x )./sqrt( v1x.*v1x+v1y.*v1y )./sqrt( v2x.*v2x+v2y.*v2y ) ;
    angsum=angsum+radtodeg(g_acos.*sign(g_sin));    
end;

1	function [angsum]=gcmfaces_polygonangle(px,py,x,y);
2	%[angsum]=gcmfaces_polygonangle(px,py,x,y);
3	%object: compute sum of interior angles for polygons (when input
4	% is px,py) or points vs polygons (when input is px,py,x,y)
5	%inputs: px,py are MxN matrices where each line specifies one polygon
6	%(optional) x,y are position vectors
7	%outputs: ang are the corresponding sums of interior angles
8
9	M=size(px,1); N=size(px,2);
10	doPointsInPolygon=0; P=1;
11	if nargin>2;
12	doPointsInPolygon=1;
13	sizxy=size(x);
14	x=reshape(x,[1 prod(sizxy)]);
15	y=reshape(y,[1 prod(sizxy)]);
16	P=length(x);
17	end;
18
19	angsum=zeros(M,P);
20	for ii=0:N-1;
21	ppx=circshift(px,[0 -ii]);
22	ppy=circshift(py,[0 -ii]);
23
24	if ~doPointsInPolygon;
25	%compute sum of corner angles
26	v1x=ppx(:,2)-ppx(:,1);
27	v1y=ppy(:,2)-ppy(:,1);
28	v2x=ppx(:,4)-ppx(:,1);
29	v2y=ppy(:,4)-ppy(:,1);
30	else;
31	%compute sum of sector angles
32	v1x=ppx(:,1)ones(1,P)-ones(M,1)x;
33	v1y=ppy(:,1)ones(1,P)-ones(M,1)y;
34	v2x=ppx(:,2)ones(1,P)-ones(M,1)x;
35	v2y=ppy(:,2)ones(1,P)-ones(M,1)y;
36	end;
37	g_acos=acos( ( v1x.v2x+v1y.v2y )./sqrt( v1x.v1x+v1y.v1y )./sqrt( v2x.v2x+v2y.v2y ) );
38	g_sin= ( v1x.v2y-v1y.v2x )./sqrt( v1x.v1x+v1y.v1y )./sqrt( v2x.v2x+v2y.v2y ) ;
39	angsum=angsum+radtodeg(g_acos.*sign(g_sin));
40	end;
41