function [Ao,Bo,Co,Do]=convsys(A1,B1,C1,D1,A2,B2,C2,D2)
%
% [Ao,Bo,Co,Do]=convsys(A1,B1,C1,D1,A2,B2,C2,D2)
%
% This algorithm gives the convolution of two state space representations
%
%        | A1  B1 |     | A2  B2 |
% u ==>  |        | ==> |        | ==> y
%        | C1  D1 |     | C2  D2 |
%
% The algorithm also accepts state space objects as inputs and gives out a
% state space object as output.

if ((nargin==2) && nargout>1) || ((nargin==8) && ~ismember(nargout,[0 4])) || ~ismember(nargin,[2 8])
    error('Incompatible Number of inputs or outputs');
end

if (nargin==2)
    if any([~isobject(A1) ~isobject(B1)])
        error('State Space Objects expected as input');
    else
        T1=A1;
        T2=B1;

        A1=T1.a;
        B1=T1.b;
        C1=T1.c;
        D1=T1.d;

        A2=T2.a;
        B2=T2.b;
        C2=T2.c;
        D2=T2.d;
    end
else
    checkcompatibility(A1,B1,C1,D1,'1');
    checkcompatibility(A2,B2,C2,D2,'2');
end

if (size(C1,1) ~= size(B2,2))
    error('Incompatible system dimensions')
end
Ao = [A1 zeros(length(A1),length(A2));B2*C1 A2];
Bo = [B1;B2*D1];
Co = [D2*C1 C2];
Do = D2*D1;
temps=cleanss(Ao,Bo,Co,Do);
if nargout==1
    Ao=temps;
end