CMSCOPE changed

21 files changed, 3265 insertions, 0 deletions

diff --git a/modules/cacsd/demos/cacsd.dem.gateway.sce b/modules/cacsd/demos/cacsd.dem.gateway.sce
new file mode 100755
index 000000000..892470049
--- /dev/null
+++ b/modules/cacsd/demos/cacsd.dem.gateway.sce
@@ -0,0 +1,23 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+// Copyright (C) 2011 - DIGITEO - Allan CORNET
+//
+// This file is released under the 3-clause BSD license. See COPYING-BSD.
+
+function subdemolist = demo_gateway()
+    demopath = get_absolute_file_path("cacsd.dem.gateway.sce");
+    add_demo(gettext("CACSD"), demopath + "cacsd.dem.gateway.sce");
+
+    subdemolist = [_("LQG")                , "lqg/lqg.dem"
+    _("Mixed-sensitivity")  , "mixed.dem"
+    _("PID")                , "pid.dem"
+    _("Inverted pendulum")  , "pendulum/pendule.dem"
+    _("Flat systems")       , "flat/flat.dem.gateway.sce"
+    _("Tracking")           , "tracking/track.dem.sce"
+    _("Robust control")     , "robust/rob.dem"]
+
+    subdemolist(:,2) = demopath + subdemolist(:,2);
+endfunction
+
+subdemolist = demo_gateway();
+clear demo_gateway;
diff --git a/modules/cacsd/demos/cont.dem b/modules/cacsd/demos/cont.dem
new file mode 100755
index 000000000..9eb473784
--- /dev/null
+++ b/modules/cacsd/demos/cont.dem
@@ -0,0 +1,28 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function demo_cont()
+
+  deff('[]=demoexc(fil)','exec(''SCI/modules/cacsd/demos/''+fil)')
+  while %t do
+    n=x_choose(['LQG','Mixed-sensitivity','PID'],'Select a demo');
+    select n
+    case 0
+      break
+    case 1
+      demoexc('lqg/lqg.dem')
+    case 2
+      demoexc('mixed.dem')
+    case 3
+      demoexc('pid.dem')
+    end
+  end
+  
+endfunction
+
+demo_cont();
+clear demo_cont;
diff --git a/modules/cacsd/demos/flat/car.dem.sce b/modules/cacsd/demos/flat/car.dem.sce
new file mode 100755
index 000000000..9889eca5d
--- /dev/null
+++ b/modules/cacsd/demos/flat/car.dem.sce
@@ -0,0 +1,20 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is released under the 3-clause BSD license. See COPYING-BSD.
+
+function demo_car()
+
+    exec("SCI/modules/cacsd/demos/flat/car.sci",-1);
+
+    initial               =[3;3;0;0];
+    final                 =[0;0;0;0];
+    my_handle             = scf(100001);
+    clf(my_handle,"reset");
+    toolbar(my_handle.figure_id,"off");
+    state=car_solve(initial,final);
+    display_car_trajectory(state);
+endfunction
+
+demo_car();
+clear demo_car;
diff --git a/modules/cacsd/demos/flat/car.sci b/modules/cacsd/demos/flat/car.sci
new file mode 100755
index 000000000..6ccf8757c
--- /dev/null
+++ b/modules/cacsd/demos/flat/car.sci
@@ -0,0 +1,164 @@
+function state=car_solve(initial,final)
+    //
+    //  CAR PACKING VIA FLATNESS AND FRENET FORMULAS
+    //
+    //    explicit computation and visualisation of the motions.
+    //
+    //    February 1993
+    //
+    // ............................................................
+    // :         pierre ROUCHON  <rouchon@cas.ensmp.fr>           :
+    // : Centre Automatique et Systemes, Ecole des Mines de Paris :
+    // : 60, bd Saint Michel -- 75272 PARIS CEDEX 06, France      :
+    // :    Telephone: (1) 40 51 91 15 --- Fax: (1) 43 54 18 93   :
+    // :..........................................................:
+    //
+    // initial: initial position [x,y,theta,phi]
+    // final  :   final position [x,y,theta,phi]
+    //        theta  : the car angle
+    //        phi    : the front wheel angle
+
+    bigT = 1 ;//basic time interval for one  smooth motion (s)
+    bigL = 1 ;// car length (m)
+
+    // computation of  intermediate configuration
+    x0 = max(initial(1),final(2))   ....
+    + bigL*abs(tan(initial(3))) ...
+    + bigL*abs(tan(final(3))) ...
+    + bigL*(abs(initial(2)-final(2))/bigL)^(1/2) ;
+    y0 = (initial(2)+final(2))/2 ;
+    intermediate=[x0,y0,0,0]'
+
+    // first polynomial curve
+    state=[matrix(initial,1,-1);
+    car_polynomial_curve(initial,intermediate,"direct")]
+    //
+    // second polynomial curve
+    state = [ state;
+    matrix(intermediate,1,-1)
+    car_polynomial_curve(final,intermediate,"reverse")
+    matrix(final,1,-1)]
+endfunction
+
+function state=car_polynomial_curve(initial,final,orient)
+
+    nbpt = 40 ; //  sampling of motion
+    theta1 = initial(3) ; phi1 = initial(4) ;
+    da = initial(1)-final(1)
+
+    M = [da^3     da^4     da^5
+    3*da^2   4*da^3   5*da^4
+    6*da    12*da^2  20*da^3 ] ;
+
+    q = [initial(2)-final(2)
+    tan(theta1)
+    tan(phi1)/(bigL*(cos(theta1)^3))] ;
+
+    p = inv(M)*q ;
+    tau=(0:nbpt)'/nbpt
+    phi=tau.*tau.*(3-2*tau) ;
+    if orient=="reverse" then
+        a = (1-phi)*final(1) + phi*initial(1) ;
+    else
+        a = (1-phi)*initial(1) + phi*final(1) ;
+    end
+    da=a-final(1)
+
+    f=  final(2)+ p(1).*da.^3 +    p(2).*da.^4 +    p(3).*da.^5 ;
+    df  =       3*p(1).*da.^2 +  4*p(2).*da.^3 +  5*p(3).*da.^4 ;
+    ddf =       6*p(1).*da   + 12*p(2).*da.^2 + 20*p(3).*da.^3 ;
+
+    k = ddf ./ ((1+df.*df).^(3/2)) ;
+    state=[ a   f  atan(df) atan(k*bigL)]
+endfunction
+
+function display_car_trajectory(state)
+    bigL=1
+    set figure_style new;clf();show_window()
+    a=gca()
+    drawlater()
+    a.isoview="on"
+    a.data_bounds=[min(state(:,1))-0.5*bigL, min(state(:,2))-1.5*bigL
+    max(state(:,1))+1.5*bigL, max(state(:,2))+1.5*bigL]
+    rect=matrix(a.data_bounds',-1,1)
+    xpoly(rect([1 3 3 1]),rect([2,2,4,4]),"lines",1)
+    C=build_car()
+    Cinit=[];Cend=[];Cinter=[];
+    for k=1:size(C,"*")
+        Cinit=[Cinit copy(C(k))];
+        Cinter=[Cinter,copy(C(k))];
+        Cend=[Cend,copy(C(k))]
+    end
+    // starting configuration
+    draw_car(Cinit,state(1,:))
+    // end configuration
+    draw_car(Cend,state($,:))
+    // intermediate configuration (inversion of velocity)
+    draw_car(Cinter,state(ceil(size(state,1)/2),:)) ;
+    // trajectory of the linearizing output
+    t1=polyline([state(1,1) state(1,2);state(1,1) state(1,2)]) ;
+    t1.line_style=2;
+    realtimeinit(0.1)
+    for i=1:size(state,1)
+        realtime(i)
+        drawlater()
+        draw_car(C, state(i,:))
+        t1.data=[t1.data;state(i,1) state(i,2)];
+        drawnow()
+    end
+    for i=(1:30)+size(state,1),realtime(i),end
+endfunction
+
+
+function C=build_car()
+    //build the graphic object for the car
+    //
+    //the car
+    hcar=polyline([-2,7,8,8,7,-2,-2;-2,-2,-1,1,2,2,-2]'/6)
+    hcar.foreground=2
+
+    // rear wheels
+    hwheel1=polyline([[-1 1]/8; [1 1]/6]')
+    hwheel1.thickness=2
+
+    hwheel2=polyline([[-1 1]/8; -[1 1]/6]')
+    hwheel2.thickness=2
+
+    // front wheels
+    hwheel3=polyline([[7 9]/8;[1 1]/6]')
+    hwheel3.thickness=2
+    hwheel4=polyline([[7 9]/8;-[1 1]/6]')
+    hwheel4.thickness=2
+    //return vector of handle on the objects
+    C=[hcar,hwheel1,hwheel2,hwheel3,hwheel4]
+endfunction
+
+function draw_car(C,pos)
+    if is_handle_valid(C) then
+        drawlater()
+        [x,y,theta,phi]=(pos(1),pos(2),pos(3),pos(4))
+        bigL=1
+        Rc=[cos(theta) sin(theta);-sin(theta) cos(theta)]
+        // the car
+        xy = [-2,-2;7,-2;8,-1;8,1;7,2;-2,2;-2,-2]/6
+        C(1).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+        // rear wheels
+        xy=[[-1 1]/8; [1 1]/6]'
+        C(2).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+        xy=[[-1 1]/8; -[1 1]/6]'
+        C(3).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+        // front wheels
+        xy=[(1-cos(phi)/8) (1/6-sin(phi)/8)
+        (1+cos(phi)/8) (1/6+sin(phi)/8)]
+        C(4).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+        xy=[(1-cos(phi)/8) (-1/6-sin(phi)/8)
+        (1+cos(phi)/8) (-1/6+sin(phi)/8)]
+        C(5).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+        drawnow()
+    end
+endfunction
+
+function h=polyline(xy)
+    xpoly(xy(:,1),xy(:,2),"lines")
+    h=gce()
+endfunction
diff --git a/modules/cacsd/demos/flat/fcts.sci b/modules/cacsd/demos/flat/fcts.sci
new file mode 100755
index 000000000..5a27fcacb
--- /dev/null
+++ b/modules/cacsd/demos/flat/fcts.sci
@@ -0,0 +1,287 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function [xdot]=car(t,x)
+    //
+    //
+    xdot=zeros(1,4) ;
+    // car length for the control computation
+    LCpct = bigL * 1.;
+    // calcul de u1 et u2
+    //
+    if (t<0)
+        u1=0 ; u2=0 ;
+    elseif (t < bigT)
+        //      motion 1
+        tau = t/bigT ;
+        phi=tau*tau*(3-2*tau) ;
+        dphi = 6*tau*(1-tau) ;
+        a = (1-phi)*a1 + phi*a0 ;
+        da = ((a0-a1)/bigT) * dphi ;
+        f   =   p(1).*(a-a0)^3 +    p(2).*(a-a0)^4 +    p(3).*(a-a0)^5 ;
+        df  = 3*p(1).*(a-a0)^2 +  4*p(2).*(a-a0)^3 +  5*p(3).*(a-a0)^4 ;
+        ddf = 6*p(1).*(a-a0)   + 12*p(2).*(a-a0)^2 + 20*p(3).*(a-a0)^3 ;
+        dddf= 6*p(1).*1        + 24*p(2).*(a-a0)   + 60*p(3).*(a-a0)^2 ;
+        k = ddf / ((1+df*df)^(3/2)) ;
+        dk = ( dddf - 3*df*ddf*ddf/(1+df*df)) / ((1+df*df)^(3/2)) ;
+        u1 = ((1+df*df)^(1/2)) * da ;
+        u2 =  (bigL*dk / (1+(bigL*k)^2) ) * da ;
+    elseif (t < (2*bigT) )
+        //     motion 2
+        tau = t/bigT -1 ;
+        phi=tau*tau*(3-2*tau) ;
+        dphi = 6*tau*(1-tau) ;
+        a = (1-phi)*a0 + phi*a1 ;
+        da = ((a1-a0)/bigT) * dphi ;
+        f   =   p(1).*(a-a0)^3 +    p(2).*(a-a0)^4 +    p(3).*(a-a0)^5 ;
+        df  = 3*p(1).*(a-a0)^2 +  4*p(2).*(a-a0)^3 +  5*p(3).*(a-a0)^4 ;
+        ddf = 6*p(1).*(a-a0)   + 12*p(2).*(a-a0)^2 + 20*p(3).*(a-a0)^3 ;
+        dddf= 6*p(1).*1        + 24*p(2).*(a-a0)   + 60*p(3).*(a-a0)^2 ;
+        k = ddf / ((1+df*df)^(3/2)) ;
+        dk = ( dddf - 3*df*ddf*ddf/(1+df*df)) / ((1+df*df)^(3/2)) ;
+        u1 = ((1+df*df)^(1/2)) * da ;
+        u2 =  (LC*dk / (1+(LC*k)^2) ) * da ;
+    else
+        u1=0; u2=0;
+    end
+    xdot(1)= u1 * cos(x(1,3)) ;
+    xdot(2)= u1 * sin(x(1,3)) ;
+    xdot(3)= u1 * tan(x(1,4)) / bigL ;
+    xdot(4)= u2 ;
+
+    function [xdot]=car2T(t,x)
+        //
+        xdot=zeros(1,6) ;
+        // calcul de u1 et u2
+        //
+        if (t<0)
+            u1=0 ; u2=0 ;
+        elseif (t < bigT)
+            //      motion
+            tau = t/bigT ;
+            phi=tau*tau*(3-2*tau) ;
+            dphi = 6*tau*(1-tau) ;
+            a = (1-phi)*a1 + phi*a0 ;
+            da = ((a0-a1)/bigT) * dphi ;
+            [ff df d2f d3f d4f d5f] = cr2Tfjt(a) ;
+            [k2 k1 k0 dk0]=cr2Tfk(df,d2f,d3f,d4f,d5f) ;
+            u1 =  ( (1+(d2*k2)^2)*(1+(d1*k1)^2)*(1+df^2) )^(1/2) * da ;
+            u2 =  da * bigL * dk0 / (1+(bigL*k0)^2) ;
+        else
+            u1=0; u2=0;
+        end
+        xdot(1)= u1 * cos(x(1,3)) ;
+        xdot(2)= u1 * sin(x(1,3)) ;
+        xdot(3)= u1 * tan(x(1,6)) / bigL ;
+        xdot(4)= u1 * sin(x(1,3)-x(1,4)) / d1 ;
+        xdot(5)= u1 * sin(x(1,4)-x(1,5)) * cos(x(1,3)-x(1,4)) / d2 ;
+        xdot(6)= u2 ;
+
+        function [ff,df,d2f,d3f,d4f,d5f]=cr2Tfjt(a)
+            //
+            //
+            M= [
+            (a-a0)^5      (a-a0)^6     (a-a0)^7     (a-a0)^8       (a-a0)^9
+            5*(a-a0)^4    6*(a-a0)^5   7*(a-a0)^6    8*(a-a0)^7    9*(a-a0)^8
+            20*(a-a0)^3   30*(a-a0)^4  42*(a-a0)^5   56*(a-a0)^6   72*(a-a0)^7
+            60*(a-a0)^2  120*(a-a0)^3 210*(a-a0)^4  336*(a-a0)^5  504*(a-a0)^6
+            120*(a-a0)^1  360*(a-a0)^2 840*(a-a0)^3 1680*(a-a0)^4 3024*(a-a0)^5
+            120  720*(a-a0)^1 2520*(a-a0)^2 6720*(a-a0)^3 15120*(a-a0)^4
+            ]*p ;
+            ff  = M(1) + b0 ;
+            df  = M(2) ;
+            d2f = M(3) ;
+            d3f = M(4) ;
+            d4f = M(5) ;
+            d5f = M(6) ;
+
+
+            function [k2,k1,k0,dk0]=cr2Tfk(df,d2f,d3f,d4f,d5f)
+                //
+                //
+                // computation of curvatures from derivatives of b=f(a)
+                //
+                g = (1+df^2)^(-1/2);
+                dg = - df*d2f*g^3 ;
+                d2g = - g*g*(d2f^2*g+df*d3f*g+3*df*d2f*dg) ;
+                d3g = ....
+                - 2*g*dg*(d2f^2*g+df*d3f*g+3*df*d2f*dg) ....
+                - g^2*(3*d2f*d3f*g+df*d4f*g ....
+                + 4*d2f^2*dg+4*df*d3f*dg+3*df*d2f*d2g) ;
+                //
+                k2  = d2f * g^3  ;
+                dk2 = d3f*g^3 + 3*d2f*g^2*dg ;
+                d2k2 = g^2*(d4f*g+6*d3f*dg+3*d2f*d2g) + 6*g*dg^2*d2f ;
+                d3k2 = 2*g*dg*(d4f*g+6*d3f*dg+3*d2f*d2g) ....
+                + g^2*(d5f*g+7*d4f*dg+9*d3f*d2g+3*d2f*d3g) ....
+                +6*dg^3*d2f+12*g*dg*d2g*d2f+6*g*dg^2*d3f ;
+                //
+                g2 = (1+(d2*k2)^2)^(-1/2) ;
+                dg2 = -d2^2*k2*dk2*g2^3 ;
+                d2g2 = -d2^2*g2^2*(dk2^2*g2+k2*d2k2*g2+3*k2*dk2*dg2) ;
+                //
+                h2 = g2^3*g ;
+                dh2 = g2^2*(3*dg2*g+g2*dg);
+                d2h2 = 2*g2*dg2*(3*dg2*g+g2*dg) ....
+                + g2^2*(3*d2g2*g+4*dg2*dg+g2*d2g) ;
+                //
+                k1 = g2*k2 + d2*h2*dk2 ;
+                dk1 = dg2*k2 + g2*dk2 + d2 * (dh2*dk2+h2*d2k2) ;
+                d2k1 = d2g2*k2 + 2*dg2*dk2 + g2*d2k2 ....
+                + d2 * (d2h2*dk2+2*dh2*d2k2+h2*d3k2) ;
+                //
+                g1 = (1+(d1*k1)^2)^(-1/2) ;
+                dg1 = - d1^2*k1*dk1*g1^3 ;
+                //
+                k0 = g1*k1 + d1*g1^3*g2*g*dk1 ;
+                dk0 = dg1*k1 + g1*dk1 ....
+                + d1*g1^2*(3*dg1*g2*g*dk1+g1*dg2*g*dk1 ....
+                +  g1*g2*dg*dk1+g1*g2*g*d2k1) ;
+
+
+
+                function coef=cr2Tkf(b,theta2,theta12,theta01,phi)
+                    //
+                    //
+                    M = [
+                    1*(a1-a0)^5    1*(a1-a0)^6   1*(a1-a0)^7    1*(a1-a0)^8    1*(a1-a0)^9
+                    5*(a1-a0)^4    6*(a1-a0)^5   7*(a1-a0)^6    8*(a1-a0)^7    9*(a1-a0)^8
+                    20*(a1-a0)^3   30*(a1-a0)^4  42*(a1-a0)^5   56*(a1-a0)^6   72*(a1-a0)^7
+                    60*(a1-a0)^2  120*(a1-a0)^3 210*(a1-a0)^4  336*(a1-a0)^5  504*(a1-a0)^6
+                    120*(a1-a0)^1  360*(a1-a0)^2 840*(a1-a0)^3 1680*(a1-a0)^4 3024*(a1-a0)^5
+                    ] ;
+                    //
+                    //
+                    df = tan(theta2) ;
+                    //
+                    // curvatures
+                    k2=tan(theta12)/d2;k1=tan(theta01)/d1;k0=tan(phi)/bigL;
+                    //
+                    ddf = k2*((1+df*df)^(3/2)) ;
+                    //
+                    // derivative of k2
+                    dk2ds2 = ( (1+(d2*k2)^2)/d2 )*( (1+(d2*k2)^2)^(1/2)*k1 - k2 ) ;
+                    dk2 = (1+df*df)^(1/2) * dk2ds2 ;
+                    //
+                    dddf = dk2 * (1+df*df)^(3/2)  + 3*k2*df*ddf*(1+df*df)^(1/2) ;
+                    //
+                    // second derivative of k2
+                    dk1ds1 = ( (1+(d1*k1)^2)/d1 )*( (1+(d1*k1)^2)^(1/2)*k0 - k1 ) ;
+                    dk1ds2 = dk1ds1 *  (1+(d2*k2)^2)^(1/2) ;
+                    //
+                    ddk2ds2 =  ....
+                    3*d2*k2*dk2ds2*(1+(d2*k2)^2)^(1/2)*k1 ....
+                    + (1+(d2*k2)^2)^(3/2)*dk1ds2/d2 ...
+                    - (1/d2+3*d2*k2*k2)*dk2ds2 ;
+                    //
+                    ddk2 = ....
+                    df*ddf*((1+df*df)^(-1/2))*dk2ds2 ....
+                    + (1+df*df)*ddk2ds2 ;
+                    //
+                    ddddf =  ....
+                    ddk2 * (1+df*df)^(3/2) ....
+                    + 6*dk2*df*ddf*(1+df*df)^(1/2) ....
+                    + 3*k2*ddf*ddf*(1+df*df)^(1/2) ....
+                    + 3*k2*df*dddf*(1+df*df)^(1/2) ....
+                    + 3*k2*(df*ddf)^2*(1+df*df)^(-1/2) ;
+                    //
+                    //
+                    //
+                    //
+                    coef =inv(M)* [b;df;ddf;dddf;ddddf] ;
+
+                    function []=ptcr(state)
+                        //
+                        x=state(1);y=state(2);theta=state(3);phi=state(4);
+                        //
+                        // plot the car
+                        //
+                        // rear wheels
+                        wheel1 =  [ -1/8 1/8 ;  1/6  1/6 ] ;
+                        wheel2 =  [ -1/8 1/8 ; -1/6 -1/6 ] ;
+                        // front wheels
+                        wheel3 = [(1-cos(phi)/8)   (1+cos(phi)/8);
+                        (1/6-sin(phi)/8) (1/6+sin(phi)/8)   ] ;
+                        wheel4 = [(1-cos(phi)/8)    (1+cos(phi)/8);
+                        (-1/6-sin(phi)/8) (-1/6+sin(phi)/8)  ] ;
+                        //
+                        // car + wheels
+                        xy_car = [ [  -1/3   7/6   4/3  4/3  7/6  -1/3  -1/3;
+                        -1/3  -1/3  -1/6  1/6  1/3   1/3  -1/3;
+                        ] wheel1 wheel2 wheel3 wheel4 ] ;
+                        xy_car = ....
+                        [x x x x x x x x x x x x x x x;
+                        y y y y y y y y y y y y y y y ] + ....
+                        bigL*[cos(theta) -sin(theta);sin(theta) cos(theta)]*xy_car  ;
+                        //
+                        xpoly(xy_car(1,1:7),xy_car(2,1:7),"lines") ;
+                        // The 4 wheels
+                        xpolys(matrix(xy_car(1,8:15),2,4),matrix(xy_car(2,8:15),2,4),[1 1 1 1])
+
+
+                        function ptcr2T(state)
+                            // plot the car with two trailers
+                            //
+                            x0=state(1);y0=state(2);theta0=state(3);theta1=state(4);
+                            theta2=state(5);phi=state(6);
+                            theta = theta0 ; x=x0 ; y=y0 ;
+                            // rear wheels
+                            wheel1 =  [ -1/8 1/8 ;  1/6  1/6 ] ;
+                            wheel2 =  [ -1/8 1/8 ; -1/6 -1/6 ] ;
+                            // front wheels
+                            wheel3 = [
+                            (1-cos(phi)/8)   (1+cos(phi)/8)
+                            (1/6-sin(phi)/8) (1/6+sin(phi)/8)
+                            ] ;
+                            wheel4 = [
+                            (1-cos(phi)/8)    (1+cos(phi)/8)
+                            (-1/6-sin(phi)/8) (-1/6+sin(phi)/8)
+                            ] ;
+                            //
+                            // car + wheels
+                            xy_car = [ [ -1/3   7/6   4/3  4/3  7/6  -1/3  -1/3
+                            -1/3  -1/3  -1/6  1/6  1/3   1/3  -1/3 ] wheel1 wheel2 wheel3 wheel4 ] ;
+                            xy_car = diag([x,y])*ones(2,15) + ...
+                            bigL*[cos(theta0) -sin(theta0);sin(theta0) cos(theta0)]*xy_car  ;
+                            //
+                            //
+                            //     Trailer 1
+                            x = x - d1*cos(theta1) ;
+                            y = y - d1*sin(theta1) ;
+
+                            //
+                            // shift + rotation
+                            xy_T1 =diag([x,y])*ones(2,11) + ...
+                            [cos(theta1) -sin(theta1);sin(theta1) cos(theta1)]*xy_T1  ;
+
+                            //
+                            //     Trailer 2
+                            x = x - d2*cos(theta2) ;
+                            y = y - d2*sin(theta2) ;
+
+                            //
+                            // shift + rotation
+                            xy_T2 = diag([x,y])*ones(2,11)+ ...
+                            [cos(theta2) -sin(theta2);sin(theta2) cos(theta2)]*xy_T2  ;
+                            //
+                            // plots
+                            //
+                            //  Car
+                            xpoly(xy_car(1,1:7),xy_car(2,1:7),"lines") ;
+                            //   The 4 wheels
+                            xpolys(matrix(xy_car(1,8:15),2,4),matrix(xy_car(2,8:15),2,4),[1 1 1 1])
+
+                            // trailer 1
+                            xpoly(xy_T1(1,1:5),xy_T1(2,1:5),"lines") ;
+                            // hitch and wheels
+                            xpolys(matrix(xy_T1(1,6:11),2,3),matrix(xy_T1(2,6:11),2,3),[1 1 1 1]);
+
+                            //
+                            // trailer 2
+                            xpoly(xy_T2(1,1:5),xy_T2(2,1:5),"lines") ;
+                            // hitch and wheels
+                            xpolys(matrix(xy_T2(1,6:11),2,3),matrix(xy_T2(2,6:11),2,3),[1 1 1 1]);
diff --git a/modules/cacsd/demos/flat/flat.dem.gateway.sce b/modules/cacsd/demos/flat/flat.dem.gateway.sce
new file mode 100755
index 000000000..b61b5f072
--- /dev/null
+++ b/modules/cacsd/demos/flat/flat.dem.gateway.sce
@@ -0,0 +1,12 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is released under the 3-clause BSD license. See COPYING-BSD.
+
+demopath = get_absolute_file_path("flat.dem.gateway.sce");
+
+subdemolist = ["Car"                ,"car.dem.sce"   ; ..
+"Two trailers truck" ,"truck.dem.sce" ]
+
+subdemolist(:,2) = demopath + subdemolist(:,2);
+clear demopath;
\ No newline at end of file
diff --git a/modules/cacsd/demos/flat/truck.dem.sce b/modules/cacsd/demos/flat/truck.dem.sce
new file mode 100755
index 000000000..dedc7b4c4
--- /dev/null
+++ b/modules/cacsd/demos/flat/truck.dem.sce
@@ -0,0 +1,23 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is released under the 3-clause BSD license. See COPYING-BSD.
+
+// Clear the functions defined in truck.sci to avoid warnings
+
+
+function demo_truck()
+    exec("SCI/modules/cacsd/demos/flat/truck.sci",-1);
+
+    initial   = [-2;3;0.5235988;0;0;1];
+    final     = [0;0;0;0;0;0];
+    state     = truck_solve(initial,final);
+    my_handle = scf(100001);
+    clf(my_handle,"reset");
+    toolbar(my_handle.figure_id,"off") ;
+    display_truck_trajectory(state);
+
+endfunction
+
+demo_truck();
+clear demo_truck;
diff --git a/modules/cacsd/demos/flat/truck.sci b/modules/cacsd/demos/flat/truck.sci
new file mode 100755
index 000000000..43f948302
--- /dev/null
+++ b/modules/cacsd/demos/flat/truck.sci
@@ -0,0 +1,350 @@
+
+
+function state=truck_solve(initial,final)
+    //
+    //  CAR  WITH 2 TRAILERS PACKING VIA FLATNESS AND FRENET FORMULAS
+    //
+    //    explicit computation and visualisation of the motions.
+    //
+    //    February 1993
+    //
+    // ............................................................
+    // :         pierre ROUCHON  <rouchon@cas.ensmp.fr>           :
+    // : Centre Automatique et Systemes, Ecole des Mines de Paris :
+    // : 60, bd Saint Michel -- 75272 PARIS CEDEX 06, France      :
+    // :    Telephone: (1) 40 51 91 15 --- Fax: (1) 43 54 18 93   :
+    // :..........................................................:
+    //
+    // initial: initial position [x,y,theta1,theta2,theta3,phi]
+    // final  :   final position [x,y,theta1,theta2,theta3,phi]
+    //        theta1,theta2,theta3  : the car and the trailers angles
+    //        phi    : the front wheel angle
+
+    bigT = 1 ;//basic time interval for one  smooth motion (s)
+    bigL = 1 ;// car length (m)
+    d1 = 1.5 ; d2 = 1 ; //trailers length
+
+    // computation of  intermediate configuration
+    LL=bigL+d1+d2
+    x0 = max(initial(1),final(2))   ....
+    + LL*abs(tan(initial(3))) ...
+    + LL*abs(tan(initial(4))) ...
+    + LL*abs(tan(initial(5))) ...
+    + LL*(abs(initial(2)-final(2))/(d1+d2+bigL))^(1/2) ;
+    y0 = (initial(2)+final(2))/2 ;
+    intermediate=[x0,y0,0,0,0,0]'
+
+    // first polynomial curve
+    state=truck_polynomial_curve(initial,intermediate,"direct")
+    //
+    // second polynomial curve
+    state = [ state;
+    truck_polynomial_curve(final,intermediate,"reverse")	    ]
+
+endfunction
+
+
+
+function state=truck_polynomial_curve(initial,final,orient)
+
+    nbpt = 40 ; //  sampling of motion
+    phi = initial(6) ;
+
+    theta2 = initial(3) ;
+    theta1 = initial(4)+theta2
+    theta0 = initial(5)+theta1
+
+    x0=initial(1)+d2*cos(theta2)+d1*cos(theta1) ;
+    y0=initial(2)+d2*sin(theta2)+d1*sin(theta1) ;
+    if orient<>"reverse" then
+        state = [x0 y0 theta0 theta1 theta2 phi] ;
+    else
+        state=[]
+    end
+    a0=final(1);a1=initial(1);b0=final(2)
+    db = initial(2)-final(2)
+    p = cr2Tkf(db,initial(3),initial(4),initial(5),phi) ;
+
+    tau=(0:nbpt)'/nbpt
+    phi=tau.*tau.*(3-2*tau) ;
+    if orient=="reverse" then
+        aa = (1-phi)*final(1) + phi*initial(1) ;
+    else
+        aa = (1-phi)*initial(1) + phi*final(1) ;
+    end
+    for i=1:(nbpt+1)
+        [bb,df,d2f,d3f,d4f,d5f] = cr2Tfjt(aa(i)) ;
+        [k2,k1,k0,dk0]=cr2Tfk(df,d2f,d3f,d4f,d5f) ;
+        theta2 = atan(df);
+        theta1 = atan(k2*d2)+theta2;
+        theta0 = atan(k1*d1) + theta1 ;
+        phi = atan(k0*bigL) ;
+        x0=aa(i)+d2*cos(theta2)+d1*cos(theta1) ;
+        y0=bb+d2*sin(theta2)+d1*sin(theta1) ;
+        state= [ state;
+        x0 y0 theta0 theta1 theta2 phi] ;
+    end
+endfunction
+
+function [ff,df,d2f,d3f,d4f,d5f]=cr2Tfjt(a)
+    //
+    //
+    da = a-a0
+    M= [  da^5      da^6      da^7      da^8       da^9
+    5*da^4    6*da^5    7*da^6    8*da^7     9*da^8
+    20*da^3   30*da^4   42*da^5   56*da^6    72*da^7
+    60*da^2  120*da^3  210*da^4  336*da^5   504*da^6
+    120*da^1  360*da^2  840*da^3 1680*da^4  3024*da^5
+    120       720*da^1 2520*da^2 6720*da^3 15120*da^4]*p ;
+    ff  = M(1) + b0 ;
+    df  = M(2) ;
+    d2f = M(3) ;
+    d3f = M(4) ;
+    d4f = M(5) ;
+    d5f = M(6) ;
+endfunction
+
+function coef=cr2Tkf(b,theta2,theta12,theta01,phi)
+    //
+    //
+    da = a1-a0
+    M = [1*da^5    1*da^6   1*da^7    1*da^8    1*da^9
+    5*da^4    6*da^5   7*da^6    8*da^7    9*da^8
+    20*da^3   30*da^4  42*da^5   56*da^6   72*da^7
+    60*da^2  120*da^3 210*da^4  336*da^5  504*da^6
+    120*da^1  360*da^2 840*da^3 1680*da^4 3024*da^5] ;
+    //
+    //
+    df = tan(theta2) ;
+    //
+    // curvatures
+    k2=tan(theta12)/d2;k1=tan(theta01)/d1;k0=tan(phi)/bigL;
+    //
+    ddf = k2*((1+df*df)^(3/2)) ;
+    //
+    // derivative of k2
+    dk2ds2 = ( (1+(d2*k2)^2)/d2 )*( (1+(d2*k2)^2)^(1/2)*k1 - k2 ) ;
+    dk2 = (1+df*df)^(1/2) * dk2ds2 ;
+    //
+    dddf = dk2 * (1+df*df)^(3/2)  + 3*k2*df*ddf*(1+df*df)^(1/2) ;
+    //
+    // second derivative of k2
+    dk1ds1 = ( (1+(d1*k1)^2)/d1 )*( (1+(d1*k1)^2)^(1/2)*k0 - k1 ) ;
+    dk1ds2 = dk1ds1 *  (1+(d2*k2)^2)^(1/2) ;
+    //
+    ddk2ds2 =  ....
+    3*d2*k2*dk2ds2*(1+(d2*k2)^2)^(1/2)*k1 ....
+    + (1+(d2*k2)^2)^(3/2)*dk1ds2/d2 ...
+    - (1/d2+3*d2*k2*k2)*dk2ds2 ;
+    //
+    ddk2 = ....
+    df*ddf*((1+df*df)^(-1/2))*dk2ds2 ....
+    + (1+df*df)*ddk2ds2 ;
+    //
+    ddddf =  ....
+    ddk2 * (1+df*df)^(3/2) ....
+    + 6*dk2*df*ddf*(1+df*df)^(1/2) ....
+    + 3*k2*ddf*ddf*(1+df*df)^(1/2) ....
+    + 3*k2*df*dddf*(1+df*df)^(1/2) ....
+    + 3*k2*(df*ddf)^2*(1+df*df)^(-1/2) ;
+    //
+    //
+    //
+    //
+    coef =inv(M)* [b;df;ddf;dddf;ddddf] ;
+endfunction
+
+function [k2,k1,k0,dk0]=cr2Tfk(df,d2f,d3f,d4f,d5f)
+    //
+    //
+    // computation of curvatures from derivatives of b=f(a)
+    //
+    g = (1+df^2)^(-1/2);
+    dg = - df*d2f*g^3 ;
+    d2g = - g*g*(d2f^2*g+df*d3f*g+3*df*d2f*dg) ;
+    d3g = ....
+    - 2*g*dg*(d2f^2*g+df*d3f*g+3*df*d2f*dg) ....
+    - g^2*(3*d2f*d3f*g+df*d4f*g ....
+    + 4*d2f^2*dg+4*df*d3f*dg+3*df*d2f*d2g) ;
+    //
+    k2  = d2f * g^3  ;
+    dk2 = d3f*g^3 + 3*d2f*g^2*dg ;
+    d2k2 = g^2*(d4f*g+6*d3f*dg+3*d2f*d2g) + 6*g*dg^2*d2f ;
+    d3k2 = 2*g*dg*(d4f*g+6*d3f*dg+3*d2f*d2g) ....
+    + g^2*(d5f*g+7*d4f*dg+9*d3f*d2g+3*d2f*d3g) ....
+    +6*dg^3*d2f+12*g*dg*d2g*d2f+6*g*dg^2*d3f ;
+    //
+    g2 = (1+(d2*k2)^2)^(-1/2) ;
+    dg2 = -d2^2*k2*dk2*g2^3 ;
+    d2g2 = -d2^2*g2^2*(dk2^2*g2+k2*d2k2*g2+3*k2*dk2*dg2) ;
+    //
+    h2 = g2^3*g ;
+    dh2 = g2^2*(3*dg2*g+g2*dg);
+    d2h2 = 2*g2*dg2*(3*dg2*g+g2*dg) ....
+    + g2^2*(3*d2g2*g+4*dg2*dg+g2*d2g) ;
+    //
+    k1 = g2*k2 + d2*h2*dk2 ;
+    dk1 = dg2*k2 + g2*dk2 + d2 * (dh2*dk2+h2*d2k2) ;
+    d2k1 = d2g2*k2 + 2*dg2*dk2 + g2*d2k2 ....
+    + d2 * (d2h2*dk2+2*dh2*d2k2+h2*d3k2) ;
+    //
+    g1 = (1+(d1*k1)^2)^(-1/2) ;
+    dg1 = - d1^2*k1*dk1*g1^3 ;
+    //
+    k0 = g1*k1 + d1*g1^3*g2*g*dk1 ;
+    dk0 = dg1*k1 + g1*dk1 ....
+    + d1*g1^2*(3*dg1*g2*g*dk1+g1*dg2*g*dk1 ....
+    +  g1*g2*dg*dk1+g1*g2*g*d2k1) ;
+
+endfunction
+
+function display_truck_trajectory(state)
+    bigL = 1 ; d1 = 1.5 ; d2 = 1;
+    a=gca();
+    drawlater();
+    a.isoview="on"
+    a.data_bounds=[min(state(:,1))-1.5*(d1+d2), min(state(:,2))-bigL
+    max(state(:,1))+1.5*bigL, max(state(:,2))+bigL]
+    rect=matrix(a.data_bounds',-1,1)
+    xpoly(rect([1 3 3 1]),rect([2,2,4,4]),"lines",1)
+    C=build_truck()
+    Cinit=[];
+    Cend=[];
+    Cinter=[];
+    for k=1:size(C,"*")
+        Cinit=[Cinit copy(C(k))];
+        Cinter=[Cinter,copy(C(k))];
+        Cend=[Cend,copy(C(k))]
+    end
+    // starting configuration
+    draw_truck(Cinit,state(1,:))
+    // end configuration
+    draw_truck(Cend,state($,:))
+    // intermediate configuration (inversion of velocity)
+    draw_truck(Cinter,state(ceil(size(state,1)/2),:)) ;
+    // trajectory of the linearizing output
+    x_lin = state(:,1)-d1*cos(state(:,4))-d2*cos(state(:,5)) ;
+    y_lin = state(:,2)-d1*sin(state(:,4))-d2*sin(state(:,5)) ;
+
+    t1=polyline([x_lin(1) y_lin(1);x_lin(1) y_lin(1)]) ;
+    t2=polyline([state(1,1) state(1,2);state(1,1) state(1,2)]) ;
+
+    t1.line_style=2;
+    t2.line_style=2;t2.foreground=5
+    realtimeinit(0.2)
+    for i=1:size(state,1)
+        realtime(i)
+        drawlater()
+        draw_truck(C, state(i,:))
+        t1.data=[t1.data;x_lin(i), y_lin(i)];
+        t2.data=[t2.data;state(i,1), state(i,2)];
+        drawnow()
+    end
+    for i=(1:30)+size(state,1),realtime(i),end
+endfunction
+
+
+function C=build_truck()
+    //build the graphic object for the truck
+    //
+    //the car
+    hcar=polyline([-2,7,8,8,7,-2,-2;-2,-2,-1,1,2,2,-2]'/6)
+    hcar.foreground=2
+    // rear wheels
+    hwheel1=polyline([[-1 1]/8; [1 1]/6]')
+    hwheel1.thickness=2
+
+    hwheel2=polyline([[-1 1]/8; -[1 1]/6]')
+    hwheel2.thickness=2
+
+    // front wheels
+    hwheel3=polyline([[7 9]/8;[1 1]/6]')
+    hwheel3.thickness=2
+    hwheel4=polyline([[7 9]/8;-[1 1]/6]')
+    hwheel4.thickness=2
+
+    //Trailer 1
+
+    ht1=polyline([-1,1,1,-1,-1;-1,-1,1,1,-1]'*bigL/3)
+    ht1.foreground=2
+    hwheel5=polyline([[-1 1]/8;[1 1]/6]'*bigL)
+    hwheel5.thickness=2
+    hwheel6=polyline([[-1 1]/8;-[1 1]/6]'*bigL)
+    hwheel6.thickness=2
+    hhitch1=polyline([bigL/3  d1;0 0]')
+    hhitch1.foreground=2
+    //Trailer 2
+    ht2=polyline([-1,1,1,-1,-1;-1,-1,1,1,-1]'*bigL/3)
+    ht2.foreground=2
+    hwheel7=polyline([[-1 1]/8;[1 1]/6]'*bigL)
+    hwheel7.thickness=2
+    hwheel8=polyline([[-1 1]/8;-[1 1]/6]'*bigL)
+    hwheel8.thickness=2
+    hhitch2=polyline([bigL/3  d2;0 0]')
+    hhitch2.foreground=2
+
+
+
+    //return vector of handle on the objects
+    C=[hcar,hwheel1,hwheel2,hwheel3,hwheel4,..
+    ht1,hwheel5, hwheel6,hhitch1,..
+    ht2, hwheel7,hwheel8,hhitch2]
+endfunction
+
+function draw_truck(C,pos)
+    drawlater()
+    [x,y,theta1,theta2,theta3,phi]=(pos(1),pos(2),pos(3),pos(4),pos(5),pos(6))
+    bigL = 1 ; d1 = 1.5 ; d2 = 1;
+    Rc=[cos(theta1) sin(theta1);-sin(theta1) cos(theta1)]
+    // the car
+    xy = [-2,-2;7,-2;8,-1;8,1;7,2;-2,2;-2,-2]/6
+    C(1).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+    // rear wheels
+    xy=[[-1 1]/8; [1 1]/6]'
+    C(2).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+    xy=[[-1 1]/8; -[1 1]/6]'
+    C(3).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+    // front wheels
+    xy=[(1-cos(phi)/8) (1/6-sin(phi)/8)
+    (1+cos(phi)/8) (1/6+sin(phi)/8)]
+    C(4).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+    xy=[(1-cos(phi)/8) (-1/6-sin(phi)/8)
+    (1+cos(phi)/8) (-1/6+sin(phi)/8)]
+    C(5).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+
+    //Trailer 1
+    Rc=[cos(theta2) sin(theta2);-sin(theta2) cos(theta2)]
+    x = x - d1*cos(theta2) ;
+    y = y - d1*sin(theta2) ;
+    xy = [-1,1,1,-1,-1;-1,-1,1,1,-1]'*bigL/3;
+    C(6).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+    //wheels
+    xy=[[-1 1]/8; [1 1]/6]'
+    C(7).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+    xy=[[-1 1]/8; -[1 1]/6]'
+    C(8).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+    //hitch
+    xy=[bigL/3  d1;0 0]'
+    C(9).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+
+    //Trailer 2
+    Rc=[cos(theta3) sin(theta3);-sin(theta3) cos(theta3)]
+    x = x - d2*cos(theta3) ;
+    y = y - d2*sin(theta3) ;
+    xy = [-1,1,1,-1,-1;-1,-1,1,1,-1]'*bigL/3;
+    C(10).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+    //wheels
+    xy=[[-1 1]/8; [1 1]/6]'
+    C(11).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+    xy=[[-1 1]/8; -[1 1]/6]'
+    C(12).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+    //hitch
+    xy=[bigL/3  d2;0 0]'
+    C(13).data=ones(xy)*diag([x;y])+bigL*xy*Rc
+    drawnow()
+endfunction
+
+function h=polyline(xy)
+    xpoly(xy(:,1),xy(:,2),"lines")
+    h=gce()
+endfunction
diff --git a/modules/cacsd/demos/lqg/lqg.dem b/modules/cacsd/demos/lqg/lqg.dem
new file mode 100755
index 000000000..52efc765b
--- /dev/null
+++ b/modules/cacsd/demos/lqg/lqg.dem
@@ -0,0 +1,168 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function demo_lqg()
+
+  mode(-1);
+  oldln=lines();
+  lines(0);
+  //display the diagram
+  x=[5,10,20,40,50,70,80,90];xmin=-10;xmax=100;
+  y=[22,28,30,32];ymin=12;ymax=40;
+
+  xx=[x([1 2 2 7 4 6 3 4 5 6 3 3 5 5]);x([3,2,7,7,5,8,3,4,5,6,4,4,6,6])];
+  yy=[y([3,1,1,1,3,3,2,2,2,2,2,4,2,4]);y([3,3,1,3,3,3,4,4,4,4,2,4,2,4])];
+
+  scf(100001);
+  clf();
+  show_window();
+  a=gca();
+  a.data_bounds=[xmin ymin;xmax ymax];
+  a.box='on';
+  xpolys(xx,yy)
+  xstring(28,30,'K');
+  xstring(56,30,'Plant');
+  xstring(12,28.80,'-');
+  xtitle('PLANT   and   CONTROLLER')
+
+
+  path=get_absolute_file_path('lqg.dem');
+  s=poly(0,'s');z=poly(0,'z');
+  messagebox(['Simple example of SISO LQG Design';
+           'Demo is in file '+path+'lqg.dem';
+           'Computes the LQG compensator and plots response'],"modal");
+
+  n=x_choose(['Continuous time';'Discrete time'],'Select time domain');
+
+  select n
+  case 0
+    return
+  case 1
+    mode(1)
+    s=poly(0,'s');
+    str='(s-1)/(s^2-5*s+1)';
+    rep=x_dialog('Nominal plant?',str)
+    if rep==[] then return,end
+    Plant=evstr(rep);
+    Plant=syslin('c',Plant);
+    mode(-1)
+  case 2
+    mode(1)
+    z=poly(0,'z');
+    str='(z+1)/(z^2-5*z+2)'
+    rep=x_dialog('Nominal plant?',str)
+    if rep==[] then return,end
+    Plant=evstr(rep);
+    Plant=syslin('d',Plant);
+    mode(-1)
+  end
+  
+  mode(1)
+
+  //Nominal Plant
+
+  P22=tf2ss(Plant);    //...in state-space form
+  [ny,nu,nx]=size(P22);
+  messagebox('Now enter weighting matrices',"modal");
+  rep=x_matrix('x-weighting matrix',eye(nx,nx))
+  if rep==[] then return,end
+  Qx=evstr(rep);
+  rep=x_matrix('u-weighting matrix',eye(nu,nu));
+  if rep==[] then return,end
+  Qu=evstr(rep);
+  bigQ=sysdiag(Qx,Qu);
+  rep=x_matrix('x-noise covariance matrix',eye(nx,nx))
+  if rep==[] then return,end
+  Rx=evstr(rep);
+  rep=x_matrix('y-noise covariance matrix',eye(ny,ny))
+  if rep==[] then return,end
+  Ry=evstr(rep);
+  bigR=sysdiag(Rx,Ry);
+  [Plqg,r]=lqg2stan(P22,bigQ,bigR);     //LQG pb as a standard problem
+  Klqg=lqg(Plqg,r);                     //LQG compensator
+  
+  disp(spec(h_cl(Plqg,r,Klqg)),'closed loop eigenvalues:');    //Check internal stability
+  [Slqg,Rlqg,Tlqg]=sensi(P22,Klqg);  //Sensitivity functions
+
+  disp(clean(ss2tf(Slqg)),'Sensitivity function');
+  disp(clean(ss2tf(Tlqg)),'Complementary sensitivity function');
+  
+  mode(-1);
+
+  messagebox('Closed-loop response',"modal");
+
+  resp=['Frequency response';'Time response'];
+  while %t do
+    n=x_choose(resp,'Select response(s)');
+    select n
+    case 0
+      disp("LQG demo stops!");break;
+    case 1 then
+      mode(1)
+      scf(100001);
+      clf();
+      show_window();
+      bode(Tlqg);
+      mode(-1)
+    case 2
+      if Plant(4)=='c' then
+        mode(1)
+        defv=['0.1';'20'];
+        _title='Enter Sampling period and Tmax';
+        rep=x_mdialog(_title,['Sampling period?';'Tmax?'],defv)
+        if rep==[] then break,end
+        dttmax=evstr(rep)
+        dt=evstr(dttmax(1));
+        tmax=evstr(dttmax(2));
+        t=0:dt/5:tmax;
+        n1=x_choose(['Step response?';'Impulse response?'],'Simulation:');
+        select n1
+        case 1 then
+          scf(100002);
+          clf();
+          show_window();
+          plot2d([t',t'],[(csim('step',t,Tlqg))',ones(t')]);
+        case 2 then
+          scf(100002);
+          clf();
+          show_window();
+          plot2d([t',t'],[(csim('impul',t,Tlqg))',0*t']);
+        end
+        mode(-1)
+      elseif Plant(4)=='d' then
+        mode(1)
+        defv=['30'];
+        _title='Tmax?';
+        rep=x_mdialog(_title,['Tmax='],defv)
+        if rep==[] then break,end
+        Tmax=evstr(rep);
+        mode(-1)
+        n2=x_choose(['Step response?';'Impulse response?'],'Simulation:');
+        select n2
+        case 0 then
+        break;
+        case 1 then
+          mode(1)
+          u=ones(1,Tmax);u(1)=0;
+          scf(100002);clf();show_window();
+          plot2d([(1:Tmax)',(1:Tmax)'],[(dsimul(Tlqg,u))',(ones(1:Tmax)')])
+          mode(-1)
+        case 2 then
+          mode(1)
+          u=zeros(1,Tmax);u(1)=1;
+          scf(100002);clf();show_window();
+          plot2d((1:Tmax)',(dsimul(Tlqg,u))')
+          mode(-1)
+        end
+      end
+    end
+  end
+  lines(oldln(2))
+endfunction
+
+demo_lqg()
+clear demo_lqg;
diff --git a/modules/cacsd/demos/lqg/lqg2.dem b/modules/cacsd/demos/lqg/lqg2.dem
new file mode 100755
index 000000000..1f13001f1
--- /dev/null
+++ b/modules/cacsd/demos/lqg/lqg2.dem
@@ -0,0 +1,44 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function demo_lqg2()
+
+  s=poly(0,'s');
+  Plant=[1/(s+1),s/(s-1)^2;(s-1)*s/(s^2-3*s+2),2/s];
+  Plant=syslin('c',[1/(s+1)*s/(s-1)^2]);     //Nominal Plant
+  P22=tf2ss(Plant);    //...in state-space form
+  [ny,nu,nx]=size(P22);
+  rand('seed',0);rand('normal');
+  bigQ=rand(nx+nu,nx+nu);bigQ=bigQ*bigQ';
+  bigR=rand(nx+ny,nx+ny);bigR=bigR*bigR';  //random weighting matrices
+  [Plqg,r]=lqg2stan(P22,bigQ,bigR);     //LQG pb as a standard problem
+  Klqg=lqg(Plqg,r);          //Controller
+  spec(h_cl(Plqg,r,Klqg))    //Check internal stability
+  [Slqg,Rlqg,Tlqg]=sensi(P22,Klqg);  //Sensitivity functions
+  frq=logspace(-3,3);     //10^-3  to 10^3 
+  y=svplot(Slqg);       //Computes singular values;
+  clf();
+  scf(100001);
+  gainplot(frq,y)   //Plot sing. values
+  w1=1/(s+1);
+  w2=100;
+  [Ptmp,r]=augment(P22,'SR');         //"S/KS" problem
+  Pinf=sysdiag(w1,w2,1)*Ptmp;        //Weighting functions
+  [Kinf,ro]=h_inf(Pinf,r,0,0.1,50);
+  [Sinf,Rinf,Tinf]=sensi(P22,Kinf);  //Sensitivity functions
+  y=svplot(Sinf);       //Computes singular values;
+  scf(100002);
+  show_window();
+  gainplot(frq,y)   //Plot sing. values
+
+  clf();t=0:0.01:30;u=sin(t);
+  plot2d([t',t'],[u',(flts(u,dscr(Tlqg,0.1))')])
+  
+endfunction
+  
+demo_lqg2();
+clear demo_lqg2  
\ No newline at end of file
diff --git a/modules/cacsd/demos/lqg/scheme.dem b/modules/cacsd/demos/lqg/scheme.dem
new file mode 100755
index 000000000..aa6481927
--- /dev/null
+++ b/modules/cacsd/demos/lqg/scheme.dem
@@ -0,0 +1,28 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function demo_scheme()
+
+  x=[5,10,20,40,50,70,80,90];xmin=-10;xmax=100;
+  y=[22,28,30,32];ymin=12;ymax=40;
+
+  xx=[xmin,xmin,x([1 2 2 7 4 6 3 4 5 6 3 3 5 5]);xmax,xmax,x([3,2,7,7,5,8,3,4,5,6,4,4,6,6])];
+  yy=[ymin,ymax,y([3,1,1,1,3,3,2,2,2,2,2,4,2,4]);ymin,ymax,y([3,3,1,3,3,3,4,4,4,4,2,4,2,4])];
+    
+  scf(100001);
+  clf();
+  show_window();
+  a=gca();
+  a.data_bounds=[xmin ymin;xmax ymax];
+  xpolys(xx,yy)
+  xstring(28,30,'K');xstring(56,30,'Plant');xstring(12,28.80,'-');
+  xtitle('PLANT   and   CONTROLLER')
+  
+endfunction
+
+demo_scheme();
+clear demo_scheme;
diff --git a/modules/cacsd/demos/mixed.dem b/modules/cacsd/demos/mixed.dem
new file mode 100755
index 000000000..26db4f7f9
--- /dev/null
+++ b/modules/cacsd/demos/mixed.dem
@@ -0,0 +1,87 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function demo_mixed()
+  mode(-1)
+  lines(0);
+  //display the diagram
+  x=[5,10,20,40,50,70,80,90];xmin=-10;xmax=100;
+  y=[22,28,30,32];ymin=12;ymax=40;
+
+  xx=[xmin,xmin,x([1 2 2 7 4 6 3 4 5 6 3 3 5 5]);xmax,xmax,x([3,2,7,7,5,8,3,4,5,6,4,4,6,6])];
+  yy=[ymin,ymax,y([3,1,1,1,3,3,2,2,2,2,2,4,2,4]);ymin,ymax,y([3,3,1,3,3,3,4,4,4,4,2,4,2,4])];
+    
+  scf(100001);clf();show_window();
+  plot2d(xx,yy,ones(1,16),'022');
+  xstring(28,30,'K');xstring(56,30,'Plant');xstring(12,28.80,'-');
+  xtitle('PLANT   and   CONTROLLER')
+  mode(2)
+
+  path=get_absolute_file_path('mixed.dem');
+  messagebox(['Mixed Sensitivity Controller Design';
+           'file: '+path+'mixed.dem'],"modal");
+  mode(1)
+  s=poly(0,'s');
+  str='[(s-1)/((s-1)^2*(s+2))]';
+  rep=x_dialog('Nominal plant?',str)
+  if rep==[] then return,end
+  Plant=evstr(rep);Plant=syslin('c',Plant);
+   //Nominal Plant
+  P22=tf2ss(Plant);    //...in state-space form
+  [ny,nu,nx]=size(P22);
+  [Pms,r]=augment(P22);
+  txt=['W1 (sensitivity function S)';'W2 (K*S)';'W3 (complementary sensitivity T)'];
+  default=['s';'1';'1'];
+  while %t
+    rep=x_mdialog('Set weighting functions',txt,default);
+    if rep==[] then return;end
+    W1=evstr(rep(1));W2=evstr(rep(2));W3=evstr(rep(3));
+    default=rep;
+    Pms=sysdiag(W1,W2,W3,eye(P22))*Pms;
+    gms=['gamma min= ';'gamma max= ';'# iterations'];
+    vals=['0.01';'1000';'50'];
+    reps=x_mdialog('Set interval for gamma and #iterations',gms,vals);
+    if reps==[] then return;end
+    mumin=1/evstr(reps(2))^2;
+    mumax=1/evstr(reps(1))^2;
+    iter=evstr(reps(3));
+    [K,mu]=h_inf(Pms,r,mumin,mumax,iter);
+    if K~=[] then break;end
+  end
+
+  disp(spec(h_cl(Pms,r,K)),'closed loop eigenvalues')    //Check internal stability
+  [Ssens,Rsens,Tsens]=sensi(P22,K);  //Sensitivity functions
+  messagebox('Singular values plot',"modal");
+  fcts=['S function';'R (=K*S) function';'T function'];
+  www=x_mdialog('Select sensitivity function',fcts,['Yes';'No';'Yes']);
+  if www==[] then return,end
+  
+  ww1=part(www(1),1)=='Y';
+  if ww1 then
+    scf(100002);clf();show_window();
+    gainplot(Ssens);
+    xtitle('S = Sensitivity function');
+  end
+
+  ww2=part(www(2),1)=='Y';
+  if ww2 then
+    scf(100003);clf();show_window();
+    gainplot(Rsens);
+    xtitle('R (=G*S) Sensitivity function');
+  end
+
+  ww3=part(www(3),1)=='Y';
+  if ww3 then
+    scf(100004);clf();show_window();
+    gainplot(Tsens);
+    xtitle('T = Complementary Sensitivity function');
+  end
+
+endfunction
+
+demo_mixed();
+clear demo_mixed;
diff --git a/modules/cacsd/demos/pendulum/graphics.sci b/modules/cacsd/demos/pendulum/graphics.sci
new file mode 100755
index 000000000..cf0bab7ef
--- /dev/null
+++ b/modules/cacsd/demos/pendulum/graphics.sci
@@ -0,0 +1,211 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function P=initialize_display(xg,teta)
+    clf();a=gca();a.isoview="on";
+    drawlater();//f=gcf();
+    a.data_bounds=[-0.4 -0.2;0.4 0.3];
+    a.margins=zeros(1,4);
+    y1=0;lb=l;hc=0.05;lc=0.1;teta=.25;r=0.013
+    P=build_pendulum([xg,y1],[lc,hc,lb,teta,r])
+    xsegs([-0.4 0.4],[y1-2*r y1-2*r]);
+    drawnow();
+endfunction
+
+
+function [P]=dpnd()
+    //dpnd() scheme of experiment
+    //!
+    clf();
+    drawlater()
+    a=gca();
+    a.isoview="on";
+    f = gcf() ;
+    f.axes_size = [640,480];
+    a.data_bounds=[0 0;100 100];
+    a.margins(3:4)=[0 0.2];
+    xg=40;
+    y1=25;
+    lb=40;
+    hc=10;
+    lc=20;
+    teta=.25;
+    r=2.5;
+    P=build_pendulum([xg,y1],[lc,hc,lb,teta,r])
+
+    //the floor
+    xarrows([10 90],[y1-5 y1-5],0);
+    xstring(90,y1,"x")
+
+    // the force
+    yg=y1+hc/2,
+    x2=xg+lc/2;
+    xarrows([x2 x2+10],[yg yg],2);
+    xstring(x2+20,yg,"u (force)",0,0);
+
+    // the vertical
+    y2=y1+hc;
+    xsegs([xg xg],[y2 y2+lb]);e=gce();e.line_style=2;e.segs_color=-1;
+
+    // the angle teta
+    xstring(xg+lb*sin(teta)/2,y2+lb*cos(teta),"a",0,0);
+    e=gce();e.font_size=3;
+
+    //the differential equations
+    xstring(5,-9,["a'''' = (-sin(a)*cos(a)*(m/(m+M))*a''^2 + 2/(mb*l)*(sin(a)*m*g - qm*cos(a)*u))/d"
+    "x'''' =  (u+m*(l/2)*(sin(a)a''^2-cos(a)*a''''))/(m+M);"
+    "m: weight of the pendulum"
+    "M: weight of the cart"
+    "l: length of the pendulum"])
+    drawnow()
+endfunction
+
+function P=build_pendulum(o,params)
+    xg=o(1)
+    y1=o(2)
+    lc=params(1) //width of the cart
+    hc=params(2) //height of the cart
+    lb=params(3) //length of the pendulum
+    teta=params(4) //angle  of the pendulum
+    r=params(5) //radius of wires
+    y2=y1+hc;
+    x2=xg+lc/2;
+    x1=xg-lc/2;
+
+    //cart
+    xrect([xg-lc/2,y1+hc,lc,hc]);e1=gce();
+    xfarcs([x1+lc/10-r;y1;2*r;2*r;0;360*64],1);e2=gce();
+    xfarcs([x2-2*r+lc/10-r;y1;2*r;2*r;0;360*64],1);e3=gce();
+
+    //pendulum
+    xsegs([xg,xg+lb*sin(teta)],[y2,y2+lb*cos(teta)]),
+    e4=gce();e4.thickness=2;e4.segs_color=-1;
+    P=glue([e4 e3 e2 e1])
+    P.user_data=[xg,lb]
+endfunction
+
+function P=set_pendulum(P,x,theta)
+    p=P.user_data
+    xg=p(1);lb=p(2);
+    //translation
+    drawlater();
+    e=P.children(1);e.data(1)=e.data(1)+x-xg;
+    e=P.children(2).children;e.data(1)=e.data(1)+x-xg;
+    e=P.children(3).children;e.data(1)=e.data(1)+x-xg;
+    e=P.children(4);e.data(:,1)=e.data(:,1)+x-xg;
+
+    //change the pendulum angle
+    e.data(2,:)=e.data(1,:)+[lb*sin(theta) lb*cos(theta)];
+    P.user_data(1)=x
+    drawnow();
+endfunction
+
+
+function draw1()
+    f=gcf();f.figure_position=[10 10];show_window()
+    clf();
+    drawlater();
+    f.background=color("gray");
+    f.figure_size=[850,650];
+    y=y(:,1:70);  n=size(y,2);
+    a1=gca();sca(a1);
+    a1.axes_bounds=[0 0 0.5 0.5];
+    a1.data_bounds=[1,min(y(1,:));n max(y(1,:))];
+    a1.margins(1)=0.2;
+    a1.axes_visible="on";
+    a1.x_label.text="time";
+    a1.y_label.text="position";
+    a1.box = "on";
+    p1=xpoly(1,y(1,1));p1=gce();
+
+    a2=newaxes();sca(a2);
+    a2.axes_bounds=[0.5,0,0.5,0.5];
+    a2.data_bounds=[1,min(y(2,:));n max(y(2,:))];
+    a2.margins(1)=0.2;
+    a2.axes_visible="on";
+    a2.x_label.text="time";
+    a2.y_label.text="theta";
+    a2.box = "on";
+    xpoly(1,y(2,1));;p2=gce();
+
+    a3=newaxes();
+    a3.axes_bounds= [0,0.5,1,0.5];
+    a3.isoview="on";
+    a3.data_bounds=[-0.4 -0.1;0.4 0.4];
+    a3.box = "on";
+    y1=0;lb=l;hc=0.05;lc=0.1;teta=100*y(2,1);r=0.013;xg=100*y(1,1);
+    sca(a3);
+    P=build_pendulum([xg,y1],[lc,hc,lb,teta,r])
+    xsegs([-0.4 0.4],[y1-2*r y1-2*r]);
+
+    drawnow()
+
+    for k=1:size(y,2)
+        drawlater();
+        xx=100*y(1,k);tt=100*y(2,k);
+        p1.data=[p1.data;k,y(1,k)];
+        p2.data=[p2.data;k,y(2,k)];
+        // contains draw now
+        P=set_pendulum(P,xx,tt);
+    end
+endfunction
+
+function draw2()
+    f=gcf();f.figure_position=[10 10];show_window()
+    clf();
+    drawlater()
+    f.figure_size=[850,650];
+    f.background=color("gray");
+    yd=yd(:,1:100);  n=size(yd,2);
+    c = kr*yd(5:8,:) //control
+    theta = yd(3,:) //angle
+    thetaE= yd(7,:) // estimated angle
+    x = yd(1,:)
+
+    a1=gca();sca(a1);
+    a1.axes_bounds=[0 0 0.5 0.5];
+    a1.data_bounds=[min(t1),min(c);t1(n) max(c)];
+    a1.margins(1)=0.2;
+    a1.axes_visible="on";
+    a1.x_label.text="time";
+    a1.y_label.text="Control (u)";
+    a1.box = "on";
+    p1=xpoly(t1(1),c(1));p1=gce();
+
+    a2=newaxes();sca(a2);
+    a2.axes_bounds=[0.5,0,0.5,0.5];
+    a2.data_bounds=[t1(1),min([theta thetaE]);t1(n) max([theta thetaE])];
+    a2.margins(1)=0.2;
+    a2.axes_visible="on";
+    a2.x_label.text="time";
+    a2.y_label.text="theta";
+    a2.box = "on";
+    xpoly(t1(1),theta(1));p2=gce();
+    xpoly(t1(1),thetaE(1));p3=gce();p3.line_style=2;
+
+    a3=newaxes();
+    a3.axes_bounds= [0,0.5,1,0.5];
+    a3.isoview="on";
+    a3.data_bounds=[-0.4 -0.1;0.4 0.4];
+    a3.box = "on";
+    y1=0;lb=l;hc=0.05;lc=0.1;;r=0.013;
+    sca(a3);
+    P=build_pendulum([100*x(1),y1],[lc,hc,lb,100*theta(1),r])
+    xsegs([-0.4 0.4],[y1-2*r y1-2*r]);
+    drawnow()
+
+    for k=1:n
+        drawlater();
+        xx=x(k);tt=theta(k);
+        p1.data=[p1.data;t1(k),c(k)];
+        p2.data=[p2.data;t1(k),theta(k)];
+        p3.data=[p3.data;t1(k),thetaE(k)];
+        // contains draw now
+        P=set_pendulum(P,xx,tt);
+    end
+
+endfunction
diff --git a/modules/cacsd/demos/pendulum/pendule.dem b/modules/cacsd/demos/pendulum/pendule.dem
new file mode 100755
index 000000000..ddec6c606
--- /dev/null
+++ b/modules/cacsd/demos/pendulum/pendule.dem
@@ -0,0 +1,235 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function demo_pendule()
+
+    mode(-1)
+
+
+    exec("SCI/modules/cacsd/demos/pendulum/simulation.sci", -1);
+    exec("SCI/modules/cacsd/demos/pendulum/graphics.sci", -1);
+
+    //disable displayed lines control
+    display_props=lines(); lines(0)
+    my_handle             = scf(100001);
+    clf(my_handle,"reset");
+
+    show_window();
+    wdim=xget('wdim')
+    //mode(1)
+
+    dpnd()
+
+    //
+    // equations 
+    //----------
+    //state =[x x' theta theta']  
+    //
+    mb=0.1;mc=1;l=0.3;m=4*mc+mb; //constants
+    //
+    messagebox(['Open loop simulation'
+    ' '
+    ' y0 = [0;0;0.1;0]; //initial state [x x'' theta theta'']'
+    ' t = 0.03*(1:180); //observation dates'
+    ' y = ode(y0,0,0.03*(1:180),ivpd); //differential equation integration'
+    ' //Display'
+    ' P = initialize_display(y0(1),y0(3))'
+    ' for k=1:size(y,2), set_pendulum(P,y(1,k),y(3,k));end'],"modal");
+    //
+    y0=[0;0;0.1;0];
+    y=ode(y0,0,0.03*(1:180),ivpd);
+
+    P=initialize_display(y0(1),y0(3));
+    for k=1:size(y,2), set_pendulum(P,y(1,k),y(3,k));end
+    //
+    messagebox(['Linearization'
+    ' '
+    ' x0=[0;0;0;0];u0=0;'
+    ' [f,g,h,j]=lin(pendu,x0,u0);'
+    ' pe=syslin(''c'',f,g,h,j);'
+    ' // display of the linear system'
+    ' ssprint(pe)'],"modal")
+
+    //
+    mode(1)
+    x0=[0;0;0;0];u0=0;
+    [f,g,h,j]=lin(pendu,x0,u0);
+    pe=syslin('c',f,g,h,j);
+    ssprint(pe)
+    mode(-1)
+
+    //
+    messagebox(['Checking the result'
+    ' //comparison with linear model computed by hand';
+    ''
+    ' f1=[0 1        0                0'
+    '     0 0    -3*mb*9.81/m         0'
+    '     0 0        0                1'
+    '     0 0  6*(mb+mc)*9.81/(m*l)   0];'
+    ' '
+    ' g1=[0 ; 4/m ; 0 ; -6/(m*l)];'
+    ' '
+    ' h1=[1 0 0 0'
+    '     0 0 1 0];'
+    ' '
+    ' err=norm(f-f1,1)+norm(g-g1,1)+norm(h-h1,1)+norm(j,1)'],"modal")
+
+
+    //
+    mode(1)
+    f1=[0 1        0             0
+    0 0    -3*mb*9.81/m         0
+    0 0        0             1
+    0 0  6*(mb+mc)*9.81/(m*l)   0];
+    g1=[0 ; 4/m ; 0 ; -6/(m*l)];
+    h1=[1 0 0 0
+    0 0 1 0];
+    err=norm(f-f1,1)+norm(g-g1,1)+norm(h-h1,1)+norm(j,1)
+    mode(-1)
+
+    messagebox(['Linear system properties analysis'
+    ' spec(f) //stability (unstable system !)'
+    ' n=contr(f,g) //controlability'
+    ' '
+    ' //observability'
+    ' m1=contr(f'',h(1,:)'') '
+    ' [m2,t]=contr(f'',h(2,:)'')'],"modal");
+    //---------
+    mode(1)
+    spec(f) //stability (unstable system !)
+    n=contr(f,g) //controlability
+
+    //observability
+    m1=contr(f',h(1,:)') 
+    [m2,t]=contr(f',h(2,:)')
+    mode(-1)
+    //
+    messagebox(['Synthesis of a stabilizing controller using poles placement'
+    ' '
+    '// only x and theta are observed  : contruction of an observer'
+    '// to estimate the state : z''=(f-k*h)*z+k*y+g*u'
+    ' to=0.1; ' 
+    ' k=ppol(f'',h'',-ones(4,1)/to)''  //observer gain'
+    '// compute the compensator gain'
+    ' kr=ppol(f,g,-ones(4,1)/to)'],"modal");
+
+    //-------------------------------------------------
+    //
+    //pole placement technique
+    //only x and theta are observed  : contruction of an observer
+    //to estimate the state : z'=(f-k*h)*z+k*y+g*u
+    //
+    to=0.1;  //
+    k=ppol(f',h',-ones(4,1)/to)'  //observer gain
+    //
+    //verification
+    //
+    // norm( poly(f-k*h,'z')-poly(-ones(4,1)/to,'z'))
+    //
+    kr=ppol(f,g,-ones(4,1)/to)  //compensator gain
+
+    //
+    messagebox(['Build full linear system  pendulum-observer-compensator'
+    ' '
+    ' ft=[f-g*kr            -g*kr'
+    '     0*f               f-k*h]'
+    ' '
+    ' gt=[g;0*g];'
+    ' ht=[h,0*h];'
+    ''
+    ' pr=syslin(''c'',ft,gt,ht);'
+    ''
+    '//Check the closed loop dynamic'
+    ' spec(pr.A)'
+    ' '
+    '//transfer matrix representation'
+    ' hr=clean(ss2tf(pr),1.d-10)'
+    ' '
+    '//frequency analysis'
+    ' black(pr,0.01,100,[''position'',''theta''])'
+    ' g_margin(pr(1,1))'],"modal");
+
+    //---------------------------------------------
+    //
+    //state: [x x-z]
+    //
+    mode(1)
+    ft=[f-g*kr            -g*kr
+    0*f               f-k*h];
+    gt=[g;0*g];
+    ht=[h,0*h];
+    pr=syslin('c',ft,gt,ht);
+
+    // closed loop dynamics:
+    spec(pr(2))
+    //transfer matrix representation
+    hr=clean(ss2tf(pr),1.d-10)
+
+
+    //frequency analysis
+    clf()
+    black(pr,0.01,100,['position','theta'])
+    g_margin(pr(1,1))
+    mode(-1)
+
+
+    //
+    messagebox(['Sampled system'
+    ' '
+    ' t=to/5; //sampling period'
+    ' prd=dscr(pr,t); //discrete model'
+    ' spec(prd.A) //poles of the discrete model'],"modal");
+
+    //---------------
+    //
+    mode(1)
+    t=to/5;
+    prd=dscr(pr,t);
+    spec(prd(2))
+    mode(-1)
+    //
+    messagebox(['Impulse response'
+    ' '
+    ' x0=[0;0;0;0;0;0;0;0]; //initial state'
+    ' u(1,180)=0;u(1,1)=1; //impulse'
+    ' y=flts(u,prd,x0);    //siscrete system simulation'
+    ' draw1()'],"modal");
+
+    //-----------------
+    //
+    mode(1)
+    x0=[0;0;0;0;0;0;0;0];
+    u(1,180)=0;u(1,1)=1;
+    y=flts(u,prd,x0);
+
+    draw1()
+    mode(-1)
+    //
+    messagebox(['Compensation of the non linear system with'
+    'linear regulator'
+    ''
+    ' t0=0;t1=t*(1:125);'
+    ' x0=[0 0 0.4 0   0 0 0 0]'';  // initial state'
+    ' yd=ode(x0,t0,t1,regu); //integration of differential equation'
+    ' draw2()'],"modal");
+    ;
+    //--------------------------------------
+    //
+    //simulation
+    //
+    mode(1)
+    t0=0;t1=t*(1:125);
+    x0=[0 0 0.4 0   0 0 0 0]';   //
+    yd=ode(x0,t0,t1,regu);
+    draw2()
+    mode(-1) 
+    lines(display_props(2))
+    messagebox('The end',"modal");
+endfunction
+
+demo_pendule();
+clear demo_pendule;
\ No newline at end of file
diff --git a/modules/cacsd/demos/pendulum/simulation.sci b/modules/cacsd/demos/pendulum/simulation.sci
new file mode 100755
index 000000000..c4b07e016
--- /dev/null
+++ b/modules/cacsd/demos/pendulum/simulation.sci
@@ -0,0 +1,61 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+
+
+function [xdot]=ivpd(t,x)
+    //ydot=ivpd(t,y) non linear equations of the pendulum
+    // y=[x;d(x)/dt,teta,d(teta)/dt].
+    // mb, mc, l must be predefined
+    //!
+    g=9.81;
+    u=0
+    qm=mb/(mb+mc)
+    cx3=cos(x(3))
+    sx3=sin(x(3))
+    d=4/3-qm*cx3*cx3
+    xd4=(-sx3*cx3*qm*x(4)**2+2/(mb*l)*(sx3*mb*g-qm*cx3*u))/d
+    //
+    xdot=[x(2);
+    (u+mb*(l/2)*(sx3*x(4)**2-cx3*xd4))/(mb+mc);
+    x(4);
+    xd4]
+endfunction
+
+
+function [y,xdot]=pendu(x,u)
+    //[y,xdot]=pendu(x,u) input (u) - output (y) function of the pendulum
+    //!
+    g=9.81;
+    qm=mb/(mb+mc)
+    cx3=cos(x(3))
+    sx3=sin(x(3))
+    d=4/3-qm*cx3*cx3
+    xd4=(-sx3*cx3*qm*x(4)**2+2/(mb*l)*(sx3*mb*g-qm*cx3*u))/d
+    //
+    xdot=[x(2);
+    (u+mb*(l/2)*(sx3*x(4)**2-cx3*xd4))/(mb+mc);
+    x(4);
+    xd4]
+    //
+    y=[x(1);
+    x(3)]
+endfunction
+
+
+function [xdot]=regu(t,x)
+    //xdot=regu(t,x) simulation of the compound system
+    // pendulum - observer - compensator
+    //!
+    xp=x(1:4);xo=x(5:8);
+    u =-kr*xo   //control
+    [y,xpd]=pendu(xp,u)
+    xod=(f-k*h)*xo+k*y+g*u
+    xdot=[xpd;xod]
+endfunction
+
+
diff --git a/modules/cacsd/demos/pendulum/yy b/modules/cacsd/demos/pendulum/yy
new file mode 100755
index 000000000..2b06014de
--- /dev/null
+++ b/modules/cacsd/demos/pendulum/yy
@@ -0,0 +1,720 @@
+-.321848E-04
+-.130207E-03
+-.298512E-03
+-.544646E-03
+-.879417E-03
+-.131703E-02
+-.187503E-02
+-.257380E-02
+-.343499E-02
+-.447806E-02
+-.571329E-02
+-.712915E-02
+-.867168E-02
+-.102150E-01
+-.115268E-01
+-.122434E-01
+-.118836E-01
+-.994784E-02
+-.615677E-02
+-.808250E-03
+0.504670E-02
+0.100574E-01
+0.133576E-01
+0.148354E-01
+0.148603E-01
+0.139438E-01
+0.125411E-01
+0.109804E-01
+0.946265E-02
+0.809142E-02
+0.690698E-02
+0.591346E-02
+0.509721E-02
+0.443769E-02
+0.391335E-02
+0.350446E-02
+0.319432E-02
+0.296965E-02
+0.282050E-02
+0.274016E-02
+0.272499E-02
+0.277429E-02
+0.289031E-02
+0.307829E-02
+0.334663E-02
+0.370702E-02
+0.417456E-02
+0.476764E-02
+0.550717E-02
+0.641469E-02
+0.750808E-02
+0.879339E-02
+0.102503E-01
+0.118090E-01
+0.133189E-01
+0.145145E-01
+0.149973E-01
+0.142667E-01
+0.118515E-01
+0.758906E-02
+0.197866E-02
+-.377710E-02
+-.838836E-02
+-.111920E-01
+-.122396E-01
+-.119665E-01
+-.108822E-01
+-.941152E-02
+-.784722E-02
+-.636121E-02
+-.503724E-02
+-.390376E-02
+-.295876E-02
+-.218597E-02
+-.156417E-02
+-.107209E-02
+-.690739E-03
+-.404301E-03
+-.200379E-03
+-.699076E-04
+-.699544E-05
+-.877912E-05
+-.753402E-04
+-.209708E-03
+-.417945E-03
+-.709298E-03
+-.109636E-02
+-.159512E-02
+-.222474E-02
+-.300656E-02
+-.396169E-02
+-.510590E-02
+-.644000E-02
+-.793322E-02
+-.949780E-02
+-.109559E-01
+-.120074E-01
+-.122207E-01
+-.110849E-01
+-.817674E-02
+-.347690E-02
+0.230785E-02
+0.786999E-02
+0.120354E-01
+0.143481E-01
+0.149976E-01
+0.144626E-01
+0.132403E-01
+0.117219E-01
+0.101658E-01
+0.871723E-02
+0.744240E-02
+0.635967E-02
+0.546203E-02
+0.473123E-02
+0.414567E-02
+0.368454E-02
+0.332964E-02
+0.306606E-02
+0.288227E-02
+0.277009E-02
+0.272441E-02
+0.274317E-02
+0.282722E-02
+0.298037E-02
+0.320952E-02
+0.352477E-02
+0.393963E-02
+0.447095E-02
+0.513862E-02
+0.596424E-02
+0.696813E-02
+0.816332E-02
+0.954421E-02
+0.110679E-01
+0.126261E-01
+0.140116E-01
+0.148887E-01
+0.147962E-01
+0.132240E-01
+0.981935E-02
+0.473181E-02
+-.113146E-02
+-.641460E-02
+-.101033E-01
+-.119405E-01
+-.122266E-01
+-.114650E-01
+-.101321E-01
+-.858384E-02
+-.704586E-02
+-.563915E-02
+-.441458E-02
+-.338202E-02
+-.253042E-02
+-.184005E-02
+-.128926E-02
+-.857801E-03
+-.528301E-03
+-.286746E-03
+-.122516E-03
+-.282410E-04
+0.356796E-06
+-.354184E-04
+-.137198E-03
+-.309597E-03
+-.560336E-03
+-.900411E-03
+-.134421E-02
+-.190947E-02
+-.261669E-02
+-.348751E-02
+-.454111E-02
+-.578695E-02
+-.721177E-02
+-.875847E-02
+-.102961E-01
+-.115857E-01
+-.122560E-01
+-.118211E-01
+-.215378E-02
+-.440565E-02
+-.685654E-02
+-.961237E-02
+-.127846E-01
+-.164872E-01
+-.208269E-01
+-.258797E-01
+-.316449E-01
+-.379573E-01
+-.443387E-01
+-.497688E-01
+-.523874E-01
+-.492245E-01
+-.362035E-01
+-.884836E-02
+0.357399E-01
+0.950836E-01
+0.156152E+00
+0.194303E+00
+0.187985E+00
+0.141109E+00
+0.784839E-01
+0.222909E-01
+-.176819E-01
+-.408561E-01
+-.508338E-01
+-.521148E-01
+-.485144E-01
+-.426925E-01
+-.362573E-01
+-.300591E-01
+-.244747E-01
+-.196143E-01
+-.154520E-01
+-.119003E-01
+-.884911E-02
+-.618440E-02
+-.379612E-02
+-.157980E-02
+0.564484E-03
+0.273453E-02
+0.502893E-02
+0.754974E-02
+0.104044E-01
+0.137056E-01
+0.175666E-01
+0.220891E-01
+0.273348E-01
+0.332700E-01
+0.396645E-01
+0.459227E-01
+0.508331E-01
+0.522638E-01
+0.469325E-01
+0.305473E-01
+-.120966E-02
+-.502702E-01
+-.111873E+00
+-.169589E+00
+-.197356E+00
+-.178716E+00
+-.125006E+00
+-.624139E-01
+-.100258E-01
+0.253216E-01
+0.445845E-01
+0.518415E-01
+0.515092E-01
+0.471118E-01
+0.410085E-01
+0.345771E-01
+0.285182E-01
+0.231214E-01
+0.184513E-01
+0.144600E-01
+0.110515E-01
+0.811319E-02
+0.553187E-02
+0.319888E-02
+0.101113E-02
+-.113049E-02
+-.332369E-02
+-.566774E-02
+-.826601E-02
+-.112274E-01
+-.146655E-01
+-.186922E-01
+-.234023E-01
+-.288390E-01
+-.349290E-01
+-.413655E-01
+-.474178E-01
+-.516615E-01
+-.516749E-01
+-.438709E-01
+-.238074E-01
+0.125104E-01
+0.657443E-01
+0.128464E+00
+0.180896E+00
+0.196994E+00
+0.166917E+00
+0.108336E+00
+0.471234E-01
+-.101354E-02
+-.318270E-01
+-.474765E-01
+-.523259E-01
+-.506266E-01
+-.455950E-01
+-.393044E-01
+-.329242E-01
+-.270238E-01
+-.218187E-01
+-.173351E-01
+-.135081E-01
+-.102348E-01
+-.740157E-02
+-.489611E-02
+-.261126E-02
+-.445210E-03
+0.170051E-02
+0.392375E-02
+0.632463E-02
+0.900792E-02
+0.120840E-01
+0.156669E-01
+0.198660E-01
+0.247669E-01
+0.303901E-01
+0.366147E-01
+0.430437E-01
+0.487928E-01
+0.522023E-01
+0.505474E-01
+0.399548E-01
+0.159200E-01
+-.250309E-01
+-.819504E-01
+-.144388E+00
+-.189576E+00
+-.193234E+00
+-.153095E+00
+-.915548E-01
+-.328096E-01
+0.108124E-01
+0.372657E-01
+0.496181E-01
+0.523609E-01
+0.495187E-01
+0.439960E-01
+0.375979E-01
+0.313072E-01
+0.255794E-01
+0.205672E-01
+0.162655E-01
+0.125953E-01
+0.944926E-02
+0.671330E-02
+0.427624E-02
+0.203246E-02
+-.118721E-03
+-.227533E-02
+-.453554E-02
+-.700051E-02
+-.977647E-02
+-.129752E-01
+-.167105E-01
+-.210882E-01
+-.261816E-01
+-.319836E-01
+-.383162E-01
+-.446778E-01
+-.500098E-01
+-.523964E-01
+-.488018E-01
+-.350996E-01
+-.683490E-02
+0.387109E-01
+0.102372E+00
+0.109601E+00
+0.122027E+00
+0.140233E+00
+0.165073E+00
+0.197701E+00
+0.239622E+00
+0.292742E+00
+0.359422E+00
+0.442528E+00
+0.545457E+00
+0.672110E+00
+0.826784E+00
+0.101396E+01
+0.123799E+01
+0.150271E+01
+0.181096E+01
+0.216366E+01
+0.255764E+01
+0.298181E+01
+0.341525E+01
+0.383311E+01
+0.421687E+01
+0.455792E+01
+0.485449E+01
+0.510814E+01
+0.532204E+01
+0.550020E+01
+0.564705E+01
+0.576703E+01
+0.586437E+01
+0.594283E+01
+0.600568E+01
+0.605563E+01
+0.609493E+01
+0.612537E+01
+0.614836E+01
+0.616497E+01
+0.617598E+01
+0.618190E+01
+0.618302E+01
+0.617939E+01
+0.617083E+01
+0.615695E+01
+0.613708E+01
+0.611031E+01
+0.607540E+01
+0.603072E+01
+0.597427E+01
+0.590356E+01
+0.581557E+01
+0.570678E+01
+0.557317E+01
+0.541035E+01
+0.521386E+01
+0.497943E+01
+0.470345E+01
+0.438344E+01
+0.401929E+01
+0.361587E+01
+0.318687E+01
+0.275522E+01
+0.234499E+01
+0.197202E+01
+0.164271E+01
+0.135771E+01
+0.111490E+01
+0.910841E+00
+0.741378E+00
+0.602047E+00
+0.488430E+00
+0.396420E+00
+0.322373E+00
+0.263175E+00
+0.216230E+00
+0.179425E+00
+0.151072E+00
+0.129857E+00
+0.114791E+00
+0.105165E+00
+0.100527E+00
+0.100659E+00
+0.105566E+00
+0.115481E+00
+0.130868E+00
+0.152451E+00
+0.181236E+00
+0.218558E+00
+0.266124E+00
+0.326076E+00
+0.401034E+00
+0.494146E+00
+0.609082E+00
+0.749971E+00
+0.921242E+00
+0.112735E+01
+0.137243E+01
+0.165986E+01
+0.199164E+01
+0.236693E+01
+0.277887E+01
+0.321104E+01
+0.363915E+01
+0.404066E+01
+0.440242E+01
+0.471994E+01
+0.499354E+01
+0.522575E+01
+0.542025E+01
+0.558132E+01
+0.571345E+01
+0.582098E+01
+0.590792E+01
+0.597776E+01
+0.603350E+01
+0.607758E+01
+0.611200E+01
+0.613836E+01
+0.615787E+01
+0.617144E+01
+0.617972E+01
+0.618308E+01
+0.618170E+01
+0.617550E+01
+0.616420E+01
+0.614726E+01
+0.612389E+01
+0.609299E+01
+0.605315E+01
+0.600254E+01
+0.593890E+01
+0.585948E+01
+0.576099E+01
+0.563962E+01
+0.549115E+01
+0.531111E+01
+0.509510E+01
+0.483913E+01
+0.454011E+01
+0.419659E+01
+0.381062E+01
+0.339132E+01
+0.295772E+01
+0.253476E+01
+0.214289E+01
+0.179263E+01
+0.148685E+01
+0.122448E+01
+0.100261E+01
+0.817362E+00
+0.664365E+00
+0.539142E+00
+0.437415E+00
+0.355307E+00
+0.289453E+00
+0.237014E+00
+0.195656E+00
+0.163497E+00
+0.139054E+00
+0.121189E+00
+0.109065E+00
+0.102113E+00
+0.100006E+00
+0.102643E+00
+0.110150E+00
+0.122879E+00
+0.141429E+00
+0.166667E+00
+0.199768E+00
+0.242258E+00
+0.296064E+00
+0.363575E+00
+0.447688E+00
+0.551828E+00
+0.679921E+00
+0.836282E+00
+0.102539E+01
+0.125159E+01
+0.151866E+01
+0.182938E+01
+0.158775E+00
+0.325031E+00
+0.506580E+00
+0.711896E+00
+0.950458E+00
+0.123308E+01
+0.157223E+01
+0.198218E+01
+0.247904E+01
+0.308025E+01
+0.380342E+01
+0.466411E+01
+0.567261E+01
+0.683025E+01
+0.812660E+01
+0.953753E+01
+0.110187E+02
+0.124795E+02
+0.137257E+02
+0.144327E+02
+0.143184E+02
+0.134341E+02
+0.121043E+02
+0.106247E+02
+0.915596E+01
+0.777203E+01
+0.651055E+01
+0.539178E+01
+0.442285E+01
+0.359975E+01
+0.291042E+01
+0.233848E+01
+0.186621E+01
+0.147645E+01
+0.115353E+01
+0.883678E+00
+0.654883E+00
+0.456715E+00
+0.279994E+00
+0.116454E+00
+-.415991E-01
+-.201613E+00
+-.371123E+00
+-.558084E+00
+-.771202E+00
+-.102028E+01
+-.131653E+01
+-.167290E+01
+-.210415E+01
+-.262682E+01
+-.325853E+01
+-.401663E+01
+-.491562E+01
+-.596379E+01
+-.715957E+01
+-.848912E+01
+-.992408E+01
+-.114117E+02
+-.128404E+02
+-.139795E+02
+-.144876E+02
+-.141488E+02
+-.131120E+02
+-.117197E+02
+-.102321E+02
+-.878060E+01
+-.742620E+01
+-.620097E+01
+-.512154E+01
+-.419188E+01
+-.340548E+01
+-.274882E+01
+-.220491E+01
+-.175603E+01
+-.138534E+01
+-.107769E+01
+-.819769E+00
+-.600011E+00
+-.408350E+00
+-.235880E+00
+-.745176E-01
+0.833323E-01
+0.245110E+00
+0.418429E+00
+0.611410E+00
+0.833014E+00
+0.109339E+01
+0.140418E+01
+0.177879E+01
+0.223250E+01
+0.278222E+01
+0.344567E+01
+0.423973E+01
+0.517764E+01
+0.626539E+01
+0.749836E+01
+0.885918E+01
+0.103147E+02
+0.118013E+02
+0.131821E+02
+0.141889E+02
+0.144811E+02
+0.139295E+02
+0.127658E+02
+0.113291E+02
+0.984222E+01
+0.841206E+01
+0.708936E+01
+0.590154E+01
+0.486174E+01
+0.397087E+01
+0.322023E+01
+0.259504E+01
+0.207792E+01
+0.165125E+01
+0.129860E+01
+0.100529E+01
+0.758498E+00
+0.547081E+00
+0.361310E+00
+0.192531E+00
+0.328219E-01
+-.125340E+00
+-.289409E+00
+-.467100E+00
+-.666721E+00
+-.897515E+00
+-.116999E+01
+-.149625E+01
+-.189018E+01
+-.236753E+01
+-.294554E+01
+-.364191E+01
+-.447287E+01
+-.545012E+01
+-.657714E+01
+-.784610E+01
+-.923594E+01
+-.107078E+02
+-.121844E+02
+-.134984E+02
+-.143473E+02
+-.144134E+02
+-.136671E+02
+-.124014E+02
+-.109356E+02
+-.945661E+01
+-.805114E+01
+-.676202E+01
+-.561252E+01
+-.461238E+01
+-.375968E+01
+-.304374E+01
+-.244881E+01
+-.195723E+01
+-.155163E+01
+-.121599E+01
+-.936132E+00
+-.699692E+00
+-.495938E+00
+-.315458E+00
+-.149820E+00
+0.875745E-02
+0.167748E+00
+0.334642E+00
+0.517279E+00
+0.724179E+00
+0.964889E+00
+0.125031E+01
+0.159299E+01
+0.200733E+01
+0.250952E+01
+0.311705E+01
+0.384747E+01
+0.471616E+01
+0.573300E+01
+0.689874E+01
+0.820222E+01
+0.961847E+01
+0.111016E+02
diff --git a/modules/cacsd/demos/pid.dem b/modules/cacsd/demos/pid.dem
new file mode 100755
index 000000000..7e99943a2
--- /dev/null
+++ b/modules/cacsd/demos/pid.dem
@@ -0,0 +1,167 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+function demo_pid()
+
+  mode(-1);
+  lines(0);
+  //display the diagram
+  x=[5,10,20,40,50,70,80,90];xmin=-10;xmax=100;
+  y=[22,28,30,32];ymin=12;ymax=40;
+
+  xx=[xmin,xmin,x([1 2 2 7 4 6 3 4 5 6 3 3 5 5]);xmax,xmax,x([3,2,7,7,5,8,3,4,5,6,4,4,6,6])];
+  yy=[ymin,ymax,y([3,1,1,1,3,3,2,2,2,2,2,4,2,4]);ymin,ymax,y([3,3,1,3,3,3,4,4,4,4,2,4,2,4])];
+
+  scf(100001);clf();show_window();
+  plot2d(xx,yy,ones(1,16),'022');
+  xstring(28,30,'K');xstring(56,30,'Plant');xstring(12,28.80,'-');
+  xtitle('PLANT   and   CONTROLLER')
+  mode(2);
+
+  path=get_absolute_file_path('pid.dem');
+  s=poly(0,'s');z=poly(0,'z');
+  messagebox(['Example of PID Design '
+              'file: '+path+'pid.dem'],"modal");
+
+  n=x_choose(['Continuous time';'Discrete time'],'Select time domain');
+  select n
+   case 0
+    warning('Demo stops!');return;
+   case 1
+    mode(1)
+    dom='c';
+    s=poly(0,'s');
+    str='[(s-1)/(s^2+5*s+1)]';
+    rep=x_dialog('Nominal plant?',str);
+    if rep==[] then return,end
+    Plant=evstr(rep);
+    Plant=syslin('c',Plant);
+    mode(-1)
+  case 2
+    mode(1)
+    dom='d'
+    z=poly(0,'z');
+    str='(z+1)/(z^2-5*z+2)'
+    rep=x_dialog('Nominal plant?',str);
+    if rep==[] then return,end
+    Plant=evstr(rep)
+    Plant=syslin('d',Plant);
+    mode(-1)
+  end
+  //Nominal Plant
+  P22=tf2ss(Plant);    //...in state-space form
+  [ny,nu,nx]=size(P22);
+  defv=['-1.2';'1';'0.1'];
+  if dom=='d' then defv=['-10';'1';'0.1'];end
+  while %t
+    mode(1)
+    if dom=='c' then
+      _title='Enter your PID controller K(s)=Kp*(1+T0/s+T1*s)';
+    end
+    if dom=='d' then
+      _title='Enter your PID controller K(z)=Kp*(1+T0/z+T1*z)';
+    end
+    defv=x_mdialog(_title,['Kp=';'T0=';'T1='],defv);
+    if defv==[] then warning('Demo stops!');return;end
+    Kp=evstr(defv(1));T0=evstr(defv(2));T1=evstr(defv(3));
+    if dom=='c' then
+      Kpid=tf2ss(Kp*(1+T0/s+T1*s));
+    end
+    if dom=='d' then
+      Kpid=tf2ss(Kp*(1+T0/z+T1*z));
+    end
+    W=[1, -P22;
+        Kpid,1];Winv=inv(W);
+
+    disp(spec(Winv(2)),'closed loop eigenvalues');//Check internal stability
+    if max(real(spec(Winv(2)))) > 0 then
+      messagebox('You lose: closed-loop is UNSTABLE!!!',"modal");
+    else
+      messagebox('Congratulations: closed-loop is STABLE !!!',"modal");
+      break;
+    end
+    mode(-1)
+  end
+  mode(1)
+  [Spid,Rpid,Tpid]=sensi(P22,Kpid);  //Sensitivity functions
+  Tpid(5)=clean(Tpid(5));
+
+  disp(clean(ss2tf(Spid)),'Sensitivity function');
+  disp(clean(ss2tf(Tpid)),'Complementary sensitivity function');
+
+  resp=['Frequency response';'Time response'];
+  while %t do
+    n=x_choose(resp,'Select response(s)');
+    if degree(Tpid(5))>0 then
+      warning('Improper transfer function! D(s) set to D(0)')
+      Tpid(5)=coeff(Tpid(5),0);
+    end
+    Tpid(5)=coeff(Tpid(5));
+    select n
+    case 0
+      break
+    case 1
+      mode(1)
+      scf(100002);clf(100002);show_window();bode(Tpid);
+      mode(-1)
+    case 2
+      if Plant(4)=='c' then
+        mode(1)
+        defv=['0.1';'50'];
+        _title='Enter Sampling period and Tmax';
+        rep=x_mdialog(_title,['Sampling period?';'Tmax?'],defv);
+        if rep==[] then break,end
+        dttmax=evstr(rep);
+        dt=evstr(dttmax(1));tmax=evstr(dttmax(2));
+        t=0:dt/5:tmax;
+        n1=x_choose(['Step response?';'Impulse response?'],'Simulation:');
+        if n1==0 then
+          warning('Demo stops!');return;
+        end
+        if n1==1 then
+          scf(100002);clf(100002);show_window();
+          plot2d([t',t'],[(csim('step',t,Tpid))',ones(t')])
+        end
+        if n1==2 then
+          scf(100002);clf(100002);show_window();
+          plot2d([t',t'],[(csim('impul',t,Tpid))',0*t'])
+        end
+        mode(-1)
+      elseif Plant(4)=='d' then
+        mode(1)
+        defv=['100'];
+        _title='Tmax?'
+        rep=x_mdialog(_title,['Tmax='],defv);
+        if rep==[] then break,end
+        Tmax=evstr(rep);
+        mode(-1)
+        while %t do
+          n=x_choose(['Step response?';'Impulse response?'],'Simulation:');
+          select n
+          case 0 then
+            break
+          case 1 then
+            mode(1)
+            u=ones(1,Tmax);u(1)=0;
+            scf(100002);clf(100002);show_window();
+            plot2d([(1:Tmax)',(1:Tmax)'],[(dsimul(Tpid,u))',(ones(1:Tmax)')])
+            mode(-1)
+          case 2 then
+            mode(1)
+            u=zeros(1,Tmax);u(1)=1;
+            scf(100002);clf(100002);show_window();
+            plot2d((1:Tmax)',(dsimul(Tpid,u))')
+            mode(-1)
+          end
+        end
+      end
+    end
+  end
+endfunction
+
+demo_pid();
+clear demo_pid;
diff --git a/modules/cacsd/demos/robust/como.dem b/modules/cacsd/demos/robust/como.dem
new file mode 100755
index 000000000..cd83d96eb
--- /dev/null
+++ b/modules/cacsd/demos/robust/como.dem
@@ -0,0 +1,471 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+// Examples of Standard plants (mostly from Kwakernak)
+// Define the polynomial s by typing "s=poly(0,'s')"
+// Most examples are in transfer representation but computations
+// are done in state space form
+
+//2.4.1 (Kwakernaak)
+
+H=[1 (s-1)/((s-2)*(s-3));
+   1 (s-1)/((s-2)*(s-3))];
+r=[1 1];
+
+//Transform to get d12 full rank
+
+H(:,2)=H(:,2)*(s+1);
+
+[sk,mu]=h_inf(H,r,0,1,20);
+
+//transform back
+
+sk=sk*(s+1);
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//3.9.1 
+c=0.1;r=0;
+
+H=[(1+sqrt(2)*s+s**2)/(s*s) 1/(s*s);
+0 c+c*r*s;
+(1+sqrt(2)*s+s**2)/(s*s) 1/(s*s)];
+
+r=[1 1];
+
+//[sk,mu]=h_inf(H,r,0,1,20);
+
+//compare (185) in como notes
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+c=0.1;r=1;
+
+H=[(1+sqrt(2)*s+s**2)/(s*s) 1/(s*s);...
+   0 c+c*r*s;..
+   (1+sqrt(2)*s+s**2)/(s*s) 1/(s*s)];
+
+r=[1 1];
+
+// divide 2nd column of h by (s+2) for properness
+H(:,2)=H(:,2)/(s+2);
+
+//[sk,mu]=h_inf(H,r,0,1,30);
+
+// Remove parasitic root introduced above :
+//vv=roots(sk(2))   //vv(2) = -2 approximately 
+
+//sk=sk/real((s-vv(2)))
+
+//compare (186) como notes
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//3.12.1
+
+c=0.1;
+H=[(s+1)/s 1/s;
+  0  c;
+ (s+1)/s 1/s];
+
+r=[1 1];
+
+//[sk,mu]=h_inf(H,r,0,1,20);
+
+//sk=clean(sk,0,0.00001);
+
+//compare (228)
+
+[sk,mu]=h_inf(H,r,0,10,20);
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//4.12.1
+
+H=[(s*s+sqrt(2)*s+1)/(s*s) 1/(s*s);
+ 0 1;
+(s*s+sqrt(2)*s+1)/(s*s) 1/(s*s)];
+
+r=[1 1];
+
+//[sk,mu]=h_inf(H,r,0,10,20);
+
+//compare (422) which is misprinted !
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//1990 IFAC (Kwakernaak)
+
+H=[-(1+s/4)/(1+2*s) (s+1)/(s-2) 0 1/(s-2) 0;
+   0 0 0 1 0;
+   0 0 0 0 1;
+(1+s/4)/(1+2*s) 0 0 0 0;
+ 0 (s+1)/(s-2) 1 1/(s-2) 0];
+
+r=[2 2];
+
+//[Sk,mu]=h_inf(H,r,0,1,30);
+
+//compare (33) in Matlab report 
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//Mc-Farlane Glover for 1/s^2
+
+G=[(1+s*sqrt(2)+s*s)/(s*s) 1/(s*s);
+   0 1;
+   (1+s*sqrT(2)+s*s)/(s*s) 1/(s*s)]
+
+r=[1,1]
+
+//[sk,mu]=h_inf(G,r,0,1,20)
+
+//compare (25) in Matlab report
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//Mc Farlane Glover for 1/s^3
+
+
+
+G=[(1+2*s+2*s*s+s**3)/(s**3) 1/(s**3);
+  0 1;
+  (1+2*s+2*s*s+s**3)/(s**3) 1/(s**3)]
+
+r=[1 1];
+//[sk,mu]=h_inf(G,r,0,1,20)
+
+//compare (27) in Matlab report
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//SISO Model matching
+
+G=[6*(s-1)/((s-2)*(s-3)) 1;
+   -1 0];
+r=[1 1];
+
+//multiply first row by all-pass transfer to have stabilizability
+//and detectability
+
+k=((s-2)*(s-3))/((s+2)*(s+3))
+
+G(1,:)=k*G(1,:)
+
+[sk,mu]=h_inf(G,r,0,5,20)
+
+//compare (42) in Matlab report
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//Mimo Model matching
+
+G=[0.5/(s-1) 0 1 0;
+   1/(s*s-s+1) 2/(s-1) 0 1;
+    -1 0 0 0;
+     0 -1 0 0];
+
+r=[2 2]
+
+p=(s*s-s+1)*(s-1);
+
+q=p/horner(p,-s);
+
+H=[q 0;
+   0 q]
+
+//multiplication from the left by inner H
+
+G(1:2,:)=H*G(1:2,:);
+
+//Stabilizability and detectability assumptions are now verified
+
+[sk,mu]=h_inf(G,r,0,1,30);
+
+// the parasitic (?) root -83.60 (Kwakernaak) or -84.53 (Francis) is not
+// in the poles of sk :
+// one finds the roots -0.6920857 + or -0.7932064i and -1.1030557
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//Mixed sensitivity (ex 11 in matlab notes)
+
+G=[(s+1)/s 0 1/s 0;
+   0 (s+0.5)/(s-1) 0 1/(s-1);
+   0 0 1 0;
+   0 0 0 1;
+(s+1)/s 0 1/s 0;
+0 (s+0.5)/(s-1) 0 1/(s-1)];
+
+r=[2 2]
+
+[sk,mu]=h_inf(G,r,0,1,30);
+
+sk=clean(sk,1.d-6)
+
+//compare (40) in Matlab notes
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//Laub's ex (CDC 90)
+
+A=[0.0904 0.8465;
+   0.6888 0.7152];
+
+B=[0.7092 0.3017 0.7001;
+   0.1814 0.9525 0.1593];
+
+C=[0.3088 0.7350;
+   0.5735 0.9820;
+   0.6644 0.5627];
+
+D=[0.7357 0.4156 0.0;
+   0.2588 0.1544 1.0;
+    0.0    1.0   0.0];
+
+r=[1 1];
+
+p=syslin('c',a,b,c,d);
+
+//[sk,mu]=h_inf(P,r,0,10,20);
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//Mustafa's ex
+
+A=[20 -100;1 0];
+B=[1;0];
+C=[1,-0.1];
+Gss=syslin('c',A,B,C);
+[Pss,r]=normlqg(Gss);    //The standard plant
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//Maciejowski's example
+
+At=[-0.001 0 1.1320 0 -1;
+ 0 -0.0538 -0.1712 0 0.0705;
+ 0 0 0 1 0;
+ 0 0.0485 0 -0.8556 -1.013;
+ 0 -0.2909 0 1.0532 -0.6859];
+
+B=[0 0 0;
+   -0.120 1 0;
+ 0  0 0;
+4.4190 0 -1.665;
+1.575 0 -0.0732];
+
+C=[1 0 0 0 0;
+   0 1 0 0 0;
+   0 0 1 0 0];
+
+Dt=1.d-5*eye(3,3);
+D=0*Dt;
+
+Gss=syslin('c',At,b,c,D);
+
+//gtf=ss2tf(Gss);gtf=clean(gtf,1.d-10);
+
+w1=syslin('c',((s+6)**2)/((s+0.00006)*(s+0.6)));
+
+w2=syslin('c',2000*(s+10)*(s+50)/((s+1000)**2));
+
+W1ss=tf2ss(w1);
+W2ss=tf2ss(w2);
+
+W1=[W1ss,0,0;0,W1ss,0;0,0,W1ss];
+W2=[W2ss,0,0;0,W2ss,0;0,0,W2ss];
+
+A1=W1(2);B1=W1(3);C1=W1(4);D1=W1(5);
+A2=W2(2);B2=W2(3);C2=W2(4);D2=W2(5);
+
+//Ptf=[W1, W1*Gtf;
+//     0*eye(3,3), W2*Gtf;
+//     eye(3,3), -Gtf];
+
+//Pss=tf2ss(Ptf);
+
+r=[3,3];
+
+Ap=[At,0*eye(5,6),0*eye(5,6);
+   B1*C,A1,0*eye(6,6);
+   B2*C,0*eye(6,6),A2];
+
+Bp=[0*eye(5,3),B;
+   -B1,B1*D;
+    0*eye(6,3),B2*D];
+
+Cp=[-D1*C,-C1,0*eye(3,6);
+    D2*C,0*eye(3,6),C2;
+    -C,0*eye(3,6),0*eye(3,6)];
+
+Dp=[D1,-D*Dt;
+    0*eye(3,3),D2*Dt;
+   eye(3,3),-D];
+
+P=syslin('c',Ap,Bp,Cp,Dp);
+
+[sk,mu]=h_inf(P,r,0,1,30);
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+//Safonov-Chiang
+
+A=[-2.2567e-02  -36.617  -18.897  -32.090  3.2509  -7.6257e-1;
+  9.2572e-05  -1.8997  9.8312e-01  -7.2562e-04  -1.7080e-01  -4.9652e-03;
+    1.2338e-02  11.720  -2.6316  8.7582e-04  -31.604  22.396;
+    0           0        1       0            0       0  ;
+    0           0        0       0           -30      0  ;
+    0           0        0       0            0      -30  ];
+
+B=0*ones(6,2);B(5,1)=30;B(6,2)=30;
+C=0*ones(2,6);C(1,2)=1;C(2,4)=1;
+D=0*ones(2,2);
+G=syslin('c',a,b,c);r=syssize(G);
+W1=[(s+0.01)/(1+0.01*s)  0;0    (s+0.01)/(0.01*s+1)], 
+tau=0.0005;
+W2=[1000/(s*s)   0;0    1000/(s*s*(s*tau+1))],
+
+
+//
+
+s=poly(0,'s');
+G=(s+1)/(s*s+1);
+w1=1;w3=(s+2);
+
+P=[w1, -w1*G;
+   0 , w3*G
+   1 , -G];
+
+r=[1,1];
+Pss=tf2ss(P);
+P12=Pss(1:2,2);P21=Pss(3,1);P22=Pss(3,2);
+sl=p21;
+sm=[s*eye()-sl(2),sl(3);sl(4),-sl(5)];
+detr(sm)       //Poles of G --> Transmission zeros of P21
+[sk,mu]=h_inf(P,r,0,1,20);   //Fails
+
+Ptmp=P;fil=(1+s*s)/(1+s)/(1+s);Ptmp(:,2)=Ptmp(:,2)*fil;
+[sk,mu]=h_inf(Ptmp,r,0,1,10);    //---> sk=0 mu=1 OK
+
+// Random example...
+
+nx=5;
+nu=2;ny=3;
+nw=4;nz=3;
+A =[1,2,3,0,1;4,5,6,1,0;7,8,9,1,1;0,1,0,1,1;5,4,3,2,1];
+B2=[1,2;2,1;1,1;1,0;2,0];  //(nx,nu)
+C2=[1,1,0,1,0;0,1,1,2,1;2,3,1,0,1];   //(ny,nx)
+r=[ny,nu];
+
+//%Q=rand(nx+nu,3)*rand(3,nx+nu);Q=Q*Q';
+//%R=rand(nx+ny,2)*rand(2,nx+ny);R=R*R';
+
+
+D12=[0,1;2,1;3,0];    //%(nz,nu)
+D21=[1,2,1,1;2,1,1,2;0,2,1,1];    //%(ny,nw);
+C1=[1,1,0,1,1;2,1,1,0,1;1,0,0,1,1];   //%(nz,nx)
+B1=[1,0,1,0;2,3,3,1;3,1,1,0;1,2,1,0;0,1,1,0];  //%(nx,nw)
+//%Q#=[C1 D12]'*[C1 D12];
+//%R#=[B1;D21]*[B1;D21]';
+D11=2*[1,2,1,0;2,2,1,1;0,0,2,1]; //%(nz,nw)
+D11=D11;
+D22=0*[1,2;2,1;3,1];       //%(ny,nu)
+gama=64;
+B1=B1/sqrt(gama);
+D11=D11/sqrt(gama);
+D21=D21/sqrt(gama);
+B2=B2*sqrt(gama);
+D12=D12*sqrt(gama);
+C1=C1/sqrt(gama);
+D11=D11/sqrt(gama);
+D12=D12/sqrt(gama);
+C2=C2*sqrt(gama);
+D21=D21*sqrt(gama);
+B=[B1,B2];C=[C1;C2];
+D=[D11,D12;D21,D22];
+
+P=syslin('c',A,B,C,D);
+
+[K,mu]=h_inf(P,r,0,10,20);
+
+//    Difficult example gamma between 1.6 and 1.8
+// h_inf finds 1.671244...  and
+// K=[(9.7574+0.2790*s)/(0.023419+s), 0.1346]
+// with P1=minreal(P) (state dimension 5 instead of 8 for P)
+// get same gamma but a controller of order (4,3)
+w=[-1.05033595281261e+00 1 0  -2.72165256144227e+03 0 0 0 0 0 0 0 0 0 0;
+   -1.39973553704802e+02 0 0  -2.94468368703638e+06 0 0 0 0 0 0 0 0 0 0;
+   0 0  -2.100e+02  -2.250e+04 0 0 0 0 0   3.80748204520352e-01 0 0 0 1;
+   0 0 1 0 0 0 0 0 0 0 0 0 0 0;
+   0,0,0,0, -4.89284164859002e+03, -7.98004338394794e+06,...
+ -4.33839479392625e+09,0 , 4.49639998582548e-02,0,0,0,0,0;
+   0 0 0 0 1 0 0 0 0 0 0 0 0 0;
+   0 0 0 0 0 1 0 0 0 0 0 0 0 0;
+  1.80466256786985e+00 0 0   4.58687395748091e+03 0 0 0,...
+  -2.34192037470726e-02 0 0 1  -1 0 0;
+   0 0 0 1.25100080458419e+04 0 0 0 0 0 0 0 0 0 0;
+   0 0 0 0 0 0 0 0 0   6.000e-01 0 0 0 1.57584459460774e+00;
+   6.31631898754447e-01 0 0 1.60540588511832e+03 0 0 0,...
+   4.67564402810304e+00 0 0   3.500e-01  -3.500e-01 0 0;
+ 1.80466256786985e+00 0 0   4.58687395748091e+03 0 0 0 0 0 0 1  -1 0 0;
+   0 1 0 0 0 0 0 0 1.79855999433019e+00 0 0 0 1 0];
+
+A=W(1:8,1:8);B=W(1:8,9:14);C=W(9:13,1:8);D=W(9:13,9:14);
+P=syslin('c',A,B,C,D);
+r=[2,1];
+
+[K,mu]=h_inf(P,r,0,1,40);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//Automotive fuel control
+eps=0.001;
+s=poly(0,'s');
+P=(5.498*s*s+4.007d+2*s-4.444d+5)/((s+eps)*(s**3+93.72*s*s+9.520d+3*s+1.214d+5));
+P=P*(18000/5.498);
+ro=2;
+w2=(s+30)*(s+60)/30/60/10;
+w1=20*ro/((s+eps)*(s+20));
+G=[w1,-w1*P;
+   0 ,w2*P;
+   1, -P];
+W1ss=tf2ss(w1);
+Pss=tf2ss(P);
+W2P=W2*P;
+//calcul de cw2;
+www=clean((s*eye()-Pss(2))**(-1)*bp);
+mat=[horner(www,-1),horner(www,1),horner(www,10),horner(www,20)];
+ww=clean(w2p-1);
+f=[horner(ww,-1),horner(ww,1),horner(ww,10),horner(ww,20)];
+cw2=f/mat;
+W2pss=syslin('c',Pss(2),Pss(3),cw2,1);
+
+AGss=[W1ss(2),W1ss(3)*Pss(4);0*eye(4,2),Pss(2)];
+BGss=[-W1ss(3),0*ones(2,1);0*ones(4,1),Pss(3)];
+CGss=[-W1ss(4),0*ones(1,4);0*ones(1,2),w2pss(4);0*ones(1,2),-pss(4)];
+DGss=[0,0;0,1;1,0];
+Gss=syslin('c',AGss,BGss,CGss,DGss);
+r=[1,1];
+[a,b1,b2,c1,c2,d11,d12,d21,d22]=smga(Gss,r);
+
+[K,mu]=h_inf(Gss,r,0,100,20);
+
+////////////SERE
+s=poly(0,'s')
+G=-4.714841d-3*(s**2+0.0218517*s+1.793965)*(s**2-0.3057194*s+4.752967);
+g=g*(s**2+0.4002371*s+4.812033)*(s-4.911808);
+g=g/((s**2)*(s**2+0.0221061*s+1.833885)*(s**2+7.392301d-3*s+3.18453));
+g=g/(s**2+2.716316d-2*s+4.624171);
+w1=0.33333333*s**3+0.3494322*s**2+9.157714d-2*s+0.012;
+w1=w1/((s**3)*(0.01*s+1));
+w2=2*s/(s+100);
+P=[w1,-w1*g;
+   0 , w2;
+   1, -g];
+[sk,mu]=h_inf(P,[1,1],1,1,1);
diff --git a/modules/cacsd/demos/robust/mu.dem b/modules/cacsd/demos/robust/mu.dem
new file mode 100755
index 000000000..e65ab5aa3
--- /dev/null
+++ b/modules/cacsd/demos/robust/mu.dem
@@ -0,0 +1,47 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+//test of musolve
+mode(-1)
+M1 = [ ..
+   5.2829       5.7683      -2.4004       1.2205      -6.4148    
+   9.7769e-01   2.9786      -3.0408       5.0257e-01  -2.6504    
+   7.0819       9.6324      -3.5750       3.3016      -6.7030    
+  -1.6261      -2.9763       1.6870      -1.0603       1.2211    
+   2.3056       4.3712      -2.4785       2.6152      -1.9832    ];
+
+M2 = [ ..
+  -1.1308      -1.7785       8.7974e-01  -7.5206e-01   1.2089    
+  -3.5255e-01  -5.7002e-01   2.9305e-01  -2.5442e-01   3.7691e-01
+  -1.3724      -2.1501       1.0741      -9.1188e-01   1.4669   
+   3.5839e-01   5.5101e-01  -2.7290e-01   2.3565e-01  -3.7663e-01
+  -4.9015e-01  -7.8706e-01   4.0215e-01  -3.3617e-01   5.3261e-01];
+//*******************************************************
+M=M1 +%i*M2;
+// Let the structure be all scalar blocks
+K = [1 1 1 1 1]'; 
+// Let the first, the third and the fifth blocks be real,
+// and let the rest of blocks be complex
+T = [1 2 1 2 1]';
+[D,g,mu] = musolve(M,K,T);
+spec(M'*D*M+%i*(G*M-M'*G')-mu^2*D)
+
+
+// Now, we compute it again with respect to all complex blocks
+T = [2 2 2 2 2]';
+[D,g,mu] = musolve(M,K,T);
+spec(M'*D*M+%i*(G*M-M'*G')-mu^2*D)
+
+
+T = 3*[1 1 1 1 1]';
+[D,g,mu] = musolve(M,K,T;
+//
+K = [2 3]';
+T = [2 2]';
+[D,g,mu] = musolve(M,K,T);
+
+
diff --git a/modules/cacsd/demos/robust/rob.dem b/modules/cacsd/demos/robust/rob.dem
new file mode 100755
index 000000000..6c43adae5
--- /dev/null
+++ b/modules/cacsd/demos/robust/rob.dem
@@ -0,0 +1,51 @@
+//
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is distributed under the same license as the Scilab package.
+//
+
+mode(1);
+lines(0);
+s=poly(0,'s');
+
+//MAC-FARLANE PROBLEM for G=1/s^3;
+G=syslin("c",1/s^3);
+
+[P,r]=macglov(G);clean(P)
+
+//Optimal controller K:
+
+halt(_("Press Return to continue ... \n"));
+//K Optimal controller , ro = gamaopt^-2;
+[K,ro]=h_inf(P,r,0,1,30);
+K,gamaopt=1/sqrt(ro)
+
+// Check internal stability:
+
+halt(_("Press Return to continue ... \n"));
+
+Tzw=lft(tf2ss(P),tf2ss(K));
+
+[Acl,Bcl,Ccl,Dcl]=abcd(Tzw); spec(Acl)
+
+//Optimal gain:
+
+halt(_("Press Return to continue ... \n"));
+
+ga=h_norm(Tzw)
+
+//Compare with gamaopt
+
+ga-gamaopt
+
+//Compare with theory
+
+halt(_("Press Return to continue ... \n"));
+
+[N,M]=lcf(tf2ss(G)); //Left coprime factorization of G.
+
+nk=hankelsv([N,M]);
+
+ro-( 1-nk(1) )
+
diff --git a/modules/cacsd/demos/tracking/track.dem.sce b/modules/cacsd/demos/tracking/track.dem.sce
new file mode 100755
index 000000000..5d4535450
--- /dev/null
+++ b/modules/cacsd/demos/tracking/track.dem.sce
@@ -0,0 +1,68 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+// This file is released under the 3-clause BSD license. See COPYING-BSD.
+
+mode(1)
+rand("seed", 0)
+rand("normal")
+// Define the plant
+nx=5;ny=1;nu=3;
+Plant=ssrand(ny,nu,nx);
+J=rand(Plant(5));Plant(5)=0*J;
+[F,G,H,J]=Plant(2:5);
+
+// define state space model of the Signal to track
+nw=4;nuu=2;
+A=rand(nw,nw);
+A=A-max(real(spec(A)))*eye(A);
+B=rand(nw,nuu);
+C=2*rand(ny,nw);
+D=0*rand(C*B);
+xx0=0*ones(nw,1);
+Model=syslin("c",A,B,C,D,xx0);
+
+// Input to Model (t is a vector), nuu components
+deff("[ut]=uu(t)","ut=[sin(3*t);cos(0.5*t)]");
+
+// Compute Signal to track
+dt=0.05;tmax=60;
+instants=0:dt:tmax;
+totrack=flts(uu(instants),dscr(Model,dt));    //Signal
+
+my_handle             = scf(100001);
+clf(my_handle,"reset");
+plot2d(instants',totrack',axesflag=1);
+
+halt(_("Press Return to continue ... \n"));
+
+clf(my_handle,"reset");
+
+//
+[L,M,T]=gfrancis(Plant,Model);
+
+//    Stabilizing the plant
+K=-ppol(F,G,-0.3*ones(1,nx));
+
+// Bigsyst= closed loop system: um --> [yplant;ymodel].
+// full state gain is [K, L - K*T] * (xplant, xmodel) + M * umodel
+BigA=[F+G*K,G*(L-K*T);
+0*ones(nw,nx),A];
+BigC=[H+J*K,J*(L-K*T);
+0*ones(ny,nx),C];
+BigB=[G*M;
+B];
+BigD=[J*M;
+D];
+
+x0=ones(nx,1);
+BigX0=[x0;xx0];
+Bigsyst=dscr(syslin("c",BigA,BigB,BigC,BigD,BigX0),dt);
+
+
+z=flts(uu(instants),Bigsyst);
+plot2d([instants',instants'],..
+[totrack(1,:)',z(1,:)'], axesflag=1);
+curves = gce();
+captions(curves.children,["Signal to track","Computed signal"],"upper_caption");
+xtitle("tracking");