// 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 [slm]=arhnk(a,ordre,tol)

    [lhs,rhs]=argn(0),
    if lhs<>1 then error(41),end
    if typeof(a)<>"state-space" then error(91,1),end;
    if a.dt<>"c" then error(93,1),end
    select rhs
    case 2 then istol=0;
    case 3 then istol=1;
    end;

    [a,b,c,d,x0,dom]=a(2:7);
    if(max(real(spec(a)))) > 0 then
        error(msprintf(_("%s: Wrong values for input argument #%d: Stable system expected.\n"),"arhnk",1));
    end
    domaine="c"
    wc=lyap(a',-b*b',domaine)
    wo=lyap(a,-c'*c,domaine)
    if istol==0 then [t,nn]=equil1(wc,wo);
    else [t,nn]=equil1(wc,wo,tol);
    end;
    n1=nn(1);
    ti=inv(t);a=t*a*ti;b=t*b;c=c*ti
    wc=t*wc*t';wo=ti'*wo*ti;
    if ordre>n1 then
        ordre=n1
    end;
    if ordre==n1 then
        a=a(1:n1,1:n1);b=b(1:n1,:);c=c(:,1:n1);
        if lhs==1 then a=syslin("c",a,b,c,d,0*ones(n1,1)),end
        return,
    end;
    sigma=wc(ordre+1,ordre+1)

    r=max(n1-ordre-1,1)

    n=n1
    sel=[1:ordre ordre+r+1:n];seleq=ordre+1:ordre+r
    b2=b(seleq,:);c2=c(:,seleq);
    u=-c(:,seleq)*pinv(b(seleq,:)')
    a=a(sel,sel);b=b(sel,:);c=c(:,sel);
    wo=wo(sel,sel);wc=wc(sel,sel);

    Gamma=wc*wo-sigma*sigma*eye()
    a=Gamma\(sigma*sigma*a'+wo*a*wc-sigma*c'*u*b')
    b1=Gamma\(wo*b+sigma*c'*u)
    c=c*wc+sigma*u*b';b=b1;
    d=-sigma*u+d
    //
    [n,n]=size(a)
    [u,m]=schur(a,"c")
    a=u'*a*u;b=u'*b;c=c*u;
    if m<n then
        t=sylv(a(1:m,1:m),-a(m+1:n,m+1:n),-a(1:m,m+1:n),"c")
        a=a(1:m,1:m)
        b=b(1:m,:)-t*b(m+1:n,:)
        c=c(:,1:m)
    end;
    //
    slm=syslin("c",a,b,c,d,0*ones(m,1));
endfunction