function sys = drss(n,p,m)
%DRSS  Generate random stable discrete-time state-space models.
%
%    SYS = DRSS(N)  generates an Nth-order SISO state-space model SYS.
%
%    SYS = DRSS(N,P)  generates a single-input Nth-order model with 
%    P outputs.
%
%    SYS = DRSS(N,P,M)  generates an Nth-order model with M inputs and 
%    P outputs.
%
%    In all cases, the sample time is left unspecified (set to -1).
%    To generate random discrete TF or ZPK models, convert the result SYS 
%    to the appropriate model type with the functions TF or ZPK.
%
%    See also:  RSS, TF, ZPK.

%   Clay M. Thompson 1-22-91
%   Copyright (c) 1986-97 by The MathWorks, Inc.
%   $Revision: 1.2 $  $Date: 1997/04/09 00:36:02 $

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

% Prob of an integrator is 0.10 for the first and 0.01 for all others
nint = (rand(1,1)<0.10)+sum(rand(n-1,1)<0.01);

% Prob of repeated roots is 0.05
nrepeated = floor(sum(rand(n-nint,1)<0.05)/2);

% Prob of complex roots is 0.5
ncomplex = floor(sum(rand(n-nint-2*nrepeated,1)<0.5)/2);
nreal = n-nint-2*nrepeated-2*ncomplex;

% Determine random poles
rep = 2*rand(nrepeated,1)-1;
mag = rand(ncomplex,1);
ang = pi*rand(ncomplex,1);
jay = sqrt(-1);
complex = mag.*exp(jay*ang);
re = real(complex);
im = imag(complex);

% Generate random state space model
a = zeros(n,n);
for i=1:ncomplex,
   ndx = [2*i-1,2*i];
   a(ndx,ndx) = [re(i),im(i);-im(i),re(i)];
end
ndx = [2*ncomplex+1:n];
if ~isempty(ndx),
   a(ndx,ndx) = diag([ones(nint,1);rep;rep;2*rand(nreal,1)-1]);
end
T = orth(rand(n,n));
a = T\a*T;
b = rand(n,m);
c = rand(p,n);
d = rand(p,m);
b = b .* (rand(n,m)<0.75);
c = c .* (rand(p,n)<0.75);
d = d .* (rand(p,m)<0.5);

sys = ss(a,b,c,d,-1);

% end drss
