function DAT = variogram(M, b)

%VARIOGRAM computes the empirical (semi)variogram for a given
%set of spatial data.
%
%INPUTS: (i)  M = (nx3)-matrix of locations and values,
%                 M(i) = (Xi,Yi,Zi), i = 1:n
%        (ii) b = number of distance bins to use
%
%OUTPUTS: DAT = (bX3)-matrix with DAT(i)=(Di,Vi,Ni),i=1:b
%                Di = distance value for bin i
%                Vi = variogram value at distance Di
%                Ni = number of distinct pairs in bin i

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

%Compute matrix, D, of pairwise distances

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

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

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

XX = X*U' - U*X' ;

YY = Y*U' - U*Y' ;

D = sqrt(XX.^2 + YY.^2) ;

D = tril(D) ;  %take lower triangular part

%Compute corresponding (Zi - Zj)^2 values

Z = M(:,3) ;

Z = (Z*U' - U*Z').^2 ;

%Now reduce to nonzero parts

Z = Z(find(D)) ;

D = nonzeros(D) ;

Dmax = max(D) ;

Dmin = min(D) ;

%Grouping into bins

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

%Cumulate Values

C = zeros(b,2) ; %cumulation matrix

N = length(D) 

U = ones(N,1) ;

% Handle 'Dmin' case

I = find(D==Dmin) ;

C(1,1) = C(1,1) + U(I)'*U(I) ;

C(:,2) = C(:,2) + U(I)'*Z(I) ;

%All other cases

for h=1:b
   
  Low = Dmin + ((h-1)/b)*(d - Dmin) ;
         
  High = Dmin + (h/b)*(d - Dmin) ;
 
  I = find((D>Low)&(D<=High)) ;
  
  C(h,1) = C(h,1) + U(I)'*U(I) ;
  
  C(h,2) = C(h,2) + U(I)'*Z(I) ;
  
end

  
%Construct variogram output

DAT = zeros(b,3) ;

for k = 1:b
   
   DAT(k,1) = Dmin +(((2*k)-1)/(2*b))*(d - Dmin);
   
   if C(k,1) > 0
      
      DAT(k,2) = C(k,2)/(2*C(k,1)) ; %(1/2) is for semivariogram.
      
   else
      
      D(k,2) = 0 ;
      
   end
   
   DAT(k,3) = C(k,1) ;
   
end