// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA - 
// 
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at    
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

function svm = svplot(Sl,w)
//svplot singular-value sigma-plot.
// svm = svplot(sl,w) computes for the linear dynamical system
// sl, the singular values of its transfer function matrix:
//                              -1
//             g(jw) = c(jw*i-a)  b+d
//
//           or
//                                        -1
//             g(exp(jw)) = c(exp(jw)*i-a)  b+d
//
// evaluated over the frequency range specified by w.
// sl is a sylin list (see syslin) representing the system
// [a,b,c,d] in state-space form.
// the i-th column of the output matrix svm contains the singular
// values of g(exp(jw)) for the i-th frequency value.
// svm = svplot(sl) is equivalent to
// svm = svplot(sl,logspace(-3,3))  (continuous) or
// svm = svplot(sl,logspace(-3,pi)) (discrete).
//!

  [a,b,c,d]=abcd(Sl);
  // Reduce a to Hessenberg form
  [q,a] = hess(a); b = q'*b; c = c*q;
  // Compute the singular values of the frequency response
  select Sl.dt
  case []
    warning(msprintf(gettext("%s: Input argument %d is assumed continuous time.\n"),"svplot",1)); 
    if argn(2) == 1
      w = logspace(-3,3);
    end
    nf = maxi(size(w)); nsv = mini(size(d)); j = sqrt(-1);
    svm(nsv,nf) = 0;
    for i = 1:nf
      svm(:,i) = svd(c*((j*w(i)*eye()-a)\b)+d);
    end
  case 'c'
    if argn(2) == 1
      w = logspace(-3,3);
    end
    nf = maxi(size(w)); nsv = mini(size(d)); j = sqrt(-1);
    svm(nsv,nf) = 0;
    for i = 1:nf
      svm(:,i) = svd(c*((j*w(i)*eye()-a)\b)+d);
    end
  case 'd'
    if argn(2) == 1
      w = logspace(-3,%pi);
    end
    nf = maxi(size(w)); nsv = mini(size(d)); j = sqrt(-1);
    svm(nsv,nf) = 0;
    for i = 1:nf
      svm(:,i) = svd(c*((exp(j*w(i))*eye()-a)\b)+d);
    end
  else T=Sl('dt');
    if argn(2) == 1
      w = logspace(-3,%pi);
    end
    nf = maxi(size(w)); nsv = mini(size(d)); j = sqrt(-1);
    svm(nsv,nf) = 0;
    for i = 1:nf
      svm(:,i) = svd(c*((exp(j*w(i)*T)*eye()-a)\b)+d);
    end
    
  end
endfunction
