--- a/scripts/control/base/__bodquist__.m
+++ b/scripts/control/base/__bodquist__.m
@@ -61,36 +61,36 @@
   endif
 
   ## check each argument to see if it's in the correct form
-  if (!isstruct(sys))
-    error("sys must be a system data structure");
+  if (! isstruct (sys))
+    error ("sys must be a system data structure");
   endif
 
   ## let __freqresp__ determine w if it's not already given
-  USEW = freqchkw(w);
+  USEW = freqchkw (w);
 
   ## get initial dimensions (revised below if sysprune is called)
-  [nn,nz,mm,pp ] = sysdimensions(sys);
+  [nn, nz, mm, pp] = sysdimensions (sys);
 
   ## check for an output vector and to see whether it`s correct
-  if (!isempty(outputs))
-    if (isempty(inputs))
+  if (! isempty (outputs))
+    if (isempty (inputs))
       inputs = 1:mm;                    # use all inputs
-      warning([rname,": outputs specified but not inputs"]);
-    elseif(is_signal_list(inputs) | ischar(inputs))
-      inputs = sysidx(sys,"in",inputs);
+      warning ("%s: outputs specified but not inputs", rname);
+    elseif (is_signal_list (inputs) || ischar (inputs))
+      inputs = sysidx (sys, "in", inputs);
     endif
-    if(is_signal_list(outputs) | ischar(outputs))
-      outputs = sysidx(sys,"out",outputs);
+    if (is_signal_list (outputs) || ischar (outputs))
+      outputs = sysidx (sys, "out", outputs);
     end
-    sys = sysprune(sys,outputs,inputs);
-    [nn,nz,mm,pp ] = sysdimensions(sys);
+    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 = zp(zer,pol,k,tsam,inname,outname);
+  if (is_siso (sys) && strcmp (sysgettype (sys), "ss"))
+    [zer, pol, k, tsam, inname, outname] = sys2zp (sys);
+    sys = zp (zer, pol, k, tsam, inname, outname);
   endif
 
   ## get system frequency response
@@ -98,46 +98,49 @@
 
   phase = arg(f)*180.0/pi;
 
-  if(!USEW)
+  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 = find(pd > dphase);
+    while (pcnt)
+      pd = abs (diff (phase));                    # phase variation
+      pdbig = find (pd > dphase);
+
+      ## relative variation
+      lp = length (f);
+      lp1 = lp-1;
 
-      lp = length(f);  lp1 = lp-1;              # relative variation
-      fd = abs(diff(f));
-      fm = max(abs([f(1:lp1); f(2:lp)]));
-      fdbig = find(fd > fm/10);
+      fd = abs (diff (f));
+      fm = max (abs ([f(1:lp1); f(2:lp)]));
+      fdbig = find (fd > fm/10);
 
-      bigpts = union(fdbig, pdbig);
+      bigpts = union (fdbig, pdbig);
 
-      if(isempty(bigpts) )
+      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));
+        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));
+        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
+        wnew = wnew(find (wnew != 0));
+        wnew = sort (wnew);
+        wnew = create_set (wnew);
+        if (isempty (wnew))   # all small crossovers
           pcnt = 0;
         else
 	  ## get new freq resp points, combine with old, and sort.
@@ -153,12 +156,12 @@
   endif
 
   ## ensure unique frequency values
-  [w,idx] = sort(w);
+  [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_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);
 

--- a/scripts/control/base/__freqresp__.m
+++ b/scripts/control/base/__freqresp__.m
@@ -50,82 +50,80 @@
   ## SYS_INTERNAL accesses members of system data structure
 
   ## Check Args
-  if ((nargin < 2) || (nargin > 4))
+  if (nargin < 2 || nargin > 4)
     print_usage ();
-  elseif (USEW & (nargin < 3) )
+  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");
+  elseif (USEW && isempty (w))
+    warning ("USEW = 1 but w is empty; setting USEW=0");
     USEW = 0;
   endif
 
-  DIGITAL = is_digital(sys);
+  DIGITAL = is_digital (sys);
 
   ## compute default w if needed
-  if(!USEW)
-    if(is_siso(sys))
-	sys = sysupdate(sys,"zp");
-	[zer,pol] = sys2zp(sys);
+  if (! USEW)
+    if (is_siso (sys))
+      sys = sysupdate (sys, "zp");
+      [zer, pol] = sys2zp (sys);
     else
-	zer = tzero(sys);
-	pol = eig(sys2ss(sys));
+      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);
+    [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
+    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));
+  if (DIGITAL)
+    jw = exp (i*w*sysgettsam(sys));
   else
     jw = i*w;
   endif
 
-  [nn,nz,mm,pp] = sysdimensions(sys);
+  [nn, nz, mm, pp] = sysdimensions (sys);
 
   ## now compute the frequency response - divide by zero yields a warning
-  if (strcmp(sysgettype(sys),"zp"))
+  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));
+    [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)
+    for ii = (l1+1):length(pol)
 	ff = ff ./ (jw - pol(ii));
     endfor
     ff = ff*sysk;
-
-  elseif (strcmp(sysgettype(sys),"tf"))
+  elseif (strcmp (sysgettype (sys), "tf"))
     ## transfer function form
-    [num,den] = sys2tf(sys);
-    ff = polyval(num,jw)./polyval(den,jw);
+    [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);
+    [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;
+    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);
+    [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);
+    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));
+  w = reshape (w, 1, length (w));
+  ff = reshape (ff, 1, length (ff));
 
 endfunction
 

--- a/scripts/control/base/__stepimp__.m
+++ b/scripts/control/base/__stepimp__.m
@@ -269,5 +269,5 @@
     y = [];
     t = [];
   endif
-  ## printf("##STEPIMP-DEBUG: gratulations, successfull completion.\n");
+
 endfunction  

--- a/scripts/control/base/are.m
+++ b/scripts/control/base/are.m
@@ -70,10 +70,11 @@
 
   if (nargin == 3 || nargin == 4)
     if (nargin == 4)
-      if (! (strcmp (opt, "N") || strcmp (opt, "P") ...
-             || strcmp (opt, "S") || strcmp (opt, "B") ...
-             || strcmp (opt, "n") || strcmp (opt, "p") ...
-             || strcmp (opt, "s") || strcmp (opt, "b")))
+      if (! (ischar (opt)
+	     && (strcmp (opt, "N") || strcmp (opt, "P")
+		 || strcmp (opt, "S") || strcmp (opt, "B")
+		 || strcmp (opt, "n") || strcmp (opt, "p")
+		 || strcmp (opt, "s") || strcmp (opt, "b"))))
         warning ("are: opt has an invalid value; setting to B");
         opt = "B";
       endif

--- a/scripts/control/base/ctrb.m
+++ b/scripts/control/base/ctrb.m
@@ -47,23 +47,23 @@
 
   if (nargin == 2)
     a = sys;
-  elseif (nargin == 1 && isstruct(sys))
-    sysupdate(sys,"ss");
-    [a,b] = sys2ss(sys);
+  elseif (nargin == 1 && isstruct (sys))
+    sysupdate (sys, "ss");
+    [a, b] = sys2ss (sys);
   else
     print_usage ();
   endif
 
-  if (!is_abcd(a,b))
+  if (! is_abcd (a, b))
     Qs = [];
   else
     ## no need to check dimensions, we trust is_abcd().
-    [na, ma] = size(a);
+    [na, ma] = size (a);
     ## using imb avoids name conflict with the "mb" function
-    [inb, imb] = size(b);
-    Qs = zeros(na, ma*imb);
+    [inb, imb] = size (b);
+    Qs = zeros (na, ma*imb);
     for i = 1:na
-      Qs(:, (i-1)*imb+1:i*imb) = b;
+      Qs(:,(i-1)*imb+1:i*imb) = b;
       b = a * b;
     endfor
   endif

--- a/scripts/control/base/damp.m
+++ b/scripts/control/base/damp.m
@@ -35,51 +35,53 @@
 
   ## assume a continuous system
   DIGITAL = 0;
-  if(nargin < 1 || nargin > 2)
+
+  if (nargin < 1 || nargin > 2)
     print_usage ();
   endif
-  if(isstruct(p))
+
+  if (isstruct (p))
     if (nargin != 1)
       error("damp: when p is a system, tsamp parameter is not allowed.");
     endif
-    [aa, b, c, d, t_samp] = sys2ss(p);
-    DIGITAL = is_digital(p);
+    [aa, b, c, d, t_samp] = sys2ss (p);
+    DIGITAL = is_digital (p);
   else
     aa = p;
     if (nargin == 2)
-        DIGITAL = 1;
-        t_samp = tsam;
+      DIGITAL = 1;
+      t_samp = tsam;
     endif
   endif
-  if (!issquare(aa))
-    error("damp: Matrix p is not square.")
+  if (! issquare (aa))
+    error ("damp: Matrix p is not square.")
   endif
   if (DIGITAL && t_samp <= 0.0)
-    error("damp: Sampling time tsam must not be <= 0.")
+    error ("damp: Sampling time tsam must not be <= 0.")
   endif
 
   ## all checks done.
-  e = eig(aa);
-  [n, m] = size(aa);
+  e = eig (aa);
+  [n, m] = size (aa);
   if (DIGITAL)
-    printf("  (discrete system with sampling time %f)\n", t_samp);
+    printf ("  (discrete system with sampling time %f)\n", t_samp);
   endif
-  printf("............... Eigenvalue ...........     Damping     Frequency\n");
-  printf("--------[re]---------[im]--------[abs]----------------------[Hz]\n");
+  printf ("............... Eigenvalue ...........     Damping     Frequency\n");
+  printf ("--------[re]---------[im]--------[abs]----------------------[Hz]\n");
   for i = 1:n
     pole = e(i);
     cpole = pole;
     if (DIGITAL)
-      cpole = log(pole) / t_samp;
+      cpole = log (pole) / t_samp;
     endif
-    d0 = -cos(atan2(imag(cpole), real(cpole)));
-    f0 = 0.5 / pi * abs(cpole);
-    if (abs(imag(cpole)) < eps)
-      printf("%12f         ---  %12f  %10f  %12f\n",
-             real(pole), abs(pole), d0, f0);
+    d0 = -cos (atan2 (imag (cpole), real (cpole)));
+    f0 = 0.5 / pi * abs (cpole);
+    if (abs (imag (cpole)) < eps)
+      printf ("%12f         ---  %12f  %10f  %12f\n",
+             real (pole), abs (pole), d0, f0);
     else
-      printf("%12f %12f %12f  %10f  %12f\n",
-             real(pole), imag(pole), abs(pole), d0, f0);
+      printf ("%12f %12f %12f  %10f  %12f\n",
+             real (pole), imag (pole), abs (pole), d0, f0);
     endif
   endfor
 

--- a/scripts/control/base/dare.m
+++ b/scripts/control/base/dare.m
@@ -78,9 +78,11 @@
 
 function x = dare (a, b, q, r, opt)
 
-  if (nargin == 4 | nargin == 5)
+  if (nargin == 4 || nargin == 5)
     if (nargin == 5)
-      if (opt != "N" || opt != "P" || opt != "S" || opt != "B")
+      if (! (ischar (opt)
+	     && (strcmp (opt, "N") || strcmp (opt, "P")
+		 || strcmp (opt, "S") || strcmp (opt, "B"))))
         warning ("dare: opt has an invalid value -- setting to B");
         opt = "B";
       endif
@@ -88,22 +90,20 @@
       opt = "B";
     endif
 
-    
     if ((p = issquare (q)) == 0)
       q = q'*q;
     endif
 
     ##Checking positive definiteness
-    if (isdefinite(r)<=0)
-      error("dare: r not positive definite");
+    if (isdefinite (r) <= 0)
+      error ("dare: r not positive definite");
     end
-    if (isdefinite(q)<0)
-      error("dare: q not positive semidefinite");
+    if (isdefinite (q) < 0)
+      error ("dare: q not positive semidefinite");
     end
 
-
     ## Check r dimensions.
-    [n,m] = size(b);
+    [n, m] = size (b);
     if ((m1 = issquare (r)) == 0)
       error ("dare: r is not square");
     elseif (m1 != m)
@@ -112,8 +112,11 @@
 
     s1 = [a, zeros(n) ; -q, eye(n)];
     s2 = [eye(n), (b/r)*b' ; zeros(n), a'];
-    [c,d,s1,s2] = balance(s1,s2,opt);
-    [aa,bb,u,lam] = qz(s1,s2,"S");
+
+    [c, d, s1, s2] = balance (s1, s2, opt);
+
+    [aa, bb, u, lam] = qz (s1, s2, "S");
+
     u = d*u;
     n1 = n+1;
     n2 = 2*n;

--- a/scripts/control/base/dcgain.m
+++ b/scripts/control/base/dcgain.m
@@ -30,26 +30,30 @@
 
 function gm = dcgain (sys, tol)
 
-  if((nargin < 1) || (nargin > 2) || (nargout > 1))
+  if (nargin < 1 || nargin > 2 || nargout > 1)
     print_usage ();
   endif
-  if(!isstruct(sys))
-    error("dcgain: first argument is not a system data structure.")
+  if (! isstruct (sys))
+    error ("dcgain: first argument is not a system data structure.")
   endif
-  sys = sysupdate(sys, "ss");
-  [aa,bb,cc,dd] = sys2ss(sys);
-  if (is_digital(sys))  aa = aa - eye(size(aa));  endif
-  if (nargin == 1)  tol = 1.0e-10;  endif
-  r = rank(aa, tol);
-  if (r < rows(aa))
+  sys = sysupdate (sys, "ss");
+  [aa, bb, cc, dd] = sys2ss (sys);
+  if (is_digital (sys))
+    aa = aa - eye (size (aa));
+  endif
+  if (nargin == 1)
+    tol = 1.0e-10;
+  endif
+  r = rank (aa, tol);
+  if (r < rows (aa))
     gm = [];
   else
     gm = -cc / aa * bb + dd;
   endif
-  if(!is_stable(sys))
-    [nn,nz,mm,pp] = sysdimensions(sys);
-    warning("dcgain: unstable system; dimensions [nc=%d,nz=%d,mm=%d,pp=%d]", ...
-      nn,nz,mm,pp);
+  if (! is_stable (sys))
+    [nn, nz, mm, pp] = sysdimensions (sys);
+    warning ("dcgain: unstable system; dimensions: nc=%d, nz=%d, mm=%d, pp=%d",
+	     nn, nz, mm, pp);
   endif
 
 endfunction

--- a/scripts/control/base/dgram.m
+++ b/scripts/control/base/dgram.m
@@ -62,7 +62,7 @@
 function m = dgram (a, b)
 
   if (nargin != 2)
-    print_usage();
+    print_usage ();
   endif
 
   ## let dlyap do the error checking...

--- a/scripts/control/base/dlqr.m
+++ b/scripts/control/base/dlqr.m
@@ -128,16 +128,17 @@
     qo = q;
   endif
 
-  ## Checking stabilizability and detectability (dimensions are checked inside these calls)
+  ## Checking stabilizability and detectability (dimensions are checked
+  ## inside these calls).
   tol = 200*eps;
-  if (is_stabilizable (ao, b,tol,1) == 0)
+  if (is_stabilizable (ao, b, tol, 1) == 0)
     error ("dlqr: (a,b) not stabilizable");
   endif
-  dflag = is_detectable (ao, qo, tol,1);
-  if ( dflag == 0)
+  dflag = is_detectable (ao, qo, tol, 1);
+  if (dflag == 0)
     warning ("dlqr: (a,q) not detectable");
-  elseif ( dflag == -1)
-    error("dlqr: (a,q) has non minimal modes near unit circle");
+  elseif (dflag == -1)
+    error ("dlqr: (a,q) has non minimal modes near unit circle");
   end
 
   ## Compute the Riccati solution
@@ -145,6 +146,4 @@
   k = (r+b'*p*b)\(b'*p*a + s');
   e = eig (a - b*k);
 
-
 endfunction
-

--- a/scripts/control/base/dre.m
+++ b/scripts/control/base/dre.m
@@ -97,72 +97,79 @@
 
 function [tvals, Plist] = dre (sys, Q, R, Qf, t0, tf, Ptol, maxits)
 
-  if(nargin < 6 | nargin > 8 | nargout != 2)
+  if (nargin < 6 || nargin > 8 || nargout != 2)
     print_usage ();
-  elseif(!isstruct(sys))
-    error("sys must be a system data structure")
-  elseif(is_digital(sys))
-    error("sys must be a continuous time system")
-  elseif(!ismatrix(Q) | !ismatrix(R) | !ismatrix(Qf))
-    error("Q, R, and Qf must be matrices.");
-  elseif(!isscalar(t0) | !isscalar(tf))
-    error("t0 and tf must be scalars")
-  elseif(t0 >= tf)              error("t0=%e >= tf=%e",t0,tf);
-  elseif(nargin == 6)           Ptol = 0.1;
-  elseif(!isscalar(Ptol))      error("Ptol must be a scalar");
-  elseif(Ptol <= 0)             error("Ptol must be positive");
+  elseif (! isstruct (sys))
+    error ("sys must be a system data structure")
+  elseif (is_digital (sys))
+    error ("sys must be a continuous time system")
+  elseif (! ismatrix (Q) || ! ismatrix (R) || ! ismatrix (Qf))
+    error ("Q, R, and Qf must be matrices");
+  elseif (! isscalar (t0) || ! isscalar (tf))
+    error ("t0 and tf must be scalars")
+  elseif (t0 >= tf)
+    error ("t0=%e >= tf=%e", t0, tf);
+  elseif (nargin < 7)
+    Ptol = 0.1;
+  elseif (! isscalar (Ptol))
+    error ("Ptol must be a scalar");
+  elseif (Ptol <= 0)
+    error ("Ptol must be positive");
   endif
 
-  if(nargin < 8) maxits = 10;
-  elseif(!isscalar(maxits))    error("maxits must be a scalar");
-  elseif(maxits <= 0)           error("maxits must be positive");
+  if (nargin < 8)
+    maxits = 10;
+  elseif (! isscalar (maxits))
+    error ("maxits must be a scalar");
+  elseif (maxits <= 0)
+    error ("maxits must be positive");
   endif
-  maxits = ceil(maxits);
+  maxits = ceil (maxits);
 
-  [aa,bb] = sys2ss(sys);
-  nn = sysdimensions(sys,"cst");
-  mm = sysdimensions(sys,"in");
-  pp = sysdimensions(sys,"out");
+  [aa, bb] = sys2ss (sys);
+  nn = sysdimensions (sys, "cst");
+  mm = sysdimensions (sys, "in");
+  pp = sysdimensions (sys, "out");
 
-  if(size(Q) != [nn, nn])
-    error("Q(%dx%d); sys has %d states",rows(Q),columns(Q),nn);
-  elseif(size(Qf) != [nn, nn])
-    error("Qf(%dx%d); sys has %d states",rows(Qf),columns(Qf),nn);
-  elseif(size(R) != [mm, mm])
-    error("R(%dx%d); sys has %d inputs",rows(R),columns(R),mm);
+  if (size (Q) != [nn, nn])
+    error ("Q(%dx%d); sys has %d states", rows (Q), columns (Q), nn);
+  elseif (size (Qf) != [nn, nn])
+    error ("Qf(%dx%d); sys has %d states", rows (Qf), columns (Qf), nn);
+  elseif (size (R) != [mm, mm])
+    error ("R(%dx%d); sys has %d inputs", rows (R), columns (R), mm);
   endif
 
   ## construct Hamiltonian matrix
   H = [aa , -(bb/R)*bb' ; -Q, -aa'];
 
   ## select time step to avoid numerical overflow
-  fast_eig = max(abs(eig(H)));
-  tc = log(10)/fast_eig;
-  nst = ceil((tf-t0)/tc);
-  tvals = -linspace(-tf,-t0,nst);
-  Plist = list(Qf);
-  In = eye(nn);
+  fast_eig = max (abs (eig (H)));
+  tc = log (10) / fast_eig;
+  nst = ceil ((tf-t0)/tc);
+  tvals = -linspace (-tf, -t0, nst);
+  Plist = list (Qf);
+  In = eye (nn);
   n1 = nn+1;
   n2 = nn+nn;
   done = 0;
-  while(!done)
+  while (! done)
     done = 1;      # assume this pass will do the job
     ## sort time values in reverse order
-    tvals = -sort(-tvals);
-    tvlen = length(tvals);
+    tvals = -sort (-tvals);
+    tvlen = length (tvals);
     maxerr = 0;
     ## compute new values of P(t); recompute old values just in case
-    for ii=2:tvlen
-      uv_i_minus_1 = [ In ; Plist{ii-1} ];
+    for ii = 2:tvlen
+      uv_i_minus_1 = [In; Plist{ii-1}];
       delta_t = tvals(ii-1) - tvals(ii);
-      uv = expm(-H*delta_t)*uv_i_minus_1;
+      uv = expm (-H*delta_t)*uv_i_minus_1;
       Qi = uv(n1:n2,1:nn)/uv(1:nn,1:nn);
       Plist(ii) = (Qi+Qi')/2;
       ## check error
-      Perr = norm(Plist{ii} - Plist{ii-1})/norm(Plist{ii});
-      maxerr = max(maxerr,Perr);
-      if(Perr > Ptol)
-        new_t = mean(tvals([ii,ii-1]));
+      Perr = norm (Plist{ii} - Plist{ii-1})/norm(Plist{ii});
+      maxerr = max (maxerr,Perr);
+      if (Perr > Ptol)
+        new_t = mean (tvals([ii,ii-1]));
         tvals = [tvals, new_t];
         done = 0;
       endif
@@ -170,11 +177,11 @@
 
     ## check number of iterations
     maxits = maxits - 1;
-    done = done+(maxits==0);
+    done = done + (maxits == 0);
   endwhile
-  if(maxerr > Ptol)
-    warning("dre: \n\texiting with%4d points, max rel chg. =%e, Ptol=%e\n", ...
-          tvlen,maxerr,Ptol);
+  if (maxerr > Ptol)
+    warning ("dre: exiting with %d points, max rel chg. = %e, Ptol = %e",
+             tvlen, maxerr, Ptol);
     tvals = tvals(1:length(Plist));
   endif
 

--- a/scripts/control/base/impulse.m
+++ b/scripts/control/base/impulse.m
@@ -56,11 +56,7 @@
 
 function [y, t] = impulse (sys, inp, tstop, n)
 
-  if ((nargin < 1) || (nargin > 4))
-    print_usage ();
-  endif
-
-  if (nargout > 2)
+  if (nargin < 1 || nargin > 4 || nargout > 2)
     print_usage ();
   endif
 

--- a/scripts/control/base/lqe.m
+++ b/scripts/control/base/lqe.m
@@ -92,7 +92,7 @@
 
 function [k, p, e] = lqe (a, g, c, sigw, sigv, zz)
 
-  if ( (nargin != 5) && (nargin != 6))
+  if (nargin != 5 && nargin != 6)
     error ("lqe: invalid number of arguments");
   endif
 

--- a/scripts/control/base/lqg.m
+++ b/scripts/control/base/lqg.m
@@ -70,89 +70,98 @@
 
 function [K, Q1, P1, Ee, Er] = lqg (sys, Sigw, Sigv, Q, R, input_list)
 
-  if ( (nargin < 5) | (nargin > 6))
+  if (nargin < 5 || nargin > 6)
     print_usage ();
-
-  elseif(!isstruct(sys) )
-    error("sys must be in system data structure");
+  elseif (! isstruct (sys))
+    error ("sys must be in system data structure");
   endif
 
-  DIG = is_digital(sys);
-  [A,B,C,D,tsam,n,nz,stname,inname,outname] = sys2ss(sys);
-  [n,nz,nin,nout] = sysdimensions(sys);
-  if(nargin == 5)
+  DIG = is_digital (sys);
+  [A, B, C, D, tsam, n, nz, stname, inname, outname] = sys2ss (sys);
+  [n, nz, nin, nout] = sysdimensions (sys);
+  if (nargin == 5)
     ## construct default input_list
     input_list = (columns(Sigw)+1):nin;
   endif
 
-  if( !(n+nz) )
-      error(["lqg: 0 states in system"]);
+  if (! (n+nz))
+    error("lqg: 0 states in system");
 
-  elseif(nin != columns(Sigw)+ columns(R))
-    error(["lqg: sys has ",num2str(nin)," inputs, dim(Sigw)=", ...
-          num2str(columns(Sigw)),", dim(u)=",num2str(columns(R))])
+  elseif (nin != columns (Sigw) + columns (R))
+    error ("lqg: sys has %d inputs, dim(Sigw)=%d, dim(u)=%d",
+	   nin, columns (Sigw), columns (R));
 
-  elseif(nout != columns(Sigv))
-    error(["lqg: sys has ",num2str(nout)," outputs, dim(Sigv)=", ...
-          num2str(columns(Sigv)),")"])
+  elseif (nout != columns (Sigv))
+    error ("lqg: sys has %d outputs, dim(Sigv)=%d", nout, columns (Sigv));
   endif
 
   ## check for names of signals
-  if(is_signal_list(input_list) | ischar(input_list))
-    input_list = sysidx(sys,"in",input_list);
+  if (is_signal_list (input_list) || ischar (input_list))
+    input_list = sysidx (sys, "in", input_list);
   endif
 
-  if(length(input_list) != columns(R))
-    error(["lqg: length(input_list)=",num2str(length(input_list)), ...
-          ", columns(R)=", num2str(columns(R))]);
+  if (length (input_list) != columns (R))
+    error ("lqg: length(input_list)=%d, columns(R)=%d",
+	   length (input_list), columns (R));
+  endif
+
+  if (! issquare (Sigw))
+    error ("lqg: Sigw is not square");
   endif
 
-  varname = {"Sigw","Sigv","Q","R"};
-  for kk=1:length(varname);
-    eval(sprintf("chk = issquare(%s);",varname{kk}));
-    if(! chk ) error("lqg: %s is not square",varname{kk}); endif
-  endfor
+  if (! issquare (Sigv))
+    error ("lqg: Sigv is not square");
+  endif
+
+  if (! issquare (Q))
+    error ("lqg: Q is not square");
+  endif
+
+  if (! issquare (R))
+    error ("lqg: Q is not square");
+  endif
 
   ## permute (if need be)
-  if(nargin == 6)
-    all_inputs = sysreorder(nin,input_list);
+  if (nargin == 6)
+    all_inputs = sysreorder (nin, input_list);
     B = B(:,all_inputs);
-    inname = inname(all_inputs);
+    inname = inname (all_inputs);
   endif
 
   ## put parameters into correct variables
-  m1 = columns(Sigw);
+  m1 = columns (Sigw);
   m2 = m1+1;
   G = B(:,1:m1);
   B = B(:,m2:nin);
 
   ## now we can just do the design; call dlqr and dlqe, since all matrices
   ## are not given in Cholesky factor form (as in h2syn case)
-  if(DIG)
-    [Ks, P1, Er] = dlqr(A,B,Q,R);
-    [Ke, Q1, jnk, Ee] = dlqe(A,G,C,Sigw,Sigv);
+  if (DIG)
+    [Ks, P1, Er] = dlqr (A, B, Q, R);
+    [Ke, Q1, jnk, Ee] = dlqe (A, G, C, Sigw, Sigv);
   else
-    [Ks, P1, Er] = lqr(A,B,Q,R);
-    [Ke, Q1, Ee] = lqe(A,G,C,Sigw,Sigv);
+    [Ks, P1, Er] = lqr (A, B, Q, R);
+    [Ke, Q1, Ee] = lqe (A, G, C, Sigw, Sigv);
   endif
   Ac = A - Ke*C - B*Ks;
   Bc = Ke;
   Cc = -Ks;
-  Dc = zeros(rows(Cc),columns(Bc));
+  Dc = zeros (rows (Cc), columns (Bc));
 
   ## fix state names
-  stname1 = strappend(stname,"_e");
+  stname1 = strappend (stname, "_e");
 
   ## fix controller output names
-  outname1 = strappend(inname(m2:nin),"_K");
+  outname1 = strappend (inname(m2:nin), "_K");
 
   ## fix controller input names
-  inname1 = strappend(outname,"_K");
+  inname1 = strappend (outname, "_K");
 
-  if(DIG)
-    K = ss(Ac,Bc,Cc,Dc,tsam,n,nz,stname1,inname1,outname1,1:rows(Cc));
+  if (DIG)
+    K = ss (Ac, Bc, Cc, Dc, tsam, n, nz, stname1, inname1, outname1,
+	    1:rows(Cc));
   else
-    K = ss(Ac,Bc,Cc,Dc,tsam,n,nz,stname,inname1,outname1);
+    K = ss (Ac, Bc, Cc, Dc, tsam, n, nz, stname, inname1, outname1);
   endif
 
 endfunction

--- a/scripts/control/base/lqr.m
+++ b/scripts/control/base/lqr.m
@@ -119,7 +119,7 @@
 
   ## disp("lqr: entry");
 
-  if ((nargin != 4) && (nargin != 5))
+  if (nargin != 4 && nargin != 5)
     error ("lqr: invalid number of arguments");
   endif
 
@@ -135,19 +135,19 @@
   endif
 
   ## Check q.
-  if ( ((n1 = issquare (q)) == 0) || (n1 != n))
+  if ((n1 = issquare (q)) == 0 || n1 != n)
     error ("lqr: q must be square and conformal with a");
   endif
 
   ## Check r.
-  if ( ((m1 = issquare(r)) == 0) || (m1 != m))
+  if ((m1 = issquare(r)) == 0 || m1 != m)
     error ("lqr: r must be square and conformal with column dimension of b");
   endif
 
   ## Check if n is there.
   if (nargin == 5)
     [n1, m1] = size (s);
-    if ( (n1 != n) || (m1 != m))
+    if (n1 != n || m1 != m)
       error ("lqr: z must be identically dimensioned with b");
     endif
 
@@ -162,7 +162,7 @@
 
   ## Check that q, (r) are symmetric, positive (semi)definite
 
-  if (issymmetric (q) && issymmetric (r) ...
+  if (issymmetric (q) && issymmetric (r)
       && all (eig (q) >= 0) && all (eig (r) > 0))
     p = are (ao, (b/r)*b', qo);
     k = r\(b'*p + s');
@@ -171,5 +171,4 @@
     error ("lqr: q (r) must be symmetric positive (semi) definite");
   endif
 
-  ## disp("lqr: exit");
 endfunction

--- a/scripts/control/base/lsim.m
+++ b/scripts/control/base/lsim.m
@@ -40,38 +40,40 @@
 
 function [y, x] = lsim (sys, u, t, x0)
 
-  if((nargin < 3)||(nargin > 4))
+  if (nargin < 3 || nargin > 4)
     print_usage ();
   endif
 
-  if(!isstruct(sys))
-    error("sys must be in system data structure");
+  if (! isstruct (sys))
+    error ("sys must be in system data structure");
   endif
 
-  sys = sysupdate(sys,"ss");
+  sys = sysupdate (sys,"ss");
 
-  [ncstates, ndstates, nin, nout] = sysdimensions(sys);
-  [a,b,c,d] = sys2ss(sys);
+  [ncstates, ndstates, nin, nout] = sysdimensions (sys);
+  [a, b, c, d] = sys2ss (sys);
 
-  if (nargin == 3)     x0 = zeros(columns(a),1);        endif
+  if (nargin == 3)
+    x0 = zeros (columns (a), 1);
+  endif
 
-  if(rows(u) ~= length(t))
-    error("lsim: There should be an input value (row) for each time instant");
+  if (rows (u) != length (t))
+    error ("lsim: There should be an input value (row) for each time instant");
   endif
-  if(columns(u) ~= columns(d))
-    error("lsim: U and d should have the same number of inputs");
+  if (columns (u) != columns (d))
+    error ("lsim: U and d should have the same number of inputs");
   endif
-  if(columns(x0) > 1)
-    error("lsim: Initial condition vector should have only one column");
+  if (columns (x0) > 1)
+    error ("lsim: Initial condition vector should have only one column");
   endif
-  if(rows(x0) > rows(a))
-    error("lsim: Initial condition vector is too large");
+  if (rows (x0) > rows p(a))
+    error ("lsim: Initial condition vector is too large");
   endif
 
   Ts = 0;
   t(2)-t(1);
   u=u';
-  n = max(size(t));
+  n = max (size (t));
 
   for ii = 1:(n-1)
 
@@ -81,8 +83,8 @@
     if (abs (t(ii+1) - t(ii) - Ts) > 10 * eps)
       Ts = t(ii+1) - t(ii);
       ## [F,G] = c2d(a,b,Ts);
-      dsys = c2d(sys, Ts);
-      [F,G] = sys2ss(dsys);
+      dsys = c2d (sys, Ts);
+      [F, G] = sys2ss (dsys);
     endif
 
     x(:,ii) = x0;
@@ -93,9 +95,10 @@
   x(:,n) = x0;
 
   y = c*x + d*u;
-  if(nargout == 0)
-   plot(t,y);
-   y=[];
-   x=[];
+  if (nargout == 0)
+   plot (t, y);
+   y = [];
+   x = [];
   endif
+
 endfunction

--- a/scripts/control/base/ltifr.m
+++ b/scripts/control/base/ltifr.m
@@ -59,53 +59,56 @@
 
 function out = ltifr (a, b, w)
 
-  if ((nargin < 2) || (nargin > 3))
+  if (nargin < 2 || nargin > 3)
     error("incorrect number of input arguments");
   endif
 
   if (nargin == 2)
     sys = a;
     w = b;
-    if(!isstruct(sys))
-      error("two arguments: 1st must be a system data structure");
+    if(! isstruct (sys))
+      error ("two arguments: 1st must be a system data structure");
     endif
 
-    if (!isvector(w))
-      error("w must be a vector");
+    if (! isvector (w))
+      error ("w must be a vector");
     endif
 
-    [nn,nz,mm,pp] = sysdimensions(sys);
-    if(mm != 1)       error("sys has %d > 1 inputs",mm); endif
+    [nn, nz, mm, pp] = sysdimensions (sys);
+    if (mm != 1)
+      error("sys has %d > 1 inputs", mm);
+    endif
 
-    [a,b] = sys2ss(sys);
+    [a, b] = sys2ss (sys);
 
   else
 
-    if (columns(a) != rows(b)),
-      error("ltifr:  A(%dx%d), B(%dx%d) not compatibly dimensioned", ...
-        rows(a), columns(a), rows(b), columns(b));
+    if (columns (a) != rows (b)),
+      error ("ltifr:  A(%dx%d), B(%dx%d) not compatibly dimensioned",
+             rows (a), columns(a), rows(b), columns(b));
     endif
 
-    if(columns(b) != 1)
-      error("ltifr: b(%dx%d) must be a single column vector", ...
-        rows(b),columns(b));
+    if (columns (b) != 1)
+      error ("ltifr: b(%dx%d) must be a single column vector",
+             rows(b), columns(b));
     endif
 
-    if (!issquare(a))
-      error("ltifr:  A(%dx$d) must be square.",rows(a),columns(a))
+    if (! issquare (a))
+      error ("ltifr:  A(%dx$d) must be square", rows(a), columns(a))
     endif
 
   endif
 
-  if (!isvector(w))
-    error("w must be a vector");
+  if (! isvector (w))
+    error ("w must be a vector");
   endif
 
-  ey = eye(size(a));
-  lw = length(w);
-  out = ones(columns(a),lw);
+  ey = eye (size (a));
+  lw = length (w);
+  out = ones (columns (a), lw);
 
-  for ii=1:lw,
+  for ii = 1:lw,
     out(:,ii) = (w(ii)*ey-a)\b;
   endfor
+
 endfunction

--- a/scripts/control/base/nichols.m
+++ b/scripts/control/base/nichols.m
@@ -100,7 +100,7 @@
 
   [f, w, sys] = __bodquist__ (sys, w, outputs, inputs, "nichols");
 
-  [stname,inname,outname] = sysgetsignals (sys);
+  [stname, inname, outname] = sysgetsignals (sys);
   systsam = sysgettsam (sys);
 
   ## Get the magnitude and phase of f.

--- a/scripts/control/base/nyquist.m
+++ b/scripts/control/base/nyquist.m
@@ -168,7 +168,7 @@
         th = atan2 (real (df), imag (df)) * 180 / pi;
 
         ## get angle at fmax
-        if (fi == length(f))
+        if (fi == length (f))
 	  fi = fi-1;
 	endif
         thm = th(fi);

--- a/scripts/control/base/obsv.m
+++ b/scripts/control/base/obsv.m
@@ -54,19 +54,19 @@
   if (nargin == 2)
     a = sys;
   elseif (nargin == 1 && isstruct(sys))
-    sysupdate(sys,"ss");
-    [a,b,c] = sys2ss(sys);
+    sysupdate (sys, "ss");
+    [a, b, c] = sys2ss (sys);
   else
     print_usage ();
   endif
 
-  if (!is_abcd(a,c'))
+  if (! is_abcd (a, c'))
     Qb = [];
   else
     ## no need to check dimensions, we trust is_abcd().
-    [na, ma] = size(a);
-    [nc, mc] = size(c);
-    Qb = zeros(na*nc, ma);
+    [na, ma] = size (a);
+    [nc, mc] = size (c);
+    Qb = zeros (na*nc, ma);
     for i = 1:na
       Qb((i-1)*nc+1:i*nc, :) = c;
       c = c * a;

--- a/scripts/control/base/place.m
+++ b/scripts/control/base/place.m
@@ -51,7 +51,7 @@
 
   ## check arguments
 
-  if(! isstruct (sys))
+  if (! isstruct (sys))
     error ("sys must be in system data structure format (see ss)");
   endif
   sys = sysupdate (sys, "ss");    # make sure it has state space form up to date
@@ -66,7 +66,7 @@
   is_digital (sys);
   [n, nz, m, p] = sysdimensions (sys);
   nx = n+nz;    # already checked that it's not a mixed system.
-  if(m != 1)
+  if (m != 1)
     error ("sys has %d inputs; need only 1", m);
   endif
 

--- a/scripts/control/base/rlocus.m
+++ b/scripts/control/base/rlocus.m
@@ -66,22 +66,22 @@
   endif
 
   ## Convert the input to a transfer function if necessary
-  [num,den] = sys2tf(sys);               # extract numerator/denom polyomials
-  lnum = length(num);      
-  lden = length(den);
-  # equalize length of num, den polynomials
-  if(lden < 2)
-    error("system has no poles");
-  elseif(lnum < lden)
+  [num, den] = sys2tf (sys);               # extract numerator/denom polyomials
+  lnum = length (num);
+  lden = length (den);
+  ## equalize length of num, den polynomials
+  if (lden < 2)
+    error ("system has no poles");
+  elseif (lnum < lden)
     num = [zeros(1,lden-lnum), num];  # so that derivative is shortened by one
   endif
 
-  olpol = roots(den);
-  olzer = roots(num);
-  nas = lden -lnum; # number of asymptotes
+  olpol = roots (den);
+  olzer = roots (num);
+  nas = lden - lnum; # number of asymptotes
   maxk = 0;
-  if(nas > 0)
-    cas = ( sum(olpol) - sum(olzer) )/nas;
+  if (nas > 0)
+    cas = (sum (olpol) - sum (olzer)) / nas;
     angles = (2*[1:nas]-1)*pi/nas;
     # printf("rlocus: there are %d asymptotes centered at %f\n", nas, cas);
   else
@@ -91,33 +91,33 @@
 
 
   # compute real axis break points and corresponding gains
-  dnum=polyderiv(num);
-  dden=polyderiv(den);
-  brkp = conv(den, dnum) - conv(num, dden);
-  real_ax_pts = roots(brkp);
-  real_ax_pts = real_ax_pts(find(imag(real_ax_pts) == 0));
-  k_break = -polyval(den,real_ax_pts) ./ polyval(num, real_ax_pts);
-  idx = find(k_break >= 0);
+  dnum = polyderiv (num);
+  dden = polyderiv (den);
+  brkp = conv (den, dnum) - conv (num, dden);
+  real_ax_pts = roots (brkp);
+  real_ax_pts = real_ax_pts(find (imag (real_ax_pts) == 0));
+  k_break = -polyval (den, real_ax_pts) ./ polyval (num, real_ax_pts);
+  idx = find (k_break >= 0);
   k_break = k_break(idx);
   real_ax_pts = real_ax_pts(idx);
-  if(!isempty(k_break))
-    maxk = max(max(k_break),maxk);
+  if (! isempty (k_break))
+    maxk = max (max (k_break), maxk);
   endif
   
-  if(nas == 0)
-    maxk = max(1, 2*maxk);  % get at least some root locus
+  if (nas == 0)
+    maxk = max (1, 2*maxk);  # get at least some root locus
   else
-    # get distance from breakpoints, poles, and zeros to center of asymptotes
-    dmax = 3*max(abs( [vec(olzer); vec(olpol); vec(real_ax_pts)] - cas ));
-    if(dmax == 0)
+    ## get distance from breakpoints, poles, and zeros to center of asymptotes
+    dmax = 3*max (abs ([vec(olzer); vec(olpol); vec(real_ax_pts)] - cas));
+    if (dmax == 0)
       dmax = 1;
     endif
  
     # get gain for dmax along each asymptote, adjust maxk if necessary
-    svals = cas + dmax*exp(j*angles);
-    kvals = -polyval(den,svals) ./ polyval(num, svals);
-    maxk = max(maxk, max(real(kvals)));
-  end
+    svals = cas + dmax * exp (j*angles);
+    kvals = -polyval (den, svals) ./ polyval (num, svals);
+    maxk = max (maxk, max (real (kvals)));
+  endif
   
   ## check for input arguments:
   if (nargin > 2)
@@ -129,8 +129,8 @@
     maxk = max_k;
   endif
   if (nargin > 1)
-    if(increment <= 0)
-      error("increment must be positive");
+    if (increment <= 0)
+      error ("increment must be positive");
     else
       ngain = (maxk-mink)/increment;
     endif
@@ -139,48 +139,49 @@
   endif
 
   ## vector of gains
-  ngain = max(30,ngain);
-  gvec = linspace(mink,maxk,ngain);
-  if(length(k_break))
-    gvec = sort([gvec, vec(k_break)']);
+  ngain = max (30, ngain);
+  gvec = linspace (mink, maxk, ngain);
+  if (length (k_break))
+    gvec = sort ([gvec, vec(k_break)']);
   endif
 
   ## Find the open loop zeros and the initial poles
-  rlzer = roots(num);
+  rlzer = roots (num);
 
   ## update num to be the same length as den
-  lnum = length(num);  
-  if(lnum < lden)
+  lnum = length (num);  
+  if (lnum < lden)
     num = [zeros(1,lden - lnum),num];
   endif
 
   ## compute preliminary pole sets
-  nroots = lden-1;
-  for ii=1:ngain
+  nroots = lden - 1;
+  for ii = 1:ngain
    gain = gvec(ii);
-   rlpol(1:nroots,ii)  = vec(sortcom(roots(den + gain*num)));
+   rlpol(1:nroots,ii) = vec(sortcom (roots (den + gain*num)));
   endfor
 
   ## set smoothing tolerance 
-  smtolx = 0.01*( max(max(real(rlpol))) - min(min(real(rlpol))));
-  smtoly = 0.01*( max(max(imag(rlpol))) - min(min(imag(rlpol))));
-  smtol = max(smtolx, smtoly);
-  rlpol = sort_roots(rlpol,smtolx, smtoly);   % sort according to nearest-neighbor
+  smtolx = 0.01*(max (max (real (rlpol))) - min (min (real (rlpol))));
+  smtoly = 0.01*(max (max (imag (rlpol))) - min (min (imag (rlpol))));
+  smtol = max (smtolx, smtoly);
+  ## sort according to nearest-neighbor
+  rlpol = sort_roots (rlpol, smtolx, smtoly);
 
-  done=(nargin == 4);    # perform a smoothness check
-  while((!done) & ngain < 1000)
+  done = (nargin == 4);    # perform a smoothness check
+  while (! done && ngain < 1000)
     done = 1 ;      # assume done
-    dp = abs(diff(rlpol'))';
-    maxdp = max(dp);
+    dp = abs (diff (rlpol'))';
+    maxdp = max (dp);
     
     ## search for poles whose neighbors are distant
-    if(lden == 2)
-      idx = find(dp > smtol);
+    if (lden == 2)
+      idx = find (dp > smtol);
     else
-      idx = find(maxdp > smtol);
+      idx = find (maxdp > smtol);
     endif
 
-    for ii=1:length(idx)
+    for ii = 1:length(idx)
       i1 = idx(ii);
       g1 = gvec(i1);
       p1 = rlpol(:,i1);
@@ -190,39 +191,40 @@
       p2 = rlpol(:,i2);
 
       ## isolate poles in p1, p2 
-      if( max(abs(p2-p1)) > smtol)
-        newg = linspace(g1,g2,5);
+      if (max (abs (p2-p1)) > smtol)
+        newg = linspace (g1, g2, 5);
         newg = newg(2:4);
-        gvec =  [gvec,newg];
+        gvec = [gvec,newg];
         done = 0;             # need to process new gains
       endif
     endfor
 
     ## process new gain values
-    ngain1 = length(gvec);
-    for ii=(ngain+1):ngain1
+    ngain1 = length (gvec);
+    for ii = (ngain+1):ngain1
       gain = gvec(ii);
-      rlpol(1:nroots,ii)  = vec(sortcom(roots(den + gain*num)));
+      rlpol(1:nroots,ii) = vec(sortcom (roots (den + gain*num)));
     endfor
 
-    [gvec,idx] = sort(gvec);
+    [gvec, idx] = sort (gvec);
     rlpol = rlpol(:,idx);
-    ngain = length(gvec);
-    rlpol = sort_roots(rlpol,smtolx, smtoly);   % sort according to nearest-neighbor
+    ngain = length (gvec);
+    ## sort according to nearest-neighbor
+    rlpol = sort_roots (rlpol, smtolx, smtoly);
   endwhile
   rldata = rlpol;
 
   ## Plot the data
-  if(nargout  == 0)
+  if (nargout  == 0)
     rlpolv = vec(rlpol);
-    axdata = [real(rlpolv),imag(rlpolv); real(olzer), imag(olzer)];
-    axlim = axis2dlim(axdata);
+    axdata = [real(rlpolv), imag(rlpolv); real(olzer), imag(olzer)];
+    axlim = axis2dlim (axdata);
     rldata = [real(rlpolv), imag(rlpolv) ];
-    [stn,inname,outname] = sysgetsignals(sys);
+    [stn, inname, outname] = sysgetsignals (sys);
 
     ## build plot command args pole by pole
 
-    n_rlpol = rows(rlpol);
+    n_rlpol = rows (rlpol);
     nelts = n_rlpol+1;
     if (! isempty (rlzer))
       nelts++;
@@ -236,7 +238,7 @@
     kk = 0;
     # asymptotes first
     if (n_A > 0)
-      len_A = 2*max(abs(axlim));
+      len_A = 2*max (abs (axlim));
       sigma_A = (sum(olpol) - sum(olzer))/n_A;
       for i_A=0:n_A-1
         phi_A = pi*(2*i_A + 1)/n_A;
@@ -250,7 +252,7 @@
       endfor
     endif
     # locus next
-    for ii=1:rows(rlpol)
+    for ii = 1:rows(rlpol)
       args{1,++kk} = real (rlpol (ii,:));
       args{2,kk} = imag (rlpol (ii,:));
       if (ii == 1)
@@ -260,22 +262,22 @@
       endif
     endfor
     # poles and zeros last
-    args{1,++kk} = real(olpol);
-    args{2,kk} = imag(olpol);
+    args{1,++kk} = real (olpol);
+    args{2,kk} = imag (olpol);
     args{3,kk} = "rx;open loop poles;";
-    if (! isempty(rlzer))
-      args{1,++kk} = real(rlzer);
-      args{2,kk} = imag(rlzer);
+    if (! isempty (rlzer))
+      args{1,++kk} = real (rlzer);
+      args{2,kk} = imag (rlzer);
       args{3,kk} = "go;zeros;";
     endif
 
     set (gcf,"visible","off");
     hplt = plot (args{:});
     set (hplt(kk--), "markersize", 2);
-    if (! isempty(rlzer))
-      set(hplt(kk--), "markersize", 2);
+    if (! isempty (rlzer))
+      set (hplt(kk--), "markersize", 2);
     endif
-    for ii=1:rows(rlpol)
+    for ii = 1:rows(rlpol)
       set (hplt(kk--), "linewidth", 2);
     endfor
     legend ("boxon", 2);
@@ -284,42 +286,42 @@
     xlabel (sprintf ("Root locus from %s to %s, gain=[%f,%f]: Real axis",
 		     inname{1}, outname{1}, gvec(1), gvec(ngain)));
     ylabel ("Imag. axis");
-    set (gcf,"visible","on");
+    set (gcf (), "visible","on");
     rldata = [];
   endif
 endfunction
 
 function rlpol = sort_roots (rlpol,tolx, toly)
   # no point sorting of you've only got one pole!
-  if(rows(rlpol) == 1)
-    return
+  if (rows (rlpol) == 1)
+    return;
   endif
 
   # reorder entries in each column of rlpol to be by their nearest-neighbors
-  dp = diff(rlpol')';
-  drp = max(real(dp));
-  dip = max(imag(dp));
-  idx = find( drp > tolx | dip > toly );
-  if(isempty(idx) )
-    return
+  dp = diff (rlpol')';
+  drp = max (real (dp));
+  dip = max (imag (dp));
+  idx = find (drp > tolx | dip > toly);
+  if (isempty (idx))
+    return;
   endif
 
-  [np,ng] = size(rlpol);  # num poles, num gains
+  [np, ng] = size (rlpol);  # num poles, num gains
   for jj = idx
     vals = rlpol(:,[jj,jj+1]);
     jdx = (jj+1):ng;
-    for ii=1:rows(rlpol-1)
+    for ii = 1:rows(rlpol-1)
       rdx = ii:np;
-      dval = abs(rlpol(rdx,jj+1)-rlpol(ii,jj));
-      mindist = min(dval);
-      sidx = min( find ( dval == mindist)) + ii - 1;
-      if( sidx != ii)
-        c1 = norm(diff(vals'));
-        [vals(ii,2), vals(sidx,2)] = swap( vals(ii,2), vals(sidx,2));
-        c2 = norm(diff(vals'));
-        if(c1 > c2 )
-          # perform the swap
-          [rlpol(ii,jdx), rlpol(sidx,jdx)] = swap( rlpol(ii,jdx), rlpol(sidx,jdx));
+      dval = abs (rlpol(rdx,jj+1)-rlpol(ii,jj));
+      mindist = min (dval);
+      sidx = min (find (dval == mindist)) + ii - 1;
+      if (sidx != ii)
+        c1 = norm (diff(vals'));
+        [vals(ii,2), vals(sidx,2)] = swap (vals(ii,2), vals(sidx,2));
+        c2 = norm (diff (vals'));
+        if (c1 > c2)
+          ## perform the swap
+          [rlpol(ii,jdx), rlpol(sidx,jdx)] = swap (rlpol(ii,jdx), rlpol(sidx,jdx));
           vals = rlpol(:,[jj,jj+1]);
         endif
       endif

--- a/scripts/control/base/step.m
+++ b/scripts/control/base/step.m
@@ -57,11 +57,7 @@
 
 function [y, t] = step (sys, inp, tstop, n)
 
-  if ((nargin < 1) || (nargin > 4))
-    print_usage ();
-  endif
-
-  if (nargout > 2)
+  if (nargin < 1 || nargin > 4 || nargout > 2)
     print_usage ();
   endif
 

--- a/scripts/control/base/tzero.m
+++ b/scripts/control/base/tzero.m
@@ -70,23 +70,23 @@
 function [zer, gain] = tzero (A, B, C, D)
 
   ## get A,B,C,D and Asys variables, regardless of initial form
-  if(nargin == 4)
-    Asys = ss(A,B,C,D);
-  elseif( (nargin == 1) && (! isstruct(A)))
-    print_usage ();
-  elseif(nargin != 1)
+  if (nargin == 4)
+    Asys = ss (A, B, C, D);
+  elseif (nargin == 1 && ! isstruct (A))
+    error ("tzero: expecting argument to be system structure");
+  elseif (nargin != 1)
     print_usage ();
   else
     Asys = A;
-    [A,B,C,D] = sys2ss(Asys);
+    [A, B, C, D] = sys2ss (Asys);
   endif
 
   Ao = Asys;                    # save for leading coefficient
-  siso = is_siso(Asys);
-  digital = is_digital(Asys);   # check if it's mixed or not
+  siso = is_siso (Asys);
+  digital = is_digital (Asys);   # check if it's mixed or not
 
   ## see if it's a gain block
-  if(isempty(A))
+  if (isempty (A))
     zer = [];
     gain = D;
     return;
@@ -95,44 +95,50 @@
   ## First, balance the system via the zero computation generalized eigenvalue
   ## problem balancing method (Hodel and Tiller, Linear Alg. Appl., 1992)
 
-  Asys = __zgpbal__ (Asys); [A,B,C,D] = sys2ss(Asys);   # balance coefficients
+  ## balance coefficients
+  Asys = __zgpbal__ (Asys);
+  [A, B, C, D] = sys2ss (Asys);
   meps = 2*eps*norm ([A, B; C, D], "fro");
-  Asys = zgreduce(Asys,meps);  [A, B, C, D] = sys2ss(Asys); # ENVD algorithm
-  if(!isempty(A))
+  ## ENVD algorithm
+  Asys = zgreduce (Asys, meps);
+  [A, B, C, D] = sys2ss (Asys);
+  if (! isempty (A))
     ## repeat with dual system
-    Asys = ss(A', C', B', D');   Asys = zgreduce(Asys,meps);
+    Asys = ss (A', C', B', D');
+    Asys = zgreduce (Asys, meps);
 
     ## transform back
-    [A,B,C,D] = sys2ss(Asys);    Asys = ss(A', C', B', D');
+    [A, B, C, D] = sys2ss (Asys);
+    Asys = ss (A', C', B', D');
   endif
 
   zer = [];                     # assume none
-  [A,B,C,D] = sys2ss(Asys);
-  if( !isempty(C) )
-    [W,r,Pi] = qr([C, D]');
-    [nonz,ztmp] = zgrownorm(r,meps);
-    if(nonz)
+  [A, B, C, D] = sys2ss (Asys);
+  if (! isempty (C))
+    [W, r, Pi] = qr ([C, D]');
+    [nonz, ztmp] = zgrownorm (r, meps);
+    if (nonz)
       ## We can now solve the generalized eigenvalue problem.
-      [pp,mm] = size(D);
-      nn = rows(A);
+      [pp, mm] = size (D);
+      nn = rows (A);
       Afm = [A , B ; C, D] * W';
       Bfm = [eye(nn), zeros(nn,mm); zeros(pp,nn+mm)]*W';
 
       jdx = (mm+1):(mm+nn);
       Af = Afm(1:nn,jdx);
       Bf = Bfm(1:nn,jdx);
-      zer = qz(Af,Bf);
+      zer = qz (Af, Bf);
     endif
   endif
 
-  mz = length(zer);
-  [A,B,C,D] = sys2ss(Ao);               # recover original system
+  mz = length (zer);
+  [A, B, C, D] = sys2ss (Ao);               # recover original system
   ## compute leading coefficient
-  if ( (nargout == 2) && siso)
-    n = rows(A);
-    if ( mz == n)
+  if (nargout == 2 && siso)
+    n = rows (A);
+    if (mz == n)
       gain = D;
-    elseif ( mz < n )
+    elseif (mz < n)
       gain = C*(A^(n-1-mz))*B;
     endif
   else