function [h,covar]=cov_ssrf_exp(dd,lag,Param,flag)
%  COV_SSRF_EXP calculates the explicit SSRF covariance function at the
%   limit of an infinite spectral cutoff.
%
% This routine was tested under Matlab R2011b 7.13.0.564 32-bit (win32)
%
%==========================================================================
%   Input parameters
%==========================================================================
%
%   dd:       Dimensionality (dd=1,2,3,)
%	lag:       Lag vector
%	Param:     Spartan parameters [eta0,eta1,xi]
%   flag:     If flag=0 the plots are suppressed
%
%==========================================================================
%   Output parameters
%==========================================================================
%
%	h       Normalized (by xi) lag vector
%	covar     Covariance vector
%
%	EXAMPLE
%==========================================================================
%
%	[h,covar]=cov_ssrf_exp(2,0:1:10,[1,2,3],1);
%
%  Copyright (C) 2013 Dionisis Hristopulos, 
%   email: dionisi@mred.tuc.gr
%
%  This program is free software; you can redistribute it and/or modify it
%  under the terms of the GNU General Public License as published by the
%  Free Software Foundation (version 2 of the License or  later version).
%
%  This program is distributed WITHOUT ANY WARRANTY.
%  See the GNU General Public License for more details.
%==========================================================================
%	REFERENCES
%==========================================================================
%
%   D. T. Hristopulos, Spartan Gibbs random field models for 
%   geostatistical applications, SIAM Journal of Scientific Computing,
%   24 (6), pp. 2125-2162 (2003).
%
%   D. T. Hristopulos and S. Elogne, Analytic properties and covariance 
%   functions for a new class of generalized Gibbs random fields, 
%   IEEE Transactions on Information Theory, 53(12), 4467-4679 (2007).
%
%   D. T. Hristopulos, Radial Covariance Functions Motivated by Spatial 
%   Random Field Models with Local Interactions, 2013, submitted to: 
%   Stochastic Environmental Research and Risk Assessment
%
%==========================================================================

eta0=Param(1);
eta1=Param(2);
xi=Param(3);

%   Calculate the normalized lag
h=lag/xi;

%    Replace zero lag to calculate the De l'Hospital limit
h(lag==0)=eps;

switch dd
    
    case 1
        
        if abs(eta1)<2
            beta1=0.5*(sqrt(2-eta1));  beta2=0.5*(sqrt(2+eta1));
            Gr=(exp(-h*beta2)).*( (cos(h*beta1))/beta2 +(sin(h*beta1))/beta1)/4;
        elseif  eta1==2
            Gr=((1+h).*exp(-h))/4;
        elseif eta1>2
            delta=sqrt(eta1^2-4); omega2=sqrt((eta1+delta)/2);
            omega1=sqrt((eta1-delta)/2);
            Gr=( (exp(-h*omega1))/omega1 - (exp(-h*omega2))/omega2 )/(2*delta);
        else
            disp('Permissibility conditions fall');
        end
        
    case 2
        
        if eta1==2
            Gr= h.* besselk(-1,h)/(4*pi);
            vv=1/(4*pi);
        elseif eta1 >2
            Dn=sqrt(eta1^2-4);
            zp= sqrt(- (- eta1 + Dn)/2);
            zm= sqrt(- (- eta1 - Dn)/2);
            Gr= ( besselk(0,h*zp)-besselk(0,h*zm))/(2*pi*Dn );
            vv=  log( (eta1+Dn)/(eta1-Dn))/ (4 *pi *Dn);
        elseif abs(eta1)<2
            Dn=sqrt(4-eta1^2);
            zp=sqrt(- (- eta1 + sqrt(-1)*Dn)/2);
            zm=sqrt( - (- eta1 - sqrt(-1)*Dn)/2);
            Gr=imag(besselk(0,h*zp))/ ( pi * Dn);
            vv= (pi/2 - atan(eta1/Dn)) / (2 *pi *Dn);
        else
            disp('Permissibility conditions fall');
            
        end
        
    case 3
        
        if abs(eta1)<2
            delta=sqrt(4-eta1^2);
            beta1=0.5*(sqrt(2-eta1));  beta2=0.5*(sqrt(2+eta1));
            Gr=(exp(-h*beta2)).*(sin(h*beta1)./(2*pi*h*delta));
        elseif  eta1==2
            Gr=exp(-h)/(8*pi);
        elseif eta1>2
            delta=sqrt(eta1^2-4); omega2=sqrt((eta1+delta)/2);
            omega1=sqrt((eta1-delta)/2);
            Gr=(exp(-h*omega1)-exp(-h*omega2))./(4*pi*h*delta);
        else
            disp('Permissibility conditions fall');
            
        end
        
    otherwise
        disp('Not defined for d>3');
        
end

if eta1 <= -2; covar=nan(size(h));
    
else
    covar=eta0*Gr;
    
    if nargin >3 && flag>0
        figure(1);
        plot(h,covar,'b-');
        if length(h)<20
            hold on; plot(h,covar,'r.');
        end;
        xlabel('h'); ylabel('SRRF Covariance');
        title(['d=',num2str(dd)]);
        axis tight; 
        figbolden
    end;
end
