function [DAT,maxdist,bin_size,bin_last] = variogram_plot(M,opts);

%VARIOGRAM computes the empirical (semi)variogram for a given
%set of spatial data. Bins are determined internally.[If you want 
%to increase bin size, go to line 133 and reset 'binMax'].

%Written by: TONY E. SMITH, 2/17/98
%
%INPUTS:  M = (nx3)-matrix of locations and values,
%              M(i) = (Xi,Yi,Zi), i = 1:n
%      opts = (optional) structure of additional settings
%              opts.maxdist = (approximate) maximum distance
%                             to be considered. (Default = Dmax/2)
%              opts.binmax = maximum number of bins to be used.
%                            (Default = 100)
%
%OUTPUTS: DAT = (bX2)-matrix with OUT(i)=(Di,Vi),i=1:b
%                Li = lag distance value for bin i
%                Vi = variogram value at distance Di
%              maxdist = Maximum Distance used for variogram
%              bin_size = (constant) number of distinct pairs 
%                         in all but last bin
%              bin_last = number of points in last bin

%See also VARIOGRAM_PLOT

dbstop if error;

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) ;

clear('X','Y','XX','YY') 

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

%Now reduce to nonzero distances and determine range

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

clear('D');

% Eliminate zero-distance pairs

Rows = find(S(:,3) > 0);

S = S(Rows,:);


% The matrix S has rows (i,j,d) containing the row,
% column, and values of all positive distances 

Dmin = min(S(:,3)) ;

Dmax = max(S(:,3)) ;  

dist_flag = 0;

bin_flag = 0;

if nargin == 2 % parse user input options
    
    fields = fieldnames(opts);
    
    nf = length(fields);
    
    for i=1:nf
        
        if strcmp(fields{i},'maxdist')            
            
            dist_flag = 1;
            
        elseif strcmp(fields{i},'binmax')
            
            bin_flag = 1;           
        
        end;
        
     end;
     
 end; %End Parsing


% Define maxdist


if dist_flag == 1
    
    maxdist = opts.maxdist;
    
else
    
    maxdist = Dmax/2;
    
end

%Define binmax

if bin_flag == 1
    
    binmax = opts.binmax;
    
else
    
    binmax = 100;
    
end

% Check maxdist
   
if Dmin >= maxdist 

  error('Dmin >= maxdist') ; 
  
end

% Now decide on number of bins required, bin_num

% There are two cases to be considered, depending on whether 
% the binmax constraint is binding. First let

S = sortrows(S,3); %sort by distance

n_S = size(S,1);

Rows = find(S(:,3) <= maxdist);

n0_S = sum(S(:,3) <= maxdist); %Number of distances <= maxdist 

bin_size = ceil(n0_S/binmax);

if bin_size >= 30 % Binmax constraint is in effect
    
    N = bin_size*binmax;
    
    bin_num = binmax;
    
else % must use 30 in each bin
    
    bin_size = 30;
    
    bin_num = ceil(n0_S/binsize);
    
    N = bin_num*bin_size;
    
end

% Now    bin_num = number of bins
%       bin_size = number of point-pairs in each bin
%              N = maximum number of point-pairs to be used
    
DAT = zeros(bin_num,2);

for i = 1:bin_num-1 %make all but last bin
    
    n0 = (i-1)*bin_size + 1;
    
    n1 = i*bin_size;
    
    S1 = S([n0:n1],1);
    
    S2 = S([n0:n1],2);
    
    Di = S([n0:n1],3);
    
    Li = mean(Di); %Lag distance i = average bin distance
    
    Z = M(:,3);
    
    Vi = sum( (Z(S1) - Z(S2)).^2 )/(2*bin_size); %Variogram estimate
    
    DAT(i,1) = Li;
    
    DAT(i,2) = Vi;
    
end

% Now make last bin

n0 = (bin_num -1)*bin_size +1;

n1 = bin_num*bin_size;

if n1 > n_S  %Last bin will be smaller than bin_size
    
    n1 = n_S;
    
end
    
S1 = S([n0:n1],1);
    
S2 = S([n0:n1],2);

D_last = S([n0:n1],3);

bin_last = length(D_last); % Number of points in last bin

L_last = mean(D_last); %Lag distance i = average bin distance

Z = M(:,3);

V_last = sum( (Z(S1) - Z(S2)).^2 )/(2*bin_last); %Variogram estimate

DAT(end,1) = L_last;

DAT(end,2) = V_last;

Lmax = L_last;

disp(['MAXDIST = ', num2str(maxdist)]);

Dist = DAT(:,1);  % distance values for bins
 
Val = DAT(:,2); %variogram values for bins
 
plot(Dist,Val,'.');