function [An,bn,cn,Ap,bp,cp,d] = AdditiveDecomposition(A,B,C,D)
%Additive decomposition of Rational function
%Inputs required:
%A = an nxn matrix
%c = the "c" row vector, also of lenght n
%b = the "b" column vector of lenght n
%d = scaler
%Returns two triples, possibly of different McMillan degree.
%(An,bn,cn) with eigenvalues of An located in open right plane
%(Ap,bp,cp) with eigenvalues of Ap located in open left plane
% d remains unchanged
% Written originally by Wolfgang Scherrer with minor changes by Conor
% Sexton. Any questions should be directed towards Conor Sexton or Bernard
% Hanzon at hcsexton@hotmail.com and b.hanzon@ucc.ie
N = length(C);
[V,T,m] = cschur (A,1);
Eigenvalues = sort(eig(T));
P = find(Eigenvalues>=0);
Q = length(P);
S = N - Q;
% Because CSCHUR does not zero the lower triangle
% of T (unlike SCHUR), A11 and A22 end up being artificially badly scaled (very
% small entries below the diagonal) and get rescaled in LYAP with bad
% consequences for accuracy of the original equation. A simple remedy is to
% make T upper triangular as it should be
T = triu (T);
A11h = T (1:S,1:S);
A12h = T (1:S,S+1:end);
A22h = T (S+1:end,S+1:end);
Bh = V'*B;
Ch = C*V;
M = lyap (-A11h,A22h,A12h);
K = eye(N);
K(1:S,S+1:end) = M;
Bh = K*Bh;
Ch = Ch*inv(K);
bp = Bh(1:S,1);
bn = Bh(S+1:end,1);
cp = Ch(1,1:S);
cn = Ch(1,S+1:end);
Ap = A11h;
An = A22h;
d = D;
end