function C = covar_fcns(m, z1, z2, nug)
% COV_BC : return covariance models
% for use in kriging with no-flow boundary conditions
 
% Copyright (c) 2010 Kristopher L. Kuhlman (klkuhlm at sandia dot gov)
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
% 
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
%
% this function is part of the set of MATLAB scripts 
% comprising the implementation of a universal kriging
% program that computes estimates of potential via
% cokriging, taking into account boundary conditions.
%
% Kuhlman, K.L., and E. Pardo-Igúzquiza, 2010. Universal
% cokriging of hydraulic heads accounting for boundary 
% conditions, Journal of Hydrology, 384(1–2), 14–25.
% http://dx.doi.org/10.1016/j.jhydrol.2010.01.002

% m       : a structure containing details of model
% z1,z2     : row vector {x,y} locations of points
% m.a     : range parameter 
% m.sigsq : variance or sill 
% m.par   : a parameter (depends on model)
% m.type  : model type
% nug     : appropriate nugget to use

a = m.a;
ss = m.sigsq;

if ~(m.type == 1 || m.type == 2 || m.type == 3)
    error('incorrect choice of model (1=Gaussian, 2=gen Cauchy, or 3=Matern)')
end

if (a <= 0 || ss <= 0 || nug < 0)
    error('range, variance, and nugget must all be positive')
end

% vector pointing from z2 to z1
h = z1-z2;
% its length 
lag = sqrt( squeeze(dot(h,h,2)) );
% normalized by variogram range
ha = lag./a;

if m.type == 1
    
    % NB: par not used in Gaussian
    
    % Gaussian model
    C = ss.*exp(-ha.^2);
    
elseif m.type == 2

    % par is exponent
    if m.alpha <= 0.0
        error('Cauchy exponent must be positive')
    end

    % Generalized Cauchy model
    C = ss.*(1 + ha.^2).^(-m.alpha);

elseif m.type == 3
    
    % par is smoothness parameter
    if m.nu < 2
        error('Matern smoothness parameter too small, must be >=2')
    end

    % Matern model
    C = (ss*2^(1-m.nu)/gamma(m.nu)).*ha.^m.nu.* ...
          besselk(m.nu,ha);
end

% apply nugget when h is close to 0 (numerically equal)
C(lag < eps) = ss + nug;
