function h2d = diff2(h2,N,a1)
%DIFF2 Difference of triangular 2nd order kernel in vector form
%
% H2D = DIFF2(H2,N,A) returns the triangular kernel associated to
% the generalized transfer function H2D(z1,z2)=(1-A/z1)(1-A/z2)H2(z1,z2).
%
% N is the memory length (in samples) of H2 and H2D, which are vectorized

%   CBVK Toolbox
%   JHMG, PMSB 2016


h2d = zeros(size(h2));

a2 = a1*a1;

ind = 0;

ind = ind+1;
% (1,1)
h2d(ind) = h2(ind);

ind = ind+1;
% (1,2)
h2d(ind) =  - a1*h2( ind + -1 ) + (1/2)*h2(ind);
h2d(ind) = 2*h2d(ind);

for k2=1+2:N
    ind = ind+1;
    % (1,k2)
    h2d(ind) =  - a1*(1/2)*h2( ind + -1 ) + (1/2)*h2(ind);
    h2d(ind) = 2*h2d(ind);
    
end

for k1=2:N
    ind = ind+1;
    % (k1,k1)
    h2d(ind) = a2*h2( ind + k1 - N - 2 ) - a1*2*(1/2)*h2( ind + k1 - N - 1 ) + h2(ind);
    
    if (k1 <= N-1)
        ind = ind+1;
        % (k1,k1+1)
        h2d(ind) = a2*(1/2)*h2( ind + k1 - N - 2 ) - a1*h2( ind + -1 ) - a1*(1/2)*h2( ind + k1 - N - 1 ) + (1/2)*h2(ind);
        h2d(ind) = 2*h2d(ind);
        
    end
    
    for k2=k1+0+2:N
        ind = ind+1;
        % (k1,k2)
        h2d(ind) = a2*(1/2)*h2( ind + k1 - N - 2 ) - a1*(1/2)*h2( ind + -1 ) - a1*(1/2)*h2( ind + k1 - N - 1 ) + (1/2)*h2(ind);
        h2d(ind) = 2*h2d(ind);
        
    end
end

