function C = covar_derivs(m, deriv, z1, z2, e1, e2)
% COV_BC : return covariance or its first 2 directional
% derivatives, 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     : structure containing model parameters
% m.type: model (1=gaussian, 2=gen Cauchy, 3=Matern)
% deriv  : integer flag {0,1,2} - zeroth, first, or second derivative
% z1,z2 : row vector {x,y} locations of points
% e1,e2 : unit row vectors indicating orientation of no-flow boundary
% m.a     : range parameter in Gaussian model
% m.sigsq : variance or sill 
% m.nug   : appropriate nugget term
% m.par   : a model-specific parameter

if deriv == 0
    % call functions from other m-file 
    % error-checking and setup done there
    C = covar_fcns(m,z1,z2,m.nug);
else
    
    % ensure required unit vectors are unit length
    l = sqrt(dot(e1,e1,2));
    e1 = e1./l(:,[1,1]);
    if deriv == 2
        l = sqrt(dot(e2,e2,2));
        e2 = e2./l(:,[1,1]);
    end
    
    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 ~(deriv == 0 || deriv == 1 || deriv == 2)
        error('incorrect model type (0=head covar, 1=head/flux x-covar, 2=flux covar)')
    end

    if (a <= 0 || ss <= 0 || m.nug < 0)
        error('range, variance, and nugget must all be positive')
    end 
    
    % a row vector pointing from z2 to z1
    h = z1-z2;  
    % just the length (a scalar)  
    lagsq = squeeze(dot(h,h,2));
    lag = sqrt( lagsq );
    
    % normalized by variogram range
    ha = lag./a;
    % create a unit vector version of h
    hunit = h./lag(:,[1,1]);
    
    % projection of lag vector onto first unit vector
    % used by most first and second derivatives
    hproj = dot(hunit,e1,2);
    
    if deriv == 2
        % magnitude of projection onto e1 * same for e2
        hprojsq = hproj .* dot(hunit,e2,2);
    end
    
    % ******************************************
    % Gaussian
    if m.type == 1
        
        % NB: par not used for Gaussian
        
        % negative first directional derivative
        % cross-covariance between head & gradient measurements
        % equation 34 in Pardo-Iguzquiza & Chica-Olmo WRR 2004
        if deriv == 1
            % when h->0 C* = 0 :: no problem evaluating for zero lag
                   
            C = (2*lag*hproj*ss)./a^2 .* exp(-ha.^2);
        
        % negative second directional derivative
        % covariance between gradient measurements
        elseif deriv == 2
            % when h->0, C** = -2*sigsq/a^2 :: no problem for zero lag
           
            C = 2*ss./a^2 .*(1 - 2*hprojsq.*ha.^2) .* exp(-ha.^2);
        end

    % ******************************************
    % generalized Cauchy model
    elseif m.type == 2

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

        % negative first directional derivative
        % cross-covariance between head & gradient measurements
        if deriv == 1
            % when h->0 C* = 0 :: no problem for zero lag
        
        C = (2*ss*m.alpha/a^2).* lag .* hproj.*(1 + ha.^2).^(-m.alpha - 1);

        % negative second directional derivative 
        % covariance between gradient measurements
        elseif deriv == 2
            % when h->0, C** = -2*p*sigsq/a^2 :: no problem for zero lag
         
            C = -(2*ss*m.alpha/a^2).*((m.alpha + 1).* 2.*hprojsq.*ha.^2* ...
                (1 + ha.^2).^(-m.alpha - 2) - (1 + ha.^2).^(-m.alpha - 1));
        end
    
    % ******************************************
    % Matern model: need to deal with divide-by-zero errors gracefully
    % Bessel function -> +ifinity at r=0
    elseif m.type == 3
    
        if m.nu < 2
            error('Matern smoothness parameter too small, must be >=2')
        end
        
        hag = ha(lag > eps);
        common = ss.*2^(1 - m.nu)./gamma(m.nu).*hag.^m.nu;
        
        % negative first directional derivative
        % cross-covariance between head & gradient measurements
        if deriv == 1
            % when h->0 C* = 0 :: numerical problem evaluating at zero lag
            
            hpg = hproj(lag > eps);
        
            C(lag > eps) = hpg.*common./a.*besselk(m.nu-1,hag);
            C(lag <= eps) = 0.0;  % specify manually
        
        % negative second directional derivative 
        % covariance between gradient measurements
        elseif deriv == 2
            
            hpsqg = hprojsq(lag > eps);
            cc = ss.*2^(1 - m.nu)./gamma(m.nu).*hag.^(m.nu - 2);
            
            C(lag > eps) = cc./a^2.* ...
                (besselk(m.nu-1,hag).*hag - besselk(m.nu-2,hag).*hpsqg.*hag.^2);

            % when h->0 C** = sigsq*Gamma[n+1]/(2*a^2*(n-1)*Gamma[n])
            % found using Mathematica 
            
            % the two Bessel functions in the solution become singular
            % at the origin; use this limit instead
            C(lag <= eps) = ss*gamma(1-m.nu)/ ...
                (2*a^2*(m.nu - 1)*gamma(m.nu));
        end
    end
end


