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' ;                        % 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); % Finds all nonzero components of YY

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

V = uu ./ V ;

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

ZZ = XX .* YYinv ;

% Now proceed with ZZ replacing XX./YY 

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


% Now add 360 to all negative values of theta

Q = min(theta, zeros(n)) ;  % put zeros in all positive positions

[I,J,V] = find(Q~=0);  

QQ = full(sparse(I,J,V,n,n)) ; % 0-1 matrix with 1's in negative cells of theta

theta = theta + 360*QQ ; %add 360 to all negative positions

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