*
*     --------------
*     PROGRAMME SOHO   ref : BDL-GF9504
*     --------------
*
      subroutine sohoephem(tjd,hec,ic,slong,slat,l0,b0,orient,nrot,
     +     trot,dsol,plong,plat,xphase,fillum,dpla,ienter)

*
*     -------------------------------
*     Declarations et initialisations
*     -------------------------------
*
      implicit double precision (a-h,o-z)
      character*40 fich
      double precision l0
      double precision slong(3),slat(3),orient(2),dsol(2)
      double precision plong(5,3),plat(5,3),xphase(5),
     +     fillum(5),dpla(5,2),hec(6)
      
      common/IORIG/iorig
      common/RSOHO/rsoho(6)

      common/IFLAG/iflag

*
*     --------------------------------------------
*     Ouverture du fichier des ephemerides sources
*     --------------------------------------------
*
      do i=1,6
         rsoho(i)=hec(i)
      enddo

      iflag=ienter

      iorig=ic
      nul=10
      lrec=3640
      fich='bdl82.ad'
      open (nul,file=fich,access='direct',recl=3640,status='old')
*
*      call DATEDJ(ij,im,ia,nh,nm,ns,t)
*      t0=t-2400000.5d0

*
*     -----------------------------------
*     Calcul et affichage des ephemerides
*     -----------------------------------
*
      call EPHSOL1(tjd,nul,ierr,slong,slat,l0,b0,orient,nrot,
     +     trot,dsol)
      if (ierr.ne.0) return

      call EPHPLA1(tjd,nul,ierr,plong,plat,xphase,fillum,dpla)
      if (ierr.ne.0) return
*
*
      
      close (nul)

      ienter=iflag

      end
*
*
*
      subroutine EPHSOL1 (tjd,nul,ierr,slong,slat,l0,b0,orient,nrot,
     +     trot,dsol)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul et affichage des ephemerides du Soleil (projet SOHO).
*
*     En entree :
*     t0          date julienne modifiee initiale TDB (reel dp).
*     dt          pas du calcul en jours (reel dp).
*     nul         numero d'unite logique du fichier source (entier).
*
*     En sortie :
*     ierr        indice erreur (entier).
*                 ierr=0 pas d'erreur.
*                 ierr=10 erreur de fichier.
*                 ierr=11 erreur de date.
*                 ierr=12 erreur de corps.
*                 ierr=13 erreur de centre.
*                 ierr=14 erreur de vitesse.
*                 ierr=15 erreur de plan.
*
*     Appels subroutines et fonctions :
*     POSAST, POSMOY, POSAPP, SOLPHY, SOLSYN, DIASOL,
*     DJDATE, RDCONV, CLRSCR.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      double precision l0
      double precision slong(3),slat(3),orient(2),dsol(2)

      common/IORIG/IORIG
*
      sl=0.d0
      sb=0.d0
      sp1=0.d0
      sp2=0.d0
      ds=0.d0
*
*      t=t0+2400000.5d0
*
100   continue
*
      call POSAST (tjd,11,2,nul,sl1,sb1,x,ierr)
      if (ierr.ne.0) return
*
      call POSMOY (tjd,11,2,nul,sl2,sb2,x,ierr)
      if (ierr.ne.0) return
*
      call POSAPP (tjd,11,2,nul,sl3,sb3,x,ierr)
      if (ierr.ne.0) return
*
      call SOLPHY (tjd,nul,sl0,sb0,sp1,sp2,ierr)
      if (ierr.ne.0) return
*
      call SOLSYN (tjd,nul,trot,nrot,ierr)
      if (ierr.ne.0) return
*
      ds=DIASOL(tjd,nul,ierr)
      if (ierr.ne.0) return
*
*      call DJDATE (t,3,ij,im,ia,nh,nm,ns)
*      call DJDATE (tdrs,2,ijr,imr,iar,nhr,nmr,nsr)
      call RDCONV (sl1,1114,isisl1,idsl1,imsl1,issl1,ifsl1)
      call RDCONV (sb1,1114,isisb1,idsb1,imsb1,issb1,ifsb1)
      call RDCONV (sl2,1114,isisl2,idsl2,imsl2,issl2,ifsl2)
      call RDCONV (sb2,1114,isisb2,idsb2,imsb2,issb2,ifsb2)
      call RDCONV (sl3,1114,isisl3,idsl3,imsl3,issl3,ifsl3)
      call RDCONV (sb3,1114,isisb3,idsb3,imsb3,issb3,ifsb3)
      call RDCONV (sl0,1112,isisl0,idsl0,imsl0,issl0,ifsl0)
      call RDCONV (sb0,1112,isisb0,idsb0,imsb0,issb0,ifsb0)
      call RDCONV (sp1,1112,isisp1,idsp1,imsp1,issp1,ifsp1)
      call RDCONV (sp2,1112,isisp2,idsp2,imsp2,issp2,ifsp2)
      call RDCONV (ds,1132,isids,idds,imds,isds,ifds)

      slong(1)=idsl1+(dble(ifsl1)/1d4)
      slong(2)=idsl2+(dble(ifsl2)/1d4)
      slong(3)=idsl3+(dble(ifsl3)/1d4)
      slat(1)=idsb1+(dble(ifsb1)/1d4)
      if (isisb1.eq.2) slat(1)=-slat(1)
      slat(2)=idsb2+(dble(ifsb2)/1d4)
      if (isisb2.eq.2) slat(2)=-slat(2)
      slat(3)=idsb3+(dble(ifsb3)/1d4)
      if (isisb3.eq.2) slat(3)=-slat(3)
      l0=idsl0+(dble(ifsl0)/dble(100))
      b0=idsb0+(dble(ifsb0)/dble(100))
      if (isisb0.eq.2) b0=-b0

      orient(1)=idsp1+(dble(ifsp1)/dble(100))
      if (isisp1.eq.2) orient(1)=-orient(1)
      orient(2)=idsp2+(dble(ifsp2)/dble(100))
      if (isisp2.eq.2) orient(2)=-orient(2)

      dsol(1)=imds
      dsol(2)=isds+(dble(ifds)/dble(100))
*
      end
*
*
*
      subroutine EPHPLA1 (tjd,nul,ierr,plong,plat,xphase,fillum,dpla)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul et affichage des ephemerides des planetes (projet SOHO).
*
*     En entree :
*     t0          date julienne modifiee initiale TDB (reel dp).
*     dt          pas du calcul en jours (reel dp).
*     ipla        indice de la planete (entier).
*     nul         numero d'unite logique du fichier source (entier).
*
*     En sortie :
*     ierr        indice erreur (entier).
*                 ierr=0 pas d'erreur.
*                 ierr=10 erreur de fichier.
*                 ierr=11 erreur de date.
*                 ierr=12 erreur de corps.
*                 ierr=13 erreur de centre.
*                 ierr=14 erreur de vitesse.
*                 ierr=15 erreur de plan.
*
*     Appels subroutines et fonctions :
*     POSAST, POSMOY, POSAPP, PHASE, DIAPLA,
*     DJDATE, RDCONV, REEL, CLRSCR.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      double precision plong(5,3),plat(5,3),xphase(5),
     +     fillum(5),dpla(5,2)
      double precision i,k
*      character*16 zk

      common/IORIG/IORIG
*
      pl=0.d0
      pb=0.d0
      i=0.d0
      k=0.d0
      dp1=0.d0
      dp2=0.d0
*
*      t=t0+2400000.5d0
*
100   continue
*

      do 500, ipla=1,5
         if (ipla.lt.3) then
            icorps=ipla
         else
            icorps=ipla+1
         endif
*               
         call POSAST (tjd,icorps,2,nul,pl1,pb1,x,ierr)
         if (ierr.ne.0) return
*
         call POSMOY (tjd,icorps,2,nul,pl2,pb2,x,ierr)
         if (ierr.ne.0) return
*
         call POSAPP (tjd,icorps,2,nul,pl3,pb3,x,ierr)
         if (ierr.ne.0) return
*
         call PHASE (tjd,icorps,nul,i,k,ierr)
         if (ierr.ne.0) return
*
         dp1=DIAPLA(tjd,icorps,1,nul,ierr)
         if (ierr.ne.0) return
*
         dp2=DIAPLA(tjd,icorps,2,nul,ierr)
         if (ierr.ne.0) return
*
*         call DJDATE (t,3,ij,im,ia,nh,nm,ns)
         call RDCONV (pl1,1114,isipl1,idpl1,impl1,ispl1,ifpl1)
         call RDCONV (pb1,1114,isipb1,idpb1,impb1,ispb1,ifpb1)
         call RDCONV (pl2,1114,isipl2,idpl2,impl2,ispl2,ifpl2)
         call RDCONV (pb2,1114,isipb2,idpb2,impb2,ispb2,ifpb2)
         call RDCONV (pl3,1114,isipl3,idpl3,impl3,ispl3,ifpl3)
         call RDCONV (pb3,1114,isipb3,idpb3,impb3,ispb3,ifpb3)
         call RDCONV (i,1111,isii,idi,imi,isi,ifi)
*         call REEL (k,1,3,zk)
         call RDCONV (dp1,1132,isidp1,iddp1,imdp1,isdp1,ifdp1)
         call RDCONV (dp2,1132,isidp2,iddp2,imdp2,isdp2,ifdp2)

         plong(ipla,1)=idpl1+(dble(ifpl1)/1d4)
         plong(ipla,2)=idpl2+(dble(ifpl2)/1d4)
         plong(ipla,3)=idpl3+(dble(ifpl3)/1d4)
         plat(ipla,1)=idpb1+(dble(ifpb1)/1d4)
         if (isipb1.eq.2) plat(ipla,1)=-plat(ipla,1)
         plat(ipla,2)=idpb2+(dble(ifpb2)/1d4)
         if (isipb2.eq.2) plat(ipla,2)=-plat(ipla,2)
         plat(ipla,3)=idpb3+(dble(ifpb3)/1d4)
         if (isipb3.eq.2) plat(ipla,3)=-plat(ipla,3)

         xphase(ipla)=idi+(dble(ifi)/dble(10))
         fillum(ipla)=k

         dpla(ipla,1)=isdp1+(dble(ifdp1)/dble(100))
         dpla(ipla,2)=isdp2+(dble(ifdp2)/dble(100))
500   continue
*
      end
*
*
*
      subroutine RDCONV (r,icode,isig,ix1,ix2,ix3,ix4)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Conversion de radians en degres ou en heure.
*
*     En entree :
*     r           valeur en radians (reel dp).
*     icode       code de conversion sous la forme abcd (entier).
*                 a: code de reduction.
*                    a=1 pas de reduction.
*                    a=2 reduction entre 0 et 2pi.
*                 b: code d'unite d'angle.
*                    b=1 conversion en degre.
*                    b=2 conversion en heure.
*                 c: code de conversion.
*                    c=1 degres (ou heures).
*                    c=2 degres (ou heures) et minutes.
*                    c=3 degres (ou heures), minutes et secondes.
*                 d: nombre de decimales de la derniere unite choisie
*                    dans le code de conversion (entre 0 et 8).
*
*     En sortie :
*     isig        signe du resultat (entier).
*                 isig=0  nul.
*                 isig=1  positif.
*                 isig=2  negatif.
*     ix1         partie entiere de degre ou d'heure (entier).
*     ix2         partie entiere de minute de degre ou d'heure (entier).
*     ix3         partie entiere de seconde de degre ou d'heure (entier)
*     ix4         partie decimale de la derniere unite choisie (entier).
*                 si c=1 decimales de degre ou d'heure.
*                 si c=2 decimales de minute de degre ou d'heure.
*                 si c=3 decimales de seconde de degre ou d'heure.
*
*     Remarque :
*     Si icode ne respecte pas les regles ci-dessus le resultats est :
*     isig=0, ix1=0, ix2=0, ix3=2, ix4=0.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension ks(2),t(2)
      common/herad/herad
      common/dgrad/dgrad
      data ks/+1,-1/
      data t/360.d0,24.d0/
      isig=0
      ix1=0
      ix2=0
      ix3=0
      ix4=0
*
      if (abs(r).lt.1.d-20) return
      ia=icode/1000
      ib=mod(icode,1000)/100
      ic=mod(icode,100)/10
      id=mod(icode,10)
      if (ia.lt.1.or.ia.gt.2.or.ib.lt.1.or.ib.gt.2.or.ic.lt.1.or.
     .     ic.gt.3.or.id.lt.0.or.id.gt.8) return
*
      if (r.gt.0.d0) isig=1
      if (r.lt.0.d0) isig=2
      if (ib.eq.1) c1=1.d0/dgrad
      if (ib.eq.2) c1=1.d0/herad
      if (ic.eq.1) c2=1.d0
      if (ic.eq.2) c2=60.d0
      if (ic.eq.3) c2=3600.d0
      c3=10.d0**(id)
      x=r*c1*c2*c3+0.5d0*ks(isig)
      y=t(ib)*c2*c3
*
      if (ia.eq.2) then
         x=mod(x,y)
         if (x.lt.0.d0) x=x+y
      endif
*
      if (x.eq.0.d0) then
         isig=0
         return
      endif
      if (x.gt.0.d0) isig=1
      if (x.lt.0.d0) isig=2
*
      x=abs(x)
      x=x/c3
*
      if (ic.eq.1) then
         k1=int(x)
         k4=(x-k1)*c3
      endif
*
      if (ic.eq.2) then
         x=x/60.d0
         k1=int(x)
         x=(x-k1)*60.d0
         k2=int(x)
         x=(x-k2)*c3
         k3=int(x)
      endif
*
      if (ic.eq.3) then
         x=x/3600.d0
         k1=int(x)
         x=(x-k1)*60.d0
         k2=int(x)
         x=(x-k2)*60.d0
         k3=int(x)
         x=(x-k3)*c3
         k4=int(x)
      endif
*
      ix1=k1
      ix2=k2
      ix3=k3
      ix4=k4
*
      return
      end
