Procédures | Cast3M

* CORMAN    PROCEDUR  CHAT      13/02/05    21:15:00     7691           
 
 
 
'DEBPROC' corman rit*'RIGIDITE' MO*'MMODEL' MA*'MCHAML'
  ORDRE*'ENTIER' U0*'CHPOINT' S0*'MCHAML' ffi*CHPOINT 
  ETAB12*'TABLE';
* quan,d on arrive ici on se trouve sur la config en debut
* de pas. ffi est le résidu et s0 est la contrainte en début de pas
* 
*
* attention dans cette version on n'utilise pas les approximants de Padé.
*
sgdlamda='SIGN' 1 ;
*** mess ' entree dans corman  ' ; 
IDYN=ETAB12.'DYNAMIQUE';
'SI' IDYN;
  MMA= ETAB12. 'MASSE';
  QUADTT2=4.D0 / ETAB12.'DT' / ETAB12.'DT';  
  MMAPLU=MMA * QUADTT2;
*  mess ' quadtt2 ' quadtt2;
'FINSI';
N=ORDRE ;
DEP1=U0;
ncontmax=250;
isol1=faux;
prec = 1.D-8; 
ttdep=table; ttdep.0 = U0; ttlam=table ; ttlam.0=0.;
ttsig=table; ttsig.0=s0;NA=0 ;
depreel0 = U0 'ENLE' 'LX' ;
depini=U0 'ENLE' 'LX' ;
sigreel0 = S0 ;
ffs0= 'BSIGMA' S0 MO MA;
lamreel0 = 0.d0;
IOUT=1;
ritblo='EXTR' rit 'RIGI'   'MULT';
* on veut lambda=1.d0
llld = mots 'UX' 'UY' 'UZ' 'RX' 'RY' 'RZ' ; 
IMAS=FAUX; 
IDIMV= 'VALE' 'DIME';
IMODD= 'VALE' 'MODE';
'SI' ( 'EXISTE' MO 'ELEM' 'TRI3' 'QUA4' 'TRI6' 'QUA8'
     'CUB8' 'TET4' 'PRI6' 'PYR5' 'CU20' 'TE10' 'PR15' 'PY13' )  ;
IMAS=VRAI;MOMAS= 'EXTR' MO 'ELEM''TRI3' 'QUA4' 'TRI6' 'QUA8'
'CUB8' 'TET4' 'PRI6' 'PYR5' 'CU20' 'TE10' 'PR15' 'PY13'; 
'FINSI';
ICOQ=FAUX;
'SI' ( 'EXISTE' MO 'ELEM' 'COQ2' 'DKT' );ICOQ=VRAI;
MOCOQ= 'EXTR' MO 'ELEM'  'COQ2' 'DKT'; 'FINSI';
'SI'  ('EXIS' MO ELEM 'POUT'); IOUT=0;'QUITTER' CORMAN; 'FINSI';
ff = ffi ;
lam0 = 0. ;
geoini = form; 
sigcour=s0;
*m0
deltatot= u0*0.;
'REPE' IT ;
'SI'  ('EXIS' MO ELEM 'POUT'); IOUT=0;'QUITTER' IT; 'FINSI';
'SI'  ('EXIS' MO ELEM 'COQ2'); IOUT=0;'QUITTER' IT; 'FINSI';
*'SI' (&it > 1) ; mess ' continuation ' ( &it - 1) lam0;'FINSI';
'SI' (&it > ncontmax);MESS ' Pas de solution apres ' ncontmax 
'continuations'; IOUT=0;
'QUITTER' IT; 'FINSI';
**************************Initialisation**************************
 
*
* quelques initialisation de table servant aux stockages intermediaire
*
tadep=table; talam=table;tasig=table;taksi = table;
tagrad = table; tagraf=table;
tabase = table;tabalpha = table; tabvald = table; tabpoly = table;
*
************************** calcul du premier ordre******************
*
deltau0= u0*0.;
kspp = 'KSIGMA' sigreel0 mo ma;
rit= 'RIGI' MO MA;
rit=rit 'ET' ritblo;
ritot = rit 'ET' kspp;
'SI'  IDYN;
  Ritot= ritot 'ET' MMAPLU;
'FINSI';
DE1= 'RESOU' ritot ff  'NOUNIL';
psc = 'XTX' de1 ; psc2 = psc  +  1.D0 ** 0.5D0;
lambda1 = 1.d0 /  psc2 ;
dep1 = de1 * lambda1;
s1 = 'SIGMA' 'LINE' dep1 mo ma;
*
*  calcul d'une norme de convergence
*
* remplissage des tables pour ordre 1
tadep . 1 = dep1;
talam . 1 = lambda1 ;
tasig . 1 = s1;
taksi . 1 = 'KSIGMA' s1 mo ma;
tagrad . 1 = 'GRAD' dep1 mo ma;
'SI' ICOQ ;tagraf . 1 = 'GRAF' dep1 mo ma ;'FINSI';
tabalpha. 1 = table;
tabalpha . 1 . 1 = ('XTX' tadep. 1)**0.5D0  ;
tabase. 1 = (tadep. 1)/(tabalpha. 1 . 1);
p=1;
*mess ' ordre ' p ' lambda ' lambda1;
*temps place;
*
******************** calcul des ordres superieurs*****************
*
*mess ' debut de boucle';
'REPETER'  bordre (N - 1);
  p = &bordre + 1;
*
*calcul des ksig* (si) * u(p - i + 1);
  nfois = p - 1 ;
  'REPETER'  ks nfois;
    fksi = taksi . &ks * tadep .( p - &ks);
    'SI' ( &ks 'EGA' 1);
      fks = -1.D0 * fksi;
   'SINON'; 
*menagement de l'espace memoire
      XXX = fks;
      fks = fks - fksi;
    'FINSI';
  'FIN' ks;
  fks = fks  ;
*
* calcul de bsigma ( D *0.5( gradui * graduj))
*
  'REPETER' bs nfois;
    gradi = tagrad . &bs;
    gradj = tagrad . ( p - &bs );
    'SI' ICOQ ;
      grafi = tagraf . &bs;
      grafj = tagraf . ( p - &bs );
    'FINSI';
 
* les cartes suivantes sont pour les elements coques
    'SI' ICOQ ;
      gradi = 'REDU' gradi  MOCOQ; gradj = 'REDU' gradj MOCOQ;
      grafi = 'REDU' grafi  MOCOQ; grafj = 'REDU' grafj MOCOQ;
 
      Uxxi='EXCO' gradi 'UX,X'  'SCAL'  chan TYPE 'SCAL';
      Uxyi='EXCO' gradi 'UX,Y'  'SCAL'  chan TYPE 'SCAL'  ;
      Uxzi='EXCO' gradi 'UX,Z'  'SCAL'  chan TYPE 'SCAL' ;
      Uyxi='EXCO' gradi 'UY,X'  'SCAL'  chan TYPE 'SCAL' ;
      Uyyi='EXCO' gradi 'UY,Y'  'SCAL'  chan TYPE 'SCAL' ;
      Uyzi='EXCO' gradi 'UY,Z'  'SCAL'  chan TYPE 'SCAL' ;
      Uzxi='EXCO' gradi 'UZ,X'  'SCAL'  chan TYPE 'SCAL' ;
      Uzyi='EXCO' gradi 'UZ,Y'  'SCAL'  chan TYPE 'SCAL' ;
 
      Uxxj='EXCO' gradj 'UX,X'  'SCAL'  chan TYPE 'SCAL' ;
      Uxyj='EXCO' gradj 'UX,Y'  'SCAL'  chan TYPE 'SCAL' ;
      Uxzj='EXCO' gradj 'UX,Z'  'SCAL'  chan TYPE 'SCAL' ; 
      Uyxj='EXCO' gradj 'UY,X'  'SCAL'  chan TYPE 'SCAL' ;
      Uyyj='EXCO' gradj 'UY,Y'  'SCAL'  chan TYPE 'SCAL' ;
      Uyzj='EXCO' gradj 'UY,Z'  'SCAL'  chan TYPE 'SCAL' ;
      Uzxj='EXCO' gradj 'UZ,X'  'SCAL'  chan TYPE 'SCAL' ;
      Uzyj='EXCO' gradj 'UZ,Y'  'SCAL'  chan TYPE 'SCAL' ;
      Bxxi='EXCO' grafi 'BX,X'  'SCAL'  chan TYPE 'SCAL' ;
      Bxyi='EXCO' grafi 'BX,Y' 'SCAL'  chan TYPE 'SCAL'  ;
      Byxi='EXCO' grafi 'BY,X' 'SCAL'  chan TYPE 'SCAL'  ;
      Byyi='EXCO' grafi 'BY,Y' 'SCAL'  chan TYPE 'SCAL'  ;
      Bxxj='EXCO' grafj 'BX,X' 'SCAL'  chan TYPE 'SCAL'  ;
      Bxyj='EXCO' grafj 'BX,Y' 'SCAL'  chan TYPE 'SCAL' ;
      Byxj='EXCO' grafj 'BY,X'  'SCAL'  chan TYPE 'SCAL';
* Nouvelle version
      epss=( 0.5 * ( Uxxi*Uxxj + (Uyxi*Uyxj) + (Uzxi*Uzxj ) ) ) 
         'NOMC' 'EPSS';
      eptt=( 0.5 * ( Uxyi*Uxyj + (Uyyi*Uyyj) + (Uzyi*Uzyj ) ) ) 
        'NOMC' 'EPTT';
      gast=(0.5 * (Uxyi*Uxxj + (Uyyi*Uyxj) + (Uzyi*Uzxj) + (Uxyj*Uxxi) + 
            (Uyyj*Uyxi) + (Uzxi*Uzyj) ) )'NOMC' 'GAST';
      rtss=(epss *0.) 'NOMC' 'RTSS';
      rttt=rtss 'NOMC' 'RTTT';
      rtst=rtss 'NOMC' 'RTST';
 
      ep2 = epss + eptt + gast + rtss + rttt + rtst;
 
    'FINSI';
 
* les cartes suivantes sont pour les elements massifs
 
    'SI' IMAS ;
       gradi = 'REDU' gradi MOMAS;gradj = 'REDU' gradj MOMAS;
       gru2 =(( 'TAGR' gradi) * gradj mo) * 0.5D0 ;
       'SI' ('NEG' imodd 'AXIS');
          epxyz= 'EXCO' gru2  ('MOTS' 'UX,X' 'UY,Y' 'UZ,Z')
                              ('MOTS'  'EPXX' 'EPYY' 'EPZZ');
          gaxyz= ('EXCO' gru2  ('MOTS'  'UX,Y' 'UX,Z' 'UY,Z')
                               ('MOTS' 'GAXY' 'GAXZ' 'GAYZ'))*2.d0;
          ep2 = epxyz + gaxyz;                 
        'SINON';          
            epxx = 'EXCO' gru2 UR,R   'EPRR';
            epyy = 'EXCO' gru2 UZ,Z   'EPZZ';
            epzz = 'EXCO' gru2 UT,T   'EPTT';
            gaxy = ('EXCO' gru2 UR,Z   'GARZ')*2.D0 ; 
            gaxz = ('EXCO' gru2 UR,T   'GART')*2.D0 ;
            gayz = ('EXCO' gru2 UT,Z   'GATZ')*2.D0 ;
            ep2 = epxx + epyy + epzz + gaxy + gaxz + gayz;                      
        'FINSI';
    'FINSI';
 
* fin du traitement spécifique
    ep2 = 'CHANG'  ep2 'TYPE' 'DEFORMATIONS'; 
    'SI' ( 'EGA' &bs 1);
      ep2tot = ep2 ;
    'SINON';
      XXX = ep2;
      ep2tot = ep2tot + ep2;
    'FINSI';
  'FIN' bs;
  sig = 'ELAS' ep2tot mo ma  ;
  fbs = ('BSIGMA' sig mo ma) * -1.D0 ;
  u2 = 'RESOU' ritot ( fks + fbs)   'NOUNIL'; 
  lamdai = (u2 'XTY' tadep.1   llld llld )* -1.D0 * lambda1;
 
*
*
* remplissage resultats ordre p
*
  talam . p = lamdai ;
  tadep . p = lamdai / lambda1  * dep1 + u2;
  tasig . p = tadep . p 'SIGMA' 'LINE'  mo ma + sig;
  taksi . p = 'KSIGMA' tasig . p mo ma;
  tagrad . p = tadep . p 'GRAD' mo ma;
  'SI' ICOQ;tagraf . p = tadep . p 'GRAF' mo ma;'FINSI';
* petit calcul pour voir l'evolution du rayon de convergence
*  norU1 = ('XTX' tadep. 1)**0.5D0 ;
*  norUn = ('XTX' tadep. P)**0.5D0 ;
*  aserie = ( prec *  norU1  /  norUn  ) ** ( 1.D0 / ( p - 1 ) );
*  mess ' ordre ' p ' lambda ' lamdai 'raymax' aserie; 
*
**Orthogonalisation de la base pour le calcul des approximants de Pade**
*
***  tabalpha. p = table;
*calcul somme de (Up.Wq).Wq
***  'REPE' ort (p - 1 );
***    alphapq = (tadep. p 'XTY' tabase. &ort llld llld );
***    'SI' ('EGA' &ort 1);
***      vec = tadep. p - (alphapq*tabase. &ort) ;
***    'SINON' ;
***      vec = vec - (alphapq*tabase. &ort);
***    'FINSI';
*remplissage de la table contenant les alphapq issu de la projection
***    tabalpha. p . &ort = alphapq;
***  'FIN' ort ;
 
*calcul des vecteurs de la nouvelle base
***  Norvec = ('XTX' vec)**0.5D0 ;
***  base = vec / Norvec ;
***  tabalpha. p . p = norvec ;
*remplissage de la table des vecteurs de la base orthonormee
***  tabase. p = base;
*temps place;
'FIN' bordre;
 
 
**********************************************************************
*************calcul du rayon de convergence maximal 'aserie'**********
**********************************************************************
 
norU1 = ('XTX' tadep. 1)**0.5D0 ;
norUn = ('XTX' tadep. N)**0.5D0;
aserie = ( prec *  norU1  /  norUn  ) ** ( 1.D0 / ( N - 1 ) );
*
**menage de la place memoire
*
'REPE' titi ( 'DIME' taksi );
       detr taksi. &titi;
'FIN'  titi;
'OUBLI'  taksi;
'REPE' titi ( 'DIME' tagrad );
       'DETR'  tagrad. &titi;
  'SI' ICOQ;     'DETR' tagraf. &titi;'FINSI';
'FIN' titi;
'OUBLI' tagrad; 'OUBLI'  tagraf; 'OUBLI' gradi; 'OUBLI' gradj;
 
*
****Calcul de dn (coeff de Dn) en resolvant un systeme triangulaire****
*
 
***'REPE'  dn (N - 1);
***  tabvald. &dn = table;
***  'REPE' equa &dn;
***     addi = tabalpha. (&dn + 1) . ( &dn + 1 - &equa ); 
***    'SI' ( &equa 'NEG' 1 ) ;
***      'REPE' ilig (&equa - 1);   
***        ope = tabalpha. ( &dn + 1 - &ilig ) . (&dn -  &equa + 1) ;
***        mult = ope * tabvald. &dn. &ilig ; 
***        addi = mult + addi ;
***      'FIN' ilig;
***    'FINSI';
***    den =  tabalpha. (&dn + 1 - &equa). (&dn + 1 - &equa);
***    tabvald. &dn . &equa = (-1. * addi)/ den;
***  'FIN' equa ;
***'FIN' dn;
 
******************************************************************
*****Calcul du critere pour les approximants de Pade 'apade'******
******         Recherche de la valeur apade        ***************
******************************************************************
adep = aserie*3.;
acou= adep;
 
ainf = 0.D0;
recher = faux;
epsipade =prec;
puis2=1;
**'REPE' cherchea ;
**************creation d'une table puissance de a *************
  tapuisa = table;
 
  'REPE' puis (N - 1);
     puisacou = acou**&puis;
     tapuisa. &puis = puisacou;
  'FIN' puis;
 
********
**********Calcul des polynomes Dn-1(a) pour chaque ordre***************
********
***  tabDn = table;
**la 1er valeur de cette table vaut 1**
***  tabDn. 0 = 1.D0;
***  adpoly= 1 ; 
***  'REPE'  polyDn (N - 1); 
***     adpoly = (tapuisa. &polyDn*tabvald.(N - 1). &polyDn) + adpoly;
***     tabdn.  &polyDn = adpoly;
***  'FIN' polyDn;
**
******************Expression des approximants de Pade***************
** 
************Calcul de Pn(U(a)) -  Pn-1 (U(a))*************
 
**Expression de Pn(U(a)) 
***  PnUa = U0;
*  mess ' avant boucle boupa'; temps place;
***  'REPE' boupa (N - 1);    
***    som=tapuisa.&boupa/tabDn.(N - 1)*tadep. &boupa*
***        tabDn.(N - 1 - &boupa);
***    PnUa = som + PnUa ;
***  'FIN' boupa;
*  mess ' apres boucle pnua'; temps place;
**Expression de Pn-1 (U(amp))
***  Pmoins1 = U0;
***  'REPE' boupm (N - 2);    
***     sam=tapuisa. &boupm/tabDn.(N - 2)*tadep.&boupm*
***         tabDn.(N - &boupm - 2);
***     Pmoins1 = sam + Pmoins1 ;   
***  'FIN' boupm;
*****Calcul de (Pn(a) - Pn-1 (a)) / ( Pn(a) - U0)***
*calcul norme de (Pn(a) - Pn-1 (a) )
***  nt = PnUa - Pmoins1 ;
***  nornt = (XTX nt)**0.5D0;
*calcul norme de ( Pn(a) - U0)
***  np =  PnUa - U0;
***  nornp = (XTX np)**0.5D0;
***  crit = nornt / nornp ;
**recherche de apade par dichotomie
***  'SI' (  'EGA'   recher faux);    
***     'SI' (crit < epsipade);
***        ainf = acou;
***        acou = acou + adep;
***        'ITER' cherchea;   
***     'SINON';
***        asup = acou;
***        acou = ( ainf + asup )/2;
***        recher = vrai;
***       compteur = 1;   
***     'FINSI';   
***  'SINON';
***     compteur = compteur + 1;
***     'SI' (crit '>' epsipade);
***        asup = acou;
***        acou = (ainf + asup ) /2.D0;  
***     'SINON';        
***        ainf = acou;
***        acou = (asup + ainf)/2.D0;
***     'FINSI';
***     'SI' ( compteur '>' 16) ;
***        apade = ainf;
***        'QUITTER' cherchea;
***     'SINON'; 
***        'ITERER' cherchea;
***     'FINSI';
*** 'FINSI';
***'FIN' cherchea;
*mess ' aserie ' aserie  ' apade' apade ' compteur ' compteur;
**
********************Resolution de Dn-1 = 0*********************
**
**Creation d'un listreel contenant les valeurs des coeff de Dn-1 ** 
***lll= prog 1.;
***'REPE' resol (N - 1) ;  
***    lla= prog  tabvald. (N - 1) . &resol ;   
***    lll = lll et lla;
***
***'FIN' resol;
 
****Resolution ******
***solpet = 'RACP' 'INTE' 0.D0 apade lll ;
***Recherche de la plus petite racine reelle positive***
***solpet = 1.e10;
***'SI'( ('DIME' soluu) > 0) ; 
***  'REPE' ppsol ('DIME' soluu) ;   
***     aDn = 'EXTR' soluu &ppsol ;
***     'SI' ( aDn '>' 0.D0);
***       'SI' (aDn '<' solpet) ;
***          solpet = aDn;
***       'FINSI';
***       'QUIT' ppsol;
***     'FINSI'; 
***  'FIN' ppsol ;
***'FINSI';
***
***
*mess ' solpet ' solpet;
***'SI' ( 'NEG' solpet 'VIDE');
***  'SI' ( solpet '<' apade);
***      apade =solpet ;
***'FINSI';
 
 
*mess 'aDn ' solpet ' apade' apade 'aserie' aserie ;
 
********************************************************************
********Recherche de 'a' pour que lambda vaille 1 pour les 2 cas****
********************************************************************
 
 
***xpetpa = apade;lamxpet=0.;
ifini=faux; 
*  mess 'Recherche de la valeur de a pour que lambda vaille 1: xpet' ;
prec1= prec * 0.1; 
ifiser=faux;ifipad=faux;
*'SI'( 'NON' PADE);
 
************Pour le developpement en serie**************************
 
* recherche pour savoir si on croise lambda=1 Pour cela on evalue d'abord la
* valeur obtenu pour aserie et si dlam + lam0 >1 on recherche les
*
* zeros du polynome donnant lambda
* talam contient les coeff du polynoem sauf talam .0 
  xnew = lamreel0;
  'REPE' valr N;    
     apuisn = aserie ** &valr;
     xnew=apuisn*talam.&valr + xnew;
  'FIN'  valr; 
*   mess ' xnew ' xnew;
  'SI' ('>' xnew 1.)  ;
     lll = prog -1.D0 ;
     'REPE' bobo N ; 
       lla = prog talam . &bobo; 
       lll= lll et lla;
     'FIN' bobo;
*    'LIST' lll;;
*     xsol= racp lll ;  
     xpet= racp 'INTE' 0.D0 aserie  lll;
*     mess ' c est bon  xpet ' xpet ;
     ifiser=vrai;
*     list xraci; 
* recherche de la plus petite racine relle positive   
*     'SI' ( ('DIME' xsol ) > 0);   
*       'REPE' bobo ( 'DIME' xsol);
*         ab = 'EXTR' xsol &bobo ;
*         mess ' racine ' ab;
*         'SI' ( ab '>' 0.D0 );
*           'SI' ( ab '<' autil) ;
*             xpet = ab;
*            list ( xpet - xraci);
*             ifiser=vrai;             
*           'FINSI';
**           'QUIT' bobo;
*         'FINSI';
*       'FIN' bobo;   
*     'FINSI';
  'SINON';
     xpet = aserie;
*     mess ' il faut recommencer xpet ' xpet;
  'FINSI';
*'SINON' ;
*******************Pour les approximants de Pade********************
* on ne fait ce travail que si le dev en serie n'arrive pas a la solution
***  'SI' ('NON' ifiser);
***  ttpuis = table;
***  solution = faux;
***  xpetinf = 0.;
***  lamrech=   1. ; 
***  'REPE' solua ;
***    'REPE' pui (N - 1);
***       puisxpet = xpetpa**&pui;
***       ttpuis. &pui = puisxpet;
***    'FIN' pui;
**********Calcul des polynomes Dn-1(asol)********** 
***    Dnxpet = table;
***    Dnxpet. 0 = 1.;
***    solpoly= 1 ;  
 
***    'REPE' polysol (N - 1); 
***      solpoly = (ttpuis. &polysol*tabvald.(N - 1). &polysol) + solpoly;
***      Dnxpet.  &polysol = solpoly;
***    'FIN' polysol;
 
*********Expression de Pn(lambda(asol))************* 
***    lamxpet = lam0;
 
***    'REPE' sola (N - 1);
***      sim=ttpuis.&sola/Dnxpet.(N - 1)*talam. &sola*
***          Dnxpet.(N - 1 - &sola);
***      lamxpet = sim + lamxpet ;
***    'FIN' sola;
*    'MESS'  ' solua ' &solua  ' lamxpet' lamxpet;
********Resolution de lamda egale a 1********   
***    'SI' (  'EGA' solution faux);
***      'SI' (lamxpet '<' lamrech);
***        ifpad=faux;
***        'QUITTER'solua;   
***      'SINON';
***        ifipad=vrai;
******        xpetsup = xpetpa;
***        xpetpa = ( xpetinf + xpetsup )/2;
***        solution = vrai;
***        tour = 1;   
***      'FINSI';   
***    'SINON';
***      tour = tour + 1;
***      'SI' (lamxpet > lamrech);
***        xpetsup = xpetpa;
***        xpetpa = (xpetinf + xpetsup ) / 2;
***      'SINON';
***        xpetinf = xpetpa;
***        xpetpa = (xpetsup + xpetinf)/2;
***      'FINSI';     
***      'SI' ( (lamxpet '>' (lamrech - prec1) ) 'ET' 
***             (lamxpet '<' (lamrech + prec1) )) ;
***         xpetpa = xpetinf;
***         'QUITTER' solua;
***      'SINON'; 
***         'ITERER' solua;
***      'FINSI';
***    'FINSI';   
*   mess ' nombre de tour'  tour;
***  'FIN' solua;
***  'FINSI';      
 
*'FINSI';
* mess  'serie ' xpet   'apade'  lamxpet;
*
*
* choix de la methode
*
PADE=FAUX;
'SI' ifiser ;
***  pade=faux; 
  ifini=vrai;
*** 'SINON'; 
*  'SI' ifipad;
*     ifini=vrai;
*     xpet=xpetpa;
*  'SINON';
*    'SI' ( xnew '>' lamxpet );
*       pade=faux;
*    'SINON';
*       xpet=xpetpa;
*    'FINSI';
*  'FINSI';
***    pade=faux;
'FINSI';
 
**creation de table pour les differentes valeurs de a****
valeura = table;
dray = xpet ;
valeura. 1 = dray;
 
********************************************************************
******************Calcul de lambda et de U *************************
**************   par les approximants de Pade  *********************
***********   et par le developpement en series   ****************** 
*******************************************************************
 
***********    par le developpement en series   ******************
 
***'SI' ( 'NON' pade);
 
  'REPE' valree N;      
     apuisn = valeura.1 ** &valree;
     depreel0=apuisn*tadep.&valree + depreel0;
     deltau0=apuisn*tadep.&valree + deltau0;
     sigreel0=apuisn*tasig.&valree + sigreel0;
     lamreel0=apuisn*talam.&valree + lamreel0;
  'FIN' valree;
**************   par les approximants de Pade  *********************  
***'SINON';    
***  tabamax = table;
***  ttDn = table;
***  tabPn = table;  
***  tabamax. 1 = table;  
***  'REPE' amax (N - 1);
***    puisamax = valeura. 1**&amax;
***    tabamax. &amax = puisamax;
***  'FIN' amax;       
***  ttDn. 0 = 1.; 
***  adpolyn= 1 ;
 
**Expression des Dn-1 *********
***  'REPE' Dnmax (N - 1);
***    eta = tabamax. &Dnmax * tabvald. (N - 1). &Dnmax;
***    adpolyn = eta + adpolyn;
***    ttDn. &Dnmax = adpolyn;        
***  'FIN' Dnmax;  
*  mess ' apres boucle dnmax';   temps place; 
***  'REPE' padeexp (N - 1);      
**Expression de Pn(U(a))********
***     exp = ttDn. (N - 1 - &padeexp)/ttDn. (N - 1);
***     PnUa=exp*tabamax. &padeexp*tadep. &padeexp;
***     depreel0 = PnUa + depreel0 ;
***     deltau0=PnUa + deltau0;
**Expression de lambda(U(a))****      
***     Pnlama=exp*tabamax. &padeexp*talam. &padeexp;
***     lamreel0 = Pnlama + lamreel0 ;
**Expression de sigma(U(a))****      
***     Pnsiga = exp*tabamax. &padeexp*tasig. &padeexp;
***     sigreel0 = Pnsiga + sigreel0;      
***  'FIN' padeexp; 
*  mess ' apres padeexp' ;  temps place;
***'FINSI';      
 
ttdep. 1  = depreel0;
 
ttsig. 1  = sigreel0 'PICA' ( depreel0 - U0 ) mo ;
tsetse = ttsig. 1  ;
ttlam. 1 = lamreel0;
*
*test pour savoir si on repart dans le bon sens
*
dlamda =   ttlam .  1  - ttlam .  (1  - 1)   ;
dlam = dlamda /dray ;
sgd  = signe( dlam) ; sgd= sgdlamda;
'SI' (sgdlamda 'NEG' sgd )  ;
   MESS ' On doit rediriger le pas ' ;
  'REPE' bouval1 1;
     valeura. &bouval1 = dray * &bouval1 * (-1) ;
     TAB . 'A' . ( &bouval1 +  NA) = valeura. &bouval1 +a0 ;
  'FIN' bouval1;
**********    par le developpement en series   ******************
*         list valeura ;
***   'SI' ( 'NON' pade);
      depreel0 = depreel0 'ENLE' LX;
     'REPE' valree1 N ;
       apuisn = valeura.1 ** &valree1;
       depreel0=apuisn*tadep.&valree1 + depreel0 ;
       sigreel0=apuisn*tasig.&valree1 + sigreel0 ;
       lamreel0=apuisn*talam.&valree1 + lamreel0 ;
     'FIN' valree1;
 
**************   par les approximants de Pade  *********************
***   'SINON';
***      tabamax = table;
***      ttDn = table;
***      tabPn = table;
***      tabamax. 1 = table;
***      'REPE' am1 (N - 1);
***        puisamax = valeura. 1**&am1;
***        tabamax. &am1 = puisamax;
***      'FIN' am1;
***      ttDn. 0 = 1.;
***      adpolyn= 1 ;
**Expression des Dn-1 *********
***      'REPE' Dn1 (N - 1);
***        eta = tabamax. &Dn1 * tabvald. (N - 1). &Dn1;
***        adpolyn = eta + adpolyn;
***        ttDn. &Dn1 = adpolyn;
***      'FIN' Dn1;
***      'REPE' pade1 (N - 1);
**Expression de Pn(U(a))********
***         exp = ttDn. (N - 1 - &pade1)/ttDn. (N - 1);
***         PnUa=exp*tabamax. &pade1*tadep. &pade1;
***         depreel0 = PnUa + depreel0 ;
**Expression de lambda(U(a))****
***         Pnlama=exp*tabamax. &pade1*talam. &pade1;
***         lamreel0 = Pnlama + lamreel0 ;
**Expression de sigma(U(a))****
***         Pnsiga = exp*tabamax. &pade1*tasig. &pade1;
***         sigreel0 = Pnsiga + sigreel0;
***      'FIN' pade1;
***   'FINSI';
  dlamda =    ttlam .  1  - ttlam .  (1  - 1)   ;
  dlam = dlamda /dray *  (-1 ) ;
  sgd = signe( dlam) ;
*
* MESS ' valeur de dlam aprés redirection  ' dlam  ;
* MESS ' LE SIGNE après redirection   ' SGd ;
*
  ttdep . 1   = depreel0 ;
  ttsig   . 1   = sigreel0 PICA ( depreel0 - U0 ) mo ;
  tsetse=ttsig   . 1   ;
  ttlam  . 1   = lamreel0;
*  U0= depreel0 ;
'FINSI';
sigreel0=sigreel0 'PICA'  deltau0 mo ;
deltatot=deltau0 + deltatot ;
sigpiok= ('SIGMA' deltatot mo ma) + S0;
sigcau= sigpiok 'PICA' deltatot mo;
sigreel0= sigcau;
ttdep.0 = depreel0;
ttsig.0 = S0;
ttlam.0= lamreel0;
'SI' ('NON' ifini) ; 
*mess ' lamrell0 ' lamreel0;
   ff= ff *( 1.d0 - lamreel0);
   depreel0= deltatot 'ENLE' 'LX';
  lamreel0=0.D0;
  'FORM' geoini;
  'FORM' deltatot;
  U0= deltatot;
'FINSI';
*a0 = tab . a  . (1 + na) ;
dlamda =    (ttlam .  1  - ttlam .   0   );
dlamb = dlamda /dray ;
sgdlamda = signe( dlamb) ;
*tab.'SIGNE'=sgdlamda;
'FORM' geoini;
'FORM' deltatot;
'SI' ifini; 
  'MESS' ' sortie de la MAN nombre de segments ' &it  IOUT;
  form geoini;
  'QUITTER' IT;
'FINSI';
'FIN' IT;
'FINPROC' deltatot IOUT;