subroutine orbit (day,kode,lfile,nsat,x,revnum,ier)
c*
c***********************************************************************
c*
c*          "Copyright [c] CNES 98 - tous droits reserves"
c*          **********************************************
c*
c*PRO MAGLIB
c*
c*VER 99.03.31 - V 1.0
c*VER 01.06.05 - V 2.0
c*VER 03.01.06 - V 2.1
c*
c*AUT spec. ESA - janvier 1991
c*AUT port. CISI
c*
c*
c*ROL Theme : Astronomie et calcul d'orbite
c*ROL         Calcul du vecteur d'etat a partir des donnees 
c*ROL         d'orbitographie de Cluster (compsants de Chebychev)
c*ROL         par resolution de l'equation de Kepler.
c*ROL         Un fichier de coefficients de Chebychev est necessaire.
c*
c*PAR day    (I) : date julienne ref 2000 (a partir de 01/01/2000)
c*
c*PAR kode   (I) : dimension du vecteur d'etat
c*
c*PAR lfile  (I) : numero logique associe au fichier en entree
c*
c*PAR nsat   (O) : identification du satellite 1, 2, 3 ou 4
c*
c*PAR x      (O) : vecteur d'etat : position en x, y et z du satellite
c*               : (kilometres) ; vitesse en x, y et z (kilometres/secondes)
c*
c*PAR revnum (O) : numero d'orbite
c*
c*PAR ier    (O) : code de retour
c*
c*NOT kode       : 3 = position
c*NOT kode       : 6 = position et vitesse
c*
c*NOT ier        : 0 = OK 
c*NOT ier        : 1 = date julienne erronee (trop petite)
c*NOT ier        : 2 = date julienne erronee (trop grande)
c*NOT ier        : 3 = non continuite de la date entre 2 enreg consecutifs
c*NOT ier        : 4 = valeur de kode erronee (differente de 3 ou 6)
c*NOT ier        : 5 = fichier incoherent (donnees invalides)
c*NOT ier        : 6 = erreur de lecture     
c*
c*NOT            : Le fichier doit etre ouvert au prealable.
c*NOT            : Les enregistrements doivent avoir des dates consecutives
c*NOT            : (date debut enreg i - date de fin enreg i-1 > gap 
c*NOT            : avec gap = 2.d-4)
c*
c*NOT common     : ssave
c*NOT            : Common de sauvegarde des donnees remanentes.
c*
c*INF utilise    : sans objet
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.06.05 - correction de commentaires de code
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*
c***********************************************************************
c*
      implicit none
c
c     ---------------------------------
c*FON Declaration identificateur rcs_id
c     ---------------------------------
c
      character rcs_id*100
c
c     --------------------------
c*FON Declaration des parametres
c     --------------------------
c
      double precision day
      integer kode,lfile,nsat
      double precision x(kode),revnum
      integer ier 
c
c     ----------------------------------
c*FON Declaration des variables communes
c     ----------------------------------
c
c     indicateur de lecture pour le 1er appel 
      integer mfile
c     date debut et date fin de l'enreg courant
      double precision daybeg,dayend
c     date de debut du fichier - marge 
c     date de fin du bloc courant + marge
      double precision dayfir,daylas
c
      common/ssave/daybeg,dayend,dayfir,daylas,mfile
c
c     --------------------------------- 
c*FON Declaration des variables locales
c     --------------------------------- 
c
      integer nrec
      integer iter
      integer k,i,j,l
c*LOC k,i,j,l : indices de boucles 
c
      double precision epoch
c*LOC epoch : date de reference du vecteur d'etat
c
      double precision revepo
c*LOC revepo : numero d'orbite associe
c
      double precision smaxis
c*LOC smaxis : demi grand axe associe en km
c
      double precision omotin
c*LOC omotin : mouvement moyen inverse pour une orbite Keplerienne
c
      double precision coeff(10,6)
c*LOC coeff : table des coefficients du polynome de Chebyshev pour les
c*LOC         6 composants du vecteur d'etat
c
      integer koeff
c*LOC koeff : nombre de coefficients polynomiaux (0< koeff <10) 
c
      double precision y(6)
c*LOC y : vecteur d'etat
c
      double precision rdist
c*LOC rdist : distance centre terre satellite a la date courante en km
c
      double precision s
c*LOC s : coefficient poynomial de Chebyshev
c
      double precision dmanom
      double precision arin,arm,rvwam,tam,comp,b,go,g1,bet
      double precision d,g2,g3,fx,gx,rx,ft,gt
      double precision daymid
      double precision scale
      double precision pa,p,pb
c*LOC Variables de travail intermediaires
c
      SAVE
c
c     ---------------------------------
c*FON Affectation identificateur rcs_id
c     ---------------------------------
c
      data rcs_id /"
     >$Id$"/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
      if (kode .le. 0) goto 902
      if (kode .gt. 6) goto 902
c
c     cas de la 1er lecture                   
      if (lfile .ne. mfile) goto 10
c
c    ********************************************************
c*FON Recherche du bloc contenant la date passee en parametre
c    ********************************************************
c
c     -----------------------------------------------------------------
c*FON Gestion de la lecture en fonction de la date du dernier bloc lu.
c*FON Lecture des blocs anterieurs si la date est < date debut du bloc.
c*FON Lecture des blocs suivants si la date est > date fin du bloc.
c*FON Pas de lecture dans les autres cas.
c    ------------------------------------------------------------------
c
      if (day .gt. dayend + 1.d-4) goto 20
      if (day .ge. daybeg - 1.d-4) goto 70
c
c     ------------------------------------------------------------------
c*FON Initialisation de la lecture et positionnement en debut de fichier
c     ------------------------------------------------------------------
c
10    continue
      dayfir = 99.d9
      daylas = 99.d9
      mfile  = lfile
      rewind lfile
c
c     -------------------------------------------------------------
c*FON Recherche du 1er enregistrement du bloc (numero du satellite)
c     -------------------------------------------------------------
c
20    continue
      read (lfile,1000,err=901,end=100) nsat
      if (nsat .le. 0) goto 20
      if (nsat .gt. 4) goto 20
c
c*    -----------------------------
c*FON Lecture du 2nd enregistrenent
c*    -----------------------------
c*FON  structure de l'enregistrement :
c*FON  nrec = identification de l'enreg (doit etre egal a 200 + nsat)
c*FON  daybeg = date debut de l'enreg (date julienne)
c*FON  dayend = date fin de l'enreg (date julienne)
c*FON  epoch =  epoque de reference du vecteur etat (date julienne)
c*FON  revepo = numero d'orbite
c*FON  smaxis = demi grand axe en km
c*FON  omotin = mouvement moyen inverse pour une orbite Keplerienne
c
      read (lfile,1100,err=901,end=100)
     >      nrec,daybeg,dayend,epoch,revepo,smaxis,omotin
c
c     ----------------------------------
c*FON Verification du 2nd enregistrement
c     ----------------------------------
c
      if (nrec .ne. 200 + nsat) goto 900
      if (daybeg .gt. dayend) goto 900
c
c     ------------------------------------------------------------------
c*FON Controle de la date et de la continuite des enregistrements
c*FON La date doit etre > la + petite date deja lu - marge
c*FON L'ecart entre la date fin de l'enreg precedent et la date debut de 
c*FON l'enreg courant doit etre < gap
c*FON marge = 1.d-4 et gap = 2.d-4
c     ------------------------------------------------------------------
c
      dayfir = dmin1(dayfir,daybeg - 1.d-4)
      if (day .lt. dayfir) goto 905
      if (daybeg .gt. daylas) goto 903
c
      daylas = dayend + 2.d-4
c
c     ------------------------------------------------------------------
c*FON Gestion de la lecture en fonction de la date du bloc :
c*FON date  <  date debut du bloc - marge : lecture des blocs anterieurs
c*FON date  >  date fin du bloc + marge   : lecture des blocs suivants
c     ------------------------------------------------------------------
c
      if (day .gt. dayend + 1.d-4) goto 20
      if (day .lt. daybeg - 1.d-4) goto 10
c
c     -----------------------------
c*FON Lecture du 3ieme enregistrenent
c*    -----------------------------
c*FON  structure de l'enregistrement :
c*FON  nrec = id de l`enreg (doit etre egal a 300 + nb de coeff polynomiaux)
c*FON  y(6) = vecteur d'etat pour l'orbite de Kepler (km, km/s)
c*FON  rdist = distance centre terre satellite a la date courante en km
c
      read (lfile,1200,err=901,end=900) nrec,y,rdist
c
c     ------------------------------------
c*FON verification du 3ieme enregistrement
c     ------------------------------------
c
      if (nrec .gt. 310) goto 900
      if (nrec .lt. 300) goto 900
      koeff = nrec - 300
c
      if (koeff. le. 0) goto 70
c
c     ----------------------------------------------------------
c*FON Lecture des coefficients du polynome de Chebyshev pour les
c*FON 6 composants du vecteur d'etat
c*    ----------------------------------------------------------
c
      do 60 k = 1, koeff
         read (lfile,1300,err=901,end=900) nrec,(coeff(k,i),i=1,6)
c
c        ------------------------------------------------------
c*FON    Verification de l'enregistrement
c*FON    nrec = id de l`enreg (doit etre egal a koeff + 11 * k)  
c*       ------------------------------------------------------
c
         if (11 * k + koeff .ne. nrec) goto 900
60    continue
c
c*     -----------------
c*FON  Fin de la lecture
c*     -----------------
c
70    continue
c   
c     *****************************************
c*FON Debut de resolution de l'orbite de Kepler
c     *****************************************
c
c     -------------------------------------------------------
c*FON Convertion de la difference de date en anomalie moyenne
c*FON Calcul du numero d'orbite
c     -------------------------------------------------------
c
      dmanom = (day - epoch) * 864.d2 / omotin
      revnum = revepo + dmanom / 6.2831853072d0
c
      arin = smaxis / rdist
      arm  = (rdist - smaxis) / smaxis
      rvwam = (y(1) * y(4) + y(2) * y(5) + y(3) * y(6)) * omotin
     >        / smaxis**2
c
c*    --------------------------------------------------------
c*FON Calcul de l'anomalie excentrique par iteration de Newton
c*    --------------------------------------------------------
c
      tam = dmanom - rvwam
      comp = 1.d-7 + 1.d-10 * dabs(tam)
      b = tam
c
c     ---------------------------------------------
c*FON Iterations pour resoudre l'equation de Kepler
c     ---------------------------------------------
c
      do 130 iter = 1, 15
         go  = cos(b)
         g1  = sin(b)
         bet = tam - arm * g1 + rvwam * go
         d   = (bet - b) / (1.d0 + arm * go + rvwam * g1)
         b   = b + d
c
c        -----------------
c*FON    Condition d'arret 
c        -----------------
c
         if (dabs(d) .le. comp) goto 140
130   continue
c
c     ------------------
c*FON Pas de convergence       
c     ------------------
c
      goto 900
140   continue
      go = go - d * g1
      g1 = g1 + d * go
      g2 = 1.d0 - go
      g3 = b - g1
      fx = 1.d0  - g2 * arin
      gx = (dmanom - g3) * omotin
c
      k = min0(kode,3)
      do 150 j = 1, k
         x(j) = fx * y(j) + gx * y(j+3)
150   continue
c
      if (kode .le. 3) goto 170
      rx = sqrt(x(1)**2 + x(2)**2 + x(3)**2)
      ft = -g1 * smaxis * arin / (omotin * rx)
      gt = 1.d0 - g2 * smaxis / rx
      do 160 j = 4, kode
         x(j) = ft * y(j-3) + gt * y(j)
160   continue
c
c*    ******************************************
c*FON Fin de la resolution de l'orbite de Kepler
c*    ******************************************
c
170   continue
c
c     ---------------------------------------------------------
c*FON Verification des coeff polynomiaux (1 n'est pas autorise)
c     ---------------------------------------------------------
c
      if (koeff .le. 1) goto 600
c
c     ----------------------------------------------------------------------
c*FON Date milieu et vecteur d'echelle pour le coeff polynomial de Chebyshev
c     ----------------------------------------------------------------------
c
      daymid = 0.5d0 * (daybeg + dayend)
      scale  = 4.d0 / (dayend - daybeg)
      s      = scale * (day - daymid)
      pa     = 1.d0
      p      = s * 0.5d0
c
c     --------------------------------
c*FON Initialisation du vecteur d'etat 
c     --------------------------------
c
      do 200 j = 1, kode
         x(j) = x(j) + coeff(1,j) + coeff(2,j) * p
200   continue
c
      if (koeff .le. 2) goto 600
      do 210 l = 3, koeff
         pb = pa
         pa = p
         p  = s * pa - pb
         do 220 j = 1, kode
            x(j) = x(j) + coeff(l,j)*p
220      continue
210   continue
c
      goto 600
c
c     -------------------
c*FON Gestion des erreurs
c     -------------------
c
100   continue
c
c     ---------------------------------------------------
c*FON Cas de la fin de fichier : 
c*FON Cas de date trop grande si le 1er enreg est deja lu 
c*FON Donnees invalide dans les autres cas
c     ---------------------------------------------------
c
      if (daylas .lt. 1.d9) goto 904
c
900   continue
      ier = -1
901   continue
      ier = ier + 2
902   continue
      ier = ier + 1
903   continue
      ier = ier + 1
904   continue
      ier = ier + 1
905   continue
      ier = ier + 1
c
c     -----------------------------------------------------------
c*FON Re_initialisation de mfile pour forcer une lecture conplete
c*FON au prochain appel
c     -----------------------------------------------------------
c
      mfile = -9999
600   continue
c
c     ****************
c     Fin de programme
c     ****************
c
      return
1000  format(i3)
1100  format(i3,2f12.6,f15.9,f11.3,2f13.5)
1200  format(i3,3f11.3,3f11.7,f11.3)
1300  format(i3,3f11.3,3f11.7)
      end