*      ------------------------------------------------------------
*	compute XYZ from direct access files files
*       open(10,access='direct',recl=284,file='xyzj.res')

*       version 2 (december 1998)
*      ------------------------------------------------------------

	implicit double precision (a-h,o-z)
	character*40 infile

	dimension cmx(6),cfx(4),cmy(6),cfy(4),cmz(6),cfz(4)
	dimension tbx(1000),tby(1000),tbz(1000)

***     4 sat. >> 8 sat.
	dimension freq(8),idn(8),dellt(8)
	


	write(6,'(a,$)') ' Name of the file of coefficients '
	read(5,'(a)') infile
              open(10,access='direct',recl=284,file=infile)

	write(*,*) 
     . ' date of the beginning of the ephemeris (year month day h m s)'
	read(*,*)  jan,jmois,jour,jh,jmn,s
	write(*,*) jan,jmois,jour,jh,jmn,s
	write(*,*) ' step  (day h m s)'
	read(*,*) kpj,kph,kpm,psec
	write(*,*) ' steps number '
	read(*,*) ipas

	ua=149597870.d0
	itour=0
***	nsat=4
	id=1
	read(10,rec=id) npla,nsat,(idn(i),i=1,nsat),(freq(i),i=1,nsat),
     .                 (dellt(i),i=1,nsat),ienrf,djori,jjan 

*	fj :   step in day.dddd

	fj=float(kpj)+(float(kph)/24.d0)+(float(kpm)/1440.d0)+
     .            psec/86400.d0

       id=2
       read(10,rec=id) isat,idx,ldat1,ldat2,t0,cmx,cfx,cmy,cfy,cmz,cfz

       lan1=int(ldat1/10000)
       lan2=int(ldat2/10000)
       ian=lan2
       ian1=ian+1

       call datejj(jour,jmois,jan,jh,jmn,s,djjd)
       call datejj(1,1,ian1,0,0,0.d0,djf)
       call datejj(1,1,ian,0,0,0.d0,dj0)

200    continue

       itour=itour+1
       if (itour.gt.ipas) go to 99

*      compute the date

       djc=djjd+(itour-1)*fj

	print 102

	do ksat=1,nsat

	  id=idint((djc-djori)/dellt(ksat))+idn(ksat)

	   call jjdate(djc,llj,llm,llan,llh,llmn,sss)

*	---------------------------------------------
	   print 101,ksat
101	   format(1x,'Satellite:',i2)
	  print 300, djc,llj,llm,llan,llh,llmn,sss
300       format(f15.6, 2(1x,i2),i5,2(1x,i2),f4.1)

*	---------------------------------------------
***     4 >> nsat

	if (ksat.lt.nsat.and.id.ge.idn(ksat+1)) then
	   write(*,*) 'end of file reached'
	   stop 999
	endif
	if (ksat.eq.nsat.and.id.ge.ienrf) then
	   write(*,*) 'end of file reached'
	   stop 990
	endif
*	---------------------------------------------
          read(10,rec=id) isat,idx,ldat1,ldat2,t0,cmx,cfx,cmy,cfy
     .   ,cmz,cfz

          t1=int(t0)+0.5d0
          anu=freq(ksat)
          tau=djc-t1                            
      	  AT=TAU*ANU       
      	  TBX(I)=CMX(1)+CMX(2)*TAU+CMX(3)*SIN(AT+CFX(1))                            
     .               +CMX(4)*TAU*SIN(AT+CFX(2))                                 
     .               +CMX(5)*TAU**2*SIN(AT+CFX(3))                              
     .               +CMX(6)*SIN(2*AT+CFX(4))                                   
      	  TBY(I)=CMY(1)+CMY(2)*TAU+CMY(3)*SIN(AT+CFY(1))                            
     .               +CMY(4)*TAU*SIN(AT+CFY(2))                                 
     .               +CMY(5)*TAU**2*SIN(AT+CFY(3))                              
     .               +CMY(6)*SIN(2*AT+CFY(4))                                   
      	  TBZ(I)=CMZ(1)+CMZ(2)*TAU+CMZ(3)*SIN(AT+CFZ(1))                            
     .               +CMZ(4)*TAU*SIN(AT+CFZ(2))                                 
     .               +CMZ(5)*TAU**2*SIN(AT+CFZ(3))                              
     .               +CMZ(6)*SIN(2*AT+CFZ(4))    
                            
*	conversion of the results in AU

	  tbx(i)=tbx(i)*1000.d0/ua
	  tby(i)=tby(i)*1000.d0/ua
	  tbz(i)=tbz(i)*1000.d0/ua

	   print 100,tbx(i),tby(i),tbz(i)
100	   format(1x,'X= ',d19.12,' Y= ',d19.12,' Z= ',d19.12)

	   print 102
102	   format(70(1h=))

        enddo
        go to 200

99      continue
        stop 223

        end


      SUBROUTINE JJDATE (DJJ,JR,MS,LAN,LH,MN,SEC)

*     ----------------------------------
*     from Julian date DJJ 
*     to Gregorian date 
*     day, month, year, h, m, s (JR,MS,LAN,LH,MN,SEC)
*     ----------------------------------


      IMPLICIT REAL*8 (D-D,R-S)

*  CONVERSION D'UNE DATE  DE LA PERIODE JULIENNE EN DATE  GREGORIENNE
*  DATES POSTERIEURS AU 1 JANVIER 1601   ( 2305813.5  JOURS JULIENS )

*  DJ EST INITIALISE SELON LA VALEUR DE DJJ
      INI=1
      CALL INITI(DJJ,DJ,IAN,LAN,INI)
      DDI=0.D0
      IF(DJJ-DJ) 8,8,9
  9   IAN=IAN+1
      DJ=DJ+DFLOAT(NJA(IAN-1))
      IF(DJJ-DJ) 7,8,9
   7  LAN=IAN-1
      DI=DJJ-(DJ-DFLOAT(NJA(LAN)))
      DDI=DI-DFLOAT(IDINT(DI))
      IDI=IDINT(DI)
      IM=0
      NJE=0
  12  IM=IM+1
      NJE=NJE+NJM(IM,LAN)
      IF(IDI-NJE) 10,11,12
  10  MS=IM
      IF(MS-1) 18,18,19
  18  JR=IDI+1
      GO TO 20
  19  JR=IDI- NJE+NJM(IM,LAN)+1
  20  LH=IDINT(DDI*24.D0+0.000001D0)
      RH=DDI*24.D0-DFLOAT(LH)
      MN=IDINT(RH*60.D0+0.00001D0)
      RM=RH*60.D0-DFLOAT(MN)
      SEC=RM*60.D0
      SEC=DABS(SEC)
      GO TO 3
   8  LAN=IAN
      MS=1
      GO TO 13
  11  MS=IM+1
  13  JR=1
      IF(DDI.NE.0.D0) GO TO 20
      LH=0
      MN=0
      SEC=0.
   3  RETURN
      END

      SUBROUTINE INITI(DJJ,DJ,IAN,LAN,I)

      IMPLICIT REAL *8 (D-D)
      DIMENSION INI(16)

*  DETERMINATION DE LA DATE DE DEPART POUR LE CALCUL DE LA DATE
*          GREGORIENNE ( I=1 )  OU DE LA DATE JULIENNE ( I=2 )

      DATA INI/5813,24075,42337,60599,78861,97123,115385,133647,151910,
     &  170172,188434,206696,224958,243220,261482,279744/
      GO TO (10,20),I
  10  IDI=IDINT(DJJ-2300000.5D0)
      DO 1 I=1,16
      IF(IDI.LT.INI(I)) GO TO 1
      DJ=2300000.5D0+DFLOAT(INI(I))
      IAN=1601+(I-1)*50
   1  CONTINUE
      RETURN
  20  CONTINUE
      DO 2 I=1,16
      JAN=1601+(I-1)*50
      IF(LAN.LT.JAN) GO TO 2
      DJ=2300000.5D0+DFLOAT(INI(I))
      IAN=JAN
   2  CONTINUE
      RETURN
      END
      FUNCTION NJM(IM,LAN)
      DIMENSION NJOURS(12)
      DATA NJOURS/31,28,31,30,31,30,31,31,30,31,30,31/
      NJM=NJOURS(IM)
      IF(IM.NE.2) RETURN
      A=FLOAT(LAN)
      IE=INT(A/4.+0.0001)
      IE=IE*4
      IF((LAN-IE).EQ.0) NJM=29
      IE=INT(A/100.+0.00001)
      IE=IE*100
      IF((LAN-IE).NE.0) RETURN
      IE=IE/100
      E=FLOAT(IE)
      JE=INT(E/4.+0.0001)
      JE=JE*4
      IF((IE-JE).EQ.0) RETURN
      NJM=28
      RETURN
      END
      FUNCTION NJA(LAN)
      NJA=337+NJM(2,LAN)
      RETURN
      END

      SUBROUTINE DATEJJ(JR,MS,LAN,LH,MN,SEC,DJJ)

*     ----------------------------------
*     from Gregorian date 
*     day, month, year, h, m, s (JR,MS,LAN,LH,MN,SEC)
*     to Julian date DJJ 
*     ----------------------------------

      DOUBLE PRECISION DJJ,SEC,DJ

*  CONVERSION D'UNE DATE GREGORIENNE EN DATE EN JOURS DE LA PERIODE JULIENNE
*  DATES POSTERIEURS AU 1 JANVIER 1601

*  DJ EST INITIALISE SELON LA VALEUR DE LAN
      INI=2
      CALL INITI(DJJ,DJ,IAN,LAN,INI)
      NAN=LAN-1
      DJJ=DJ
      IF(IAN.EQ.LAN) GO TO 2
      DO 6 I=IAN,NAN
      NJ=NJA(I)
      DJJ=DJJ+DFLOAT(NJ)
   6  CONTINUE
   2  CONTINUE
      IF(MS.EQ.1) GO TO 3
      MSM=MS-1
      DO 16 IM=1,MSM
      NJ=NJM(IM,LAN)
      DJJ=DJJ+DFLOAT(NJ)
  16  CONTINUE
  3   DJJ=DJJ+DFLOAT(JR-1)
      DJJ=DJJ+DFLOAT(LH)/24.D0
      DJJ=DJJ+DFLOAT(MN)/1440.D0
      DJJ=DJJ+SEC/86400.D0
      RETURN
      END
