function OUT = var_spher_plot(M,opts)

%VAR_SPHER_PLOT estimates a spherical variogram for a given data set
%using weighted nonlinear least squares.

%Written by: TONY E. SMITH, 2/9/01

% Functions called: variogram.m, var_fit.m, spher.m, var_plot.m
%
% INPUTS:  M = (nx3)-matrix of locations and data 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)
%
% SCREEN OUTPUT: Range, Sill, Nugget, plus Number of Iterations in fitting.
%
%              NOTE: If Iterations exceed 600 you will get an error
%              message telling you that the algorthm has failed to converge
%              after 600 iterations (which is the default MAX ITERATIONS)
%
% LOT OUTPUT: Plot of Variogram plus Spherical Fit
%
% OUTPUTS: OUT = cell(6,1) with
%         (i)   OUT{1}= [r,s,a] (= estimated [range,sill,nugget])
%         (ii)  OUT{2}= lag distance values, (Di, i=1:b) 
%         (iii) OUT{3}= empirical variogram values Gi at Di
%         (iv)  OUT{4}= spherical values Vi at Di
%         (v)   OUT{5}= number of iterations in FMINS
%         (vi)  OUT{6}= Max Distance used in variogram estimates



if nargin > 1
   
       
	[DAT,d] = variogram(M,opts) ;

else
   
   [DAT,d] = variogram(M) ;

end

D = DAT(:,1) ;

G = DAT(:,2) ;

%N = DAT(:,3);

y = zeros(1,3) ;

y(1) = max(D) ;      %initial estimate of range = max(D)

y(2) = std(M(:,3))^2 ;   %initial estimate of sill = sample variance

y(3) = max([0,min(G)]) ; %initial estimate of nugget

res = var_fit(y,D,G) ;

p = res.param;

iter = res.iter ;


OUT = cell(6,1) ;

OUT{1} = p ;

OUT{2} = D ;

OUT{3} = G ;

OUT{4} = spher(p,D) ;

OUT{5} = iter ;

OUT{6} = d;

var_plot(OUT);

%At this point we check to see if variogram is flat

if p(2) - p(3) < 10^(-5)
    
    p(1) = 0;
    
    disp(' ');
	disp('*************************************************');
	disp(' ');
	disp('VARIOGRAM is essentially FLAT. RANGE set to ZERO');
	disp(' ');
	disp('**************************************************');
    
    OUT{1} = p ;
    
end


%Make Screen Output


%SCREEN OUTPUT

info.fmt = '%12.3f'; %print format

   
   info.rnames = strvcat(' ','RANGE','SILL','NUGGET','ITERATIONS');
  
   mat(:,1) = [p';iter];
   
disp(' ');
disp('SPHERICAL VARIOGRAM:');
mprint_new(mat,info) %Display Regression Results

disp(['MAXDIST = ',num2str(OUT{6})]);


