1 files changed, 75 insertions, 0 deletions
diff --git a/modules/cacsd/macros/svplot.sci b/modules/cacsd/macros/svplot.sci
new file mode 100755
index 000000000..1a2cc310d
--- /dev/null
+++ b/modules/cacsd/macros/svplot.sci
@@ -0,0 +1,75 @@
+// 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.1-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 = max(size(w)); nsv = min(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 = max(size(w)); nsv = min(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 = max(size(w)); nsv = min(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 = max(size(w)); nsv = min(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