function [z, s2] = pkrige(M, v, co, c, a, mod);
% [z, s2] = pkrige(M, v, co, c, a, mod);
% returns a punctual kriging estimate of the value z
% and the error variance s2 of an unknown point whose
% xy coordinates are given by the row vector v
% M is the n by 3 matrix of x,y,z values of the
% n control points. co, c and a are the parameters
% of the fitted semivariogram, and mod is the
% model used: mod = 1: spherical; mod = 2: exponential
% mod = 3: Gaussian.
% Written by Gerry Middleton, December 1996
[n col] = size(M);
X = [M(:,1:2); v];
H = zeros(n+1);       % calculated distance matrix (see Edist)
for i = 1:2
   X1 = X(:,i)*ones(1,n+1);
   DX = X1-X1';
   H = H + (DX .^2);
end;
H = sqrt(H);
Gam = zeros(n+1);
if mod == 1           % calculate semivariance matrix
   for i = 1:n+1
      for j = 1:n+1
         if i == j 
            Gam(i,j) = 0;
         elseif H(i,j) > a
            Gam(i,j) = co + c;
         else
            h = H(i,j)/a;
            Gam(i,j) =  c*1.5*h - 0.5*h^3;
         end
      end
    end
    fprintf('Semivariances calculated by spherical model\n');
elseif mod == 2
   Gam = co + c*(1 - exp(-3*abs(H/a)));
   fprintf('Semivariances calculated by the exponential model\n');
elseif mod == 3
   Gam = co + c*(1 - exp(-(sqrt(3)*(H/a)) .^2));
   fprintf('Semivariances calculated by the Gaussian model\n');
else
   error('Wrong entry for mod: should be 1,2, or 3');
end;
Gam                   % comment out if no printout required
gamp = Gam(:,n+1);    % vector gamma(p)
gamp(n+1) = 1.0;
Gam(:,n+1) = [ones(n,1); 0];  % modify Gam 
Gam(n+1,1:n) = ones(1,n);
w = Gam\gamp;         % weights, last term is lambda
z = w(1:n)'*M(:,3);   % predicted value
s2 = w'*gamp;
plot(M(:,1), M(:,2),'o');  % comment these lines out if
hold on;                   % no plot is required
plot(v(1), v(2), '+');
axis('equal');
dx = 0.02*(max(M(:,1)) - min(M(:,1)));
for i = 1:n
   text(M(i,1)+dx, M(i,2), num2str(M(i,3)));
end
text(v(1)+dx, v(2), num2str(z));
hold off