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 LONGITUDE
%              Yi MUST BE LATITUDE
%              Zi is the value at the point (Xi,Yi)
%
%              Xi - Xj measures the horizontal distance between locations
%
%              Yi - Yi measures the vertical distance between locations
%
%              The angle between observations in degrees is given by:
%
%                   theta = [atan{(Yi-Yj)./(Xi-Xj)}]*180/pi
%
%  OUTPUT:  OUT = cell( n x (n+3) )with
%
%                                       cell(Xi,Yi,Zi,THETAi)
%

indata= M;

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

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

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

Z = M( :, 3);                             % column vector of values (prices/residuals)

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

XX = X*U' - U*X';                         % n x n matrix of longitude differences

YY = Y*U' - U*Y';                         % n x n matrix of latitude differences

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

theta = zeros(n,n);

% 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 positionstgttg

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


tempe=(theta<=45)|(theta>=315);            % identify properties within +/- 45 degrees

east=XX<0;                               % identify properties due east

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

numeast = U'*tempeast*U;                   % number of observations due east

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

symmde=disteast+disteast';                % symmetric matrix of distances

D =tril(symmde);

Dmin = min(min(D));

Dmax = max(max(D));

if(Dmin < Dmax/2)
   
   d = Dmax/2;
   
else
   
   error('Dmin>=Dmax/2');
   
end

[SS(:,1),SS(:,2),SS(:,3)]  =  find(D);

clear('D');

S(:,1) = (Z(SS(:,1))-Z(SS(:,2))).^2;

S(:,2)=SS(:,3);

clear('SS');

S = sortrows(S,2);

V = S(:,1);

D = S(:,2);

clear('S');

i = 1;

b = 1;

c = 0;

VV = 0;

d1 = D(1);

PAIRS = length(V);

NUM = max([ceil(PAIRS/50),30]);

while D(i) <= d
   
   c = c + 1;
   
   VV = VV + V(i);
   
   if (c>=NUM)&(D(i) < D(i+1));
      
      OUT(b,1) = (D(i) + d1)/2;
      
      OUT(b,2) = VV/(2*c);
      
      OUT(b,3) = c;
      
      b = b + 1;
      
      d1 = D(i);
      
      c = 0;
      
      VV = 0;
      
   end
   
   i = i + 1 ;
   
end

if (c<=20)&(b>1)
   
   OUT(b-1,1) = (D(i) + d1)/2;
   
   OUT(b-1,2) = (VV + OUT(b-1,2)*(2*OUT(b-1,3)))/(2*(c + OUT(b-1,3)));
   
   OUT(b-1,3) = OUT(b-1,3) + c;
   
else
   
   OUT(b,1) = (D(i) + d1)/2;
   
   OUT(b,2) = VV/(2*c);
   
   OUT(b,3) = c;
   
end





tempn = ((theta>=45)&(theta<=135))|((theta>=225)&(theta<=315));

                                          % identify properties within 45-135 degrees

north=YY<0;                                % identify properties due north

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

numnorth=U'*tempnort*U                    % number of properties north

tempw=(theta<=45)|(theta>=315);

west=XX>0;

tempwest=tempw.*west;

numwest=U'*tempwest*U

temps=((theta>=45)&(theta<=135))|((theta>=225)&(theta<=315));

south=YY>0;

tempsout=temps.*south;

numsouth=U'*tempsout*U