--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,36 @@
+2007-11-07  Muthiah Annamalai  <muthuspost@gmail.com>
+
+	* control/base/bode_bounds.m, control/base/dgram.m,
+	control/base/dlyap.m, control/base/freqchkw.m,
+	control/base/gram.m, control/base/place.m,
+	control/hinf/hinf_ctr.m, control/hinf/hinfsyn_chk.m,
+	control/hinf/hinfsyn_ric.m, control/system/is_sample.m,
+	control/system/is_signal_list.m, control/system/ss2tf.m,
+	control/system/sys2fir.m, control/system/sysgettsam.m,
+	control/system/sysgettype.m, control/system/sysreorder.m,
+	control/system/tf2sys.m, control/system/zp2tf.m,
+	control/util/axis2dlim.m, control/util/swap.m,
+	control/util/zgfmul.m, control/util/zgfslv.m,
+	control/util/zginit.m, control/util/zgreduce.m,
+	control/util/zgrownorm.m, control/util/zgscal.m,
+	control/util/zgsgiv.m, control/util/zgshsr.m, general/isa.m,
+	geometry/inpolygon.m, linear-algebra/housh.m,
+	miscellaneous/compare_versions.m, miscellaneous/inputname.m,
+	miscellaneous/run.m, quaternion/qconj.m,
+	quaternion/qcoordinate_plot.m, quaternion/qderiv.m,
+	quaternion/qderivmat.m, quaternion/qinv.m, quaternion/qmult.m,
+	quaternion/qtrans.m, quaternion/qtransvmat.m, signal/fractdiff.m,
+	signal/freqz_plot.m, signal/periodogram.m, signal/rectangle_lw.m,
+	signal/rectangle_sw.m, signal/sinc.m, signal/triangle_lw.m,
+	signal/triangle_sw.m, signal/yulewalker.m, sparse/colperm.m,
+	sparse/etreeplot.m, sparse/nonzeros.m, sparse/spalloc.m,
+	sparse/spones.m, sparse/spy.m, specfun/isprime.m,
+	statistics/distributions/empirical_cdf.m,
+	statistics/distributions/empirical_inv.m,
+	statistics/distributions/empirical_pdf.m,
+	statistics/models/logistic_regression_derivatives.m,
+	statistics/models/logistic_regression_likelihood.m: Check nargin.
+
 2007-11-07  David Bateman  <dbateman@free.fr>
 
 	* general/gradient.m: Correctly convert deltax and deltay scalar

--- a/scripts/control/base/bode_bounds.m
+++ b/scripts/control/base/bode_bounds.m
@@ -36,32 +36,51 @@
 
 function [wmin, wmax] = bode_bounds (zer, pol, DIGITAL, tsam)
 
+  if (nargin != 4)
+    print_usage ();
+  endif
+
   ## make sure zer,pol are row vectors
-  if(!isempty(pol)) pol = reshape(pol,1,length(pol)); endif
-  if(!isempty(zer)) zer = reshape(zer,1,length(zer)); endif
+  if (! isempty (pol))
+    pol = reshape (pol, 1, length (pol));
+  endif
+  if (! isempty (zer))
+    zer = reshape (zer, 1, length (zer));
+  endif
 
   ## check for natural frequencies away from omega = 0
   if (DIGITAL)
     ## The 2nd conditions prevents log(0) in the next log command
-    iiz = find(abs(zer - 1) > norm(zer) * eps && abs(zer) > norm(zer) * eps);
-    iip = find(abs(pol - 1) > norm(pol) * eps && abs(pol) > norm(pol) * eps);
+    iiz = find (abs(zer-1) > norm(zer)*eps && abs(zer) > norm(zer)*eps);
+    iip = find (abs(pol-1) > norm(pol)*eps && abs(pol) > norm(pol)*eps);
 
     ## avoid dividing empty matrices, it would work but looks nasty
-    if (!isempty(iiz)) czer = log(zer(iiz))/tsam;
-    else               czer = [];                 endif
+    if (! isempty (iiz))
+      czer = log (zer(iiz))/tsam;
+    else
+      czer = [];
+    endif
 
-    if (!isempty(iip)) cpol = log(pol(iip))/tsam;
-    else               cpol = [];                 endif
-
+    if (! isempty (iip))
+      cpol = log (pol(iip))/tsam;
+    else
+      cpol = [];
+    endif
   else
     ## continuous
-    iip = find((abs(pol)) > (norm(pol) * eps));
-    iiz = find((abs(zer)) > (norm(zer) * eps));
+    iip = find (abs(pol) > norm(pol)*eps));
+    iiz = find (abs(zer) > norm(zer)*eps));
 
-    if(!isempty(zer)) czer = zer(iiz);
-    else              czer = [];                endif
-    if(!isempty(pol)) cpol = pol(iip);
-    else              cpol = [];                endif
+    if (! isempty (zer))
+      czer = zer(iiz);
+    else
+      czer = [];
+    endif
+    if (! isempty (pol))
+      cpol = pol(iip);
+    else
+      cpol = [];
+    endif
   endif
   
   if (isempty (iip) && isempty (iiz))
@@ -69,8 +88,8 @@
     wmin = -1;
     wmax = 3;
   else
-    wmin = floor(log10(min(abs([cpol,czer]))));
-    wmax = ceil(log10(max(abs([cpol,czer]))));
+    wmin = floor (log10 (min (abs ([cpol, czer]))));
+    wmax = ceil (log10 (max (abs ([cpol, czer]))));
   endif
 
   ## expand to show the entirety of the "interesting" portion of the plot
@@ -78,6 +97,8 @@
   wmax++;
 
   ## run digital frequency all the way to pi
-  if (DIGITAL) wmax = log10(pi/tsam); endif
+  if (DIGITAL)
+    wmax = log10 (pi/tsam);
+  endif
 
 endfunction

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

--- a/scripts/control/base/dlyap.m
+++ b/scripts/control/base/dlyap.m
@@ -90,6 +90,10 @@
 
 function x = dlyap (a, b)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   if ((n = issquare (a)) == 0)
     warning ("dlyap: a must be square");
   endif
@@ -129,19 +133,19 @@
       blksiz = 1;
     endif
 
-    Ajj = kron (s (j:j1, j:j1), s) - eye (blksiz*n);
+    Ajj = kron (s(j:j1,j:j1), s) - eye (blksiz*n);
 
-    rhs = reshape (b (:, j:j1), blksiz*n, 1);
+    rhs = reshape (b (:,j:j1), blksiz*n, 1);
 
     if (j1 < n)
-      rhs2 = s*(x (:, (j1+1):n) * s (j:j1, (j1+1):n)');
+      rhs2 = s*(x(:,(j1+1):n) * s(j:j1,(j1+1):n)');
       rhs = rhs + reshape (rhs2, blksiz*n, 1);
     endif
 
     v = - Ajj\rhs;
-    x (:, j) = v (1:n);
+    x(:,j) = v (1:n);
 
-    if(blksiz == 2)
+    if (blksiz == 2)
       x (:, j1) = v ((n+1):blksiz*n);
     endif
 

--- a/scripts/control/base/freqchkw.m
+++ b/scripts/control/base/freqchkw.m
@@ -29,6 +29,10 @@
 
 function USEW = freqchkw (w)
 
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   if (isempty (w))
     USEW = 0;
   elseif (! isvector (w))

--- a/scripts/control/base/gram.m
+++ b/scripts/control/base/gram.m
@@ -29,6 +29,10 @@
 
 function m = gram (a, b)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   ## Let lyap do the error checking...
 
   m = lyap (a, b*b');

--- a/scripts/control/base/place.m
+++ b/scripts/control/base/place.m
@@ -45,78 +45,81 @@
 
 function K = place (sys, P)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   ## check arguments
 
-  if(!isstruct(sys))
-    error("sys must be in system data structure format (see ss)");
+  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
-  if(!is_controllable(sys))
-    error("sys is not controllable.");
-  elseif( min(size(P)) != 1)
-    error("P must be a vector")
+  sys = sysupdate (sys, "ss");    # make sure it has state space form up to date
+  if (! is_controllable (sys))
+    error ("sys is not controllable");
+  elseif (min (size (P)) != 1)
+    error ("P must be a vector")
   else
-    P = reshape(P,length(P),1); # make P a column vector
+    P = P(:); # make P a column vector
   endif
   ## system must be purely continuous or discrete
-  is_digital(sys);
-  [n,nz,m,p] = sysdimensions(sys);
+  is_digital (sys);
+  [n, nz, m, p] = sysdimensions (sys);
   nx = n+nz;    # already checked that it's not a mixed system.
   if(m != 1)
-    error(["sys has ", num2str(m)," inputs; need only 1"]);
+    error ("sys has %d inputs; need only 1", m);
   endif
 
   ## takes the A and B matrix from the system representation
-  [A,B]=sys2ss(sys);
-  sp = length(P);
-  if(nx == 0)
-    error("place: A matrix is empty (0x0)");
-  elseif(nx != length(P))
-    error(["A=(",num2str(nx),"x",num2str(nx),", P has ", num2str(length(P)), ...
-	"entries."])
+  [A, B] = sys2ss (sys);
+  sp = length (P);
+  if (nx == 0)
+    error ("place: A matrix is empty (0x0)");
+  elseif (nx != length (P))
+    error ("A=(%dx%d), P has %d entries", nx, nx, length (P))
   endif
 
   ## arguments appear to be compatible; let's give it a try!
   ## The second step is the calculation of the characteristic polynomial ofA
-  PC=poly(A);
+  PC = poly (A);
 
   ## Third step: Calculate the transformation matrix T that transforms the state
   ## equation in the controllable canonical form.
 
   ## first we must calculate the controllability matrix M:
-  M=B;
-  AA=A;
+  M = B;
+  AA = A;
   for n = 2:nx
-    M(:,n)=AA*B;
-    AA=AA*A;
+    M(:,n) = AA*B;
+    AA = AA*A;
   endfor
 
   ## second, construct the matrix W
-  PCO=PC(nx:-1:1);
-  PC1=PCO;      # Matrix to shift and create W row by row
+  PCO = PC(nx:-1:1);
+  PC1 = PCO;      # Matrix to shift and create W row by row
 
   for n = 1:nx
     W(n,:) = PC1;
-    PC1=[PCO(n+1:nx),zeros(1,n)];
+    PC1 = [PCO(n+1:nx), zeros(1,n)];
   endfor
 
-  T=M*W;
+  T = M*W;
 
   ## finaly the matrix K is calculated
-  PD = poly(P); # The desired characteristic polynomial
+  PD = poly (P); # The desired characteristic polynomial
   PD = PD(nx+1:-1:2);
   PC = PC(nx+1:-1:2);
 
   K = (PD-PC)/T;
 
   ## Check if the eigenvalues of (A-BK) are the same specified in P
-  Pcalc = eig(A-B*K);
+  Pcalc = eig (A-B*K);
 
-  Pcalc = sortcom(Pcalc);
-  P = sortcom(P);
+  Pcalc = sortcom (Pcalc);
+  P = sortcom (P);
 
-  if(max( (abs(Pcalc)-abs(P))./abs(P) ) > 0.1)
-    disp("Place: Pole placed at more than 10% relative error from specified");
+  if (max ((abs(Pcalc)-abs(P))./abs(P) ) > 0.1)
+    warning ("place: Pole placed at more than 10% relative error from specified");
   endif
 
 endfunction

--- a/scripts/control/hinf/hinf_ctr.m
+++ b/scripts/control/hinf/hinf_ctr.m
@@ -57,6 +57,10 @@
 
 function K = hinf_ctr (dgs, F, H, Z, g)
 
+  if (nargin != 5)
+    print_usage ();
+  endif
+
   nw = dgs.nw;
   nu = dgs.nu;
   nz = dgs.nz;
@@ -78,7 +82,7 @@
 
   nout = nz + ny;
   nin = nw + nu;
-  nstates = size(A, 1);
+  nstates = size (A, 1);
 
   F11 = F(1:(nw-ny),:);
   F12 = F((nw-ny+1):nw,:);
@@ -94,33 +98,33 @@
   D1122 = D11((nz-nu+1):nz,(nw-ny+1):nw);
 
   ## D11ik may be the empty matrix, don't calculate with empty matrices
-  [nd1111,md1111] = size(D1111);
-  md1112 = length(D1112);
-  md1121 = length(D1121);
+  [nd1111, md1111] = size (D1111);
+  md1112 = length (D1112);
+  md1121 = length (D1121);
 
-  if ((nd1111 == 0) || (md1112 == 0))
+  if (nd1111 == 0) || md1112 == 0)
     d11hat = -D1122;
   else
-    xx = inv(g*g*eye(nz-nu) - D1111*D1111');
+    xx = inv (g*g*eye(nz-nu) - D1111*D1111');
     d11hat = -D1121*D1111'*xx*D1112 - D1122;
   endif
   if (md1112 == 0)
-    d21hat = eye(ny);
+    d21hat = eye (ny);
   elseif (nd1111 == 0)
-    d21hat = chol(eye(ny) - D1112'*D1112/g/g);
+    d21hat = chol (eye(ny) - D1112'*D1112/g/g);
   else
-    xx = inv(g*g*eye(nz-nu) - D1111*D1111');
-    xx = eye(ny) - D1112'*xx*D1112;
-    d21hat = chol(xx);
+    xx = inv (g*g*eye(nz-nu) - D1111*D1111');
+    xx = eye (ny) - D1112'*xx*D1112;
+    d21hat = chol (xx);
   endif
   if (md1121 == 0)
-    d12hat = eye(nu);
+    d12hat = eye (nu);
   elseif (md1111 == 0)
-    d12hat = chol(eye(nu) - D1121*D1121'/g/g)';
+    d12hat = chol (eye(nu) - D1121*D1121'/g/g)';
   else
-    xx = inv(g*g*eye(nw-ny) - D1111'*D1111);
-    xx = eye(nu)-D1121*xx*D1121';
-    d12hat = chol(xx)';
+    xx = inv (g*g*eye(nw-ny) - D1111'*D1111);
+    xx = eye (nu)-D1121*xx*D1121';
+    d12hat = chol (xx)';
   endif
 
   b2hat = (B2+H12)*d12hat;
@@ -138,13 +142,13 @@
 
   ## non-zero D22 is a special case
   if (d22nz)
-    if (rank(eye(nu) + d11hat*D22) < nu)
+    if (rank (eye(nu) + d11hat*D22) < nu)
       error(" *** cannot compute controller for D22 non-zero.");
     endif
 
     d22new = [D22, zeros(ny,ny); zeros(nu,nu), 0*D22'];
-    xx = inv(eye(nu+ny) + d22new*dhat);
-    mhat = inv(eye(nu+ny) + dhat*d22new);
+    xx = inv (eye(nu+ny) + d22new*dhat);
+    mhat = inv (eye(nu+ny) + dhat*d22new);
     ahat = ahat - bhat*((eye(nu+ny)-xx)/dhat)*chat;
     bhat = bhat*xx;
     chat = mhat*chat;
@@ -152,6 +156,6 @@
 
   endif
 
-  K = ss(ahat,bhat(:,1:ny),chat(1:nu,:),dhat(1:nu,1:ny));
+  K = ss (ahat, bhat(:,1:ny), chat(1:nu,:), dhat(1:nu,1:ny));
 
 endfunction

--- a/scripts/control/hinf/hinfsyn_chk.m
+++ b/scripts/control/hinf/hinfsyn_chk.m
@@ -80,6 +80,10 @@
 
 function [retval, Pc, Pf] = hinfsyn_chk (A, B1, B2, C1, C2, D12, D21, g, ptol)
 
+  if (nargin != 9)
+    print_usage ();
+  endif
+
   Pc = Pf = [];
 
   ## Construct the two Hamiltonians
@@ -88,39 +92,39 @@
   Hf = [ A' , g2*C1'*C1 - C2'*C2; -B1*B1' , -A];
 
   ## check if Hc, Hf are in dom(Ric)
-  Hcminval = min(abs(real(eig(Hc))));
-  Hfminval = min(abs(real(eig(Hf))));
-  if(Hcminval < ptol);
-    disp("hinfsyn_chk: Hc is not in dom(Ric)");
+  Hcminval = min (abs (real (eig (Hc))));
+  Hfminval = min (abs (real (eig (Hf))));
+  if (Hcminval < ptol);
+    warning ("hinfsyn_chk: Hc is not in dom(Ric)");
     retval = 0;
     return
   endif
   if(Hfminval < ptol)
-    disp("hinfsyn_chk: Hf is not in dom(Ric)");
+    warning ("hinfsyn_chk: Hf is not in dom(Ric)");
     retval = 0;
     return
   endif
 
   ## Solve ARE's
-  Pc = are(A, B2*B2'-g2*B1*B1',C1'*C1);
-  Pf = are(A',C2'*C2-g2*C1'*C1,B1*B1');
+  Pc = are (A, B2*B2'-g2*B1*B1', C1'*C1);
+  Pf = are (A', C2'*C2-g2*C1'*C1, B1*B1');
 
-  Pceig = eig(Pc);
-  Pfeig = eig(Pf);
-  Pcfeig = eig(Pc*Pf);
+  Pceig = eig (Pc);
+  Pfeig = eig (Pf);
+  Pcfeig = eig (Pc*Pf);
 
-  if(min(Pceig) < -ptol)
-    disp("hinfsyn_chk: Pc is not >= 0");
+  if (min (Pceig) < -ptol)
+    warning ("hinfsyn_chk: Pc is not >= 0");
     retval = 0;
     return
   endif
-  if(min(Pfeig) < -ptol)
-    disp("hinfsyn_chk: Pf is not >= 0");
+  if (min (Pfeig) < -ptol)
+    warning ("hinfsyn_chk: Pf is not >= 0");
     retval = 0;
     return
   endif
-  if(max(abs(Pcfeig)) >= g*g)
-    disp("hinfsyn_chk: rho(Pf*Pc) is not < g^2");
+  if (max (abs (Pcfeig)) >= g*g)
+    warning ("hinfsyn_chk: rho(Pf*Pc) is not < g^2");
     retval = 0;
     return
   endif

--- a/scripts/control/hinf/hinfsyn_ric.m
+++ b/scripts/control/hinf/hinfsyn_ric.m
@@ -47,39 +47,44 @@
 
 function [Xinf, x_ha_err] = hinfsyn_ric (A, BB, C1, d1dot, R, ptol)
 
+  if (nargin != 6)
+    print_usage ();
+  endif
+
   x_ha_err = 0;        # assume success
   Xinf = [];                 # default return value
-  n = issquare(A);
-  nw = issquare(R);
-  if(rank(R) != nw)    x_ha_err = 6;
+  n = issquare (A);
+  nw = issquare (R);
+  if (rank (R) != nw)
+    x_ha_err = 6;
   else                 # build hamiltonian Ha for X_inf
     xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1, BB'];
     Ha = [A, 0*A; -C1'*C1, -A'] - xx;
     x_ha_err = 0;
-    [d, Ha] = balance(Ha);
-    [u, s] = schur(Ha, "A");
-    rev = real(eig(s));
+    [d, Ha] = balance (Ha);
+    [u, s] = schur (Ha, "A");
+    rev = real (eig (s));
 
-    if (any(abs(rev) <= ptol))  # eigenvalues near the imaginary axis
+    if (any (abs (rev) <= ptol))  # eigenvalues near the imaginary axis
       x_ha_err = 1;
-    elseif (sum(rev > 0) != sum(rev < 0))
+    elseif (sum (rev > 0) != sum (rev < 0))
       ## unequal number of positive and negative eigenvalues
       x_ha_err = 2;
     else
       ## compute positive Riccati equation solution
       u = d * u;
       Xinf = u(n+1:2*n,1:n) / u(1:n,1:n);
-      if (!all(all(finite(Xinf))))
+      if (! all (all (finite (Xinf))))
         x_ha_err = 3;
-      elseif (norm(Xinf-Xinf') >= 10*ptol)
+      elseif (norm (Xinf-Xinf') >= 10*ptol)
         ## solution not symmetric
         x_ha_err = 4;
       else
         ## positive semidefinite?
         ## force symmetry (faster, avoids some convergence problems)
         Xinf = (Xinf + Xinf')/2;
-        rev = eig(Xinf);
-        if (any(rev <= -ptol))
+        rev = eig (Xinf);
+        if (any (rev <= -ptol))
           x_ha_err = 5;
         endif
       endif

--- a/scripts/control/system/is_sample.m
+++ b/scripts/control/system/is_sample.m
@@ -28,6 +28,10 @@
 
 function out = is_sample (Ts)
 
+  if (nargin != 1)
+    print_usage (); 
+  endif
+
   out = (isscalar (Ts) && (Ts == abs (Ts)) && (Ts != 0));
 
 endfunction

--- a/scripts/control/system/is_signal_list.m
+++ b/scripts/control/system/is_signal_list.m
@@ -24,9 +24,13 @@
 
 function flg = is_signal_list (mylist)
 
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   flg = iscell (mylist);
-  if(flg)
-    flg = (rows(mylist) == 1 | columns(mylist) == 1);
+  if (flg)
+    flg = (rows (mylist) == 1 || columns (mylist) == 1);
   end
   if (flg)
     for ii = 1:length (mylist)

--- a/scripts/control/system/ss2tf.m
+++ b/scripts/control/system/ss2tf.m
@@ -59,16 +59,20 @@
 
 function [num, den] = ss2tf (a, b, c, d)
 
+  if (nargin != 4)
+    print_usage ();
+  endif
+
   ## Check args
-  [n,m,p] = abcddim(a,b,c,d);
+  [n, m, p] = abcddim (a, b, c, d);
   if (n == -1)
     num = [];
     den = [];
     error("ss2tf: Non compatible matrix arguments");
-  elseif ( (m != 1) | (p != 1))
+  elseif (m != 1 || p != 1)
     num = [];
     den = [];
-    error(["ss2tf: not SISO system: m=",num2str(m)," p=",num2str(p)]);
+    error ("ss2tf: not SISO system: m=%d, p=%d", m, p);
   endif
 
   if(n == 0)
@@ -77,21 +81,21 @@
     den = 1;
   else
     ## First, get the denominator coefficients
-    den = poly(a);
+    den = poly (a);
 
     ## Get the zeros of the system
-    [zz,g] = tzero(a,b,c,d);
+    [zz, g] = tzero (a, b, c, d);
 
     ## Form the Numerator (and include the gain)
-    if (!isempty(zz))
-      num = g * poly(zz);
+    if (! isempty (zz))
+      num = g * poly (zz);
     else
       num = g;
     endif
 
     ## the coefficients must be real
-    den = real(den);
-    num = real(num);
+    den = real (den);
+    num = real (num);
   endif
+
 endfunction
-

--- a/scripts/control/system/sys2fir.m
+++ b/scripts/control/system/sys2fir.m
@@ -30,6 +30,10 @@
 
 function [c, tsam, inname, outname] = sys2fir (sys)
 
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   ## let sys2tf do most of the work
 
   [num, den, tsam, inname, outname] = sys2tf (sys);

--- a/scripts/control/system/sysgettsam.m
+++ b/scripts/control/system/sysgettsam.m
@@ -24,8 +24,12 @@
 
 function T = sysgettsam (sys)
 
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   if (! isstruct (sys))
-    print_usage ();
+    error ("sysgettsam: expecting argument to be system structure");
   endif
 
   T = sys.tsam;

--- a/scripts/control/system/sysgettype.m
+++ b/scripts/control/system/sysgettype.m
@@ -39,11 +39,15 @@
 
 function systype = sysgettype (sys)
 
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   if (! isstruct (sys))
     error ("sysgettype: input sys is not a structure");
   endif
 
   typestr = {"tf", "zp", "ss"};
-  systype = typestr{ sys.sys(1) + 1};
+  systype = typestr{sys.sys(1) + 1};
 
 endfunction

--- a/scripts/control/system/sysreorder.m
+++ b/scripts/control/system/sysreorder.m
@@ -44,6 +44,10 @@
 
 function pv = sysreorder (vlen, list)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   ## disp('sysreorder: entry')
 
   pv = 1:vlen;

--- a/scripts/control/system/tf2sys.m
+++ b/scripts/control/system/tf2sys.m
@@ -63,7 +63,7 @@
 ## name changed to tf Feb 2004
 
 function outsys = tf2sys (varargin)
-  
+
   warning("tf2sys is deprecated.  Use tf() instead.");
   outsys = tf(varargin{:});
 

--- a/scripts/control/system/zp2tf.m
+++ b/scripts/control/system/zp2tf.m
@@ -37,18 +37,22 @@
 
 function [num, den] = zp2tf (zer, pol, k)
 
+  if (nargin != 3)
+    print_usage ();
+  endif  
+
   ## Find out whether data was entered as a row or a column vector and
   ## convert to a column vector if necessary.
 
-  [rp,cp] = size(pol);
-  [rz,cz] = size(zer);
+  [rp, cp] = size (pol);
+  [rz, cz] = size (zer);
 
-  if(!(isvector(zer) | isempty(zer)) )
-    error(sprintf("zer(%dx%d) must be a vector",rz,cz));
-  elseif(!(isvector(pol) | isempty(pol)) )
-    error(sprintf("pol(%dx%d) must be a vector",rp,cp));
-  elseif(length(zer) > length(pol))
-    error(sprintf("zer(%dx%d) longer than pol(%dx%d)",rz,cz,rp,cp));
+  if (! (isvector (zer) || isempty (zer)))
+    error ("zer(%dx%d) must be a vector", rz, cz);
+  elseif(! (isvector (pol) || isempty (pol)))
+    error ("pol(%dx%d) must be a vector", rp, cp);
+  elseif (length (zer) > length (pol))
+    error ("zer(%dx%d) longer than pol(%dx%d)", rz, cz, rp, cp);
   endif
 
   ## initialize converted polynomials
@@ -58,18 +62,24 @@
   ## call __zp2ssg2__ if there are complex conjugate pairs left, otherwise
   ## construct real zeros one by one.  Repeat for poles.
 
-  while(!isempty(zer))
-    if( max(abs(imag(zer))) )     [poly, zer] = __zp2ssg2__ (zer);
-    else                          poly = [1, -zer(1)];
-                                  zer = zer(2:length(zer));      endif
-    num = conv(num,poly);
+  while(! isempty (zer))
+    if (max (abs (imag (zer))))
+      [poly, zer] = __zp2ssg2__ (zer);
+    else
+      poly = [1, -zer(1)];
+      zer = zer(2:length(zer));
+    endif
+    num = conv (num, poly);
   endwhile
 
-  while(!isempty(pol))
-    if( max(abs(imag(pol))) )     [poly, pol] = __zp2ssg2__ (pol);
-    else                          poly = [1, -pol(1)];
-                                  pol = pol(2:length(pol));      endif
-    den = conv(den,poly);
+  while (! isempty (pol))
+    if (max (abs (imag (pol))))
+      [poly, pol] = __zp2ssg2__ (pol);
+    else
+      poly = [1, -pol(1)];
+      pol = pol(2:length(pol));
+    endif
+    den = conv (den, poly);
   endwhile
 
 endfunction

--- a/scripts/control/util/axis2dlim.m
+++ b/scripts/control/util/axis2dlim.m
@@ -39,13 +39,13 @@
 
 function axvec = axis2dlim (axdata)
 
-  if(isempty(axdata))
+  if (nargin < 1 || isempty (axdata))
     axdata = 0;
   endif
 
   ## compute axis limits
-  minv = min(axdata);
-  maxv = max(axdata);
+  minv = min (axdata);
+  maxv = max (axdata);
   delv = (maxv-minv)/2;      # breadth of the plot
   midv = (minv + maxv)/2;    # midpoint of the plot
   axmid = [midv(1), midv(1), midv(2), midv(2)];

--- a/scripts/control/util/swap.m
+++ b/scripts/control/util/swap.m
@@ -31,6 +31,10 @@
 
 function [a1, b1] = swap (a, b)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   a1 = b;
   b1 = a;
 

--- a/scripts/control/util/zgfmul.m
+++ b/scripts/control/util/zgfmul.m
@@ -33,6 +33,10 @@
 
 function y = zgfmul (a, b, c, d, x)
 
+  if (nargin != 5)
+    print_usage ();
+  endif 
+
   [n,m] = size(b);
   [p,m1] = size(c);
   nm = n+m;

--- a/scripts/control/util/zgfslv.m
+++ b/scripts/control/util/zgfslv.m
@@ -27,6 +27,10 @@
 
 function x = zgfslv (n, m, p, b)
 
+  if (nargin != 4)
+    print_usage ();
+  endif
+
   nmp = n+m+p;
   gam1 = (2*n)+m+p;    gam2 = n+p;     gam3 = n+m;
 

--- a/scripts/control/util/zginit.m
+++ b/scripts/control/util/zginit.m
@@ -34,6 +34,10 @@
 
 function zz = zginit (a, b, c, d)
 
+  if (nargin != 4)
+    print_usage ();
+  endif
+
   [nn,mm] = size(b);
   [pp,mm] = size(d);
 

--- a/scripts/control/util/zgreduce.m
+++ b/scripts/control/util/zgreduce.m
@@ -25,6 +25,10 @@
 
 function retsys = zgreduce (Asys, meps)
 
+  if (nargin != 2)
+    print_usage ();
+  endif  
+
   ## SYS_INTERNAL accesses members of system data structure
 
   is_digital(Asys);             # make sure it's pure digital/continuous

--- a/scripts/control/util/zgrownorm.m
+++ b/scripts/control/util/zgrownorm.m
@@ -26,6 +26,10 @@
 
 function [sig, tau] = zgrownorm (mat, meps)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   rownorm = [];
   for ii = 1:rows (mat)
     rownorm(ii) = norm (mat(ii,:));

--- a/scripts/control/util/zgscal.m
+++ b/scripts/control/util/zgscal.m
@@ -34,6 +34,10 @@
 
 function x = zgscal (a, b, c, d, z, n, m, p)
 
+  if (nargin != 8)
+    print_usage ();
+  endif
+
   ## initialize parameters:
   ## Givens rotations, diagonalized 2x2 block of F, gcg vector initialization
 

--- a/scripts/control/util/zgsgiv.m
+++ b/scripts/control/util/zgsgiv.m
@@ -30,6 +30,10 @@
 
 function [a, b] = zgsgiv (c, s, a, b)
 
+  if (nargin != 4)
+    print_usage ();
+  endif
+
   t1 = c*a + s*b;
   t2 = -s*a + c*b;
   a = t1;

--- a/scripts/control/util/zgshsr.m
+++ b/scripts/control/util/zgshsr.m
@@ -38,8 +38,12 @@
 
 function x = zgshsr (y)
 
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   if (! isvector (y))
-    error (sprintf ("y(%dx%d) must be a vector", rows (y), columns (y)));
+    error ("y(%dx%d) must be a vector", rows (y), columns (y));
   endif
   x = vec (y);
   m = length (x);

--- a/scripts/general/isa.m
+++ b/scripts/general/isa.m
@@ -25,5 +25,11 @@
 ## Adapted-by: jwe
 
 function retval = isa (x, cname)
+
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   retval = strcmp (class (x), cname);
+
 endfunction

--- a/scripts/geometry/inpolygon.m
+++ b/scripts/geometry/inpolygon.m
@@ -38,15 +38,19 @@
 
 function [IN, ON] = inpolygon (X, Y, xv, yv)
 
-  if ( !(isreal(X) && isreal(Y) && ismatrix(Y) && ...
-	 ismatrix(Y) && size_equal(X,Y)) )
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! (isreal (X) && isreal (Y) && ismatrix (Y) && ismatrix (Y)
+	 && size_equal (X, Y)))
     error ("inpolygon: first two arguments must be real matrices of same size");
-  elseif ( !(isreal(xv) && isreal(yv) && isvector(xv) && ...
-	     isvector(yv) && size_equal(xv,yv)) )
+  elseif (! (isreal (xv) && isreal (yv) && isvector (xv) && isvector (yv)
+	     && size_equal (xv, yv)))
     error ("inpolygon: last two arguments must be real vectors of same size");
   endif
 
-  npol = length(xv);
+  npol = length (xv);
   do_boundary = (nargout >= 2);
   
   IN = zeros (size(X), "logical");
@@ -76,6 +80,7 @@
     endif
     j = i;
   endfor
+
 endfunction
 
 %!demo

--- a/scripts/linear-algebra/housh.m
+++ b/scripts/linear-algebra/housh.m
@@ -35,8 +35,12 @@
 
 function [housv, beta, zer] = housh (x, j, z)
 
-  ## check for valid inputs
-  if (!isvector (x) && ! isscalar (x))
+  if (nargin != 3)
+    print_usage ();
+  endif
+
+  ## Check for valid inputs.
+  if (! isvector (x) && ! isscalar (x))
     error ("housh: first input must be a vector")
   elseif (! isscalar(j))
     error ("housh: second argment must be an integer scalar")

--- a/scripts/miscellaneous/compare_versions.m
+++ b/scripts/miscellaneous/compare_versions.m
@@ -70,6 +70,10 @@
 
 function out = compare_versions (v1, v2, operator)
 
+  if (nargin != 3)
+    print_usage ();
+  endif
+
   ## Make sure that the version numbers are valid.
   if (! (ischar (v1) && ischar (v2)))
     error ("compare_versions: both version numbers must be strings");
@@ -184,6 +188,7 @@
 
   ## Reverse the output if not is given.
   out = xor (not_op, out);
+
 endfunction
 
 ## tests

--- a/scripts/miscellaneous/inputname.m
+++ b/scripts/miscellaneous/inputname.m
@@ -25,7 +25,13 @@
 ## @end deftypefn
 
 function s = inputname (n)
+
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   s = evalin ("caller", sprintf ("deblank (argn(%d,:));", n));
+
 endfunction
 
 ## Warning: heap big magic in the following tests!!!

--- a/scripts/miscellaneous/run.m
+++ b/scripts/miscellaneous/run.m
@@ -29,6 +29,11 @@
 ## PKG_ADD: mark_as_command run
 
 function run (s)
+
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   [d, f, ext] = fileparts (s);
   if (! isempty (d))
     if (exist (d, "dir"))

--- a/scripts/quaternion/qconj.m
+++ b/scripts/quaternion/qconj.m
@@ -32,6 +32,10 @@
 
 function retval = qconj (q)
 
+  if (nargin != 1 )
+    print_usage ();
+  endif
+
   [a, b, c, d] = quaternion (q);
 
   retval = quaternion (-a, -b, -c, d);

--- a/scripts/quaternion/qcoordinate_plot.m
+++ b/scripts/quaternion/qcoordinate_plot.m
@@ -38,6 +38,10 @@
 
 function qcoordinate_plot (qf, qb, qv)
 
+  if (nargin != 3 )
+    print_usage ();
+  endif
+
   degrees = pi / 180;
   d180 = 180 * degrees;
 

--- a/scripts/quaternion/qderiv.m
+++ b/scripts/quaternion/qderiv.m
@@ -43,6 +43,10 @@
 
 function Dmat = qderivmat (Omega)
 
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   Omega = vec (Omega);
 
   if (length (Omega) != 3)

--- a/scripts/quaternion/qderivmat.m
+++ b/scripts/quaternion/qderivmat.m
@@ -43,6 +43,10 @@
 
 function Dmat = qderivmat (Omega)
 
+  if (nargin != 1 )
+    print_usage ();
+  endif
+
   Omega = vec (Omega);
 
   if (length (Omega) != 3)

--- a/scripts/quaternion/qinv.m
+++ b/scripts/quaternion/qinv.m
@@ -32,6 +32,10 @@
 
 function retval = qinv (q)
 
+  if (nargin !=  1)
+    print_usage ();
+  endif
+
   if (norm (q) != 0)
     retval = qconj (q) / sum (q .* q);
   else

--- a/scripts/quaternion/qmult.m
+++ b/scripts/quaternion/qmult.m
@@ -40,7 +40,11 @@
 ## Adapted-By: jwe
 
 function retval = qmult (a, b)
-  
+
+  if (nargin != 2 )
+    print_usage ();
+  endif
+
   [a1, b1, c1, d1] = quaternion (a);
   [a2, b2, c2, d2] = quaternion (b);
   

--- a/scripts/quaternion/qtrans.m
+++ b/scripts/quaternion/qtrans.m
@@ -28,6 +28,10 @@
 
 function v = qtrans (v, q)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   if (! isvector (v) || length (v) != 4)
     error ("qtrans: v(%d,%d) must be a quaternion", rows (v), columns (v));
   elseif (! isvector (q) || length (q) != 4)

--- a/scripts/quaternion/qtransvmat.m
+++ b/scripts/quaternion/qtransvmat.m
@@ -29,6 +29,10 @@
 
 function Aib = qtransvmat (qib)
 
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   if (! isvector(qib) || length (qib) != 4)
     error ("qtransvmat: q(%d,%d) must be a quaternion", rows (qib), \
 	   columns (qib));

--- a/scripts/signal/fractdiff.m
+++ b/scripts/signal/fractdiff.m
@@ -28,6 +28,10 @@
 
 function retval = fractdiff (x, d)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   N = 100;
 
   if (! isvector (x))

--- a/scripts/signal/freqz_plot.m
+++ b/scripts/signal/freqz_plot.m
@@ -25,38 +25,42 @@
 
 function freqz_plot (w, h)
 
-    n = length (w);
+  if (nargin != 2)
+    print_usage ();
+  endif
 
-    ## ## exclude zero-frequency
-    ## h = h (2 : length (h));
-    ## w = w (2 : length (w));
-    ## n = n-1;
+  n = length (w);
 
-    mag = 20 * log10 (abs (h));
-    phase = unwrap (arg (h));
-    maxmag = max (mag);
+  ## ## exclude zero-frequency
+  ## h = h (2 : length (h));
+  ## w = w (2 : length (w));
+  ## n = n-1;
 
-    subplot (3, 1, 1);
-    plot (w, mag);
-    grid ("on");
-    legend("Pass band (dB)");
-    axis ([w(1), w(n), maxmag-3, maxmag], "labely");
+  mag = 20 * log10 (abs (h));
+  phase = unwrap (arg (h));
+  maxmag = max (mag);
 
-    subplot (3, 1, 2);
-    plot (w, mag);
-    grid ("on");
-    legend ("Stop band (dB)");
-    if (maxmag - min (mag) > 100)
-      axis ([w(1), w(n), maxmag-100, maxmag], "labely");
-    else
-      axis ("autoy", "labely");
-    endif
+  subplot (3, 1, 1);
+  plot (w, mag);
+  grid ("on");
+  legend("Pass band (dB)");
+  axis ([w(1), w(n), maxmag-3, maxmag], "labely");
 
-    subplot (3, 1, 3);
-    plot (w, phase*360/(2*pi));
-    grid ("on");
-    legend ("Phase (degrees)");
-    xlabel ("Frequency");
-    axis ([w(1), w(n)], "autoy", "label");
+  subplot (3, 1, 2);
+  plot (w, mag);
+  grid ("on");
+  legend ("Stop band (dB)");
+  if (maxmag - min (mag) > 100)
+    axis ([w(1), w(n), maxmag-100, maxmag], "labely");
+  else
+    axis ("autoy", "labely");
+  endif
+
+  subplot (3, 1, 3);
+  plot (w, phase*360/(2*pi));
+  grid ("on");
+  legend ("Phase (degrees)");
+  xlabel ("Frequency");
+  axis ([w(1), w(n)], "autoy", "label");
 
 endfunction

--- a/scripts/signal/periodogram.m
+++ b/scripts/signal/periodogram.m
@@ -28,6 +28,10 @@
 
 function retval = periodogram (x)
 
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   [r, c] = size(x);
 
   if (r == 1)

--- a/scripts/signal/rectangle_lw.m
+++ b/scripts/signal/rectangle_lw.m
@@ -28,6 +28,10 @@
 
 function retval = rectangle_lw (n, b)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   retval = zeros (n, 1);
   t = floor (1 / b);
 

--- a/scripts/signal/rectangle_sw.m
+++ b/scripts/signal/rectangle_sw.m
@@ -28,6 +28,10 @@
 
 function retval = rectangle_sw (n, b)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   retval = zeros (n, 1);
   retval(1) = 2 / b + 1;
 

--- a/scripts/signal/sinc.m
+++ b/scripts/signal/sinc.m
@@ -34,6 +34,10 @@
 
 function result = sinc (x)
 
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   result = ones (size (x));
 
   i = (x != 0);

--- a/scripts/signal/triangle_lw.m
+++ b/scripts/signal/triangle_lw.m
@@ -28,6 +28,10 @@
 
 function retval = triangle_lw (n, b)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   retval = 1 - (0 : n-1)' * b;
   retval = max ([retval'; (zeros (1, n))])';
 

--- a/scripts/signal/triangle_sw.m
+++ b/scripts/signal/triangle_sw.m
@@ -28,6 +28,10 @@
 
 function retval = triangle_sw (n, b)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   retval = zeros(n,1);
   retval(1) = 1 / b;
 

--- a/scripts/signal/yulewalker.m
+++ b/scripts/signal/yulewalker.m
@@ -31,6 +31,10 @@
 
 function [a, v] = yulewalker (c)
 
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   p = length (c) - 1;
 
   if (columns (c) > 1)

--- a/scripts/sparse/colperm.m
+++ b/scripts/sparse/colperm.m
@@ -26,6 +26,11 @@
 ## @end deftypefn
 
 function p = colperm (s)
+
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   [i, j] = spfind (s);
   idx = find (diff ([j; Inf]) != 0);
   [dummy, p] = sort (idx - [0; idx(1:(end-1))]);

--- a/scripts/sparse/etreeplot.m
+++ b/scripts/sparse/etreeplot.m
@@ -27,5 +27,10 @@
 ## @end deftypefn
 
 function etreeplot (s, varargin)
+
+  if (nargin < 1)
+    print_usage ();
+  endif
+
   treeplot (etree (s+s'), varargin{:});
 endfunction

--- a/scripts/sparse/nonzeros.m
+++ b/scripts/sparse/nonzeros.m
@@ -22,6 +22,11 @@
 ## @end deftypefn
 
 function t = nonzeros (s)
+
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   if (issparse (s))
     [i, j, t] = spfind (s);
   else

--- a/scripts/sparse/spalloc.m
+++ b/scripts/sparse/spalloc.m
@@ -41,5 +41,10 @@
 ## @end deftypefn
 
 function s = spalloc (r, c, nz)
+
+  if (nargin < 2)
+    print_usage ();
+  endif
+
   s = sparse (r, c);
 endfunction

--- a/scripts/sparse/spones.m
+++ b/scripts/sparse/spones.m
@@ -23,6 +23,11 @@
 ## @end deftypefn
 
 function s = spones (s)
+
+  if (nargin != 1)
+    print_usage ();
+  endif
+
   if (issparse (s))
     [i, j, v, m, n] = spfind (s);
   else

--- a/scripts/sparse/spy.m
+++ b/scripts/sparse/spy.m
@@ -30,6 +30,10 @@
 
 function spy (S, varargin) 
 
+  if (nargin < 1)
+    print_usage ();
+  endif
+
   markersize = NaN;
   if (numel (i) < 1000)
     LineSpec = "*";

--- a/scripts/specfun/isprime.m
+++ b/scripts/specfun/isprime.m
@@ -35,6 +35,11 @@
 ## @end deftypefn
 
 function t = isprime (n)
+
+  if (nargin < 1)
+    print_usage ();
+  endif
+
   if (! isscalar (n))
     nel = numel (n);
     t = n;

--- a/scripts/statistics/distributions/empirical_cdf.m
+++ b/scripts/statistics/distributions/empirical_cdf.m
@@ -28,6 +28,10 @@
 
 function cdf = empirical_cdf (x, data)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   if (! isvector (data))
     error ("empirical_cdf: data must be a vector");
   endif

--- a/scripts/statistics/distributions/empirical_inv.m
+++ b/scripts/statistics/distributions/empirical_inv.m
@@ -28,6 +28,10 @@
 
 function inv = empirical_inv (x, data)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   if (! isvector (data))
     error ("empirical_inv: data must be a vector");
   endif

--- a/scripts/statistics/distributions/empirical_pdf.m
+++ b/scripts/statistics/distributions/empirical_pdf.m
@@ -28,6 +28,10 @@
 
 function pdf = empirical_pdf (x, data)
 
+  if (nargin != 2)
+    print_usage ();
+  endif
+
   if (! isvector (data))
     error ("empirical_pdf: data must be a vector");
   endif

--- a/scripts/statistics/models/logistic_regression_derivatives.m
+++ b/scripts/statistics/models/logistic_regression_derivatives.m
@@ -29,6 +29,10 @@
 
 function [dl, d2l] = logistic_regression_derivatives (x, z, z1, g, g1, p)
 
+  if (nargin != 6)
+    print_usage ();
+  endif
+
   ## first derivative
   v = g .* (1 - g) ./ p; v1 = g1 .* (1 - g1) ./ p;
   dlogp = [(dmult (v, z) - dmult (v1, z1)), (dmult (v - v1, x))];
@@ -39,4 +43,4 @@
   d2l = [z, x]' * dmult (w, [z, x]) - [z1, x]' * dmult (w1, [z1, x]) ...
       - dlogp' * dlogp;
 
-endfunction
\ No newline at end of file
+endfunction

--- a/scripts/statistics/models/logistic_regression_likelihood.m
+++ b/scripts/statistics/models/logistic_regression_likelihood.m
@@ -29,6 +29,10 @@
 
 function [g, g1, p, dev] = logistic_regression_likelihood (y, x, beta, z, z1)
 
+  if (nargin != 5)
+    print_usage ();
+  endif
+
   e = exp ([z, x] * beta); e1 = exp ([z1, x] * beta);
   g = e ./ (1 + e); g1 = e1 ./ (1 + e1);
   g = max (y == max (y), g); g1 = min (y > min(y), g1);