function [a,b,c,d] = mvtf2ss(num,comden,iu)
%MVTF2SS        Multivariable Transfer function to state-space conversion.
%       [A,B,C,D] = MVTF2SS(NUM,Comden)  calculates the state-space representation:
%               .
%               x = Ax + Bu
%               y = Cx + Du
%
%       of the system:
%                         NUM(s)
%               H(s) = ----------
%                       Comden(s)
%
%       The vector Comden contains the coefficients of the
%       denominator in descending powers of s.  Matrix NUM contains the
%       numerator coefficients with as many rows as there are outputs y.
%       The number of columns of NUM equals
%       length(Comden) * number of system inputs.
%
%       [A,B,C,D] = MVTF2SS(NUM,Comden,iu) calculates the state-space
%       representation for input iu and all the outputs.

%       Dr M.P. Ford 4th September 1987
% Copyright (c) 1987 by GEC Engineering Research Centre & Cambridge Control Ltd
%
% History:
%       Made compatible with `parallel' in Control T/B ver.3, 9.7.93,JMM.
%       MRN0012
%       MRN0019

nargs = nargin;
error(nargchk(2,3,nargs));
[mnum,nnum] = size(num);   % mnum is always the number of outputs
                           % NB nnum is NOT the number of inputs
ld = length(comden);
if nnum > ld
   if rem(nnum,ld)~=0
      error('Numerator matrix not consistent with Common Denominator.');
   end
   ni=nnum/ld;         % number of inputs of system
else
   ni=1;
   num=[zeros(mnum,ld-nnum),num]; % pad out num to length comden
end
k=1:ld;

if nargs == 3    % just one input iu
  if max(size(iu))~=1
     error('Input specified is not a scalar');
  end
  if (iu<=0)|(iu>ni)
     error('Input specified does not exist in this system');
  end
  num=num(:,(iu-1)*ld+k);
  ni=1;
end

% Check for leading zeros in Common denominator
iz = find(comden ~= 0);
if iz(1)~=1
   error('Order of Numerator greater than Common denominator.')
end
% The above check assumes the Numerator has not also been padded out
% with zeros

%  normalize the numerator matrix
num = num./comden(1);

% Handle constant case:
[mnum,nnum]=size(num);
if ni == nnum
   a = 0;
   b = zeros(1,nnum);
   c = zeros(mnum,1);
   d = num;
return
end

% Normalize the rest of comden to comden(1),
comden = comden(2:ld) ./ comden(1);

outputs = 1:mnum;    % This line added 9.7.93, JMM.
for j=1:ni        % for the rest of the columns
    numj=num(:,k+(j-1)*ld);
    dj = numj(:,1);
    aj = [-comden; eye(ld-2,ld-1)];
    bj = eye(ld-1,1);
    cj = numj(:,2:ld) - numj(:,1) * comden;
    if j==1
       a=aj;
       b=bj;
       c=cj;
       d=dj;
    else
       % Correction made here 9.7.93, JMM:
       [a,b,c,d]=parallel(a,b,c,d,aj,bj,cj,dj,[],[],outputs,outputs);
                 % `parallel' connection, but of outputs only.
    end
end

if ni ~= 1   % find observable part
   [a,b,c,t,k]=obsvf(a,b,c);
   sk=sum(k);
   %if sk~=length(comden)
   %   disp(['MVTF2SS  Warning order of system greater',...
   %   ' than order of common denominator.'])
   %end
   [m,n]=size(a);
   i=m-sk+1:m;
   a=a(i,i);
   b=b(i,:);
   c=c(:,i);
end