PROGRAM URANUS_EPH

  IMPLICIT NONE

  REAL(8):: DJJ
  REAL(8), DIMENSION(3):: XYZ
  REAL(8), DIMENSION(300,3,73,5):: coef
  REAL(8), DIMENSION(73,5):: a, b

  real(8), dimension(6):: XV
  integer:: ns

  INTEGER:: nsat, i, j, loop, isat

  !Mise en memoire des coef de CHEB

  OPEN(31,file='NOE-coefcheb1.res'); OPEN(32,file='NOE-coefcheb2.res')
  OPEN(33,file='NOE-coefcheb3.res'); OPEN(34,file='NOE-coefcheb4.res')
  OPEN(35,file='NOE-coefcheb5.res')

  DO isat=1,5 !nsat  
     DO loop=1,73
        READ(30+isat,*) a(loop,isat), b(loop,isat)
        DO i=1,3
           READ(30+isat,*) coef(:,i,loop,isat)
        ENDDO !i
     ENDDO !loop
  ENDDO !isat

  CLOSE(31); CLOSE(32); CLOSE(33); CLOSE(34); CLOSE(35)
  write(*,*)"fin de mise en memoire!"

  !Fin de mise en memoire

    WRITE(*,*)"Entrez la date desiree (en jours juliens TDB)"
    READ(*,*) DJJ
    WRITE(*,*)"Entrez le numero du satellite (1=Ariel, 2=Umbriel, 3=Titania, 4=Oberon, 5=miranda)"
    READ(*,*) isat

    CALL EPHM(DJJ, isat, XYZ)

    WRITE(*,*)"Voici les positions en km en repere Uranocentrique (attention ca n'est pas le barycentre du systeme!)  & 
         equateur moyen terrestre J2000. "
    WRITE(*,*) XYZ

CONTAINS

  SUBROUTINE EPHM(DJJ, isat, XYZ)

    !********************************************************************
    !*          Ephemerides preliminaires des satellites d'Uranus       *
    !*                (ajustement entre 1954 et 2003)                   *
    !*                                                                  *
    !*                        (Version 1.3)                             *
    !*                                                                  *
    !*------------------------------------------------------------------*
    !*                   V.Lainey (Lainey@oma.be)                       *
    !*                                                                  *
    !*    ROB- 3, Avenue Circulaire, B-1180 Bruxelles (Belgium)         *
    !*------------------------------------------------------------------*
    !*    IMCCE - 77, Avenue Denfert-Rochereau 75014 Paris (France)     *
    !*__________________________________________________________________*
    !*                                                                  *
    !*  DJJ: jour julien (TDB)                                          *
    !*                                                                  *
    !* isat: numero du satellite                                        *
    !*       (1=Ariel, 2=Umbriel, 3=Titania, 4=Oberon, 5=Miranda)       *
    !*                                                                  *
    !*  XYZ: (x,y,z) Uranocentrique J2000 (equateur terrestre moyen)    *
    !*------------------------------------------------------------------*
    !* references: V.Lainey (priv. com.)                                *
    !********************************************************************
    !*                          21/02/2006                              *
    !********************************************************************

    IMPLICIT NONE

    REAL(8), INTENT(IN):: DJJ
    INTEGER, INTENT(IN):: isat

    REAL(8), DIMENSION(3), INTENT(OUT):: XYZ

    REAL(8):: date, tmp, x, vals
    REAL(8), DIMENSION(0:299):: TN
    INTEGER:: loop, i, k, ns

    date=DJJ-2453736.5d0 ! 01/01/2006 0H00

    IF((date<0.d0).OR.(date>2189.90875d0))THEN; WRITE(*,*)"Erreur de date: ", date; STOP
    ENDIF !dateMax: 30 decembre 2011 a 21H48 36 sec (TDB)

    loop=INT(date/29.99875d0)+1

    x=(date-0.5d0*(b(loop,isat)+a(loop,isat)))/(0.5d0*(b(loop,isat)-a(loop,isat)))
    TN(0)=1.d0
    TN(1)=x

    DO i=1,3

       vals=0.d0

       DO k=2,299 !ndegre-1
          TN(k)=2.d0*x*TN(k-1)-TN(k-2)
       ENDDO !k

       DO k=0,299 !ndegre-1
          vals=vals+coef(k+1,i,loop,isat)*TN(k)
       ENDDO !k
       XYZ(i)=vals-0.5d0*coef(1,i,loop,isat)

    ENDDO !i

  END SUBROUTINE EPHM

END PROGRAM URANUS_EPH
