function kernels = voltkernels(f,g,h,x,names,values,P,N,T,vb,df,vd,tp,wtylr,D,Dg,K)
%VOLTKERNELS computes the triangular Volterra kernels of a linear-analytic
%system, using Carleman Bilinearizations
%
% kernels = VOLTKERNELS(f,g,H,X,NAMES,VALUES,P,N,T,VB) computes samples
% spaced by T seconds of the first P <=5 triangular Volterra kernels up to
% order P of the linear-analytic system x'= f(x) + g(x)u, y = hx. The
% kernels are returned in vector form in a cell.
% 
% NAMES is a structure with the names of the parameters of f, g and h and
% VALUES is a structure (or a cell of structures with different sets) of
% parameter values. 
% 
% VB ~= 0 sets verbose mode on (default off)
%
% kernels = VOLTKERNELS(f,g,H,X,NAMES,VALUES,P,N,T,VB,DF,VD,TP) with DF=1
% (default 0) differentiates results to obtain ZOH kernels when integrator
% has been added at the input of the system (presently for P<5); with VD=1
% (default 0), kernels are returned in diagonal form, to be used in
% simulations with kvdiags; with TP=1 (default 0) applies T^p factors to
% kernels
%
% kernels = VOLTKERNELS(f,g,H,X,NAMES,VALUES,P,N,T,VB,DF,VD,TP,WTYLR) uses
% pre-calculated Taylor expansion coefficients
%
% kernels = VOLTKERNELS(f,g,H,X,NAMES,VALUES,P,N,T,VB,DF,VD,TP,WTYLR,D,Dg,K)
% uses also pre-calculated matrices in the bilinearization
%
% Modify this function accordingly if only a certain operation on the
% kernels for each set of parameters is required, and not returning the
% kernels themselves

% CBVK Toolbox
% JHMG, PMSB 2016

% ----------------------------------------------------

if (df==1 & P>4)
    error('cbvk:voltkernelsP','Diagonalization not yet implemented for P=5');
end

% Number of parameter sets
if iscell(values)
    NS = length(values);
else
    NS = 1;
end

kernels=cell(NS,P);

% Carleman Bilinearization
% Loop of parameter sets

if nargin<10
    vb=0;
elseif nargin <11
    vd=0;df=0;tp=0
end
        

for ns=1:NS
    if nargin<14
        [F,G,b,c] = carlbl(f,g,h,x,names,values,P,vb);
    elseif nargin<15
        [F,G,b,c] = carlbl(f,g,h,x,names,values,P,vb,wtylr);
    else
        [F,G,b,c] = carlbl(f,g,h,x,names,values,P,vb,wtylr,D,Dg,K);
    end
end

% Dimension
m=length(x);

% Kernels
% Loop of parameter sets
for ns=1:NS
    if vb
        if NS==1
            disp('Kernel computation')
        else
            disp(['Kernel computation, parameter set ' num2str(ns)])
        end
    end
    
    if NS==1
        kernels = vkrnsvdx(F,G,b,c,N,T,m,P);
    else
        kernels(ns,:)= vkrnsvdx(F{ns},G{ns},b{ns},c,N,T,m,P);
    end
    
    % T^p factors
    if tp
        for p=1:P
            kernels{ns,p}=kernels{ns,p}*(T^p);
        end
    end
    
    if (df & tp)
        warning('cbvk:voltkernelsD','Differentiation not usually applied together with T^p factors');
    end
    % Differentiation
    if df
        kernels{ns,1}=kernels{ns,1}-[0 kernels{ns,1}(1:N-1)];
        if P>1
            kernels{ns,2}=diff2(kernels{ns,2},N,1);
            if P>2
                kernels{ns,3}=diff3(kernels{ns,3},N,1);
                if P>3
                    kernels{ns,4}=diff4(kernels{ns,4},N,1);
                    % Not yet implemented
%                     if P>4
%                         kernels{ns,5}=diff5(kernels{ns,5},N,1);
%                     end
                end
            end
        end
    end
    
    % Diagonalization
    if vd
        for p=2:P
            kernels{ns,p}=kvdiags(kernels{ns,p},p,N);
        end
    end
end


        
    




