% This script computes empirical variograms.
% It requires AREA, a , to be given along with 
% a data file containing locations: loc(i) (Xi, Yi) 

a = 325844 ;  % area input

[m,n] = size(loc); %loc = data matrix

i = 1;

D = zeros(m*(m-1)/2, 3) ;

% Start procedure.

i = 1 ;

k = 1 ;  %counter for D.

Dmax = 0 ;

while i < m
   
   j = i + 1 ;
   
   while j <= m
      
      D(k,1) = i ;
      
      D(k,2) = j ;
      
      D(k,3) = norm(loc(i,:) - loc(j,:));
      
      if D(k,3) > Dmax 
         
         Dmax = D(k,3);
         
      end
      
      j = j + 1 ;
      
      k = k + 1 ;
      
   end
   
   i = i + 1 ;
   
end

Dist = D(:,3) ;

i = 1 ;

ind = 1 ;

K1 = [1:10] ;

K2 = K1 - K1 ;

while i <= length(Dist)
   
   if Dist(i) == 0   %avoids identical locations
      
      ind = 1 ;
      
   else
      
      ind  = ceil(10 * (Dist(i)/Dmax));
      
   end
   
   K2(ind) = K2(ind) + 2 ;  %Here we add both ij and ji
   
   i = i + 1 ;
   
end

K3 = K1 - K1 ;      

K3(1) = K2(1) ;

nn = 2 ;

while nn <= length(K3)
   
   K3(nn) = K3(nn-1) + K2(nn) ;
   
   nn = nn + 1 ;
   
end

K = K1 - K1 ;

i = 1 ;

while i <= length(K)
   
   K(i) = K3(i) .* (a/(m.^2));
   
   i = i + 1 ;
end

i = 1;

L = K1 - K1 ;

while i <= length(K)

   L(i) = sqrt(K(i)/pi) - Dmax *(i/10) ;
   
   i = i + 1 ;
   
end