--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/base/__bodquist__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,156 @@
+## Copyright (C) 1996, 1998 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{f}, @var{w}] =} bodquist (@var{sys}, @var{w}, @var{out_idx}, @var{in_idx})
+## used internally by bode, nyquist; compute system frequency response.
+##
+## @strong{Inputs}
+## @table @var
+## @item sys
+## input system structure
+## @item w
+## range of frequencies; empty if user wants default
+## @item out_idx
+## list of outputs; empty if user wants all
+## @item in_idx
+## list of inputs; empty if user wants all
+## @item rname
+## name of routine that called bodquist ("bode" or "nyquist")
+## @end table
+## @strong{Outputs}
+## @table @var
+## @item w
+## list of frequencies
+## @item f
+## frequency response of sys; @math{f(ii) = f(omega(ii))}
+## @end table
+## @strong{Note} bodquist could easily be incorporated into a Nichols
+## plot function; this is in a "to do" list.
+##
+## Both bode and nyquist share the same introduction, so the common parts are
+## in bodquist.  It contains the part that finds the number of arguments,
+## determines whether or not the system is SISO, and computes the frequency
+## response.  Only the way the response is plotted is different between the
+## two functions.
+## @end deftypefn
+
+function [f, w] = bodquist (sys, w, outputs, inputs, rname)
+
+  ## check number of input arguments given
+  if (nargin != 5)
+    usage("[f,w] = bodquist(sys,w,outputs,inputs,rname)");
+  endif
+
+  ## check each argument to see if it's in the correct form
+  if (!is_struct(sys))
+    error("sys must be a system data structure");
+  endif
+
+  ## let freqresp determine w if it's not already given
+  USEW = freqchkw(w);
+
+  ## get initial dimensions (revised below if sysprune is called)
+  [nn,nz,mm,pp ] = sysdimensions(sys);
+
+  ## check for an output vector and to see whether it`s correct
+  if (!isempty(outputs))
+    if (isempty(inputs))
+      inputs = 1:mm;                    # use all inputs
+      warning([rname,": outputs specified but not inputs"]);
+    endif
+    sys = sysprune(sys,outputs,inputs);
+    [nn,nz,mm,pp ] = sysdimensions(sys);
+  endif
+
+  ## for speed in computation, convert local copy of
+  ## SISO state space systems to zero-pole  form
+  if( is_siso(sys) & strcmp( sysgettype(sys), "ss") )
+    [zer,pol,k,tsam,inname,outname] = sys2zp(sys);
+    sys = zp2sys(zer,pol,k,tsam,inname,outname);
+  endif
+
+  ## get system frequency response
+  [f,w] = freqresp(sys,USEW,w);
+
+  phase = arg(f)*180.0/pi;
+
+  if(!USEW)
+    ## smooth plots
+    pcnt = 5;           # max number of refinement steps
+    dphase = 5;         # desired max change in phase
+    dmag = 0.2;         # desired max change in magnitude
+    while(pcnt)
+      pd = abs(diff(phase));                    # phase variation
+      pdbig = vec(find(pd > dphase));
+
+      lp = length(f);  lp1 = lp-1;              # relative variation
+      fd = abs(diff(f));
+      fm = max(abs([f(1:lp1); f(2:lp)]));
+      fdbig = vec(find(fd > fm/10));
+
+      bigpts = union(fdbig, pdbig);
+
+      if(isempty(bigpts) )
+        pcnt = 0;
+      else
+        pcnt = pcnt - 1;
+        wnew = [];
+        crossover_points = find ( phase(1:lp1).*phase(2:lp) < 0);
+        pd(crossover_points) = abs(359.99+dphase - pd(crossover_points));
+        np_pts = max(3,ceil(pd/dphase)+2);              # phase points
+        nm_pts = max(3,ceil(log(fd./fm)/log(dmag))+2);  # magnitude points
+        npts = min(5,max(np_pts, nm_pts));
+
+        w1 = log10(w(1:lp1));
+        w2 = log10(w(2:lp));
+        for ii=bigpts
+          if(npts(ii))
+            wtmp = logspace(w1(ii),w2(ii),npts(ii));
+            wseg(ii,1:(npts(ii)-2)) = wtmp(2:(npts(ii)-1));
+          endif
+        endfor
+        wnew = vec(wseg)'; # make a row vector
+        wnew = wnew(find(wnew != 0));
+        wnew = sort(wnew);
+        wnew = create_set(wnew);
+        if(isempty(wnew))   # all small crossovers
+          pcnt = 0;
+        else
+          [fnew,wnew] = freqresp(sys,1,wnew);    # get new freq resp points
+          w = [w,wnew];                 # combine with old freq resp
+          f = [f,fnew];
+          [w,idx] = sort(w);            # sort into order
+          f = f(idx);
+          phase = arg(f)*180.0/pi;
+        endif
+      endif
+    endwhile
+  endif
+
+  ## ensure unique frequency values
+  [w,idx] = sort(w);
+  f = f(idx);
+
+  w_diff = diff(w);
+  w_dup = find(w_diff == 0);
+  w_idx = complement(w_dup,1:length(w));
+  w = w(w_idx);
+  f = f(w_idx);
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/base/__freqresp__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,133 @@
+## Copyright (C) 1996, 1998 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{out} =} freqresp (@var{sys}, @var{USEW}@{,@var{w}@});
+## Frequency response function - used internally by @code{bode}, @code{nyquist}.
+## minimal argument checking; "do not attempt to do this at home"
+##
+## @strong{Inputs}
+## @table @var
+## @item sys
+## system data structure
+## @item USEW
+## returned by @code{freqchkw}
+## @item optional
+## must be present if @var{USEW} is true (nonzero)
+## @end table
+## @strong{Outputs}
+## @table @var
+## @item @var{out}
+## vector of finite @math{G(j*w)} entries (or @math{||G(j*w)||} for MIMO)
+## @item w
+## vector of corresponding frequencies
+## @end table
+## @end deftypefn
+
+## Author: R. Bruce Tenison <btenison@eng.auburn.edu>
+## Created: July 11, 1994
+
+function [ff, w] = freqresp (sys, USEW, w);
+
+  ## SYS_INTERNAL accesses members of system data structure
+
+  save_val = empty_list_elements_ok;
+  empty_list_elements_ok = 1;
+
+  ## Check Args
+  if( (nargin < 2) || (nargin > 4) )
+    usage ("[ff,w] = freqresp(sys,USEW{,w})");
+  elseif( USEW & (nargin < 3) )
+    error("USEW=1 but w was not passed.");
+  elseif( USEW & isempty(w))
+    warning("USEW=1 but w is empty; setting USEW=0");
+    USEW=0;
+  endif
+
+  DIGITAL = is_digital(sys);
+
+  ## compute default w if needed
+  if(!USEW)
+    if(is_siso(sys))
+      sys = sysupdate(sys,"zp");
+      [zer,pol] = sys2zp(sys);
+    else
+      zer = tzero(sys);
+      pol = eig(sys2ss(sys));
+    endif
+
+    ## get default frequency range
+    [wmin,wmax] = bode_bounds(zer,pol,DIGITAL,sysgettsam(sys));
+    w = logspace(wmin,wmax,50);
+  else
+    w = reshape(w,1,length(w));         # make sure it's a row vector
+  endif
+
+  ## now get complex values of s or z
+  if(DIGITAL)
+    jw = exp(i*w*sysgettsam(sys));
+  else
+    jw = i*w;
+  endif
+
+  [nn,nz,mm,pp] = sysdimensions(sys);
+
+  ## now compute the frequency response - divide by zero yields a warning
+  if (strcmp(sysgettype(sys),"zp"))
+    ## zero-pole form (preferred)
+    [zer,pol,sysk] = sys2zp(sys);
+    ff = ones(size(jw));
+    l1 = min(length(zer)*(1-isempty(zer)),length(pol)*(1-isempty(pol)));
+    for ii=1:l1
+      ff = ff .* (jw - zer(ii)) ./ (jw - pol(ii));
+    endfor
+
+    ## require proper  transfer function, so now just get poles.
+    for ii=(l1+1):length(pol)
+      ff = ff ./ (jw - pol(ii));
+    endfor
+    ff = ff*sysk;
+
+  elseif (strcmp(sysgettype(sys),"tf"))
+    ## transfer function form
+    [num,den] = sys2tf(sys);
+    ff = polyval(num,jw)./polyval(den,jw);
+  elseif (mm==pp)
+    ## The system is square; do state-space form bode plot
+    [sysa,sysb,sysc,sysd,tsam,sysn,sysnz] = sys2ss(sys);
+    n = sysn + sysnz;
+    for ii=1:length(jw);
+      ff(ii) = det(sysc*((jw(ii).*eye(n)-sysa)\sysb)+sysd);
+    endfor;
+  else
+    ## Must be state space... bode
+    [sysa,sysb,sysc,sysd,tsam,sysn,sysnz] = sys2ss(sys);
+    n = sysn + sysnz;
+    for ii=1:length(jw);
+      ff(ii) = norm(sysc*((jw(ii)*eye(n)-sysa)\sysb)+sysd);
+    endfor
+
+  endif
+
+  w = reshape(w,1,length(w));
+  ff = reshape(ff,1,length(ff));
+
+  ## restore global variable
+  empty_list_elements_ok = save_val;
+endfunction
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/base/__stepimp__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,279 @@
+## Copyright (C) 1996, 1998 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{y}, @var{t}] = } stepimp (@var{sitype}, @var{sys} [, @var{inp}, @var{tstop}, @var{n}])
+## Impulse or step response for a linear system.
+## The system can be discrete or multivariable (or both).
+## This m-file contains the "common code" of step and impulse.
+##
+## Produces a plot or the response data for system sys.
+##
+## Limited argument checking; "do not attempt to do this at home".
+## Used internally in @code{impulse}, @code{step}. Use @code{step}
+## or @code{impulse} instead.
+## @end deftypefn
+## @seealso{step and impulse}
+
+## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de>
+## Created: October 2, 1997
+## based on lsim.m of Scottedward Hodel
+
+function [y, t] = stepimp (sitype, sys, inp, tstop, n)
+
+  if (sitype == 1)         IMPULSE = 0;
+  elseif (sitype == 2)     IMPULSE = 1;
+  else                     error("stepimp: invalid sitype argument.")
+  endif
+  sys = sysupdate(sys,"ss");
+
+  USE_DEF = 0;   # default tstop and n if we have to give up
+  N_MIN = 50;    # minimum number of points
+  N_MAX = 2000;  # maximum number of points
+  T_DEF = 10.0;  # default simulation time
+
+  ## collect useful information about the system
+  [ncstates,ndstates,NIN,NOUT] = sysdimensions(sys);
+  TSAMPLE = sysgettsam(sys);
+
+  if (nargin < 3)                      inp = 1;
+  elseif (inp < 1 | inp > NIN)         error("Argument inp out of range")
+  endif
+
+  DIGITAL = is_digital(sys);
+  if (DIGITAL)
+    NSTATES = ndstates;
+    if (TSAMPLE < eps)
+      error("stepimp: sampling time of discrete system too small.")
+    endif
+  else        NSTATES = ncstates;       endif
+  if (NSTATES < 1)
+    error("step: pure gain block (n_states < 1), step response is trivial");
+  endif
+  if (nargin < 5)
+    ## we have to compute the time when the system reaches steady state
+    ## and the step size
+    ev = eig(sys2ss(sys));
+    if (DIGITAL)
+      ## perform bilinear transformation on poles in z
+      for i = 1:NSTATES
+        pole = ev(i);
+        if (abs(pole + 1) < 1.0e-10)
+          ev(i) = 0;
+        else
+          ev(i) = 2 / TSAMPLE * (pole - 1) / (pole + 1);
+        endif
+      endfor
+    endif
+    ## remove poles near zero from eigenvalue array ev
+    nk = NSTATES;
+    for i = 1:NSTATES
+      if (abs(ev(i)) < 1.0e-10)
+        ev(i) = 0;
+        nk = nk - 1;
+      endif
+    endfor
+    if (nk == 0)
+      USE_DEF = 1;
+      ## printf("##STEPIMP-DEBUG: using defaults.\n");
+    else
+      ev = ev(find(ev));
+      x = max(abs(ev));
+      t_step = 0.2 * pi / x;
+      x = min(abs(real(ev)));
+      t_sim = 5.0 / x;
+      ## round up
+      yy = 10^(ceil(log10(t_sim)) - 1);
+      t_sim = yy * ceil(t_sim / yy);
+      ## printf("##STEPIMP-DEBUG: nk=%d   t_step=%f  t_sim=%f\n",
+      ##   nk, t_step, t_sim);
+    endif
+  endif
+
+  if (DIGITAL)
+    ## ---- sampled system
+    if (nargin == 5)
+      n = round(n);
+      if (n < 2)
+        error("stepimp: n must not be less than 2.")
+      endif
+    else
+      if (nargin == 4)
+        ## n is unknown
+      elseif (nargin >= 1)
+        ## tstop and n are unknown
+        if (USE_DEF)
+          tstop = (N_MIN - 1) * TSAMPLE;
+        else
+          tstop = t_sim;
+        endif
+      endif
+      n = floor(tstop / TSAMPLE) + 1;
+      if (n < 2)  n = 2;  endif
+      if (n > N_MAX)
+        n = N_MAX;
+        printf("Hint: number of samples limited to %d by default.\n", \
+               N_MAX);
+        printf("  ==> increase \"n\" parameter for longer simulations.\n");
+      endif
+    endif
+    tstop = (n - 1) * TSAMPLE;
+    t_step = TSAMPLE;
+  else
+    ## ---- continuous system
+    if (nargin == 5)
+      n = round(n);
+      if (n < 2)
+        error("step: n must not be less than 2.")
+      endif
+      t_step = tstop / (n - 1);
+    else
+      if (nargin == 4)
+        ## only n in unknown
+        if (USE_DEF)
+          n = N_MIN;
+          t_step = tstop / (n - 1);
+        else
+          n = floor(tstop / t_step) + 1;
+        endif
+      else
+        ## tstop and n are unknown
+        if (USE_DEF)
+          tstop = T_DEF;
+          n = N_MIN;
+          t_step = tstop / (n - 1);
+        else
+          tstop = t_sim;
+          n = floor(tstop / t_step) + 1;
+        endif
+      endif
+      if (n < N_MIN)
+        n = N_MIN;
+        t_step = tstop / (n - 1);
+      endif
+      if (n > N_MAX)
+        tstop = (n - 1) * t_step;
+        t_step = tstop / (N_MAX - 1);
+        n = N_MAX;
+      endif
+    endif
+    tstop = (n - 1) * t_step;
+    [jnk,B] = sys2ss(sys);
+    B = B(:,inp);
+    sys = c2d(sys, t_step);
+  endif
+  ## printf("##STEPIMP-DEBUG: t_step=%f n=%d  tstop=%f\n", t_step, n, tstop);
+
+  F = sys.a;
+  G = sys.b(:,inp);
+  C = sys.c;
+  D = sys.d(:,inp);
+  y = zeros(NOUT, n);
+  t = linspace(0, tstop, n);
+
+  if (IMPULSE)
+    if (!DIGITAL && (D'*D > 0))
+      error("impulse: D matrix is nonzero, impulse response infinite.")
+    endif
+    if (DIGITAL)
+      y(:,1) = D / t_step;
+      x = G / t_step;
+    else
+      x = B;
+      y(:,1) = C * x;
+      x = F * x;
+    endif
+    for i = 2:n
+      y(:,i) = C * x;
+      x = F * x;
+    endfor
+  else
+    x = zeros(NSTATES, 1);
+    for i = 1:n
+      y(:,i) = C * x + D;
+      x = F * x + G;
+    endfor
+  endif
+
+  if(nargout == 0)
+    ## Plot the information
+    oneplot();
+    gset nogrid
+    gset nologscale
+    gset autoscale
+    gset nokey
+    clearplot();
+    if (gnuplot_has_multiplot)
+      if (IMPULSE)
+        gm = zeros(NOUT, 1);
+        tt = "impulse";
+      else
+        ssys = ss2sys(F, G, C, D, t_step);
+        gm = dcgain(ssys);
+        tt = "step";
+      endif
+      ncols = floor(sqrt(NOUT));
+      nrows = ceil(NOUT / ncols);
+      for i = 1:NOUT
+        subplot(nrows, ncols, i);
+        title(sprintf("%s: | %s -> %s", tt,sysgetsignals(sys,"in",inp,1), ...
+          sysgetsignals(sys,"out",i,1)));
+        if (DIGITAL)
+          [ts, ys] = stairs(t, y(i,:));
+          ts = ts(1:2*n-2)';  ys = ys(1:2*n-2)';
+          if (length(gm) > 0)
+            yy = [ys; gm(i)*ones(size(ts))];
+          else
+            yy = ys;
+          endif
+          grid("on");
+          xlabel("time [s]");
+          ylabel("y(t)");
+          plot(ts, yy);
+        else
+          if (length(gm) > 0)
+            yy = [y(i,:); gm(i)*ones(size(t))];
+          else
+            yy = y(i,:);
+          endif
+          grid("on");
+          xlabel("time [s]");
+          ylabel("y(t)");
+          plot(t, yy);
+        endif
+      endfor
+      ## leave gnuplot in multiplot mode is bad style
+      oneplot();
+    else
+      ## plot everything in one diagram
+      title([tt, " response | ", sysgetsignals(sys,"in",inp,1), ...
+        " -> all outputs"]);
+      if (DIGITAL)
+        stairs(t, y(i,:));
+      else
+        grid("on");
+        xlabel("time [s]");
+        ylabel("y(t)");
+        plot(t, y(i,:));
+      endif
+    endif
+    y=[];
+    t=[];
+  endif
+  ## printf("##STEPIMP-DEBUG: gratulations, successfull completion.\n");
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/system/__abcddims__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,38 @@
+## Copyright (C) 1997 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{y}, @var{my}, @var{ny}] =} abcddims (@var{x})
+##
+## Used internally in @code{abcddim}.  If @var{x} is a zero-size matrix,
+## both dimensions are set to 0 in @var{y}.
+## @var{my} and @var{ny} are the row and column dimensions of the result.
+## @end deftypefn
+
+## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
+## Created: February 1997
+
+function [y, my, ny] = abcddims (x)
+
+  y = x;
+  if(isempty(y))
+    y = [];
+  endif
+  [my,ny] = size(y);
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/system/__syschnamesl__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,116 @@
+## Copyright (C) 1996 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{old_names} =} syschnamesl (@var{olist}, @var{old_names}, @var{inames}, @var{listname})
+## used internally in syschnames
+## item olist: index list
+## old_names: original list names
+## inames: new names
+## listname: name of index list
+##
+## combines the two string lists old_names and inames
+## @end deftypefn
+
+function old_names = syschnamesl (olist, old_names, inames, listname)
+
+  probstr = [];
+  if( max(olist) > rows(old_names) )
+    probstr = ["index list value(s) exceed(s) number of signals (", ...
+      num2str(rows(old_names)),")"];
+
+  elseif( length(olist) > rows(inames) )
+    probstr = ["index list dimension exceeds number of replacement names (", ...
+      num2str(rows(inames)),")"];
+
+  elseif(isempty(olist))
+    probstr = [];    # do nothing, no changes
+
+  elseif(min(size(olist)) != 1 )
+    probstr = "index list must be either a vector or an empty matrix";
+
+  elseif(max(olist) > rows(old_names))
+    probstr = ["max(",listname,")=",num2str(max(olist))," > ", ...
+        num2str(rows(old_names)),", too big"];
+
+  elseif(min(olist) < 1)
+    probstr = ["min(",listname,")=",num2str(min(olist))," < 1, too small"];
+
+  else
+    if( length(olist)  == 1)
+        len_in = columns(inames);
+        len_out = columns(old_names);
+
+      if (len_in < len_out)
+        inames(1,(len_in+1):(len_out)) = zeros(1,(len_out - len_in));
+      endif
+
+      old_names(olist,1:length(inames)) = inames;
+    elseif(length(olist) > 1)
+      for ii=1:length(olist)
+        mystr = inames(ii,:);
+        len_my = columns(mystr);
+        len_out = columns(old_names);
+
+        if (len_my < len_out)
+          mystr(1,(len_my+1):(len_out)) = " "*ones(1,(len_out - len_my));
+          len_my = len_out;
+        endif
+
+        old_names(olist(ii),1:len_my) = mystr;
+      endfor
+    endif
+  endif
+  if(!isempty(probstr))
+    ## the following lines are NOT debugging code!
+    disp("Problem in syschnames: old names are")
+    outlist(old_names," ")
+    disp("new names are")
+    outlist(inames,"    ")
+    disp("list indices are")
+    disp(olist)
+    error(sprintf("syschnames: \"%s\" dim=(%d x %d)--\n\t%s\n", ...
+        listname, rows(olist), columns(olist),probstr));
+  endif
+
+  ## change zeros  to blanks
+  if( find(old_names == 0) )
+    ## disp("syschnamesl: old_names contains zeros ")
+    ## old_names
+    ## disp("/syschnamesl");
+
+    [ii,jj] = find(old_names == 0);
+    for idx=1:length(ii)
+      old_names(ii(idx),jj(idx)) = " ";
+    endfor
+
+    ## disp("syschnamesl: old_names fixed zeros ")
+    ## old_names
+    ## disp("/syschnamesl");
+  endif
+
+  ## just in case it's not a string anymore
+  if( !isstr(old_names) )
+    old_names = setstr(old_names);
+  endif
+
+  ## disp("syschnamesl: exit, old_names=")
+  ## old_names
+  ## disp("/syschnamesl: exiting")
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/system/__syscont_disc__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,54 @@
+## Copyright (C) 1996, 1998 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{n_tot}, @var{st_c}, @var{st_d}, @var{y_c}, @var{y_d}] =} syscont_disc(@var{sys})
+## Used internally in syscont and sysdisc.
+##
+## @strong{Inputs}
+## @var{ sys} is a system data structure.
+##
+## @strong{Outputs}
+## @table @var
+## @item n_tot
+## total number of states
+## @item st_c
+## vector of continuous state indices (empty if none)
+## @item st_d
+## vector of discrete state indices (empty if none)
+## @item y_c
+## vector of continuous output indices
+## @item y_d
+## vector of discrete output indices
+## @end table
+## @end deftypefn
+
+## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
+## Created: February 1997
+
+function [n_tot, st_c, st_d, y_c, y_d] = syscont_disc (sys)
+
+  ## get ranges for discrete/continuous states and outputs
+  [nn,nz,mm,pp,yd] = sysdimensions(sys);
+  n_tot = nn + nz;
+  st_c = 1:(nn);
+  st_d = nn + (1:nz);
+  y_c = find(yd == 0);          # y_c, y_d will be empty if there are none.
+  y_d = find(yd == 1);
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/system/__sysdefioname__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,57 @@
+## Copyright (C) 1996, 1998 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{ioname} =} sysdefioname (@var{n},@var{str} @{,@var{m}@})
+## return default input or output names given @var{n}, @var{str}, @var{m}.
+## @var{n} is the final value, @var{str} is the string prefix, and @var{m}
+## is start value
+##
+## used internally, minimal argument checking
+##
+## @strong{Example} @code{ioname = sysdefioname(5,"u",3)}
+## returns the list:
+## @example
+## ioname =
+## (
+##   [1] = u_3
+##   [2] = u_4
+##   [3] = u_5
+## )
+## @end example
+## @end deftypefn
+
+function ioname = sysdefioname (n, str, m)
+
+  if (nargin < 2 | nargin > 3)
+    usage("ioname = sysdefioname(n,str[,m])");
+  endif
+
+  if (nargin == 2)           m = min(1,n);            endif
+
+  ioname = list();
+  jj = 1;
+  if(n > 0 & m > 0 & m <= n)
+    for ii = m:n
+      ioname(ii+1-m) = sprintf("%s_%d",str,ii);
+    endfor
+  elseif(m > n)
+    error("str=%s; start value m=%d > final value n=%d",str,m,n);
+  endif
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/system/__sysdefstname__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,42 @@
+## Copyright (C) 1996 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{stname} =} sysdefstname (@var{n}, @var{nz})
+## return default state names given @var{n}, @var{nz}
+##
+## used internally, minimal argument checking
+## @end deftypefn
+
+function stname = sysdefstname (n, nz)
+
+  stname = list();
+  if(n > 0)
+    for ii = 1:n
+      stname(ii) = sprintf("x_%d",ii);
+    endfor
+  endif
+
+  ## Set default names for discrete states
+  if(nz > 0)
+    for ii = (n+1):(n+nz)
+      stname(ii) = sprintf("xd_%d",ii);
+    endfor
+  endif
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/system/__sysgroupn__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,57 @@
+## Copyright (C) 1996, 1998 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{names} =} sysgroupn (@var{names})
+## names = sysgroupn(names)
+## Locate and mark duplicate names
+## inputs:
+## names: list of signal names
+## kind: kind of signal name (used for diagnostic message purposes only)
+## outputs:
+## returns names with unique suffixes added; diagnostic warning
+## message is printed to inform the user of the new signal name
+##
+## used internally in sysgroup and elsewhere.
+## @end deftypefn
+
+function names = sysgroupn (names, kind)
+
+  ## check for duplicate names
+  l = length(names);
+  ii = 1;
+  while(ii <= l-1)
+    st1 = nth(names,ii);
+    jj = ii+1;
+    while ( jj <= l)
+      st2 = nth(names,jj);
+      if(strcmp(st1,st2))
+        suffix = ["_",num2str(jj)];
+        warning("sysgroup: %s name(%d) = %s name(%d) = %s", ...
+          kind,ii,kind,jj,st1);
+        strval = sprintf("%s%s",st2,suffix);
+        names(jj) = strval;
+        warning("sysgroup:     changed %s name %d to %s",kind,jj,strval);
+        ## restart the check (just to be sure there's no further duplications)
+        ii = 0; jj = l;
+      endif
+      jj = jj+1;
+    endwhile
+    ii = ii+1;
+  endwhile
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/system/__tf2sysl__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,35 @@
+## Copyright (C) 1996 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{vec} =} tf2sysl (@var{vec})
+## used internally in @ref{tf2sys}.
+## strip leading zero coefficients to get the true polynomial length
+## @end deftypefn
+
+function vec = tf2sysl (vec)
+
+  while( (length(vec) > 1) & (vec(1) == 0) )
+    vec = vec(2:length(vec));
+  endwhile
+
+  if(vec(1) == 0)
+    warning("tf2sys: polynomial has no nonzero coefficients!")
+  endif
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/system/__zp2ssg2__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,68 @@
+## Copyright (C) 1996, 1998 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{poly}, @var{rvals}] =} zp2ssg2 (@var{rvals})
+## Used internally in @code{zp2ss}
+## Extract 2 values from @var{rvals} (if possible) and construct
+## a polynomial with those roots.
+## @end deftypefn
+
+## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
+## Created: August 1996
+
+function [poly, rvals] = zp2ssg2 (rvals)
+
+  ## locate imaginary roots (if any)
+  cidx = find(imag(rvals));
+
+  if(!isempty(cidx))
+    ## select first complex root, omit from cidx
+    r1i = cidx(1);      r1 = rvals(r1i);     cidx = complement(r1i,cidx);
+
+    ## locate conjugate root (must be in cidx list, just in case there's
+    ## roundoff)
+    err = abs(rvals(cidx) - r1');
+    minerr = min(err);
+    c2i = find(err == minerr);
+    r2i = cidx(c2i);
+    r2 = rvals(r2i);
+    cidx = complement(r2i,cidx);
+
+    ## don't check for divide by zero, since 0 is not complex.
+    if(abs(r2 - r1')/abs(r1) > 1e-12)
+      error(sprintf("r1=(%f,%f); r2=(%f,%f), not conjugates.", ...
+        real(r1),imag(r1),real(r2),imag(r2)));
+    endif
+
+    ## complex conjugate pair
+    poly = [1, -2*real(r1), real(r1)^2+imag(r1)^2];
+  else
+    ## select two roots (they're all real)
+    r1 = rvals(1);
+    r2 = rvals(2);
+    poly = [1, -(r1+r2), (r1*r2)];
+    r1i = 1;  r2i = 2;
+  endif
+
+  ## remove roots used
+  idx = complement([r1i, r2i],1:length(rvals));
+  rvals = rvals(idx);
+
+endfunction
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/util/__outlist__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,80 @@
+## Copyright (C) 1996, 1998 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} outlist (@var{lmat}@{, @var{tabchar}, @var{yd}, @var{ilist} @})
+## Prints an enumerated list of strings.
+## internal use only; minimal argument checking performed
+##
+## @strong{Inputs}
+## @table @var
+## @item        lmat
+## list of strings
+## @item        tabchar
+## tab character (default: none)
+## @item   yd
+## indices of strings to append with the string "(discrete)"
+## (used by @var{sysout}; minimal checking of this argument)
+## @math{yd = [] } indicates all outputs are continuous
+## @item ilist
+## index numbers to print with names.
+##
+## default: @code{1:rows(lmat)}
+## @end table
+##
+## @strong{Outputs}
+## prints the list to the screen, numbering each string in order.
+## @end deftypefn
+
+## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
+## Created: December 1995
+
+function str_val = outlist (name_list, tabchar, yd, ilist)
+
+  ## save for restore later
+  save_empty = empty_list_elements_ok;
+  empty_list_elements_ok = 1;
+
+  if( nargin < 1 | nargin > 4 )
+    usage("str_val = outlist(x[,tabchar,yd,ilist])");
+  endif
+
+  m = length(name_list);
+  if(nargin < 4)           ilist = 1:m;          endif
+  if(nargin ==1)
+    empty_list_elements_ok = 1;
+    tabchar = "";
+  endif
+
+  if(nargin < 3)             yd = zeros(1,m);
+  elseif(isempty(yd))        yd = zeros(1,m);          endif
+
+  str_val = "";
+  dstr = list(""," (discrete)");
+  if((m >= 1) && (is_list(name_list)))
+    for ii=1:m
+      str_val = sprintf("%s%s%d: %s%s\n",str_val,tabchar, ilist(ii), ...
+          nth(name_list,ii),nth(dstr,yd(ii)+1));
+    endfor
+  else
+    str_val = sprintf("%sNone",tabchar);
+  endif
+
+  empty_list_elements_ok = save_empty;
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/util/__zgpbal__.m	Fri Jan 14 06:00:04 2000 +0000
@@ -0,0 +1,109 @@
+## Copyright (C) 1996 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{retsys} =} zgpbal (@var{Asys})
+##
+## used internally in @code{tzero}; minimal argument checking performed
+##
+## implementation of zero computation generalized eigenvalue problem
+## balancing method (Hodel and Tiller, Allerton Conference, 1991)
+## Based on Ward's balancing algorithm (SIAM J. Sci Stat. Comput., 1981)
+##
+## zgpbal computes a state/input/output weighting that attempts to
+## reduced the range of the magnitudes of the nonzero elements of [a,b,c,d]
+## The weighting uses scalar multiplication by powers of 2, so no roundoff
+## will occur.
+##
+## zgpbal should be followed by zgpred
+## @end deftypefn
+
+## References:
+## ZGEP: Hodel, "Computation of Zeros with Balancing," 1992, submitted to  LAA
+## Generalized CG: Golub and Van Loan, "Matrix Computations, 2nd ed" 1989
+
+## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
+## Created: July 24, 1992
+## Conversion to Octave by R. Bruce Tenison July 3, 1994
+
+function retsys = zgpbal (Asys)
+
+  if( (nargin != 1) | (!is_struct(Asys)))
+    usage("retsys = zgpbal(Asys)");
+  endif
+
+  Asys = sysupdate(Asys,"ss");
+  [a,b,c,d] = sys2ss(Asys);
+
+  [nn,mm,pp] = abcddim(a,b,c,d);
+
+  np1 = nn+1;
+  nmp = nn+mm+pp;
+
+  ## set up log vector zz, incidence matrix ff
+  zz = zginit(a,b,c,d);
+
+  ## disp("zgpbal: zginit returns")
+  ## zz
+  ## disp("/zgpbal")
+
+  if (norm(zz))
+    ## generalized conjugate gradient approach
+    xx = zgscal(a,b,c,d,zz,nn,mm,pp);
+
+    for i=1:nmp
+      xx(i) = floor(xx(i)+0.5);
+      xx(i) = 2.0^xx(i);
+    endfor
+
+    ## now scale a
+    ## block 1: a = sigma a inv(sigma)
+    for i=1:nn
+      a(i,1:nn) = a(i,1:nn)*xx(i);
+      a(1:nn,i) = a(1:nn,i)/xx(i);
+    endfor
+    ## block 2: b= sigma a phi
+    for j=1:mm
+      j1 = j+nn;
+      b(1:nn,j) = b(1:nn,j)*xx(j1);
+    endfor
+    for i=1:nn
+      b(i,1:mm) = b(i,1:mm)*xx(i);
+    endfor
+    for i=1:pp
+      i1 = i+nn+mm;
+      ## block 3: c = psi C inv(sigma)
+      c(i,1:nn) = c(i,1:nn)*xx(i1);
+    endfor
+    for j=1:nn
+      c(1:pp,j) = c(1:pp,j)/xx(j);
+    endfor
+    ## block 4: d = psi D phi
+    for j=1:mm
+      j1 = j+nn;
+      d(1:pp,j) = d(1:pp,j)*xx(j1);
+    endfor
+    for i=1:pp
+      i1 = i + nn + mm;
+      d(i,1:mm) = d(i,1:mm)*xx(i1);
+    endfor
+  endif
+
+  retsys = ss2sys(a,b,c,d);
+endfunction
+