--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/are.m	Tue Aug 10 20:56:55 1993 +0000
@@ -0,0 +1,60 @@
+function x = are(a,b,c,opt)
+
+#
+# usage: x = are(a,b,c{,opt})
+#
+# solves algebraic riccati equation
+#
+#    a' x + x a - x b x + c = 0
+#
+# for identically dimensioned square matrices a,b,c.  If b (c) is not square,
+# then the function attempts to use b*b' (c'*c) instead.
+#
+# opt is an option passed to the eigenvalue balancing routine; default is 'B'
+#
+# see also: balance
+
+if((nargin == 3) || (nargin == 4))
+  if(nargin==4)
+    if((opt ~= 'N') || (opt ~= 'P') || (opt ~= 'S') || (opt ~= 'B'))
+      printf('warning: are: opt has an illegal value; setting to B');
+      opt = 'B';
+    endif
+  else
+    opt = 'B';
+  endif
+  if( (n = is_square(a) ) == 0)
+    error('are: a is not square')
+  endif
+  if(is_controllable(a,b) == 0)
+    printf('warning: are: a,b are not controllable')
+  endif
+  if( (m = is_square(b)) == 0)
+    b = b*b';
+    m = rows(b);
+  endif
+  if(is_observable(a,c) == 0)
+    printf('warning: are: a,c are not observable')
+  endif
+  if( (p = is_square(c)) == 0)
+    c = c'*c;
+    p = rows(c);
+  endif
+  if( (n ~= m) || (n ~= p))
+    error('are: a,b,c not conformably dimensioned.')
+  endif
+
+  # should check for controllability/observability here
+  # use Boley-Golub (Syst. Contr. Letters, 1984) method, not the
+  #                    n-1
+  # rank([ B A*B ... A^   *B]) method 
+  
+  [u,s] = schur(balance([a,-b;-c,-a'],opt),'A');
+  n1=n+1;
+  n2 = 2*n;
+  x = u(n1:n2,1:n)/u(1:n,1:n);
+else
+  error('usage: x = are(a,b,c)')
+endif
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/is_controllable.m	Tue Aug 10 20:56:55 1993 +0000
@@ -0,0 +1,47 @@
+function retval = is_controllable (a,b,tol)
+
+# usage: is_controllable (a,b{,tol})
+#
+# returns 1 the pair(a,b) is controllable, then return value is the 
+# dimension of x.  0therwise, returns a value of 0
+#
+# See also: size, rows, columns, length, is_matrix, is_scalar, is_vector
+
+# This should really use the method below, but I'm being lazy for now:
+#
+# Controllability is determined by applying Arnoldi iteration with complete
+# re-orthogonalization to obtain an orthogonal basis of the Krylov subspace
+#                   n-1
+# span([b,a*b,...,a^   b]).  tol is a roundoff paramter, 
+# set to 2*eps if omitted
+
+  if ((nargin == 2) || (nargin == 3))
+    n = is_square(a);
+    [nr, nc] = size (b);
+    if((n == 0) || (n != nr) || (nc == 0))
+      retval = 0;
+    else
+
+      m = b;
+      tmp = b;
+      for ii=1:(n-1)
+        tmp = a*tmp;
+        m = [m,tmp];
+      endfor
+
+      # if n is of any significant size, m will be low rank, so be careful!
+      if(nargin == 3)
+        if(is_scalar(tol))
+          retval = (rank(m,tol) == n);
+        else
+          error('is_controllable: tol must be a scalar')
+        endif
+      else
+        retval = (rank(m) == n);
+      endif
+    endif
+  else
+    error ("usage: is_controllable (a,b)");
+  endif
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/is_observable.m	Tue Aug 10 20:56:55 1993 +0000
@@ -0,0 +1,18 @@
+function retval = is_observable (a,c,tol)
+
+# usage: is_observable (a,c{,tol})
+#
+# returns 1 the pair(a,c) is observable, then return value is the 
+# dimension of x.  0therwise, returns a value of 0
+#
+# See also: size, rows, columns, length, is_matrix, is_scalar, is_vector
+
+  if (nargin == 2)
+    retval = is_controllable(a',c');
+  elseif (nargin == 3)
+    retval = iscontrollable(a',c',tol);
+  else
+    error ("usage: is_observable (a,c{,tol})");
+  endif
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/is_square.m	Tue Aug 10 20:56:55 1993 +0000
@@ -0,0 +1,21 @@
+function retval = is_square (x)
+
+# usage: is_square (x)
+#
+# If x is square, then return value is the dimension of x.
+# otherwise, returns a value of 0
+#
+# See also: size, rows, columns, length, is_matrix, is_scalar, is_vector
+
+  if (nargin == 1)
+    [nr, nc] = size (x);
+    if( nr == nc) 
+      retval = nr;
+    else
+      retval = 0;
+    endif
+  else
+    error ("usage: is_square (x)");
+  endif
+
+endfunction