function OUT = tgt_theta(M)

%	theta computes the angle between observations in a plane
%
%	INPUT:  M = ( n x 3 ) matrix of locations
%            
%              M(i) = (Xi,Yi,Zi), i = 1:n
%
%              Xi MUST BE LATITUDE
%              Yi MUST BE LONGITUDE
%              Zi is the value at the point (Xi,Yi)
%
%              Xi - Xj measures the vertical distance between locations
%
%              Yi - Yi measures the horizontal distance between locations
%
%              The angle between observations in degrees is given by:
%
%                   theta = [atan{(Xi-Xj)./(Yi-Yj)}]*180/pi
%
%  OUTPUT:  OUT = cell( n x (n+3) )with
%
%					cell(Xi,Yi,Zi,THETAi)
%

n = length(M( :,1));                      % number of observations

X = M( :, 1);                             % column vector of latitudes

Y = M( :, 2);                             % column vector of longitudes

U = ones(n,1);                            % column vector of ones

XX = X*U' - U*X' + eye(n);                % n x n matrix of latitude differences


YY = (Y*U' - U*Y');                       % n x n matrix of longitude differences

D = sqrt(XX.^2 + YY.^2);                  % compute distances for all properties

% The following procedure computes the
% matrix XX./YY with zeros in all (ij)-components
% where YY(ij) = 0.

[I,J,V] = find(YY);

uu = ones(length(V),1);

V = uu ./ V ;

YYinv = full(sparse(I,J,V,n,n));

ZZ = XX .* YYinv ;

theta = ((atan(ZZ))*180/pi);          % computing the angle theta

OUT = [M,theta] ;                         % concatenate M and theta

temp=abs(theta)<=45;                      % identify properties within +/- 45 degrees

east=YY<=0;                               % identify properties due east

tempeast=temp.*east;                      % identify properties in both sets

disteast=tempeast.*D;                     % matrix of distances for properties east

symmde=disteast+disteast';



%temp = (theta<=135)&(theta>=45);          % identify properties within 45-135 degrees

%north=XX<0;                               % identify properties due north

%tempnort=temp.*north;                     % identify properties in both sets