      subroutine POSAST (tdj,icorps,iplan,nul,alpha,delta,dist,ierr)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul des coordonnees astrometriques J2000 d'un corps.
*
*     En entree :
*     tdj         date julienne TDB (reel dp).
*     icorps      indice du corps (entier).
*                 icorps=01 Mercure.
*                 icorps=02 Venus.
*                 icorps=04 Mars.
*                 icorps=05 Jupiter.
*                 icorps=06 Saturne.
*                 icorps=11 Soleil.
*     iplan       code du plan de reference (enier).
*                 iplan=1 coordonnees equatoriales.
*                 iplan=2 coordonnees ecliptiques.
*     nul         numero d'unite logique du fichier des sources des    f
*                 ephemerides (entier).
*
*     En sortie :
*     alpha       ascension droite ou longitude en radian (reel dp).
*     delta       declinaison ou latitude en radian (reel dp).
*     dist        distance au centre geometrique en ua (reel dp).
*     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 :
*     BDL82, C2000, OBLIMO, ECLIEQ, RECSPH.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension v1(6),v2(9),v3(3)
      common/t2000/t2000
      common/tlua/tlua
*
      alpha=0.d0
      delta=0.d0
      dist=0.d0
*
      ierr=12
      if (icorps.lt.1.or.icorps.eq.3.or.(icorps.gt.6.and.icorps.ne.11))
     .return
      ierr=15
      if (iplan.lt.1.or.iplan.gt.2) return
*
      call C2000 (tdj,12,nul,v1,ierr)
      if (ierr.ne.0) return
      if (icorps.eq.11) then
         do i=1,9
            v2(i)=0.d0
         enddo
      else
         call BDL82 (tdj,icorps,11,2,nul,v2,ierr)
         if (ierr.ne.0) return
      endif
      v3(1)=v2(1)
      v3(2)=v2(2)
      v3(3)=v2(3)
*
      do n=1,2
         x=v3(1)-v1(1)
         y=v3(2)-v1(2)
         z=v3(3)-v1(3)
         r=sqrt(x*x+y*y+z*z)
         tau=r*tlua
         if (n.eq.1) dist=r
         v3(1)=v2(1)-v2(4)*tau+v2(7)*tau*tau/2.d0
         v3(2)=v2(2)-v2(5)*tau+v2(8)*tau*tau/2.d0
         v3(3)=v2(3)-v2(6)*tau+v2(9)*tau*tau/2.d0
      enddo
*
      do i=1,3
         v2(i)=v3(i)-v1(i)
      enddo
*
      if (iplan.eq.1) then
         epsm=OBLIMO(t2000)
         call ECLIEQ (epsm,v2,1,v3)
         call RECSPH (v3,alpha,delta,dist)
      else
         call RECSPH (v2,alpha,delta,dist)
      endif
*
      return
      end
*
*
*
      subroutine POSMOY (tdj,icorps,iplan,nul,alpha,delta,dist,ierr)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul des coordonnees moyennes de la date d'un corps.
*
*     En entree :
*     tdj         date julienne TDB (reel dp).
*     icorps      indice du corps (entier).
*                 icorps=01 Mercure.
*                 icorps=02 Venus.
*                 icorps=04 Mars.
*                 icorps=05 Jupiter.
*                 icorps=06 Saturne.
*                 icorps=11 Soleil.
*     iplan       code du plan de reference (enier).
*                 iplan=1 coordonnees equatoriales.
*                 iplan=2 coordonnees ecliptiques.
*     nul         numero d'unite logique du fichier des sources des
*                 ephemerides (entier).
*
*     En sortie :
*     alpha       ascension droite ou longitude en radian (reel dp).
*     delta       declinaison ou latitude en radian (reel dp).
*     dist        distance au centre geometrique en ua (reel dp).
*     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 :
*     C2000, OBLIMO, QPREC, REQUAM, RECLIM, ECLIEQ, RECSPH.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension v1(3),v2(3),qp(3)
      common/t2000/t2000
*
      alpha=0.d0
      delta=0.d0
      dist=0.d0
*
      ierr=12
      if (icorps.lt.1.or.icorps.eq.3.or.(icorps.gt.6.and.icorps.ne.11))
     .return
      ierr=15
      if (iplan.lt.1.or.iplan.gt.2) return
*
      call C2000 (tdj,icorps,nul,v1,ierr)
      if (ierr.ne.0) return
*
      if (iplan.eq.1) then
         epsm=OBLIMO(t2000)
         call ECLIEQ (epsm,v1,1,v2)
         call QPREC (t2000,tdj,1,qp)
         call REQUAM (v2,qp,1,v1)
         call RECSPH (v1,alpha,delta,dist)
      else
         call QPREC (t2000,tdj,2,qp)
         call RECLIM (v1,qp,1,v2)
         call RECSPH (v2,alpha,delta,dist)
      endif
*
      return
      end
*
*
*
      subroutine POSAPP (tdj,icorps,iplan,nul,alpha,delta,dist,ierr)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul des coordonnees apparentes d'un corps.
*
*     En entree :
*     tdj         date julienne TDB (reel dp).
*     icorps      indice du corps (entier).
*                 icorps=01 Mercure.
*                 icorps=02 Venus.
*                 icorps=04 Mars.
*                 icorps=05 Jupiter.
*                 icorps=06 Saturne.
*                 icorps=11 Soleil.
*     iplan       code du plan de reference (enier).
*                 iplan=1 coordonnees equatoriales.
*                 iplan=2 coordonnees ecliptiques.
*     nul         numero d'unite logique du fichier des sources des
*                 ephemerides (entier).
*
*     En sortie :
*     alpha       ascension droite ou longitude en radian (reel dp).
*     delta       declinaison ou latitude en radian (reel dp).
*     dist        distance au centre geometrique en ua (reel dp).
*     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 :
*     BDL82, C2000, OBLIMO, QPREC, NUTA, REQUAM, REQUAV, RECLIM, RECLIV,
*     ECLIEQ, RECSPH.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension v1(6),v2(9),v3(3),qp(3)
      common/t2000/t2000
      common/tlua/tlua
*
      alpha=0.d0
      delta=0.d0
      dist=0.d0
*
      ierr=12
      if (icorps.lt.1.or.icorps.eq.3.or.(icorps.gt.6.and.icorps.ne.11))
     .return
      ierr=15
      if (iplan.lt.1.or.iplan.gt.2) return
*
      call C2000 (tdj,12,nul,v1,ierr)
      if (ierr.ne.0) return
      if (icorps.eq.11) then
         do i=1,9
            v2(i)=0.d0
         enddo
      else
         call BDL82 (tdj,icorps,11,2,nul,v2,ierr)
         if (ierr.ne.0) return
      endif
      v3(1)=v2(1)
      v3(2)=v2(2)
      v3(3)=v2(3)
*
      do n=1,2
         x=v3(1)-v1(1)
         y=v3(2)-v1(2)
         z=v3(3)-v1(3)
         r=sqrt(x*x+y*y+z*z)
         tau=r*tlua
         if (n.eq.1) then
            dist=r
            v3(1)=v2(1)-v2(4)*tau+v2(7)*tau*tau/2.d0
            v3(2)=v2(2)-v2(5)*tau+v2(8)*tau*tau/2.d0
            v3(3)=v2(3)-v2(6)*tau+v2(9)*tau*tau/2.d0
         endif
      enddo
*
      v1(1)=v2(1)-v1(1)-(v2(4)-v1(4))*tau+v2(7)*tau*tau/2.d0
      v1(2)=v2(2)-v1(2)-(v2(5)-v1(5))*tau+v2(8)*tau*tau/2.d0
      v1(3)=v2(3)-v1(3)-(v2(6)-v1(6))*tau+v2(9)*tau*tau/2.d0
*
*
      call RECSPH (v1,r1,r2,r3)
*
      if (iplan.eq.1) then
         epsm=OBLIMO(t2000)
         call ECLIEQ (epsm,v1,1,v2)
         call QPREC (t2000,tdj,1,qp)
         call NUTA (tdj,3,dpsi,deps)
         epsm=OBLIMO(tdj)
         epsv=epsm+deps
         call REQUAM (v2,qp,1,v1)
         call REQUAV (v1,epsm,dpsi,deps,1,v2)
         call RECSPH (v2,alpha,delta,x)
      else
         call QPREC (t2000,tdj,2,qp)
         call NUTA (tdj,1,dpsi,deps)
         call RECLIM (v1,qp,1,v2)
         call RECLIV (v2,dpsi,1,v1)
         call RECSPH (v1,alpha,delta,x)
      endif
*
      return
      end
*
*
*
      subroutine SOLPHY (tdj,nul,l0,b0,p1,p2,ierr)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul des donnees pour les observations physiques du soleil.
*
*     En entree :
*     tdj         date julienne TDB (reel dp).
*     nul         numero d'unite logique du fichier des sources des
*                 ephemerides (entier).
*
*     En sortie :
*     l0          longitude heliographique du centre du dique solaire
*                 en radian (reel dp).
*     b0          latitude heliographique du centre du disque solaire
*                 en radian (reel dp).
*     p1          angle de position de l'extremite nord de l'axe de
*                 rotation par rapport au pole de l'equateur
*                 en radian (reel dp).
*     p2          angle de position de l'extremite nord de l'axe de
*                 rotation par rapport au pole de l'ecliptique
*                 en radian (reel dp).
*     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 :
*     POSAPP, NUTA, OBLIMO.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      double precision l0
      common/dpi/dpi
      common/dgrad/dgrad
      common/t2000/t2000
      common/sjul/sjul
*
      l0=0.d0
      b0=0.d0
      p1=0.d0
      p2=0.d0
*
      gi=7.25d0*dgrad
      go=73.d0+40.d0/60.d0+50.25d0/3600.d0*(tdj-2396758.5d0)/365.25d0
      gm=292.766d0+(2430000.5d0-tdj)*14.18439716d0
      go=mod(go*dgrad,dpi)
      gm=mod(gm*dgrad,dpi)
*
      call POSAPP (tdj,11,2,nul,xlon,xlat,xdist,ierr)
      if (ierr.ne.0) return
      call NUTA (tdj,3,dpsi,deps)
      epsv=OBLIMO(tdj)+deps
*
      gl=xlon-dpsi
      c=cos(gl-go)
      s=sin(gl-go)
*
      l0=atan2(-s*cos(gi),-c)+gm
      l0=mod(l0,dpi)
      if (l0.lt.0.d0) l0=l0+dpi
*
      b0=asin(s*sin(gi))
*
      x=atan(-cos(xlon)*tan(epsv))
      y=atan(-c*tan(gi))
      p1=x+y
      p2=y
*
      return
      end
*
*
*
      subroutine SOLSYN (tdj,nul,tdrs,nrs,ierr)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul du debut et du numero d'une rotaton synodique du Soleil.
*
*     En entree :
*     tdj         date julienne TDB (reel dp).
*     nul         numero d'unite logique du fichier des sources des
*                 ephemerides (entier).
*
*     En sortie :
*     tdrs        date julienne TDB du debut de la rotation (reel dp).
*     nrs         numero de la rotation (entier).
*     ierr        indice erreur (entier).
*                 ierr=0  pas d'erreur.
*                 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.
*                 ierr=20 erreur dans l'algorithme.
*
*     Appels subroutines et fonctions :
*     POSAPP, NUTA.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      double precision l0
      common/pi/pi
      common/dpi/dpi
      common/dgrad/dgrad
*
      data rsyn/27.2753d0/
      data prec/1.d-4/
* Changed from nmax of 5 on 1/4/96 SPP
      data nmax/10/
*
      tdrs=0.d0
      nrs=0
      t1=tdj
      gi=7.25d0*dgrad
*
      do n=1,nmax
         go=73.d0+40.d0/60.d0+50.25d0/3600.d0*(t1-2396758.5d0)/365.25d0
         gm=292.766d0+(2430000.5d0-t1)*14.18439716d0
         go=mod(go*dgrad,dpi)
         gm=mod(gm*dgrad,dpi)
         call POSAPP (t1,11,2,nul,gla,xlat,xdist,ierr)
         if (ierr.ne.0) return
         call NUTA (t1,1,dpsi,deps)
         gl=gla-dpsi
         l0=atan2(-sin(gl-go)*cos(gi),-cos(gl-go))+gm
         l0=mod(l0,dpi)
         if (l0.lt.0.d0) l0=l0+dpi
         if (n.eq.1) then
            t2=t1-(dpi-l0)*rsyn/dpi
            if (abs(l0-pi).lt.0.35d0) t2=t2+5.d0
         else
            if (l0.lt.pi) t2=t1+l0*rsyn/dpi
            if (l0.gt.pi) t2=t1-(dpi-l0)*rsyn/dpi
         endif
* Added next line 1/4/96 SPP.
         t3=t2-t1
*
         if (abs(t2-t1).lt.prec) goto 100
         t1=t2
      enddo
*
      ierr=20
      return
*
100   continue
      tdrs=t2
      nrs=((tdrs+1.d0)-2398167.4d0)/rsyn+1
      ierr=0
      return
      end
*
*
*
      double precision function DIASOL (tdj,nul,ierr)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul du diametre apparent du Soleil.
*
*     En entree :
*     tdj         date julienne TDB (reel dp).
*     nul         numero d'unite logique du fichier des sources des
*                 ephemerides (entier).
*
*     En sortie :
*     DIASOL      diametre apparent du Soleil en radian (reel dp).
*     ierr        indice erreur (entier).
*                 ierr=0  pas d'erreur.
*                 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 :
*     C2000.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension v(3)
      common/sdrad/sdrad
      common/ua/ua
      common/rsol/rsol
*
      DIASOL=0.d0
*
      call C2000 (tdj,11,nul,v,ierr)
      if (ierr.ne.0) return
      d=sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
*
      DIASOL=2.d0*asin(rsol/ua/d)
*
      return
      end
*
*
*
      double precision function DIAPLA (tdj,icorps,icode,nul,ierr)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul du diametre apparent d'une planete.
*
*     En entree :
*     tdj         date julienne TDB (reel dp).
*     icorps      indice du corps (entier).
*                 icorps=01 Mercure.
*                 icorps=02 Venus.
*                 icorps=04 Mars.
*                 icorps=05 Jupiter.
*                 icorps=06 Saturne.
*     icode       indice du code de calcul (integer).
*                 icode=1 diametre equatorial.
*                 icode=2 diametre polaire.
*                 par defaut : icode=1.
*     nul         numero d'unite logique du fichier des sources des
*                 ephemerides (entier).
*
*     En sortie :
*     DIAPLA      diametre apparent en radian (reel dp).
*     ierr        indice erreur (entier).
*                 ierr=0  pas d'erreur.
*                 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 :
*     C2000.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension v(3)
      common/sdrad/sdrad
      common/ua/ua
      common/ray/ray(2,6)
*
      DIAPLA=0.d0
*
      ierr=12
      if (icorps.lt.1.or.icorps.eq.3.or.icorps.gt.6) return
*
      k=icode
      if (k.lt.1.or.k.gt.3) k=1
*
      call C2000 (tdj,icorps,nul,v,ierr)
      if (ierr.ne.0) return
      d=sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))*ua
*
      DIAPLA=2.d0*asin(ray(k,icorps)/d)
*
      return
      end
*
*
*
      subroutine PHASE (tdj,icorps,nul,i,k,ierr)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul de la phase des planetes.
*
*     En entree :
*     tdj         date julienne TDB (reel dp).
*     icorps      indice du corps (entier).
*                 icorps=01 Mercure.
*                 icorps=02 Venus.
*                 icorps=04 Mars.
*                 icorps=05 Jupiter.
*                 icorps=06 Saturne.
*     nul         numero d'unite logique du fichier des sources des
*                 ephemerides (entier).
*
*     En sortie :
*     i           angle de phase en radian (reel dp).
*     k           phase : fraction illuminee du disque (reel dp).
*     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 :
*     POSAPP.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      double precision i,k
      common/pi/pi
*
      i=0.d0
      k=0.d0
*
      ierr=12
      if (icorps.lt.1.or.icorps.eq.3.or.icorps.gt.6) return
*
      call POSAPP (tdj,11,1,nul,as,ds,rs,ierr)
      if (ierr.ne.0) return
      call POSAPP (tdj,icorps,1,nul,ac,dc,rc,ierr)
      if (ierr.ne.0) return
*
      ss=sin(ds)
      cs=cos(ds)
      sc=sin(dc)
      cc=cos(dc)
      sa=sin(as-ac)
      ca=cos(as-ac)
*
      elong=acos(ss*sc+cs*cc*ca)
*
      i=atan(rs*sin(elong)/(rc-rs*cos(elong)))
      if (i.lt.0.d0) i=i+pi
*
      k=(1.d0+cos(i))/2.d0
      return
      end
*
*
*
      subroutine C2000 (tdj,icorps,nul,r,ierr)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul de coordonnees rectangulaires ecliptiques moyennes J2000
*     (cas de la Terre ou de Soho considere comme objet ou comme centre
*     du repere).
*
*     En entree :
*     tdj         date julienne TDB (reel dp).
*     icorps      indice du corps (entier).
*                 icorps=01 Mercure.
*                 icorps=02 Venus.
*                 icorps=04 Mars.
*                 icorps=05 Jupiter.
*                 icorps=06 Saturne.
*                 icorps=11 Soleil.
*                 icorps=12 Terre ou Soho (voir indice IORIG).
*     nul         numero d'unite logique du fichier source (entier).
*
*     En sortie :
*     r(9)        table des coordonnees J2000 (reel dp).
*                 r(1): position x (ua).
*                 r(2): position y (ua).
*                 r(3): position z (ua).
*                 r(4): derivee premiere x (ua/j).
*                 r(5): derivee premiere y (ua/j).
*                 r(6): derivee premiere z (ua/j).
*     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.
*
*     En common :
*     IORIG       indice de l'origine des coordonnees (entier).
*                 1 : Terre (coordonnees geocentriques).
*                 2 : Soho  (coordonnees sohocentriques).
*
*     Appels subroutines et fonctions :
*     BDL82, SOHO.
*
*     Remarque :
*
*     Dans cette subroutine le parametre icorps=12 peut representer
*     la Terre ou Soho. C'est l'indice IORIG qui permet de differencier
*     les deux corps.
*
*     On distingue 5 cas de coordonnees ecliptiques pour lesquels la
*     Terre ou Soho est concerne :
*     1/ Position et vitesse J2000 heliocentriques de la Terre.
*     2/ Position geocentrique J2000 du Soleil ou d'une planete.
*     3/ Position et vitesse J2000 heliocentriques de Soho.
*     4/ Position sohocentrique J2000 du Soleil.
*     5/ Position sohocentrique J2000 d'une planete.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension r(6),v1(6),v2(6)
      common/IORIG/IORIG
*
*     Position et vitesse heliocentriques ecliptiques J2000 de la Terre.
*
      if (IORIG.eq.1.and.icorps.eq.12) then
         call BDL82 (tdj,12,11,1,nul,v1,ierr)
         if (ierr.ne.0) return
         do i=1,6
            r(i)=v1(i)
         enddo
         return
      endif
*
*     Position geocentrique ecliptique J2000 du Soleil ou d'une planete.
*
      if (IORIG.eq.1.and.icorps.ne.12) then
         call BDL82 (tdj,icorps,12,0,nul,v1,ierr)
         if (ierr.ne.0) return
         do i=1,3
            r(i)=v1(i)
         enddo
         return
      endif
*
*     Position et vitesse heliocentriques ecliptiques J2000 de Soho.
*
      if (IORIG.eq.2.and.icorps.eq.12) then
         tdjm=tdj-2400000.5d0
         call SOHO (tdjm,v1,ierr)
         if (ierr.ne.0) return
         do i=1,6
            r(i)=v1(i)
         enddo
         return
      endif
*
*     Position sohocentrique ecliptique J2000 du Soleil.
*
      if (IORIG.eq.2.and.icorps.eq.11) then
         tdjm=tdj-2400000.5d0
         call SOHO (tdjm,v1,ierr)
         if (ierr.ne.0) return
         do i=1,3
            r(i)=-v1(i)
         enddo
         return
      endif
*
*     Position sohocentrique ecliptique J2000 d'une planete.
*
      if (IORIG.eq.2.and.icorps.le.6) then
         tdjm=tdj-2400000.5d0
         call BDL82 (tdj,icorps,11,0,nul,v1,ierr)
         if (ierr.ne.0) return
         call SOHO (tdjm,v2,ierr)
         if (ierr.ne.0) return
         do i=1,3
            r(i)=v1(i)-v2(i)
         enddo
         return
      endif
*
      end
*
*
*
      subroutine SOHO (tdjm,r,ierr)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul des coordonnees heliocentriques ecliptiques moyennes J2000
*     de Soho (positions et vitesses rectangulaires).
*
*     En entree :
*     tdjm        date julienne modifiee TDB (reel dp).
*
*     En sortie :
*     r(6)        coordonnees heliocentriques J2000 de Soho (reel dp).
*                 r(1) : position x (ua).
*                 r(2) : position y (ua).
*                 r(3) : position z (ua).
*                 r(4) : derivee premiere x (ua/j).
*                 r(5) : derivee premiere y (ua/j).
*                 r(6) : derivee premiere z (ua/j).
*     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 :
*     BDL82.
*
*     Remarque :
*
*     Cette subroutine est une simulation de celle qui sera fournie par
*     l'ESA.
*
*     On suppose que Soho est au point de Lagrange L2 entre le Soleil et
*     la Terre a la distance Soleil-Soho = 0.98999 ua et que la vitesse
*     angulaire heliocentrique de Soho est egale a celle de la Terre.
*
*     Deux cas se presentent :
*
*     1/ L'ESA fournit les coordonnees rectangulaires heliocentriques
*        moyennes ecliptiques J2000 de Soho (positions et vitesses) :
*        -> l'ensemble des instructions ci-dessous est a remplacer.
*
*     2/ L'ESA fournit les coordonnees rectangulaires heliocentriques
*        moyennes ecliptiques J2000 de Soho (positions seulement) :
*        -> il faut conserver les instructions ci-dessous correspondant
*           au calcul des vitesses.
*
*     Attention !
*     La subroutine BDL82 appelee ci-dessous utilise les dates juliennes
*     comme systeme de datation et le TDB comme echelle de temps.
*     Il faut en tenir compte si l'ESA utilise les dates juliennes
*     modifiees comme systeme de datation et/ou une echelle de temps
*     differente du TDB.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension r(6),v(6)

      common/RSOHO/rsoho(6)

      common/IFLAG/iflag

      sum1=sqrt(rsoho(1)**2+rsoho(2)**2+rsoho(3)**2)
      sum2=sqrt(rsoho(4)**2+rsoho(5)**2+rsoho(6)**2)

      if (sum1.lt.1.0d-6) then
*
*     Calcul des positions rectangulaires.
*     
         tdj=tdjm+2400000.5d0
         d=0.98999d0
         call BDL82 (tdj,12,11,0,10,v,ierr)
         r(1)=v(1)*d
         r(2)=v(2)*d
         r(3)=v(3)*d
*
*     Calcul des vitesses rectangulaires.
*
         tdj=tdjm+2400000.5d0
         d=0.98999d0
         call BDL82 (tdj,12,11,1,10,v,ierr)
         r(4)=v(4)*d
         r(5)=v(5)*d
         r(6)=v(6)*d
*
         iflag=1
*
      elseif (sum2.lt.1.0d-6) then
         r(1)=rsoho(1)
         r(2)=rsoho(2)
         r(3)=rsoho(3)
         tdj=tdjm+2400000.5d0
         d=0.98999d0
         call BDL82 (tdj,12,11,1,10,v,ierr)
         r(4)=v(4)*d
         r(5)=v(5)*d
         r(6)=v(6)*d
*         
         iflag=2
*
      else
         r(1)=rsoho(1)
         r(2)=rsoho(2)
         r(3)=rsoho(3)
         r(4)=rsoho(4)
         r(5)=rsoho(5)
         r(6)=rsoho(6)
*
         iflag=3
         ierr=0
*
      endif
*
      return
      end

