function [Q,R,H,P] = pareq(A,B,C,D,E,J,K,L,rho)

% Constructs the Chow-Willsky matrices for the system
%  dx = Ax + Bu + Ed + Kf
%   y = Cx + Du + Jd + Lf
%
% When disturbances and faults are present:
% [Q,R,H,P] = pareq(A,B,C,D,E,J,K,L [,rho])
%
% When there are no disturbances:
% [Q,R,P] = pareq(A,B,C,D,K,L [,rho])
%
% When there are no faults:
% [Q,R,H] = pareq(A,B,C,D,E,J [,rho])
%
% The default value of 'rho' is n.
%
% (c) Mattias Nyberg 1997, 1998

if ((nargin == 8) | (nargin == 9))
  
    n = size(A, 1);
	if nargin==8
      rho = n;
	end
    k = size(B, 2);
    m = size(C, 1);
    kd = size(E, 2);
    kf = size(K, 2);
	  
    if ((size(B,1)~=n) | (size(E,1)~=n)  | (size(K,1)~=n))
      error('Wrong matrix dimensions 1!');
    end
    if (size(D,1)~=m | size(J,1)~=m | size(L,1)~=m)
      error('Wrong matrix dimensions 2!');
    end
    if (size(D,2)~=k | size(J,2)~=kd | size(L,2)~=kf)
      error('Wrong matrix dimensions 3!');
    end

	Qr = D;
	for i=0:rho-1,
	  Qr = [C*A^i*B, Qr];
	end
	Q = Qr;
	for i=1:rho,
	  Q = [[Qr(:,i*k+1:(rho+1)*k), zeros(m,i*k)]; Q];
	end
	
	R = [];
	for i=0:rho,
	  R = [R; C*A^i];
	end

	Hr = J;
	for i=0:rho-1,
	  Hr = [C*A^i*E, Hr];
	end
	H = Hr;
	for i=1:rho,
	  H = [[Hr(:,i*kd+1:(rho+1)*kd), zeros(m,i*kd)]; H];
	end
	
	Pr = L;
	for i=0:rho-1,
	  Pr = [C*A^i*K, Pr];
	end
	P = Pr;
	for i=1:rho,
	  P = [[Pr(:,i*kf+1:(rho+1)*kf), zeros(m,i*kf)]; P];
	end

else
    n = size(A, 1);
	if nargin==6
      rho = n;
	else
	  rho = K;
	end
    K = E;
    L = J;    
    k = size(B, 2);
    m = size(C, 1);
    kf = size(K, 2);

    if ((size(B,1)~=n) | (size(K,1)~=n))
      error('Wrong matrix dimensions 1!');
    end
    if (size(D,1)~=m | size(L,1)~=m)
      error('Wrong matrix dimensions 2!');
    end
    if (size(D,2)~=k | size(L,2)~=kf)
      error('Wrong matrix dimensions 3!');
    end

	Qr = D;
	for i=0:rho-1,
	  Qr = [C*A^i*B, Qr];
	end
	Q = Qr;
	for i=1:rho,
	  Q = [[Qr(:,i*k+1:(rho+1)*k), zeros(m,i*k)]; Q];
	end
	
	R = [];
	for i=0:rho,
	  R = [R; C*A^i];
	end

	Pr = L;
	for i=0:rho-1,
	  Pr = [C*A^i*K, Pr];
	end
	P = Pr;
	for i=1:rho,
	  P = [[Pr(:,i*kf+1:(rho+1)*kf), zeros(m,i*kf)]; P];
	end
	    
    H = P;
end;