function sys = rssnb(n,p,m,varargin)
%RSSNB   Generate random stable continuous-time state-space models
%with H_inf norm = 1.
%
%   SYS = RSSNB(N) generates an Nth-order SISO state-space model SYS.
%
%   SYS = RSSNB(N,P) generates a single-input Nth-order model with 
%   P outputs.
%
%   SYS = RSSNB(N,P,M) generates an Nth-order model with P outputs
%   and M inputs.
%
%   SYS = RSSNB(N,P,M,S1,...,Sk) generates a S1-by-...-by-Sk array of
%   state-space models with N states, P outputs, and M inputs.
%
%   See also RSS, DRSS, DRSSNB

%   Copyright 2003-2004 The MathWorks, Inc. 

switch nargin
case 0
   n=max([1,round(abs(10*randn(1,1)))]);
   p=max([1,round(4*randn(1,1))]);
   m=max([1,round(4*randn(1,1))]);
case 1
   m=1;
   p=1;
case 2
   m=1;
end
arraydims= [varargin{:}];

% Check all inputs are positive integers
sizes = [m n p arraydims];
if ~isequal(sizes,round(sizes)) | ~all(isfinite(sizes)) | any(sizes<0) 
   error('Input arguments must be non negative integers.')
end

% Generate random state matrices
a = randn([n,n,arraydims]);
b = randn([n,m,arraydims]);
c = randn([p,n,arraydims]);
d = randn([p,m,arraydims]);

% Scale poles to lie in unit disk and then bilinear transform
I = eye(n);
sqrt2 = sqrt(2);
for k = 1:prod(arraydims)
  a_d = a(:,:,k);
  b_d = b(:,:,k);
  c_d = c(:,:,k);
  d_d = d(:,:,k);

  e = eig(a_d);
  emax = max([abs(e); 1e-8]);
  a_d = (1-rand^2)*a_d/emax;  
  
  iad = inv(I+a_d);
  
  a(:,:,k) = iad*(a_d-I);
  b(:,:,k) = sqrt2*iad*b_d;
  c(:,:,k) = sqrt2*c_d*iad;
  d(:,:,k) = d_d-c_d*iad*b_d;
end

% Normalize each system to have H_inf gain = 1
sys = ss(a,b,c,d);
sysg = norm(sys,inf);
for k = 1:prod(arraydims)
  sys(:,:,k) = sys(:,:,k)/sysg(k);
end