function sys = drssnb(n,p,m,varargin)
%DRSSNB  Generate random stable discrete-time state-space models
%with H_inf norm = 1.
%
%   SYS = DRSSNB(N) generates an Nth-order SISO state-space model SYS.
%
%   SYS = DRSSNB(N,P) generates a single-input Nth-order model with 
%   P outputs.
%
%   SYS = DRSSNB(N,P,M) generates an Nth-order model with P outputs
%   and M inputs.
%
%   SYS = DRSSNB(N,P,M,S1,...,Sk) generates a S1-by-...-by-Sk array of
%   state-space models with N states, P outputs, and M inputs.
%
%   In all cases, the sample time is left unspecified (set to -1).
%
%   See also RSS, DRSS, RSSNB.

%   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
for k = 1:prod(arraydims)
  e = eig(a(:,:,k));
  emax = max([abs(e); 1e-8]);
  a(:,:,k) = (1-rand^2)*a(:,:,k)/emax;  
end

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