function OUT = var_spher_wtd_plot(M,varargin)

% VAR_SPHER_WTD estimates a spherical variogram parameters, p = (r,s,a), 
% for a given data set using weighted nonlinear least squares. This is an 
% approximation of GLS in which the covariance matrix between estimators
% is assumed to be diagonal and the variance of each estimator, Gi, at distance, 
% Di, with bin size, Ni, and theoretical value, Vi(p) = spher(Di;p), is estimated by
% 
% var(Gi;p) = (Vi(p)^2)/Ni
% 
% Hence the weighted least squares estimator is given by the p that minimizes:
% 
% err = sum[ ((Gi - Vi(p))^2)/var(Gi;p) ]
%
%Written by: TONY E. SMITH, 9/7/04

%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
%         varargin = optional bandwidth specification (if not 
%                    specified, the default value is Dmax/2).
%
%SCREEN OUTPUT: Range, Sill, Nugget, plus Number of Iterations in fitting.

%PLOT 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}= 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}= Bandwidth used in variogram estimates


if nargin > 1
   
   d = varargin{:} ;
      
	[DAT,d] = variogram(M,d) ;

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_wtd(y,D,G,N) ;

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


%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

VARIOGRAM_BANDWIDTH = OUT{6}