function [h,covar,err]=cov_ssrf_num(dd,lag,Param,flag)

%  COV_SSRF_NUM calculates the covariance function for Spartan spatial RFs
%   using numerical integration (adaptive Gauss-Kronrod quadrature
%   method). The spectral integral is given by Equation (10) in Hristopulos
%   and Elogne (2007).
%   The only difference is that the normalized lag (by the characeristic
%   length xi) is used herein, with the cooresponding
%   modification of the integrand due to the change of variables r-->
%   h=r/xi.
%
% This routine was tested under Matlab 6.5 R13
%
%==========================================================================
%   Input parameters
%==========================================================================
%
%   dd:        Number of spatial dimensions
%	lag:      Lag vector
%	Param:    Spartan parameters [eta0,eta1,xi,kmax]
%   flag:     If flag=0 plotting is suppressed
%
%==========================================================================
%   Output parameters
%==========================================================================
%
%	h:       Normalized (by xi) lag vector
%	covar:   Covariance matrix
%   err:     Integration error
%
%	EXAMPLE
%==========================================================================
%
%	[h,covar,err]=cov_ssrf_num(2,0:1:10,[1,2,3,2*pi],1);
%
%==========================================================================
%
%   FUNCTIONS USED
%
%   QUADGK: adaptive Gauss-Kronrod quadrature method (Matlab function)
%
%==========================================================================
%
%  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, Positive Semidefinite Functions and Karhunen-Loeve
%   Expansions Related to Spartan Random Fields, 2013, submitted to:
%   SIAM Journal on Uncertainty Quantification
%
%==========================================================================

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

if isnan(xi), disp('xi=nan'); return; end;

if eta1<=-2
    disp('Permissibility conditions fall');
    covar=nan(size(lag)); return;
elseif dd>3
    disp('Covariance not defined for d>3');
    covar=nan(size(lag)); return;
end;

bd=1/(2*pi)^dd;
dh=dd/2-1;
g0=sqrt(bd);

%   Define the wavenumber cutoff if undefined

if length(Param)==3, kmax=20/xi; flag=0; end;

%   Convert lags into normalized lag vector (used for matrix of lags)

if size(lag,1) >1 && size(lag,2) >1
    resh=1;
    h=lag(:)/xi;
    %    Replace zero lag to calculate the De l'Hospital limit
    h(h==0)=eps;
    
else
    h=lag/xi;
    resh=0;
    %    Replace zero lag to calculate the De l'Hospital limit
    h(lag==0)=eps;
end;

n=length(h);
%==========================================================================
%   Older definition of the integrand of the spectral representation using
%   inline functions
%==========================================================================
% %   Define the integrand of the spectral density in the spectral
% %   representation (Schoenberg, 1939;
%   P1(1)=d; P1(2)=dh;
%   P2=[eta1 xi];
%   P3=h*xi;
%
% % num='(x.^(P1(2)+1)).* besselj(P1(2),x*P3) * P2(2)^P1(2) ./';
% % den='( (x*P2(2)).^4 + P2(1)*(P2(2)*x).^2 + 1)';
% 
% %   Construct the integrand of the spectral representation
% fun=strcat(num,den);
% density=inline(fun,3);


density = @(x,P1,P2,P3) (x.^(P1(2)+1)).* besselj(P1(2),x*P3) * P2(2)^(P1(2)+2) ...
    ./ ( (x*P2(2)).^4 + P2(1)*(P2(2)*x).^2 + 1);

%   Define parameters used in the spectral density

%   First, dimension parameters
pp1=[dd dh];

%   Second, SSRF parameters
pp2=[eta1 xi];

%==========================================================================
%   The code segment below separates the integration domain in two
%   segments. Presumably most of the density is in the first segment. This
%   trick does not seem necessary for calculating the integral with QUADGK.
%==========================================================================
%   
% %   Define an arbitrary break point for the integration. This aims to
% %   separate the dominant part of the spectral density (up to the minimum
% %   of kmax or kup: kup *xi =5)
% 
% k1=min(5/xi,kmax);
% 
% cor1=zeros(size(h));
% cor2=zeros(size(h));

% if kmax>k1
%     for i=1:n
%         pp3=h(i)*xi;
%         cor1(i)=quadl(@(x)density(x,pp1,pp2,pp3),0,k1);
%         cor2(i)=quadl(@(x)density(x,pp1,pp2,pp3),k1,kmax);
%     end;
% elseif k1>=kmax
%     for i=1:n
%         pp3=h(i)*xi;
%         cor1(i)=quadl(@(x)density(x,pp1,pp2,pp3),0,kmax);
%         cor2(i)=0;
%     end;
%     
% end;
% 
% covar=eta0*g0*(cor1+cor2)./(h.^dh);

cor=zeros(size(h)); errbnd=zeros(size(h)); 
    for i=1:n
        pp3=h(i)*xi; 
        [cor(i),errbnd(i)]=quadgk(@(x)density(x,pp1,pp2,pp3),0,kmax);
    end;


covar=eta0*g0*cor./(h.^dh);
if resh==1
    covar=reshape(covar,size(lag));
end
err=eta0*g0*errbnd./(h.^dh);

if nargin >3 && flag>0 && resh==0
    figure(2);
    plot(h,covar,'b*',h,covar-err,'r-',h,covar+err,'r-');
%     ciplot(covar-err,covar+err,covar,'b-');
%     hold on; plot(h,covar,'b-');
    if length(h)<30
        hold on; plot(h,covar,'r.');
    end;
    xlabel('h'); ylabel('SRRF Covariance');
    title(['d=',num2str(dd)]);
    axis tight;
    figbolden
elseif nargin >3 && flag>0 && resh==1
    figure();
    surf(covar); shading flat;
    xlabel('h1'); xlabel('h2');  zlabel('SRRF Covariance');
    title(['d=',num2str(dd)]);
    axis tight;
    figbolden
end;
