c     programme bda/progf/ephem.f
c
c     Calcul des ephemerides pour les binaires spectroscopiques
c     a: On donne une date civile, on recoit la phase et la Vr
c     b: on donne un intervalle de dates civile, on recoit un
c        fragment de la courbe de Vr
c     c: on entre une phase ou un intervalle de phase et l'on
c        recoit les dates civiles correspondantes
c
      double precision pj,dobser,pinf,psup
      character noet*16,isb*1,rep*1,fres*40
      character*3 lmois(12)
c     character mois*5
      logical l
      dimension nflag(2),nj(12)
      data lmois/'jan','fev','mar','avr','mai','jui','jul',
     1            'aou','sep','oct','nov','dec'/
      data nj/31,28,31,30,31,30,31,31,30,31,30,31/
      data po2/0./,fk2/0./,impr/0/,iopen/0/,iferm/1/,mult/0/
      data fres/'/home/vaud/mermio/bda/tmp/ephem.out'/
c
      pi2=6.28318531
      degra=pi2/360.
      grade=1./degra
c
c     Entete
c

      write(6,110)
  110 format(/,' Calcul des ephemerides des etoiles SB')
c
c     Identification de l'etoile etudiee
c
    1 write(6,105)
  105 format(/,' Quelle est l''etoile etudiee?  ',$)
      read(5,106) noet
  106 format(a16)
c
      write(6,107)
  107 format(/,' Quelle est sa periode?  ',$)
      read(5,*) perio
c
      write(6,108)
  108 format(/,' Quelle est l''epoque T?  ',$)
      read(5,*) epoq
c
      write(6,116)
  116 format(/,' Quelle est la vitesse systemique?  ',$)
      read(5,*) vsys
c
      write(6,112)
  112 format(/,' Quelle est l''amplitude?  ',$)
      read(5,*) gk1
c
      write(6,109)
  109 format(/,' Quelle est l''excentricite?  ',$)
      read(5,*) exc
      if(exc.eq.0.) then
        om1=0.
        goto 2
      else
        write(6,111)
  111   format(/,' Quelle est l''orientation omega?  ',$)
        read(5,*) om1
      endif
c
    2 write(6,113) noet
  113 format(/,1x,a16' est-elle SB2 (o/n)?  ',$)
      read(5,114) isb
  114 format(a1)
      if(isb.eq.'o') then
        om2=om1+180.
        if(om2.gt.360.) om2=om2-360.
        om1=om1*degra
        om2=om2*degra
        write(6,115)
  115   format(/,' Donnez le rapport des masses (Mb/Ma):  ',$)
        read(5,*) rap
        gk2=gk1/rap
        nflag(1)=1
        nflag(2)=1
      else
        om1=om1*degra
        nflag(1)=1
        nflag(2)=0
      endif
c
c     Menu du programme
c
    5 write(6,120)
  120 format(//,' Possibilites offertes par le programme:',/
     1/,'   Calcul de la phase pour une date civile donnee ..... 1',
     2/,'   Calcul des phases pour un intervalle de dates ...... 2',
     3/,'   Calcul de la date pour une phase donnee ............ 3',
     4/,'   Calcul des dates correspondant a des phases ........ 4',
     5/,'   Plusieurs phases au cours d''une seule nuit ......... 5',/
     5/,'   Ouverture d''un fichier de resultats ................ 6',
     5/,'   Fermeture du fichier des resultats ................. 7',
     5/,'   Creation d''un fichier par etoile ................... 8',
     6/,'   Determination pour une autre etoile ................ 9',/
     7/,'   Fin de la procedure ............................... 10',
     8//,'   Que desirez-vous calculer?  ',$)
      read(5,*) ich
      goto(10,20,30,40,50,60,70,80,1,90),ich
c
c     Calcul d'une phase pour une date civile donnee
c
   10 write(6,130)
  130 format(//,'   Donnez le jour, le mois et l''annee:  ',$)
      read(5,*,err=10) jour,mois,iann
      jour=jour+1
      if(jour.gt.nj(mois)) then
        jour=1
        mois=mois+1
          if(mois.eq.13) then
            mois=1
            iann=iann+1
          endif
      endif
      if(iann.lt.99) iann=iann+1900
      goto 13
c
c     Calcul des phases pour un intervalle donne
c
   20 write(6,140)
  140 format(//,'   Donnez le jour, le mois et l''annee du debut:  '
     1,$)
      read(5,*,err=20) jour1,mois1,iann1
      jour1=jour1+1
      if(jour1.gt.nj(mois1)) then
        jour1=1
        mois1=mois1+1
          if(mois1.eq.13) then
            mois1=1
            iann1=iann1+1
          endif
      endif
      if(iann1.lt.99) iann1=iann1+1900
   21 write(6,150)
  150 format(//,'   Donnez le jour, le mois et l''annee de la fin:  '
     1,$)
      read(5,*,err=21) jour2,mois2,iann2
      jour2=jour2+1
      if(jour2.gt.nj(mois2)) then
        jour2=1
        mois2=mois2+1
          if(mois2.eq.13) then
            mois2=1
            iann2=iann2+1
          endif
      endif
      if(iann2.lt.99) iann2=iann2+1900
      jour=jour1
      mois=mois1
      iann=iann1
      goto 13
c
   22 if(iann.eq.iann2) then
        if(mois.eq.mois2) then
          if(jour.gt.jour2) then
            write(6,141)
  141       format(/,' Un autre intervalle (o/n)?  ',$)
            read(5,*) rep
            if(rep.eq.'o') then
              impr=0
              goto 20
            else
              impr=0
              goto 5
            endif
          endif
        endif
      endif
      goto 13
c
c
c     do 11 i=1,12
c     if(mois(1:3).eq.lmois(i)) then
c        month=i
c        goto 12
c     else
c        read(mois,400) month
c 400    format(i3)
c     endif
c  11 continue
c  12 continue
c
c
c     Calcul de la date pour une phase donnee
c
   30 write(6,160)
  160 format(/,' Quelle phase:  ',$)
      read(5,*) fas
   31 write(6,161)
  161 format(/,' Donnez la date limite inferieure:  ',$)
      read(5,*,err=31) jour,mois,iann
      jour=jour+1
      if(jour.gt.nj(mois)) then
        jour=1
        mois=mois+1
          if(mois.eq.13) then
            mois=1
            iann=iann+1
          endif
      endif
      if(iann.lt.99) iann=iann+1900
      tusin=0.0
      call rjul(tusin,jour,mois,iann,pinf)
c
   32 write(6,162)
  162 format(/,' Donnez la date limite superieure:  ',$)
      read(5,*,err=32) jour,mois,iann
      jour=jour+1
      if(jour.gt.nj(mois)) then
        jour=1
        mois=mois+1
          if(mois.eq.13) then
            mois=1
            iann=iann+1
          endif
      endif
      if(iann.lt.99) iann=iann+1900
      tusin=0.0
      call rjul(tusin,jour,mois,iann,psup)
c
      obser=epoq+perio*fas
   33 if(obser.ge.pinf.and.obser.le.psup) then
        osr=obser+2440000.0
        dobser=dble(osr)
        call juldat(dobser,tus,jr,ms,ia)
        phas=fas
        goto  15
      else
        if(obser.gt.psup) then
          write(6,163)
  163     format(/,' Une autre phase (o/n)?  ',$)
          read(5,*) rep
          if(rep.eq.'o') then
            impr=0
            goto 30
          else
            impr=0
            goto 5
          endif
        else
          obser=obser+perio
          goto 33
        endif
      endif
c
   40 continue
   50 continue
c
c     Ouverture d'un fichier resultat en ecriture
c
   60 if(iferm.eq.0) then
        write(6,215)
  215   format(/,'  Le fichier est deja ouvert!' )
        goto 5
      endif
c
      inquire(file=fres,exist=l)
      if(.not.l) then
        open(15,file=fres,status='new')
      else
        write(6,220)
  220   format(/,'  Ce fichier existe deja.',//,
     1'                 Voulez-vous ecrire a la suite (o/n)?  ',$)
        read(5,*) rep
        if(rep.eq.'o') then
          open(15,file=fres,status='old',fileopt='eof')
          iopen=1
        else
          open(15,file=fres,status='old')
        endif
      endif         
      iopen=1
      iferm=0
      goto 5
c
c     Fermeture du fichier
c
   70 if(iferm.eq.1) then
        write(6,216)
  216   format(/,'  Le fichier est deja ferme!' )
      else
        close(15,status='keep')
        iferm=1
      endif
      goto 5
c
c     Determination de la periode julienne
c
   13 frac=0.0
      tusin=frac*86400.
      call rjul(tusin,jour,mois,iann,pj)
c
c     Calcul de la phase
c
   14 if(pj.gt.epoq) then
        ph=(pj-epoq)/perio
        ipha=int(ph)
        phas=ph-float(ipha)
      else
        ph=(epoq-pj)/perio
        ipha=int(ph)
        pha=ph-float(ipha)
        phas=1.-pha
      endif
c
c     Creation d'un fichier par etoile
c
   80 mult=1
c     goto 1
c
c     Calcul de la vitesse radiale
c
c
   15 call calvr(vsys,phas,exc,gk1,gk2,om1,om2,nflag,vr1,vr2)
c
      if(ich.eq.2.and.impr.eq.1) goto 62
      if(ich.eq.3.and.impr.eq.1) goto 62
      impr=1
      write(6,197) noet
  197 format(//,' Ephemerides pour l''etoile ',a16,/)
      if(iopen.eq.1) write(15,197) noet
      if(isb.eq.'o') then
        write(6,198)
        if(iopen.eq.1) write(15,198) 
  198   format('     Date           P J       Phase   ',
     1  '   A    Vr    B',/)
      else
        write(6,199)
        if(iopen.eq.1) write(15,199)
  199   format('     Date           P J       Phase   ',
     1  '      Vr ',/)
      endif
c
   62 if(ich.eq.3) then
        jou=jr-1
        moi=ms
        ian=ia
        pj=obser
        phas=fas
        goto 63
      else
        jou=jour-1
        moi=mois
        ian=iann
      endif
   63 if(jou.eq.0) then
        moi=mois-1
        if(moi.eq.0) then
          moi=12
          ian=iann-1
        endif
        jou=nj(moi)
      endif
      if(isb.eq.'o') then
        write(6,200) jou,lmois(moi),ian,pj,phas,vr1,vr2
  200   format(2x,i2,1x,a3,1x,i4,3x,'244 ',f7.2,3x,f5.2,5x,f7.2,
     1  3x,f7.2,' km/s')
        if(iopen.eq.1) write(15,200) jou,lmois(moi),ian,pj,
     1    phas,vr1,vr2
      else
        write(6,210) jou,lmois(moi),ian,pj,phas,vr1
  210   format(2x,i2,1x,a3,1x,i4,3x,'244 ',f7.2,3x,f5.2,5x,f7.2,
     1  ' km/s')
        if(iopen.eq.1) write(15,210) jou,lmois(moi),ian,pj,
     1    phas,vr1
      endif
c
      if(ich.eq.1) then
        write(6,117)
  117   format(/,' Une autre date (o/n)?  ',$)
        read(5,*) rep
        if(rep.eq.'o') then
          goto 10
        else
          impr=0
          goto 5
        endif
      endif
c
      if(ich.eq.2) then
        jour=jour+1
        lim=nj(mois)+1
        if(jour.gt.lim) then
          jour=1
          mois=mois+1
            if(mois.eq.13) then
              mois=1
              iann=iann+1
            endif
        endif
        goto 22
      endif
c
      if(ich.eq.3) then
        obser=obser+perio
        goto 33
      endif
c
      if(ich.eq.4) goto 5
c
   90 stop
      end
c
c
c
      subroutine rjul(tusin,j,m,ia,pj)
c
c     calcul de la periode julienne pour une date et 
c     une heure TU donnee
c
c     tus = temps universel moyen, en seconde
c     J, M, IA = jour, mois, annee civile
c     rjul = periode julienne geocentrique
c     donne la periode julienne modulo 2440000
c     precision 1/1000 de jour
c     valable de 1970 a 1995
c
      dimension lmois(11)
      double precision pj,pjul,sjul,tus
      data lmois/31,59,90,120,151,181,212,243,273,304,334/
c
      mod(i,k)=i-(i/k)*k
c
      tus=dble(tusin)
      nbj=0
      m1=m-1
      if(m1) 1,1,10
   10 nbj=lmois(m1)
      if(m-2) 1,1,2
    2 if(mod(ia,4)) 1,4,1
    4 nbj=nbj+1
    1 nbj=nbj+j
      a1=dble(ia-1)
      pjul=dble(nbj)+(365.25D0*(a1-1287.D0))-0.25D0
      nbj=int(pjul-248513.D0)
      pjul=dble(nbj)
      sjul=(tus-43200.D0)*1.1574D-5
      pj=pjul+sjul
c
      return
      end
c
c
c
      subroutine calvr(vsyst,phase,pe,fk1,fk2,po1,po2,nflag,vrc1,vrc2)
c
c     Calcul de la courbe des vitesses
c
      real m
      dimension nflag(2)
      data pi2/6.28318531/
c
c  60 do 7 i=1,201 
c     phase=float(i-1)/200.
c
c     definition de l'anomalie moyenne:
c
      m=pi2*phase
c
c     resolution de l'equation de Kepler
c
      e=m 
      do 6 j=1,5 
    6 e=e+(m-e+pe*sin(e))/(1.-pe*cos(e))
      w=2.*atan(sqrt((1.+pe)/(1.-pe))*tan(e/2.))  
      if(nflag(1).eq.0) then
        vrc1=0.
        else
        vrc1=vsyst+fk1*(pe*cos(po1)+cos(w+po1)) 
      endif
      if(nflag(2).eq.0) then
        vrc2=0.
        else
        vrc2=vsyst+fk2*(pe*cos(po2)+cos(w+po2))     
      endif
c   7 continue
c
      return
      end
c
c
c
      subroutine juldat(dj,tus,jr,ms,ia)
c
c     sous-programme d'interet general
c     calcul de l'heure et de la date a partir de la periode julienne
c
c     dj          periode julienne donnee (en double precision)
c                 (si dj est heliocentrique, tranformer dj
c                 prealablement en utilisant hjulp)
c     tus         temps universel en secondes de la pj donnee
c     jr,ms,ia    jours, mois et annee de la date a la periode pj
c
c
c
      double precision dj,a,julien
      double precision r
      logical bisxt
      dimension lmois(0:11)
      data lmois/0,31,59,90,120,151,181,212,243,273,304,334/
c
      a=dj-2415020.0
c
      a=a/365.2425
      ia=a
c
c     recherche du jour et du mois
c
      nj=(a-ia)*365.2425 + 1
      ia=ia + 1900
      do 1 i=11,0,-1
      j=i
      lm=lmois(i)
      if(bisxt (ia)  .and.  i .ge. 2) lm=lm+1
      if(nj .gt. lm) goto 2
    1 continue
    2 ms=j+1
      jr=nj-lm
c
c     calcul de l'heure
c
      r=julien (0.,jr,ms,ia)
      tus=(dj-r)*86400.
      if(tus .lt. 0.) then
        call minus(jr,ms,ia)
        tus=tus+86400
        endif
      return
c
  100 format(' subroutine juldat: pj inferieur a 1.1.1900 a 0h tu')
      end
c
c
c
      subroutine minus (jr,ms,ia)
c
c     soustraction de 1 jour a la date jr.ms.ia
c
      logical bisxt
      dimension mois(12)
      data mois/31,28,31,30,31,30,31,31,30,31,30,31/
c
      jr=jr-1
      if(jr .eq. 0) then
        ms=ms-1
        if(ms .eq. 0) then
          ia=ia-1
          ms=12
          endif
        jr=mois(ms)
        if(bisxt (ia)  .and.  ms .eq. 2) jr=jr+1
      endif
      return
      end
c
c
c
      logical function bisxt (ia)
c
c     determination d'une annee bissextile
c
c     bisxt    .true.  ia est une annee bissextile
c     bisxt    .false.  ia n'est pas une annee bissextile
c     ia        annee donnee
c
      logical mod100
c
      mod100=mod(ia,100) .eq. 0
      bisxt=(.not. mod100  .and.  mod(ia,4) .eq. 0)   .or.
     1      (mod100  .and.  mod(ia/100,4) .eq. 0)
      return
      end
c
c
c
      double precision function julien (tus,j,m,ia)
c
c     calcul de la periode julienne pour une date et un intervalle tu
c     quelconque a partir du 0 janvier 1900 0h tu (31 decembre 1899)
c
c     tus          temps universel, en secondes
c     j, m, ia     jour, mois et annee civile
c     julien       periode julienne (en double precision)
c
      logical bisxt
      dimension lmois(12)
      data lmois/0,31,59,90,120,151,181,212,243,273,304,334/
c
      nbj=j-1 + lmois(m)
      if(m .gt. 2 .and. bisxt (ia)) nbj=nbj + 1
      ia1=ia-1
      jul= 365.25*(4713.+ia1) + 0.8
      jul=jul + nbj - 13
      julien=jul + (dble(tus)/86400.0d0) - 0.5
      return
      end