//                navet.bas : demo de icse
//  calcul trajectoire optimale de rentree d'une navette spatiale
//  *************************************************************
//
// Copyright INRIA
libn=ilib_for_link('icsenb','icsenb.o',[],'f')
nlink=link('./'+libn,['icsenb','icsenf'],'f')

//clear icse contexte
clear nu u uc uv itu dtu y0 tob binf bsup b fy fu obs don;
// les doubles slash introduisent des commentaires
//
t0=0.d0;    // instant initial
tf=1.d0;    // instant final
dtf=0;
ermx=1.d-6; // test d'arret absolu sur la valeur du second membre dans
            // la resolution de l'etat
iu=[0,0,1]; //  iu   :indications sur la structure du controle
    //    iu(1)=1 si l'etat initial depend du controle constant,0 sinon
    //    iu(2)=1 si l'etat initial depend du controle variable,0 sinon
    //    iu(3)=1 si le second membre depend du controle constant,0 sinon
nuc=1;      // nombre de parametres independants du temps
nuv=1;      // nombre de parametres dependants du temps
ilin=0;     // indicateur de linearite :
       // 0 pour un systeme non affine
       // 1 pour un systeme affine dont la partie lineaire n'est pas autonome
       // ilin=2 pour un systeme affine dont la partie lineaire est autonome
nti=150;    //nombre de pas de temps correspondant a dti (premier pas de temps)
dti=tf/nti;
ntf=00;     // nombre de pas de temps correspondant a dtf (second pas de temps)
            // si l'on utilise un seul pas de temps,on doit prendre ntf=0
ny=4;       // dimension de l'etat a un instant donne
nea=0;      // nombre d'equations algebriques (eventuellement nul)
itmx=10;    // nombre maximal d'iterations dans la resolution de
            // l'equation d'etat discrete a un pas de temps donne
nex=1;      // nombre d'experiences effectuees
nob=3;      // dimension du vecteur des mesures pour une experience donnee
            // en un instant donne
ntob=1;     // nombre d'instants de mesure pour une experience donnee
ntobi=1;    // nombre d'instants de mesure correspondant a dti (premier
            // pas de temps)
//
nu=nuc+nuv*(nti+ntf+1); // dimension du vecteur des parametres de controle
//
//  uc(1,nuc)          :controle constant
echtf=2000;
uc=[2500/echtf];
//  uv(1,nuv*(nti+ntf)):controle variable
//if nuv>0, uv(1,nuv*(nti+ntf+1))=0; end;
 alpha0=.20704/.029244; legu='alpha initial : ann. cz';  //annulation cz
 alpha0=17.391;         legu='alpha initial : finesse max';//finesse maximum
legu=' navette americaine .'+legu;
if nuv>0, uv=alpha0*ones(1,nuv*(nti+ntf+1)); end;
//  itu(1,nitu)        :tableau de travail entier reserve a
//                      l'utilisateur
itu=[0];
//  dtu(1,ndtu)        :tableau de travail double precision reserve
//                      a l'utilisateur
raddeg=%pi/180;
//dtu=[ 249.9,         ..//s      1
//      .078540,       ..//cx0    2
//     -.0061592,      ..//cx1    3
//      .621408e-3,    ..//cx2    4
//     -.207040,       ..//cz0    5
//      .029244,       ..//cz1    6
//       83388,        ..//zm     7
//     3.9860119e14,   ..//zmu    8
//     6378166,        ..//rt     9
//     1.2,            ..//ro0   10
//     6700,           ..//h     11
//     raddeg,         ..//      12
//     echtf,          ..//echtf 13
//     0,0,0,0,0,0,0,  ..// inutilises    14 a 20
//     1000,           ..//mise a echelle v    21
//     1,              ..//mise a echelle gam  22
//     1.e5,           ..//mise a echelle z    23
//     1,           ..//mise a echelle l       24
//     1.e6      ];      //cpenal              25
dtu=[ 249.9,         ..
      .078540,       ..
     -.0061592,      ..
      .621408e-3,    ..
     -.207040,       ..
      .029244,       ..
       83388,        ..
     3.9860119e14,   ..
     6378166,        ..
     1.2,            ..
     6700,           ..
     raddeg,         ..
     echtf,          ..
     0,0,0,0,0,0,0,  ..
     1000,           ..
     1,              ..
     1.e5,           ..
     1,           ..
     1.e6      ];      //cpenal              25
//
y0=[7803, -1*raddeg, 121920, 0]; // etat initial
//          (valeur arbitraire si iu(1) ou iu(2) est non nul)
y0=y0./dtu(1,21:24);   //mise a l'echelle de y0
//  tob(1,ntob)        :instants de mesure (compatibilite avec ntob
//                      et ntobi)
tob=[1];
binf=-20*ones(1,nu);//  borne inf des parametres
binf(1,1)=2500/echtf;
bsup=40*ones(1,nu);// borne sup des parametres
bsup(1,1)=4000/echtf;
//
obs(nob,ny)=0; // matrice d'observation
//
//don=[762,            ..//vfin         1
//     -5*raddeg,      ..//gamma final  2
//     24384    ];     ..//zfin         3
don=[762,            ..
     -5*raddeg,      ..
     24384    ];
don=don./dtu(1,21:23);      // mise a l'echelle
nomf='icsenf';   // noms de subroutines de dynamique
legfb=' croissant ';
// changements pour calculer en temps retrograde
retro=1;
if retro>0;
legfb=' retrograde '; ..
  don1=don; ..
  don=y0(1,1:3); ..
  y0(1,1:3)=don1; ..
  nomf='icsenb' ; end ;
//
nap=20;     // nombre d'appels du simulateur
imp=2;      // niveau de debug pour optim
large=100;  // taille de nu au dela de laquelle on choisit un optimiseur
// pour les problemes de grande taille (alg='gc' dans l'appel de optim)
//
exec('icseinit.sce');   
[co,u,g,itv,dtv]=icse(u,nomf,nap,imp);
ulink(nlink)


