function [C,diffuse] = acovf(T,R,K,Z,H,D,U,Omega,eigval,order)
% ACOVF  [Not a public function] Autocovariance function for general state space.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

%**************************************************************************

realsmall = getrealsmall();

if isempty(eigval)
    eigval = ordeig(T);
end

ny = size(Z,1);
[nx,nb] = size(T);
nf = nx - nb;
ne = size(R,2);

if ny == 0
    Z = zeros(0,nb);
    H = zeros(0,ne);
end

Tf = T(1:nf,:);
Ta = T(nf+1:end,:);
Rf = R(1:nf,:);
Ra = R(nf+1:end,:);

% Detect unit root elements of y, xf, alpha, and xb.
unitroots = abs(abs(eigval(1:nb)) - 1) <= realsmall;
dy = any(abs(Z(:,unitroots)) > realsmall,2).';
df = any(abs(Tf(:,unitroots)) > realsmall,2).';
da = unitroots;
if ~isempty(U)
    db = any(abs(U(:,unitroots)) > realsmall,2).';
else
    db = da;
end
Caa = zeros(nb);

% Solve Lyapunov equation for the contemporaneous covariance matrix of the
% stable elements of the vector alpha.
Caa(~da,~da) = covfun.lyapunov(Ta(~da,~da),Ra(~da,:)*Omega*Ra(~da,:).');

Ra_Omega_Rft = Ra*Omega*Rf.';
Cff = Tf*Caa*Tf.' + Rf*Omega*Rf.';
Cyy = Z*Caa*Z.' + H*Omega*H.';
Cyf = Z*Ta*Caa*Tf.' + Z*Ra_Omega_Rft;
Cya = Z*Caa;
Cfa = Tf*Caa*Ta.' + Ra_Omega_Rft.';

C = zeros(ny+nf+nb,ny+nf+nb,1+order);
C(:,:,1) = [ ...
    Cyy,Cyf,Cya; ...
    Cyf',Cff,Cfa; ...
    Cya',Cfa',Caa; ...
    ];
C(:,:,1) = (C(:,:,1) + C(:,:,1).')/2;
diffuse = [dy,df,db];

if order > 0
    TT = [Z*Ta;Tf;Ta];
    for i = 1 : order
        C(1:end,:,i+1) = TT*C(ny+nf+1:end,:,i);
    end
end

if ~isempty(U)
    for i = 0 : order
        C(ny+nf+1:end,:,i+1) = U*C(ny+nf+1:end,:,i+1);
        C(:,ny+nf+1:end,i+1) = C(:,ny+nf+1:end,i+1)*U.';
    end
end

C(diffuse,:,:) = Inf;
C(:,diffuse,:) = Inf;

end