function OUT = o_krige(y,L0,L,h,p)

% O_KRIGE carries out Ordinary Kriging by first fitting a spherical covariogram 
% to data (LO,y), and then using the estimated parameters to construct kriging 
% predictions at a set of locations, L. [The option to use a prespecified set of 
% range, sill and nugget parameters is also available. In this case, no variogram
% will be estimated.]    

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

%Functions called: distance_mat.m, var_spher.m, spher_mat.m, spher.m
%
%INPUTS: (i)   y  = (n,1) vector of dependent variable values
%                   at n sample locations.
%        (ii)  L0 = (n,2) matrix of sample location coordinates.
%        (iii) L  = (N,2) matrix of prediction location coordinates.
%        (iv)  h  = bandwidth for Kriging prediction.
%                   NOTE: If there are less than 3 points within distance, h,
%                   of a given point, then the three closest points are used.
%         (v)  p  = optional variogram parameters
%
%OUTPUTS: OUT = (4x1) cell with
%         OUT{1} = estimated mean (regression constant)
%         OUT{2} = (3,1) vector of estimated variogram parameters
%                  [r,s,a]
%         OUT{3} = N-vector of kriged values at locations in L
%         OUT{4} = N-vector of standard errors at locations in L

%Error Check

dbstop if error;

%Initial Settings

n = length(L0) ;

U = ones(n,1) ;

M(:,1:2) = L0 ;

D = distance_mat(L0) ;

M(:,3) = y ;

disp(' ');
disp('Beginning covariogram estimation');
disp(' ');

if nargin < 5    
       
    res = var_spher(M) ;
    
    p = res.param ; %parameter vector for variogram
    
end    

clear('M') ;

%Construction of Covariance Matrix
	
C = p(2)*U*U' - spher_mat(p,D) + p(3)*eye(n) ;

%Check rank of C

if rank(C) < size(C,1) %not of full rank
    
    C = C + .0000001*eye(n); %add to diag (like ridge regression)
    
end
   
clear('D') ;

A = chol(C)';
        
AA = A\eye(n) ; 

mu = AA*U\AA*y ; % same as (U'*inv(C)*y)/(U'*inv(C)*U) ;
   
%Note: mu and p are now defined for use in Kriging

%Compute Kriged Values and Standard Errors

disp(' ');
disp('Begin Kriging');
disp(' ');

N = length(L(:,1)) ;

KRIGE = zeros(N,1) ;

STD_ERR = zeros(N,1) ;


RES = y - U*mu ;

for i = 1:N   
    
%     if i == 384
%         
%         i = i;
%         
%     end
     
   %Identify points within distance h of location i
   
   X1 = L(i,1)* U - L0(:,1) ;
   
   X2 = L(i,2)* U - L0(:,2) ;
   
   DD = sqrt(X1.^2 + X2.^2) ;
   
   pos = max(zeros(n,1),(h*U - DD)) ; 
   
   K = find(pos) ;  %List of nonzero positions in DD
   
   if length(K) < 3
      
      %find three closest points
      
      list(:,1) = [1:n]' ;
      
      list(:,2) = DD(:) ;
      
      ngh = sortrows(list,2) ; 
      
      K = ngh(1:3,1) ;
      
      
   end
                
      
   %Compute local covariances
   
   CC = inv(C(K,K)) ; %inverse of local cov matrix
   
   d = DD(K) ; %vector of distances to i
   
   c = p(2)*U(K) - spher(p,d); %local cov vector
   
   %Correct c at any zero-distance points   
   
   KK = length(K) ;
   
   for j = 1:KK
      
      if d(j) == 0 
         
         c(j) = c(j) + p(3) ;
         
      end
      
   end
   
   Z = U(K) ;
   
   KRIGE(i) = mu + (c' + ((1 - Z'*CC*c)/(Z'*CC*Z)))*CC*RES(K) ;   
    
   Z = U(K) ;
   
   ERR_VAR = p(2) - c'*CC*c + ((1 - Z'*CC*c)^2)/(Z'*CC*Z) ;
   
   STD_ERR(i) = sqrt(ERR_VAR) ;   
        
end


%Prepare Output

OUT{1} = mu ;
OUT{2} = p ;
OUT{3} = KRIGE ;
OUT{4} = STD_ERR;