! BEGIN cluster_orbitlib.f
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
! |                                                                    |
! | Modules permettants la lecture d'un fichier d'orbite de CLUSTER    |
! | et l'extraction des la position et de la vitesse de chaque s/s     |
! | dans le repere GEI.                                                |
! | Ces modules ont ete donne par l'ESA (voir le DDID).                |
! | Ils sont completes par quelques utilitaire pour interfacer la      |
! | bibli qualifaclib.                                                 |
! |                                                                    |
! |                 P. Robert, CNRS/CETP,  1996                        |
! |                              Revu Mars 2001                        |
! |                              Revu Aout 2002 (dirpath impose)       |
! |                              Revu Janv 2004 (diag. erreurs)        |
! |                                                                    |
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
!
      subroutine openorbi(dirpath,iyear,imon,iday,ifc)

      integer ifc(4)
      character*(*) dirpath
      character*25 ficorbi1,ficorbi2,ficorbi3,ficorbi4
      character*160 filestr
!
! *** path des fichiers d'orbite
!
!
!     ---------------------------------------------------------------+--
!     calcul des noms des fichiers orbit et ouverture
!
!     P. Robert, CETP, Juin 2001
!     ---------------------------------------------------------------+--
!
!
! *** calcul du nom des fichiers
!
      write(ficorbi1,210) iyear,imon,iyear-2000,imon,iday,1
      write(ficorbi2,210) iyear,imon,iyear-2000,imon,iday,2
      write(ficorbi3,210) iyear,imon,iyear-2000,imon,iday,3
      write(ficorbi4,210) iyear,imon,iyear-2000,imon,iday,4
!
  210 format(i4.4,'_',i2.2,'/',3i2.2,'_sat',i1,'.orb')
!
      print*
      print*, 'openorbi: dir  =', trim(dirpath)
      print*, '          files=', ficorbi1
      print*, '                ', ficorbi2
      print*, '                ', ficorbi3
      print*, '                ', ficorbi4
!
! *** ouverture des fichiers d'orbite
!
      close(ifc(1))
      close(ifc(2))
      close(ifc(3))
      close(ifc(4))
!
      call longmot(dirpath,nd)
!
      filestr=dirpath(1:nd)//'/'//ficorbi1
      open(ifc(1),file=filestr,status='old',err=91)
      filestr=dirpath(1:nd)//'/'//ficorbi2
      open(ifc(2),file=filestr,status='old',err=92)
      filestr=dirpath(1:nd)//'/'//ficorbi3
      open(ifc(3),file=filestr,status='old',err=93)
      filestr=dirpath(1:nd)//'/'//ficorbi4
      open(ifc(4),file=filestr,status='old',err=94)
!
      return
!
   91 continue
      dirpath=dirpath//ficorbi1
      print*,'***/ CAN''T OPEN FILE :',trim(dirpath)
      stop 'CLUPOSVIT      *** ABORTED ! NO ORBIT 1 FILE            ***'
   92 continue
      dirpath=dirpath//ficorbi2
      print*,'***/ CAN''T OPEN FILE :',trim(dirpath)
      stop 'CLUPOSVIT      *** ABORTED ! NO ORBIT 2 FILE            ***'
   93 continue
      dirpath=dirpath//ficorbi3
      print*,'***/ CAN''T OPEN FILE :',trim(dirpath)
      stop 'CLUPOSVIT      *** ABORTED ! NO ORBIT 3 FILE            ***'
   94 continue
      dirpath=dirpath//ficorbi4
      print*,'***/ CAN''T OPEN FILE :',trim(dirpath)
      stop 'CLUPOSVIT      *** ABORTED ! NO ORBIT 4 FILE            ***'
!
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cluposvit(iyear,imon,iday,ih,im,is,ifc,pos,vit,rev)
!
! *** declaration des variables
!
      integer ifc(4)
      real*8 djday,revn,posvit(6),sec
      real   pos(3,4),vit(3,4),rev(4)
      character*27 cerror(6)
!
      data cerror /'1="DAY" TOO EARLY', '2="DAY" TOO LATE', &
                   '3=TIME GAP IN DATA', '4=WRONG VALUE OF "KODE"', &
                   '5=FILE CONTENT INCONSISTENT', &
                   '6=READ ERROR FROM DATA FILE'/
!
!     ---------------------------------------------------------------+--
!     calcul les positions et vitesse des 4 sat de CLUSTER
!
!     P. Robert, CETP, Juin 2001
!     ---------------------------------------------------------------+--
!
      sec= dble(is)
      kode= 6
!
      call jd2000(djday,iyear,imon,iday,ih,im,sec)
!
! *** boucle sur les 4 satellites
!
      do 30 isat=1,4
!
      call orbit(djday,kode,ifc(isat),ierr,numsat,posvit,revn)
!
      if(ierr.ne.0) then
                    print*, '*** cluposvit: sat#',isat,' at HH MM SS=', &
                            ih,im,is, ' error ',cerror(ierr)
                    print*, '    position, velocity and rev. number ', &
                            'set to zero'
                    do ix=1,3
                    pos(ix,isat)=0.
                    vit(ix,isat)=0.
                    enddo
                    rev(isat)=0.
                    return
                    endif
!
! *   boucle sur les composantes pour le chargement de pos
!
      do 20 ix=1,3
      pos(ix,isat)= real(posvit(ix))
      vit(ix,isat)= real(posvit(ix+3))
   20 continue
!
      rev(isat)=real(revn)
!
   30 continue
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine closeorbi(ifc)
!
      integer ifc(4)
!
!     ---------------------------------------------------------------+--
!     fermeture des fichiers d'orbite
!
!     P. Robert, CETP, Juin 2001
!     ---------------------------------------------------------------+--
!
! *** fermeture des fichiers
!
      close(ifc(1))
      close(ifc(2))
      close(ifc(3))
      close(ifc(4))
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cosatpos(dmjdc20,ic1,ic2,ic3,ic4,satpos,satvit,satrev)
!
      real*8 dgei1(6),dgei2(6),dgei3(6),dgei4(6)
      real*8 drev1   ,drev2   ,drev3   ,drev4
      real*8 dmjdc20
!
      real   satpos(3,4),satvit(3,4),satrev(4)
!
!     ******************************************************************
!     calcule la position et la vitesse des satellites dans le GEI
!     ******************************************************************
!
      data ns1,ns2,ns3,ns4 /1,2,3,4/
!
      call orbit(dmjdc20,6,ic1,ierr1,ns1,dgei1,drev1)
      call orbit(dmjdc20,6,ic2,ierr2,ns2,dgei2,drev2)
      call orbit(dmjdc20,6,ic3,ierr3,ns3,dgei3,drev3)
      call orbit(dmjdc20,6,ic4,ierr4,ns4,dgei4,drev4)
!
      if(ierr1.ne.0.or. &
         ierr2.ne.0.or. &
         ierr3.ne.0.or. &
         ierr4.ne.0)  print*, 'erreur=', ierr1,ierr2,ierr3,ierr4
!
      do 10 i=1,3
      satpos(i,1)=real(dgei1(i))
      satpos(i,2)=real(dgei2(i))
      satpos(i,3)=real(dgei3(i))
      satpos(i,4)=real(dgei4(i))
   10 continue
!
      do 20 i=4,6
      satvit(i,1)=real(dgei1(i))
      satvit(i,2)=real(dgei2(i))
      satvit(i,3)=real(dgei3(i))
      satvit(i,4)=real(dgei4(i))
   20 continue
!
      satrev(1)=real(drev1)
      satrev(2)=real(drev2)
      satrev(3)=real(drev3)
      satrev(4)=real(drev4)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine corevmoy(satrev,revmoy)
!
      real satrev(4)
!
!     ******************************************************************
!     compute le revolution number moyen
!     ******************************************************************
!
      revmoy=0.
!
      do 10 j=1,4
      revmoy=revmoy+abs(satrev(j))
   10 continue
      revmoy=revmoy/4.
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cvsatpos(satpos)
!
      real satpos(3,4)
!
!     ******************************************************************
!     convertit satellite position de km en RT
!     ******************************************************************
!
      do 10 i=1,3
      do 10 j=1,4
      satpos(i,j)=satpos(i,j)/6378.
   10 continue
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cbaryclus(s,gx,gy,gz)
!
! **********************************************************************
! compute barycentre of 4 points of cluster
! **********************************************************************
!
!
      real s(3,4)
!
      gx=(s(1,1)+s(1,2)+s(1,3)+s(1,4))/4.
      gy=(s(2,1)+s(2,2)+s(2,3)+s(2,4))/4.
      gz=(s(3,1)+s(3,2)+s(3,3)+s(3,4))/4.
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine transclus(s,gx,gy,gz)
!
! **********************************************************************
! tranlate position of the 4 points of cluster of gx,gy,gz
! **********************************************************************
!
!
      real s(3,4)
!
      do 10 j=1,4
      s(1,j)=s(1,j)-gx
      s(2,j)=s(2,j)-gy
      s(3,j)=s(3,j)-gz
   10 continue
!
      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine gei_to_gse(pos)
!
      real pos(3,4)
!
! *** transforme pos(3,4) de gei vers gse
!
      do 10 isat=1,4
!
      call tgeigse(pos(1,isat),pos(2,isat),pos(3,isat),x,y,z)
!
      pos(1,isat)= x
      pos(2,isat)= y
      pos(3,isat)= z
!
   10 continue
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine gei_to_gsm(pos)
!
      real pos(3,4)
!
! *** transforme pos(3,4) de gei vers gsm
!
      do 10 isat=1,4
!
      call tgeigsm(pos(1,isat),pos(2,isat),pos(3,isat),x,y,z)
!
      pos(1,isat)= x
      pos(2,isat)= y
      pos(3,isat)= z
!
   10 continue
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine longmot(mot,nc)
!
!     ------------------------------------------------------------------
! *   Objet  : calcul de la longueur utile d'une chaine de caracteres
! *   Classe : depouillement specifique GEOS/UBF
! *   Auteur : P. Robert, CETP, 1998
!     ------------------------------------------------------------------
!
      character mot*(*)
!
!                    *********************
!
      nbc=len(mot)
!
! *** nc= longueur au de la de laquelle il n'y a que des blancs
!         (mais la chaine 1:nc peu contenir des blancs)
!
      do 10 i=1,nbc
      ii=nbc-i+1
      if(mot(ii:ii).ne.' ') go to 12
   10 continue
      ii=0
   12 continue
      nc=ii
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      SUBROUTINE ORBIT(DAY,KODE,LFILE,IERROR,NSAT,X,REVNUM)
!P    ORBIT: RETRIEVAL ROUTINE FOR COMPRESSED CLUSTER ORBIT
!
!     INPUT:
!I    DAY (R*8) = MODIFIED JULIAN DAY, FROM 2000, FOR THE STATE VECTOR
!I    KODE (I*4) = NUMBER OF COMPONENTS OF STATE VECTOR = DIM. OF ARRAY
!              X(); = 3 FOR S/C POSITION, = 6 FOR POSITION & VELOCITY
!I    LFILE (I*4) = LOGICAL NUMBER OF INPUT DATA FILE
!     OUTPUT:
!O    IERROR (I*4) = RETURN CODE: 0=NO ERROR, 1='DAY' TOO EARLY, 2=TOO
!            LATE, 3=TIME GAP IN DATA, 4=WRONG VALUE OF 'KODE',
!            5=FILE CONTENT INCONSISTENT, 6=READ ERROR FROM DATA FILE
!O    NSAT (I*4) = SATELLITE NUMBER; 1, 2, 3, 4
!O    X(KODE) (R*8) = SPACECRAFT POSITION, KM (AND VELOCITY, KM/S)
!O    REVNUM (R*8) = REVOLUTION NUMBER
!
!F    READS A SEQUENTIAL FORMATTED FILE WITH LOGICAL NUMBER 'LFILE'
!
!lm
!lm   Re-written by Laurent Mirioni (CETP) on april, 2004 because of
!lm   "Irreductible loop" detected by Foresys.
!lm
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION Y(6),COEFF(10,6),X(KODE)
!     INITIALISE FILE NUMBER TO FORCE FILE READING AT FIRST CALL
      DATA MFILE/-9999/,DAYFIR/99.D9/,DAYLAS/99.D9/,DAYEND/0.D0/, &
           DAYBEG/0.D0/
      SAVE
!     ... EXECUTABLE STATEMENTS ...
!
!
!     INITIALISE ERROR CODES
      IERROR = 0
      IF (KODE.GT.0 .AND. KODE.LE.6) THEN
!
!        ALWAYS REWIND IF A NEW FILE NUMBER IS USED
         IF (LFILE .EQ. MFILE) THEN
!
!           CHECK IF 'DAY' IS INSIDE LAST READ RECORD BLOCK
            IF (DAY .GT. DAYEND+1.D-4) THEN
  100          CONTINUE
!                 READ 1ST RECORD IN A BLOCK
                  READ (LFILE,10000,ERR = 1500,END = 1300) NSAT
!                 IF: NSTA = A SATELLITE NUMBER; THEN THIS IS 1ST RECORD IN
!                 A BLOCK
                  IF (NSAT.GT.0 .AND. NSAT.LE.4) THEN
!
!                    READ 2ND RECORD IN THE BLOCK
!F                   NREC = RECORD IDENTIFICATION, SHALL BE = 200 + NSAT
!F                   DAYBEG = BEGIN TIME OF THE RECORD (MJD)
!F                   DAYEND = END TIME OF THE RECORD (MJD)
!F                   EPOCH = EPOCH OF REFERENCE STATE VECTOR (MJD)
!F                   REVEPO = REVOLUTION NUMBER AT EPOCH
!F                   SMAXIS = SEMIMAJOR AXIS FOR THE KEPLER ORBIT
!F                   OMOTIN = INVERSE MEAN MOTION FOR THE KEPLER ORBIT
                     READ (LFILE,10001,ERR = 1500,END = 1300) &
                     NREC, DAYBEG, DAYEND, EPOCH, REVEPO, SMAXIS, OMOTIN
!                    CHECK CONSISTENCY OF FILE
                     IF (NREC .NE. 200+NSAT) GOTO 1400
                     IF (DAYBEG .GT. DAYEND) GOTO 1400
!
!                    DAYFIR = START TIME OF 1ST RECORD ON 1ST BLOCK ON
!                    THE FILE - MARGIN
                     DAYFIR = MIN(DAYFIR,DAYBEG-1.D-4)
!                    ERROR RETURN IF 'DAY' IS BEFORE START OF FILE (WITH MARGIN)
                     IF (DAY .LT. DAYFIR) GOTO 1800
!                    ERROR RETURN IF THERE IS A GAP FROM LAST BLOCK
!                    (WITH MARGIN)
                     IF (DAYBEG .GT. DAYLAS) GOTO 1600
!                    DAYLAS = END TIME OF LAST READ RECORD BLOCK + MARGIN
                     DAYLAS = DAYEND + 2.D-4
!
!                    CONTINUE READ IF 'DAY' IS AFTER END OF THIS RECORD BLOCK
                     IF (DAY .LE. DAYEND+1.D-4) GOTO 200
                  END IF
               GOTO 100
!              REWIND WHEN 'DAY' IS EARLIER THAN START OF PRESENT RECORD BLOCK
  200          IF (DAY .GE. DAYBEG-1.D-4) GOTO 1000
            ELSEIF (DAY .GE. DAYBEG-1.D-4) THEN
               GOTO 1100
            END IF
         END IF
  300    CONTINUE
!
!           INITIALISE THE READING FROM THE FILE
            DAYFIR = 99.D9
            DAYLAS = 99.D9
            MFILE = LFILE
            REWIND (LFILE)
  400       CONTINUE
!              READ 1ST RECORD IN A BLOCK
!
               READ (LFILE,10000,ERR = 1500,END = 1300) NSAT
!              IF: NSTA = A SATELLITE NUMBER; THEN THIS IS 1ST RECORD IN A BLOCK
               IF (NSAT.GT.0 .AND. NSAT.LE.4) THEN
!
!                 READ 2ND RECORD IN THE BLOCK
!F                NREC = RECORD IDENTIFICATION, SHALL BE = 200 + NSAT
!F                DAYBEG = BEGIN TIME OF THE RECORD (MJD)
!F                DAYEND = END TIME OF THE RECORD (MJD)
!F                EPOCH = EPOCH OF REFERENCE STATE VECTOR (MJD)
!F                REVEPO = REVOLUTION NUMBER AT EPOCH
!F                SMAXIS = SEMIMAJOR AXIS FOR THE KEPLER ORBIT
!F                OMOTIN = INVERSE MEAN MOTION FOR THE KEPLER ORBIT
                  READ (LFILE,10001,ERR = 1500,END = 1300) &
                  NREC, DAYBEG, DAYEND, EPOCH, REVEPO, SMAXIS, OMOTIN
!                 CHECK CONSISTENCY OF FILE
                  IF (NREC .NE. 200+NSAT) GOTO 1400
                  IF (DAYBEG .GT. DAYEND) GOTO 1400
!
!                 DAYFIR = START TIME OF 1ST RECORD ON 1ST BLOCK ON
!                 THE FILE - MARGIN
                  DAYFIR = MIN(DAYFIR,DAYBEG-1.D-4)
!                 ERROR RETURN IF 'DAY' IS BEFORE START OF FILE (WITH MARGIN)
                  IF (DAY .LT. DAYFIR) GOTO 1800
!                 ERROR RETURN IF THERE IS A GAP FROM LAST BLOCK (WITH MARGIN)
                  IF (DAYBEG .GT. DAYLAS) GOTO 1600
!                 DAYLAS = END TIME OF LAST READ RECORD BLOCK + MARGIN
                  DAYLAS = DAYEND + 2.D-4
!
!                 CONTINUE READ IF 'DAY' IS AFTER END OF THIS RECORD BLOCK
                  IF (DAY .LE. DAYEND+1.D-4) GOTO 500
               END IF
            GOTO 400
!           REWIND WHEN 'DAY' IS EARLIER THAN START OF PRESENT RECORD BLOCK
  500       IF (DAY .GE. DAYBEG-1.D-4) GOTO 1000
         GOTO 300
!
!        READ 3RD RECORD IN THE BLOCK
!F       NREC = RECORD IDENTIFICATION, SHALL BE = 300 + NUMBER OF POL.COEFF.
!F       Y(6) = REFERENCE STATE VECTOR FOR KEPLER ORBIT (KM, KM/S)
!F       RDIST = S/C EARTH CENTRE DISTANCE AT EPOCH
 1000    READ (LFILE,10002,ERR = 1500,END = 1400) NREC, Y, RDIST
!
!        CHECK CONSISTENCY OF FILE
         IF (NREC .GT. 310) GOTO 1400
         IF (NREC .LT. 300) GOTO 1400
!        KOEFF = NUMBER OF POLYNOMIAL COEFFICIENTS, BETWEEN 0 AND 10
         KOEFF = NREC - 300
!
!        IF THERE ARE NO COEFFICIENTS IN THIS BLOCK
         IF (KOEFF .GT. 0) THEN
            DO K = 1,KOEFF
!F             NREC = RECORD IDENT. = KOEFF + 11*K
!F             COEFF(10,6) = MATRIX WITH UP TO 10 COEFFICIENTS OF THE CHEBYSHEV
!F             POLYNOMIAL FOR EACH OF THE 6 COMPONENTS OF THE STATE VECTOR
               READ (LFILE,10003,ERR=1500,END=1400) &
               NREC,(COEFF(K,I), I = 1,6)
!
!              CHECK CONSISTENCY OF FILE
               IF (11*K+KOEFF .NE. NREC) GOTO 1400
            END DO
         END IF
!
!        TIME CONVERTED TO DIFFERENCE IN MEAN ANOMALY
!        END OF BLOCK READING SEQUENCE
 1100    DMANOM = (DAY-EPOCH) * 864.D2 / OMOTIN
!        ORBIT NUMBER
         REVNUM = REVEPO + DMANOM/6.2831853072D0
!
!        START MODELLING KEPLER ORBIT
         ARIN = SMAXIS / RDIST
         ARM = (RDIST-SMAXIS) / SMAXIS
         RVWAM = (Y(1)*Y(4)+Y(2)*Y(5)+Y(3)*Y(6)) * OMOTIN / SMAXIS**2
!        CALC. OF ECC. ANOMALY BY NEWTON'S ITERATION
         TAM = DMANOM - RVWAM
         COMP = 1.D-7 + 1.D-10*ABS(TAM)
         B = TAM
!        ITERATIONS TO SOLVE KEPLER'S EQUATION:
         DO 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
!           THIS GIVES THE ACCURACY      1.D-14 IN B & THE G'S
            IF (ABS(D) .LE. COMP) GOTO 1200
         END DO
         GOTO 1400
 1200    GO = GO - D*G1
         G1 = G1 + D*GO
         G2 = 1.D0 - GO
         G3 = B - G1
         FX = 1.D0 - G2*ARIN
         GX = (DMANOM-G3) * OMOTIN
!
         K = MIN(KODE,3)
         DO J = 1,K
            X(J) = FX*Y(J) + GX*Y(J+3)
         END DO
!
         IF (KODE .GT. 3) THEN
            RX = SQRT(X(1)**2+X(2)**2+X(3)**2)
            FT = -G1*SMAXIS*ARIN/(OMOTIN*RX)
            GT = 1.D0 - G2*SMAXIS/RX
            DO J = 4,KODE
               X(J) = FT*Y(J-3) + GT*Y(J)
            END DO
         END IF
!        END OF MODELLING KEPLER ORBIT
!
!        CHECK IF POLYNOMIAL COEFFICIENTS ARE REQUIRED (1 IS NOT WORTH WHILE)
         IF (KOEFF .LE. 1) RETURN
!        MID-POINT & SCALE FACTOR FOR CHEBYSHEV POLYNOMIAL
         DAYMID = 0.5D0 * (DAYBEG+DAYEND)
         SCALE = 4.D0 / (DAYEND-DAYBEG)
!        ADD CHEBYSHEV POLYNOMIAL TO KEPLER STATE VECTOR
         S = SCALE * (DAY-DAYMID)
         PA = 1.D0
         P = S * 0.5D0
!
!        'KODE' = NUMBER OF COMPONENTS OF THE STATE VECTOR
         DO J = 1,KODE
            X(J) = X(J) + COEFF(1,J) + COEFF(2,J)*P
         END DO
!
         IF (KOEFF .LE. 2) RETURN
         DO L = 3,KOEFF
            PB = PA
            PA = P
            P = S*PA - PB
            DO J = 1,KODE
               X(J) = X(J) + COEFF(L,J)*P
            END DO
         END DO
         RETURN
!        END-OF-FILE ONLY IF AT LEAST ONE RECORD HAS BEEN READ
!        ERROR RETURNS; IERROR = 5, 6, 4, 3, 2 OR 1
 1300    IF (DAYLAS .LT. 1.D9) GOTO 1700
 1400    IERROR = -1
 1500    IERROR = IERROR + 2
      END IF
      IERROR = IERROR + 1
 1600 IERROR = IERROR + 1
 1700 IERROR = IERROR + 1
 1800 IERROR = IERROR + 1
!     FORCE A RE-INITIALISATION OF READ AT NEXT CALL AFTER AN ERROR
      MFILE = -9999
      RETURN
!
!     ... FORMAT DECLARATIONS ...
!
10000 FORMAT (I3)
10001 FORMAT (I3,2F12.6,F15.9,F11.3,2F13.5)
10002 FORMAT (I3,3F11.3,3F11.7,F11.3)
10003 FORMAT (I3,3F11.3,3F11.7)
      END

!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      SUBROUTINE PR2000(DAY,P)
!P  COMPUTES THE PRECESSION MATRIX P(3,3) FOR CONVERTING A VECTOR
! IN MEAN GEOCENTRIC EQUATORIAL SYSTEM OF 2000.0 TO MEAN-OF-DATE.
! REF: THE ASTRONOMICAL ALMANAC 1985 PAGE B18.
!
!INPUT:  DAY = MJD2000 = MOD. JULIAN DAY FOR THE MEAN-OF-DATE SYSTEM
!            = MJD(1950) - 18262.0
!
!OUTPUT: P(3,3) = PRECESSION MATRIX FOR THE TRANSFORMATION:
!     R(MEAN-OF-DATE) = P(,)*R(2000)
!
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION P(3,3)
!
! CONVERT TO STANDARD EPOCH J2000.0 = 2000 JAN 1 AT 12:00:00
      T = DAY - 0.5D0
!
!  GZ=GREEK Z(A), ZA=Z(A), TH=THETA, ACCORDING TO THE REFERENCE.
! ORIGINAL, WITH TJC = (DAY-0.5D0)/36525.D0  IN JULIAN CENTURIES:
!     GZ = RAD*TJC*(0.6406161D0 + TJC*(839.D-7 + TJC*5.D-6))
!     ZA = GZ + RAD*TJC*TJC*(2202.D-7 + TJC*1.D-7)
!     TH = RAD*TJC*(0.5567530D0 - TJC*(1185.D-7 + TJC*116.D-7))
!
      GZ = T*(0.3061153D-6 + T*(0.10976D-14 + T*0.179D-20))
      ZA = GZ + T*T*(0.2881D-14 + T*0.358D-22)
      TH = T*(0.2660417D-6 - T*(0.1550D-14 + T*0.41549D-20))
!
      CGZ=DCOS(GZ)
      SGZ=DSIN(GZ)
      CZA=DCOS(ZA)
      SZA=DSIN(ZA)
      CTH=DCOS(TH)
      STH=DSIN(TH)
      P(1,1) = CGZ*CZA*CTH - SGZ*SZA
      P(1,2) = -SGZ*CZA*CTH - CGZ*SZA
      P(1,3) = -CZA*STH
      P(2,1) = CGZ*SZA*CTH + SGZ*CZA
      P(2,2) = -SGZ*SZA*CTH + CGZ*CZA
      P(2,3) = -SZA*STH
      P(3,1) = CGZ*STH
      P(3,2) = -SGZ*STH
      P(3,3) = CTH
      RETURN
      END
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      SUBROUTINE DJ2000(DAY,I,J,K,JHR,MI,SEC)
!P  COMPUTES CALENDER DATE FROM MODIFIED JULIAN DAY 2000
!   VALID FOR DATES BETWEEN 1950/JAN/1 AND 2099/DEC/31.
!   MJD(2000) = MJD(1950) - 18262.0 IS = 0 ON 2000/01/01 AT 00:00:00.
!
!I  (REAL*8) DAY = MOD. JULIAN DAY, REFERRED TO 2000 (MAY BE NEGATIVE).
!O  (INTEGERS): I=YEAR, J=MONTH, K=DAY, JHR=HOUR, MI=MINUTE
!O  (REAL*8): SEC=SECOND.
!
      IMPLICIT REAL*8(A-H,O-Z)
!  MAKE SURE TO ROUND-OFF ONLY DOWN, ALSO FOR NEGATIVE MJD:
      JDAY = INT(DAY + 18262.D0)
      L = (4000*(JDAY + 18204))/1461001
      N = JDAY - (1461*L)/4 + 18234
      M = (80*N)/2447
      K = N - (2447*M)/80
      JJ = M/11
      J = M + 2 - 12*JJ
      I = 1900 + L + JJ
      SEC = (DAY - DFLOAT(JDAY-18262))*24.D0
      JHR = INT(SEC)
      SEC = (SEC - DFLOAT(JHR))*6.D1
      MI = INT(SEC)
      SEC = (SEC - DFLOAT(MI))*6.D1
      RETURN
      END
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      SUBROUTINE JD2000(DAY,JEAR,MONTH,KDAY,JHR,MI,SEC)
!P GIVES THE NEW MOD. JULIAN DAY (MJD=0.0 ON 2000/JAN/1 AT 0:00:00)
!P FOR INPUT CALENDAR DATES BETWEEN 1950/JAN/1 AND 2099/DEC/31.
!
!   MJD(2000) = MJD(1950) - 18262.0
!
!I  (INT*4) JEAR = YEAR WITH 2 OR 4 DIGITS; 2 DIGITS => 1950 TO 2049
!I  (INT*4) MONTH = MONTH
!I  (INT*4) KDAY = DAY
!I  (INT*4) JHR = HOUR
!I  (INT*4) MI = MINUTE
!I  (REAL*8) SEC = SECOND.
!O  (REAL*8) DAY = MOD. JUL. DAY, REFERRED TO 2000.
!
      IMPLICIT REAL*8(A-H,O-Z)
      JJ = (14 - MONTH)/12
      L = JEAR - JJ - 1900*(JEAR/1900) + 100*(2000/(JEAR+1951))
      IDAY=KDAY-36496+(1461*L)/4+(367*(MONTH-2+JJ*12))/12
      DAY = DFLOAT(IDAY)
      DAY = DAY + (DFLOAT((JHR*60 + MI)*60) + SEC)/864.D2
      RETURN
      END
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
!
