1 files changed, 56 insertions, 0 deletions
diff --git a/modules/cacsd/macros/riccati.sci b/modules/cacsd/macros/riccati.sci
new file mode 100755
index 000000000..89ce0ecb8
--- /dev/null
+++ b/modules/cacsd/macros/riccati.sci
@@ -0,0 +1,56 @@
+// 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 [x1,x2]=riccati(a,b,c,dom,typ)
+    //[x1,[x2]]=riccati(a,b,c,dom,[typ]) is a Riccati solver
+    // x=x1/x2 solves:
+    // a'*x+x*a-x*b*x+c=0 (continuous time case)
+    // a'*x*a-(a'*x*b1/(b2+b1'*x*b1))*(b1'*x*a)+c-x with b=b1/b2*b1'
+    // (discrete time case)
+    // If called with LHS=1 (one output argument) riccati returns x.
+    //
+    // -- a,b,c real matrices nxn, b and c symmetric.
+    // -- dom = 'c' or 'd' for the time domain (continuous or discrete)
+    // -- typ = 'eigen' --->block diagonalization
+    //    typ = 'schur' --->schur method
+    // See also ric_desc
+    //!
+
+    [lhs,rhs]=argn(0),
+    if rhs==4 then typ="eigen",end,
+    ham=[a -b;-c -a'],
+    [n,n]=size(a),
+    if part(dom,1)=="c" then
+        select  typ,
+        case "schur" then
+            [s,d]=schur(ham,"c"),
+            if d<>n then
+                error(msprintf(gettext("%s: Wrong dimension (%d) of stable subspace: %d expected.\n"),"riccati",d, n))
+            end
+            s=s(:,1:n),
+        case "eigen" then
+            [hb,u1]=bdiag(ham),
+            [u2,d]=schur(hb,"c"),
+            u=u1*u2,
+            if d<>n then
+                error(msprintf(gettext("%s: Wrong dimension (%d) of stable subspace: %d expected.\n"),"riccati",d, n))
+            end
+            s=u(:,1:n),
+        end,
+    else
+        aa=[eye(n,n) b;0*ones(n,n) a'],bb=[a  0*ones(n,n);-c eye(n,n)],
+        [bs,as,s,n1]=schur(bb,aa,"d");
+        if n1<>n then
+            error(msprintf(gettext("%s: Wrong dimension (%d) of stable subspace: %d expected.\n"),"riccati",n1, n))
+        end
+        s=s(:,1:n1);
+    end,
+    x1=s(n+1:2*n,:),x2=s(1:n,:),
+    if lhs==1 then x1=x1/x2,end
+endfunction