*
*     ------------------
*     PROGRAMME COMETEPH   ref : BDL-GF9610
*     ------------------
*
*     Auteurs :
*     P. Rocher et G. Francou.
*     Bureau des Longitudes.
*     Equipe : Dynamique du Systeme Solaire.
*     77, avenue Denfert Rochereau. 75014 PARIS.
*     E-mail : rocher@imcce.fr     francou@imcce.fr
*     tel : (33) 01 40 51 22 72  (33) 01 40 51 22 60
*     fax : (33) 01 46 33 28 34
*
*     Destinataire :
*     J. Bernard - CT/TI/MS/AM.
*     CNES.
*     18 avenue E. Belin 31055 Toulouse Cedex.
*     E-mail : bernard@venus.cst.cnes.fr
*     tel : (33) 05 61 25 00 00
*     fax : (33) 05 61 27 35 40
*
*     -----
*     Objet
*     -----
*
*     Ce programme est un exemple d'utilisation de la subroutine COMET
*     qui calcule les coordonnees heliocentriques rectangulaires des
*     cometes (ecliptique J2000).
*     Ce travail est realise dans le cadre de la mission Rosetta.
*     Langage utilise : Fortran 77.
*
*     -----------------------
*     Sources des Ephemerides
*     -----------------------
*
*     Les calculs sont effectues a partir des Ephemerides Tchebychev
*     des Cometes calculees par P. Rocher et installes sur les fichiers
*     a acces direct : cnesxx.ad.
*
*     --------------------------
*     Specifications des calculs
*     --------------------------
*
*     Periode de validite : elle depend de la Comete.
*     Repere : Equinoxe et Ecliptiques J2000.
*     Echelle de temps : Temps Dynamique Barycentrique.
*
*     -------------------------------
*     Declarations et Initialisations
*     -------------------------------
*
      implicit double precision (a-h,o-z)
      character rep
      character*6 ncnes
*
      dimension r(6)
*
      data nul/10/
*
*     ----------------------
*     Lecture des parametres
*     ----------------------
*
100   continue
*
      call LCOMET (nul,ncnes,num,tdeb,tfin)
      call LDATE (tdeb,tfin,t0,dt)
      t=t0
*
*     -----------------------------------
*     Calcul et Affichage des Ephemerides
*     -----------------------------------
*
200   continue
*
      call COMET (t,nul,r,ierr)
      if (ierr.ne.0) goto 300
*
      call DJDATE (t,3,ij,im,ia,nh,nm,ns)
*
      call CLRSCR
      write (*,1001) ncnes,ij,im,ia,nh,nm,ns,t,r
*
*     --------------------------------
*     Autre Calcul ou Fin du Programme
*     --------------------------------
*
      write (*,1002)
      read (*,1003) rep
*
      if (rep.eq.'n'.or.rep.eq.'N') then
         close (nul)
         goto 100
      endif
*
      if (rep.eq.'q'.or.rep.eq.'Q') then
         call CLRSCR
         write (*,1004)
         read (*,*)
         close (nul)
         stop
      endif
*
      if (rep.eq.'+'.and.t+dt.le.tfin) t=t+dt
      if (rep.eq.'-'.and.t-dt.ge.tdeb) t=t-dt
      goto 200
*
*     --------------
*     Message erreur
*     --------------
*
300   continue
      call CLRSCR
      write (*,1005) ierr
      read (*,*)
      stop
*
*     -------
*     Formats
*     -------
*
1001  format (2x,77('-')/
     .        2x,'EPHEMERIDES COMETES-CNES',
     .        t50,'Bureau des Longitudes - GF9610'/
     .        2x,77('-')//
     .        2x,'Comete',2x,a6///
     .        2x,'Coordonnees Heliocentriques'//
     .        2x,'Repere : Equinoxe et Ecliptique moyens J2000'//
     .        2x,'Date :  ',i2.2,'/',i2.2,'/',i4,
     .        3x,'Heure : ',i2.2,'h ',i2.2,'m ',i2.2,'s TDB'//
     .        2x,'Date Julienne : ',f13.5///
     .        2x,'X ',f15.11,' ua  ',
     .        6x,'Y ',f15.11,' ua  ',
     .        6x,'Z ',f15.11,' ua  '//
     .        2x,'X''',f15.11,' ua/j',
     .        6x,'Y''',f15.11,' ua/j',
     .        6x,'Z''',f15.11,' ua/j'/)
1002  format (2x,77('-')/
     .        2x,'(+) Date suiv',3x,'(-) Date prec',
     .        3x,'(N) Nouveau',3x,'(Q) Quitter',
     .        5x,'Votre choix ? ',$)
1003  format (a1)
1004  format (/2x,'Programme termine : Appuyez sur Entree',$)
1005  format (/2x,'Erreur ',i2.2,' : Appuyez sur Entree',$)
*
      end
*
*
*
      subroutine COMET (tdj,nul,r,ierr)
*-----------------------------------------------------------------------
*
*     ref : BDL-GF9610
*
*     Objet :
*     Calcul des coordonnees rectangulaires heliocentriques d'une Comete.
*     Repere : Equinoxe et Ecliptique J2000.
*
*     En entree :
*     tdj         Date Julienne TDB (reel dp).
*     nul         Numero d'Unite Logique du Fichier Comete (entier).
*
*     En sortie :
*     r(3)        Table des Coordonnees Rectangulaires (reel dp).
*                 r(1): position x (ua)  r(4) : vitesse x (ua/j).
*                 r(2): position y (ua)  r(5) : vitesse y (ua/j).
*                 r(3): position z (ua)  r(6) : vitesse z (ua/j).
*     ierr        Indice Erreur (entier).
*                 ierr=0 pas d'erreur.
*                 ierr=10 erreur de fichier.
*                 ierr=11 erreur de date.
*
*     Remarque :
*     Le Fichier Comete cnesxx.ad (acces direct) doit etre ouvert avant
*     le premier appel a cette subroutine, par exemple :
*     open (10,file='cnesxx.ad',access='direct',recl=80,status='old').
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      character*6 ncnes,niau
      character*30 nom
      integer delta
*
      parameter (ncf=10)
*
      dimension r(6)
      dimension coef(ncf,3),t(ncf),u(ncf)
*
      data nul0/0/
*
      do n=1,6
         r(n)=0.d0
      enddo
*
      if (nul.ne.nul0) then
         read (nul,rec=1,iostat=nerr) nbloc,tdeb,delta,ncnes,niau,nom
         if (nerr.ne.0) goto 100
         tfin=tdeb+nbloc*delta
         nul0=1
      endif
*
      ierr=11
      if (tdj.lt.tdeb.or.tdj.gt.tfin) return
*
      n=(tdj-tdeb)/delta
      if (tdj.eq.tfin) n=n-1
      t0=tdeb+n*delta
      nrec=3*n+1
*
      do n=1,3
         nrec=nrec+1
         read (nul,rec=nrec,iostat=nerr) (coef(i,n),i=1,ncf)
         if (nerr.ne.0) goto 100
      enddo
*
      x=2.d0*(tdj-t0)/delta-1.d0
*
      t(1)=1.d0
      t(2)=x
      u(1)=1.d0
      u(2)=2.d0*x
*
      do i=3,ncf
         t(i)=2.d0*x*t(i-1)-t(i-2)
         u(i)=2.d0*x*u(i-1)-u(i-2)
      enddo
*
      do n=1,3
         do i=1,ncf
            r(n)=r(n)+t(i)*coef(i,n)
            if (i.gt.1) r(n+3)=r(n+3)+(i-1)*u(i-1)*coef(i,n)
         enddo
         r(n+3)=r(n+3)*2.d0/delta
      enddo
*
      ierr=0
      return
*
100   continue
      ierr=10
      return
*
      end
*
*
*
      subroutine LCOMET (nul,ncnes,num,tdeb,tfin)
*-----------------------------------------------------------------------
*
*     ref : BDL-GF9610
*
*     Objet :
*     Lecture des parametres lies a une Comete.
*
*     En entree :
*     nul         Numero d'Unite Logique du Fichier Comete (entier).
*
*     En sortie :
*     ncnes       Nom cnes de la Comete (character*6).
*     num         Numero de la Comete (entier).
*     tdeb        Date Julienne Debut de periode de validite (reel dp).
*     tfin        Date Julienne Fin de periode de validite (reel dp).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      logical fexist,ok
      integer delta
      character*2 anum
      character*6 ncnes,niau
      character*30 nom
      character*40 fich
*
100   continue
      call CLRSCR
      write (*,1000)
*
      write (*,1001)
      read (*,1002,iostat=nerr) num
      if (nerr.ne.0.or.num.lt.1.or.num.gt.999) then
         write (*,2001)
         read (*,*)
         goto 100
      endif
*
      ok=.true.
      write (anum,'(i2.2)') num
      fich='cnes'//anum//'.ad'
      inquire (file=fich,exist=fexist)
*
      if (fexist) then
         open (nul,file=fich,access='direct',recl=80,status='old',
     .      iostat=nerr)
         if (nerr.eq.0) then
            read (nul,rec=1,iostat=nerr) nbloc,tdeb,delta,ncnes,niau,nom
            if (nerr.eq.0) then
               tfin=tdeb+nbloc*delta
            else
               ok=.false.
            endif
         else
            ok=.false.
         endif
      else
         ok=.false.
      endif
*
      if (.not.ok) then
         write (*,2002)
         read (*,*)
         goto 100
      else
         return
      endif
*
1000  format (2x,77('-')/
     .        2x,'EPHEMERIDES COMETES-CNES',
     .        t50,'Bureau des Longitudes - GF9610'/
     .        2x,77('-'))
1001  format (//2x,'Choix de la Comete'/
     .        /2x,'Entrez le Numero Cnes de la Comete (xxx) : ',$)
1002  format (i4)
*
2001  format (/2x,'Erreur de reponse : Appuyez sur Entree',$)
2002  format (/2x,'Erreur de Fichier : Appuyez sur Entree',$)
*
      end
*
*
*
      subroutine LDATE (tdeb,tfin,t0,dt)
*-----------------------------------------------------------------------
*
*     ref : BDL-GF9610
*
*     Objet :
*     Lecture des parametres de datation.
*
*     En entree :
*     tdeb        Date Julienne Debut de periode de validite (reel dp).
*     tfin        Date Julienne Din de periode de validite (reel dp).
*
*     En sortie :
*     t0          Date Julienne initiale TDB (reel dp).
*     dt          Pas en Jours entre deux dates (reel dp).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      character rep
*
100   continue
      call CLRSCR
      write (*,1000)
*
      call DJDATE (tdeb,0,j1,m1,ia1,nh,nm,ns)
      call DJDATE (tfin,0,j2,m2,ia2,nh,nm,ns)
      idat1=ia1*10000+m1*100+j1
      idat2=ia2*10000+m2*100+j2
*
      write (*,1001) idat1,idat2
      read (*,2001,iostat=nerr) idat
      if (nerr.ne.0) then
         write (*,2000)
         read (*,*)
         goto 100
      endif
*
      ia=idat/10000
      m=mod(idat,10000)/100
      j=mod(idat,100)
      call DATEDJ (j,m,ia,0,0,0,t0)
      if (t0.eq.0.d0.or.t0.lt.tdeb.or.t0.gt.tfin) then
         write (*,2000)
         read (*,*)
         goto 100
      endif
*
      write (*,1002)
      read (*,2002,iostat=nerr) iheur
      if (nerr.ne.0) then
         write (*,2000)
         read (*,*)
         goto 100
      endif
*
      ih=iheur/10000
      im=mod(iheur,10000)/100
      is=mod(iheur,100)
      call DATEDJ (j,m,ia,ih,im,is,t0)
      if (t0.eq.0.d0.or.t0.lt.tdeb.or.t0.gt.tfin) then
         write (*,2000)
         read (*,*)
         goto 100
      endif
*
      write (*,1003)
      read (*,2003,iostat=nerr) dt
      if (nerr.ne.0.or.dt.le.0.d0.or.dt.ge.1000.d0) then
         write (*,2000)
         read (*,*)
         goto 100
      endif
*
      write (*,1004)
      read (*,2004) rep
      if (rep.ne.' '.and.rep.ne.'O'.and.rep.ne.'o') goto 100
      return
*
1000  format (2x,77('-')/
     .        2x,'EPHEMERIDES COMETES-CNES',
     .        t50,'Bureau des Longitudes - GF9610'/
     .        2x,77('-'))
1001  format (//2x,'Entrez la date initiale (aaaammjj)',
     .        2x,'entre ',i8,' et ',i8,t66,' ?  ',$)
1002  format (//2x,'Entrez l''heure initiale (hhmmss)  ',
     .        2x,'entre 000000   et 240000',t66,' ?  ',$)
1003  format (//2x,'Entrez le pas du calcul  (jours)  ',
     .        2x,'entre 1.00     et 1000.00',t66,' ?  ',$)
1004  format (//2x,'Confirmez-vous ces choix (o/n) ? ',$)
2000  format (/2x,'Erreur de reponse : Appuyez sur Entree',$)
2001  format (i9)
2002  format (i7)
2003  format (f10.5)
2004  format (a1)
*
      end
*
*
*
      subroutine DJDATE (tdj,icode,jour,mois,kan,nh,nm,ns)
*-----------------------------------------------------------------------
*
*     ref : BDL-GF9504
*
*     Objet :
*     Passage d'une date julienne dans le calendrier julien/gregorien.
*
*     En entree :
*     tdj         date julienne (reel dp).
*     icode       indce code de la conversion (entier)
*                 icode=0 conversion tdj -> jour-mois-kan et 00-00-00.
*                 icode=1 conversion tdj -> jour-mois-kan et nh-00-00.
*                 icode=2 conversion tdj -> jour-mois-kan et nh-nm-00.
*                 icode=3 conversion tdj -> jour-mois-kan et nh-nm-ns.
*
*     En sortie :
*     jour        quantieme du mois (entier).
*     mois        numero du mois (entier).
*     kan         millesime astronomique de l'annee (entier).
*     nh          nombre d'heures (entier).
*     nm          nombre de minutes (entier).
*     ns          nombre de secondes (entier).
*
*     Remarque :
*     On utilise le formulaire de J. Meeus.
*     La date julienne tdj doit etre superieure ou egale a 0.
*     Si tdj<0, alors jour=mois=kan=nh=nm=ns=0.
*     L'annee kan est comptee en millesimes astronomiques.
*     L'annee 1 avant J.C. des historiens correspond a l'an 0 des
*     astronomes et ainsi de suite.
*     Ex : l'an 80 avant J.C. correspond a l'an astronomique -79.
*
*     Rappel :
*     Le calendier gregorien commence le 15/10/1582 a 0h (DJ2299160.5).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      jour=0
      mois=0
      kan=0
      nh=0
      nm=0
      ns=0
*
      if (tdj.lt.0.d0.or.icode.lt.0.or.icode.gt.3) return
*
      if (icode.eq.0) dt=0.0d0
      if (icode.eq.1) dt=0.5d0/24.d0
      if (icode.eq.2) dt=0.5d0/1440.d0
      if (icode.eq.3) dt=0.5d0/86400.d0
*
      t=tdj+dt+0.5d0
      z=dint(t)
      f=t-z
*
      if (z.lt.2299161.d0) then
         a=z
      else
         alpha=dint((z-1867216.25d0)/36524.25d0)
         a=z+1.d0+alpha-dint(alpha/4.d0)
      endif
*
      b=a+1524.d0
      c=dint((b-122.1d0)/365.25d0)
      d=dint(365.25d0*c)
      e=dint((b-d)/30.6001d0)
*
      jour=b-d-dint(30.6001d0*e)
*
      if (e.lt.13.5d0) then
         mois=e-1.d0
      else
         mois=e-13.d0
      endif
*
      if (mois.lt.2.5d0) then
         kan=c-4715.d0
      else
         kan=c-4716.d0
      endif
*
      f=f*24.d0
      nh=f
      f=(f-nh)*60.d0
      nm=f
      f=(f-nm)*60.d0
      ns=f
*
      return
      end
*
*
*
      subroutine DATEDJ (jour,mois,kan,nh,nm,ns,tdj)
*-----------------------------------------------------------------------
*
*     ref : BDL-GF9504
*
*     Objet :
*     Passage du calendrier gregorien/julien a une date julienne.
*
*     En entree :
*     jour        quantieme du mois (entier).
*     mois        numero du mois (entier).
*     kan         millesime astronomique de l'annee (entier).
*     nh          nombre d'heures (entier).
*     nm          nombre de minutes (entier).
*     ns          nombre de secondes (entier).
*
*     En sortie :
*     tdj         date julienne (reel dp).
*
*     Remarque :
*     On utilise le formulaire de J. Meeus.
*     Cette fonction n'est valable que pour les dates posterieures au
*     1/1/4713 12h avant J.C.
*     L'annee kan est comptee en millesimes astronomiques.
*     L'annee 1 avant J.C. des historiens correspond a l'an 0 des
*     astronomes et ainsi de suite.
*     Ex : l'an 80 avant J.C. correspond a l'an astronomique -79.
*
*     Rappel :
*     Le calendier gregorien commence le 15/10/1582 a 0h (DJ2299160.5).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      tdj=0.d0
*
      if ((kan.lt.-4713)
     .   .or.(mois.lt.1.or.mois.gt.12)
     .   .or.(jour.lt.1.or.jour.gt.31)
     .   .or.(nh.lt.0.or.nh.gt.24)
     .   .or.(nm.lt.0.or.nm.gt.59)
     .   .or.(ns.lt.0.or.ns.gt.59)) return
*
      a=0.d0
      b=0.d0
      c=0.d0
*
      if (mois.gt.2) then
         y=kan
         m=mois
      else
         y=kan-1.d0
         m=mois+12
      endif
*
      if (y.lt.0) then
         c=-0.75d0
      else
         idate=kan*10000+mois*100+jour
         if (idate.ge.15821015) then
            a=dint(y/100.d0)
            b=2.d0-a+dint(a/4.d0)
         endif
      endif
*
      tdj=dint(365.25d0*y+c)
     .   +dint(30.6001d0*(m+1))
     .   +jour+b
     .   +dfloat(nh)/24.d0
     .   +dfloat(nm)/1440.d0
     .   +dfloat(ns)/86400.d0
     .   +1720994.5d0
*
      return
      end
*
*
*
      subroutine CLRSCR
*-----------------------------------------------------------------------
*
*     ref : BDL-GF9504
*
*     Objet :
*     Effacement de l'ecran.
*
*     Remarque :
*     L'effacement d'ecran fonctionne avec les systemes d'exploitation
*     DOS et UNIX seulement.
*     Pour rendre operationnelle la subroutine CLRSCR, il faut retirer
*     les commentaires :
*     *DOS  dans le cas du systeme d'exploitation DOS,
*     *UNX  dans le cas du systeme d'exploitation UNIX.
*
*-----------------------------------------------------------------------
*
*DOS  character*7 escscr
*DOS  escscr=char(27)//char(91)//'2J'//char(27)//char(91)//'H'
*DOS  write (*,'(2x,a)') escscr
*
*UNX  call system ('clear')
*
      character*7 escscr
      escscr=char(27)//char(91)//'2J'//char(27)//char(91)//'H'
      write (*,'(2x,a)') escscr
*
      return
      end
