function [tp,ep]=krig(mod_opt,xc,yc,x,y,var,krig_para,model,model_para)
% krig function returns the kriged results at (xc,yc) using observed
% data at var(x,y). The kriging resolution is specified by res.
% Input Parameters:
% mod_opt:  determine what kind of kriging model is selected:   
%       mod_opt=1:     object analysis
%       mod_opt=2:     ordinary kriging -- second order stationary case
%       mod_opt=3;     kriging  -- universal kriging with linear drift w/o revision
% xc, yc:    coordinates vectors for specified kriging points  
%            if the kriging grid is N x M, the their length is N*M.
%    x,y:    vectors with the coordinates of the observation points
%    var:    vector with the measurements of the variable to be mapped
%
%% kriging parameters
% krig_para = [blk_opt1,10,range2,res,kmin,kmax,eplim,rot_ang,aspect];
% blk_opt:  option for whether using station to station kriging ...
%            or station - block kriging     
%       blk_opt=1:     station to station (excluding the observation itself)
%       blk_opt=2:     block to station
%     nb:    no. of elements in a block, i.e, block size = nb * nb             
%  range:    search radius for kriging
%    res:    kriging cell resolution 
%   kmin:    minimum observed data points used for kriging, > 3
%   kmax:    maximum observed data points used for kriging
%            setting a reasonable kmax can speed up computations
%  eplim:    normalized kriging error limit
% anisotropy_para =[rot_ang aspect]
%  rot_ang:  rotate angle for anisotropy
%  aspect:   take into account anisotropy  scale(y)/scale(x)
%  model:    index of variogram or covariance model
%  model_para: variogram or covariance model parameters
%
% Output parameters:
%
% tp:    the estimate of the measured variable at the grid-points
% ep:    the estimation variance for tp, divided by the total variance
% 
% By Dezhang Chu, Woods Hole Oceanographic Institution
% Last modification  10 August, 1998
% %
% Reference:
%  A. J. Journel and Ch. J. Huijbregts, Mining Geostatistics
%  Academic Press, New York, pp. 304-313, 1978. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%  Kriging Software Package  version 1.0 beta,   October 29, 1998
%%  Copyright (c) 1998, property of Dezhang Chu and Woods Hole Oceanographic
%%  Institution.  All Rights Reserved.

global Max_Memory

dbstop if error
blk_opt=krig_para(1);
nb=krig_para(2);
range=krig_para(3);
res=krig_para(4);
kmin=krig_para(5);
kmax=krig_para(6);
eplim=krig_para(7);
tic

n=length(x);		% n = number of data points
ngx=length(xc);
ngy=length(yc);
if n*n*8 > Max_Memory
   ninc= round(n/1500);
   [x,y,var]=datareduction(x,y,var,ninc,1);
   disp(sprintf('size reduction factor = %g',ninc));
   n=length(x);
elseif ngx*ngy*8 > Max_Memory
   disp(sprintf(' Grid points = %g', ngx*ngy))
end
   
if length(krig_para) == 9
  anisotropy_para=krig_para(8:9);
% dealing with anisotropy
  ang=anisotropy_para(1);
  aspect=anisotropy_para(2);
  if abs(ang) > eps & abs(1-aspect) > 1e-5
%% observations
    cz=[x(:) y(:)];
    rx=1;
    ry=aspect;
    [cz,rot]=trans(cz,ang,rx,ry);
    x=cz(:,1);y=cz(:,2);
%% predictions
    cz=[xc(:) yc(:)];
    [cz,rot]=trans(cz,ang,rx,ry);
    xc=cz(:,1);yc=cz(:,2);
  end
end

if kmin < 3 
   kmin=3;
end
   
x=reshape(x,1,n); 	% makes sure that x is row vector
y=reshape(y,1,n);
nv=length(xc);				% nv = number of grid points
xc=reshape(xc,1,nv);
yc=reshape(yc,1,nv);
   
var=reshape(var,n,1); 	% and t is column vector

if model > 0				% variogram
   model_type=2;
else
   model_type=1;			% covariance
end

if mod_opt == 3 			% for universal kriging with linear drift 
  delx=-(xc(ones(n,1),:)'-x(ones(nv,1),:));
  dely=-(yc(ones(n,1),:)'-y(ones(nv,1),:));
end
pad=zeros(1,3);		% needed to pad K in the main routine


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if blk_opt == 2 & model < 0
   model_para1=[0 1 model_para(3:5)];
% The covariance between a grid block and itself:
   CVV=blk2blk(0,0,0,0,nb,res,model,model_para1); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


NaNcnt=0;
Anomaly_cnt=0;
ep=zeros(nv,1);
tp=zeros(nv,1);
rindx=floor(range/res)+1;
val_m=max(abs(var));

% Ordinary Kriging routine
hkJob=findobj(gcf,'Tag','JobPercent');
nv_cnt=round(0.05*nv);
for j=1:nv				% kriging through gridded points
  if rem(j,nv_cnt) == 0 
     percent=floor(100*j/nv);
     JobStr=sprintf('   %g percent done in %g sec.',percent,round(toc));
     set(hkJob,'String',JobStr);
  end
  difx=xc(j)-x;
  dify=yc(j)-y;
  r=sqrt( difx.*difx+dify.*dify);  % grid to station
  if kmin >= n
     indx_sort=1:n;
     k=1:n;nk=n;
  else
     [r_sort, indx_sort]=sort(r);
     ind=find(r_sort <= range);
     nd=length(ind);
     if isempty(ind)
        k=1:min(kmin,n);
     elseif ind(nd) >= kmin & ind(nd) < kmax 
        k = 1:ind(nd);
     elseif ind(nd) >= kmax
        k = 1:kmax;
     elseif ind(nd) < kmin
        k = 1:min(kmin,n);
     end
     nk=length(k);
  end
  %% station to grid 
  if blk_opt == 2 & model < 0        
     M20=sta2blk(x(indx_sort(k)),y(indx_sort(k)),xc(j),yc(j),nb,res,model,model_para);
  else
     M20=variogrammodel(model,r(indx_sort(k)),model_para);
  end
 %%%% model selection %%%%%
  switch mod_opt
      case 1		% Objective mapping, follow Journel and Huijbregts, p. 307
        M2=M20';
      case 2		% Ordinary Kriging, follow Journel and Huijbregts, p. 307
        M2=[M20 1 ]';
      case 3		% Universal Kriging w/ Linear drift, follow Journel and Huijbregts, p. 319
        M2=[M20 1 0 0]';
  end
   
% Construct the coefficient matrix and its inverse matrix
% to avoid repeated computation in the loop
  if  j == 1 | kmin < n 
    rs=sqrt(((x(ones(length(k),1),indx_sort(k))' ...
          -x(ones(length(k),1),indx_sort(k))).^2+...
          (y(ones(length(k),1),indx_sort(k))' ...
          -y(ones(length(k),1),indx_sort(k))).^2));
    K0=variogrammodel(model,rs,model_para);
    if model_type == 1
       K0(1:nk+1:nk^2)=ones(nk,1);	   		% diagonal = 1
    else
       K0(1:nk+1:nk^2)=zeros(nk,1);	   		% diagonal = 0
    end
    K=[K0 ones(nk,1);ones(1,nk) 0];

      switch mod_opt
      case 1		% Objective mapping, follow Journel and Huijbregts, p. 307
        K=K0;
      case 2		% Ordinary Kriging, follow Journel and Huijbregts, p. 307
        K=[K0 ones(nk,1);ones(1,nk) 0];
      case 3		% Universal kriging w/ linear drift, follow Journel and Huijbregts, p. 319
        K=[K0 ones(nk,1) x(indx_sort(k))' y(indx_sort(k))'; ...
               ones(1,nk) pad; delx(j,indx_sort(k)) pad; dely(j,indx_sort(k)) pad];
    end
    K_inv=inv(K);
  end			
  lambda=K_inv*M2;
  indx=find(isnan(M2));
  if length(indx) > 0 
     disp(sprintf('nans in lambda at j=%g ',j));
  end
  lnan=find(isnan(var(indx_sort(k))));
  if length(lnan) > 0
     var(k(lnan))=zeros(length(lnan),1);
  end
  tp(j)=sum_nan(lambda(1:nk).*var(indx_sort(k)));
  if blk_opt == 2 & model < 0
     ep(j)=CVV-sum_nan(lambda.*M2);	% Error estimate for tp(j)
  else
     if model_type == 1
        ep(j)=1-sum_nan(lambda.*M2);
     else
        ep(j)=sum_nan(lambda.*M2);
     end
  end
  if abs(tp(j)) > 3.0*val_m  
    Anomaly_cnt=Anomaly_cnt+1;
    tp(j)=val_m;
  end
  if abs(ep(j)) > eplim 
    tp(j)=NaN;
    ep(j)=abs(ep(j));
  end
end

if NaNcnt > 1 | Anomaly_cnt > 1
   disp(sprintf('NaNcnt = %g, Anomaly_cnt = %g',NaNcnt,Anomaly_cnt))
end
JobStr=sprintf('   %g percent done in %g sec.',100,round(toc));
set(hkJob,'String',JobStr);

