function results = var_fit(y,D,G)

%VAR_FIT = least-squares fit of spherical variogram.

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

%MODIFIED with Truncation at max(D)/2, 2/23/04
%
%Functions called: err_fit.m
%
%INPUTS:
%
%  (i)   y = [ r0 s0 n0] = initial values of [r s n]
%  (ii)  D = vector of distance values
%  (iii) G = vector of empirical variogram values
%  
%OUTPUTS: results.param = [r,s,a] = estimated range, sill, and 
%                  nugget values of the spherical variogram.
%         results.iter  = number of iterations (max = 600)

u = ones(length(D),1) ;

% [p,opt] = fminsearch('err_fit',y,[],u,D,G);
% 
% p = abs(p) ;

timeo = clock;

[p,err,exitflag,output] = fminsearch('err_fit',y,[],u,D,G);

time_fminsearch = etime(clock,timeo); %check if too long

%*********************************
%p(3) = max([0 p(3)]); %Make nonegative
%************************************

%*********************************
%Check for nonnegative nugget

if p(3) < 0
          
    y0 = y(1:2);
    
    [p0,err,exitflag,output] = fminsearch('err_fit_zeroed',y0,[],u,D,G);
    
    p = [p0,0];
    
	disp(' ');
	disp('*******************************************************');
	disp(' ');
	disp('New estimation done with NUGGET truncated to ZERO');
	disp(' ');
	disp('*******************************************************');
    
end   
    
%************************************

%Check for failure of convergence

if exitflag == 0 
    
	fprintf(1,'VAR_FIT: convergence not obtained in %3d iterations \n',output.iterations);
	
	y0 = y(2:3);
	
	[p0,err,exitflag,output] = fminsearch('err_fit_truncated',y0,[],u,D,G); %re-solves with range = max(D)
	
	p = [max(D),p0];
	
	disp(' ');
	disp('*******************************************************');
	disp(' ');
	disp('New estimation done with RANGE truncated to MAX DISTANCE');
	disp(' ');
	disp('(Depending on data, you may wish to increase max distance)');
	disp(' ');
	disp('*******************************************************');

end

%Check for Range larger that max(D)


if p(1) > max(D) 
    
	y0 = y(2:3);
	
	[p0,err,exitflag,output] = fminsearch('err_fit_truncated',y0,[],u,D,G); %re-solves with range = max(D)
	
	p = [max(D),p0];
	
	disp(' ');
	disp('*******************************************************');
	disp(' ');
	disp('New estimation done with RANGE truncated to MAX DISTANCE');
	disp(' ');
	disp('(Depending on data, you may wish to increase max distance)');
	disp(' ');
	disp('*******************************************************');

end

%********************************
%p(3) = max([0 p(3)]); %Make nonegative
%*********************************

if p(3) < 0  %Check for negative nugget
    
    %re-solves with [r,n] =[max(D),0]
    
    r0 = max(D);
    
    s0 = y(2);
    
    [p0,err,exitflag,output] = fminsearch('err_fit_degenerate',s0,[],r0,u,D,G); 

    p = [r0,p0,0];
    
    disp(' ');
	disp('************************************************************');
	disp(' ');
	disp('Very bad fit. RANGE set to MAX DISTANCE and NUGGET set to ZERO');
	disp(' ');
	disp('************************************************************');

end;

results.iter = output.iterations;


%Finally, be sure that the nugget is not bigger than the sill.
%(This is Cressie's 'relative sill' approach)

p(2) = p(3) + abs(p(2) - p(3)) ;

results.param = p ;

results.iter = output.iterations; %Number of iterations used in FMINS