function OUT = var_sim(L,Z,d)

%VAR_SIM estimates a spherical variogram for a collection of k data sets
%using VAR_SPHER.

%Written by: TONY E. SMITH, 2/21/98

%Functions called: variogram.m, var_fit.m, spher.m

%INPUTS:  L = (nx2)-matrix of locations (Xi,Yi), i = 1:n,
%         Z = (nxk)-matrix of k data sets (Zi) for locations in L.
%         d = (optional) bandwidth
%
%OUT = (kx4)-matrix [ri,si,ai,vi] with (ri,si,ai) denoting the 
%      variogram parameter estimates, and vi the sample variance 
%      estimate for data set i.

[n,k] = size(Z);

OUT = zeros(k,3);

M(:,1:2) = L;

index = 10;

for i = 1:k    
   
   M(:,3) = Z(:,i); 
   
   if nargin == 2
       
       [DAT,bw] = variogram(M) ;
       
   else
       
       [DAT,bw] = variogram(M,d) ;
       
   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) ;
    
	OUT(i,:) = res.param;
    
      
   %screen counter
   
   if 100*(i/k) >= index
      
      percent_done = index
      
      index = index + 10;
      
   end %end simulation loop 
      
end

%compute sample variance (with normalizer = n-1)
   
v = (std(Z).^2)';
   
OUT = [OUT,v];