      BLOCK DATA
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Constantes partagees.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      common/t2000/t2000
      common/sjul/sjul
      common/ua/ua
      common/rmlt/rmlt
      common/tlua/tlua
      common/rsol/rsol
      common/ray/ray(2,6)
      common/rfk5/rfk5(3,3)
*
      common/sdrad/sdrad
      common/dgrad/dgrad
      common/herad/herad
      common/dpi/dpi
      common/pi/pi
*
      data t2000/2451545.d0/
      data sjul/36525.d0/
      data ua/149597870.d0/
      data rmlt/0.01230002d0/
      data tlua/0.577551831d-2/
      data rsol/696000.d0/
      data ray/
     . 2439.7d0,  2439.7d0,
     . 6051.9d0,  6051.9d0,
     . 6378.14d0, 6356.75d0,
     . 3397.d0,   3375.d0,
     . 71492.d0,  66854.d0,
     . 60268.d0,  54364.d0/
      data rfk5/
     .  1.000000000000d0, -0.000000440360d0, +0.000000190919d0,
     . +0.000000440360d0,  1.000000000000d0, -0.000000188592d0,
     . -0.000000190919d0, +0.000000188593d0,  1.000000000000d0/
*
      data pi/3.141592653589793d0/
      data dpi/6.283185307179586d0/
      data sdrad/0.4848136811095360d-5/
      data dgrad/0.1745329251994330d-1/
      data herad/0.2617993877991495d0/
*
      end
*
*
*
      subroutine BDL82 (tdj,icorps,icent,ider,nul,r,ierr)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul de coordonnees rectangulaires ecliptiques J2000 moyennes
*     des principaux corps du systeme solaire a partir des ephemerides
*     Tchebychev sources BDL82 (systeme de reference FK5).
*
*     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.
*     icent       indice du centre du repere (entier).
*                 icent=11  Soleil.
*                 icent=12  Terre.
*     ider        indice derivees (entier).
*                 ider=0 : positions sans derivees.
*                 ider=1 : positions et derivees premieres.
*                 ider=2 : positions et derivees premieres et secondes.
*     nul         numero d'unite logique du fichier source (entier).
*
*     En sortie :
*     r(9)        table des coordonnees rectangulaires (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).
*                 r(7): derivee seconde x (ua/j/j).
*                 r(8): derivee seconde y (ua/j/j).
*                 r(9): derivee seconde z (ua/j/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.
*
*     Remarque :
*     On utilise un fichier source d'ephemerides BDL82 qui doit etre
*     ouvert dans la procedure appelante (numero d'unite logique nul).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension r(9),v(9),resul(9,3),ic(3),ic0(3),tdj0(3),ider0(3)
      dimension bloc(638),coef(636),loc(3,16)
      dimension pc(15),dc(15),ac(15)
      common/rmlt/rmlt
      common/rfk5/rfk5(3,3)
      equivalence (coef(1),bloc(3))
*
      data nul0/0/,irec0/0/
      data ic/0,3,10/,ic0/3*0/,tdj0/3*0.d0/,ider0/3*0/
      data pc/1.d0,2.d0,0.d0,12*0.d0/,npc0/0/,pc2/2.d0/
      data dc/0.d0,1.d0,0.d0,12*0.d0/,ncd0/0/,dc3/100.d0/
      data ac/0.d0,0.d0,4.d0,12*0.d0/,nca0/0/,ac4/100.d0/
*
      if (nul.ne.nul0) then
         read (nul,rec=1) ident,tdeb,tfin,delta,nbloc,ndmot,loc
         ierr=10
         if (ident.ne.2) goto 200
         f=-rmlt/(1.d0+rmlt)
      endif
*
      ierr=11
      if (tdj.lt.tdeb.or.tdj.gt.tfin) goto 200
      ierr=12
      if (icorps.lt.1.or.icorps.eq.3.or.
     .    icorps.gt.6.and.icorps.ne.11.and.icorps.ne.12) goto 200
      ierr=13
      if (icent.lt.11.or.icent.gt.12.or.icent.eq.icorps) goto 200
      ierr=14
      if (ider.lt.0.or.ider.gt.2) goto 200
      nder=3*(ider+1)
*
      irec=(tdj-tdeb)/delta+2
      if (tdj.eq.tfin) irec=irec-1
      if (nul.ne.nul0.or.irec.ne.irec0) then
         read (nul,rec=irec) (bloc(i),i=1,ndmot)
         dj1=bloc(1)
         dj2=bloc(2)
         irec0=irec
      endif
      tf=(tdj-dj1)/delta
*
      if (icent.eq.12.and.icorps.eq.11.or.
     .    icent.eq.11.and.icorps.eq.12) then
          n1=2
          n2=3
      endif
      if (icent.eq.12.and.icorps.ne.11) then
         ic(1)=icorps
         n1=1
         n2=3
      endif
      if (icent.eq.11.and.icorps.ne.12) then
         ic(1)=icorps
         n1=1
         n2=1
      endif
*
      do 100 n=n1,n2
*
         if (nul.eq.nul0.and.ic(n).eq.ic0(n).and.tdj.eq.tdj0(n).and.
     .       ider.le.ider0(n)) goto 100
*
         ic0(n)=ic(n)
         tdj0(n)=tdj
         ider0(n)=ider
*
         k=ic(n)
         tp=tf*loc(3,k)
         ncf=loc(2,k)
         ierr=12
         if (ncf.eq.0) goto 200
         l=loc(1,k)+3*ncf*idint(tp-int(tf))
*
         cp=2.d0*(tp-int(tp)+int(tf))-1.d0
         pc(2)=cp
         if (pc(2).ne.pc2.or.ncf.gt.npc0) then
            npc0=ncf
            pc2=pc(2)
            dt=cp+cp
            do i=3,ncf
               pc(i)=dt*pc(i-1)-pc(i-2)
            enddo
         endif
         do i=1,3
            resul(i,n)=0.d0
            do j=1,ncf
               jp=ncf-j+1
               jt=l+ncf*(i-1)+jp-1
               resul(i,n)=resul(i,n)+pc(jp)*coef(jt)
            enddo
         enddo
         if (ider.eq.0) goto 100
*
         cd=2.d0*loc(3,k)/delta
         dc(3)=dt+dt
         if (dc(3).ne.dc3.or.ncf.gt.ncd0) then
            ncd0=ncf
            dc3=dc(3)
            do i=4,ncf
               dc(i)=dt*dc(i-1)+2.d0*pc(i-1)-dc(i-2)
            enddo
         endif
         do i=1,3
            resul(i+3,n)=0.d0
            do j=2,ncf
               jp=ncf-j+2
               jt=l+ncf*(i-1)+jp-1
               resul(i+3,n)=resul(i+3,n)+dc(jp)*coef(jt)
            enddo
            resul(i+3,n)=resul(i+3,n)*cd
         enddo
         if (ider.eq.1) goto 100
*
         ca=cd*cd
         ac(4)=24.d0*cp
         if (ac(4).ne.ac4.or.ncf.gt.nca0) then
            nca0=ncf
            ac4=ac(4)
            do i=5,ncf
               ac(i)=dt*ac(i-1)+4.d0*dc(i-1)-ac(i-2)
            enddo
         endif
         do i=1,3
            resul(i+6,n)=0.D0
            do j=3,ncf
               jp=ncf-j+3
               jt=l+ncf*(i-1)+jp-1
               resul(i+6,n)=resul(i+6,n)+ac(jp)*coef(jt)
            enddo
            resul(i+6,n)=resul(i+6,n)*ca
         enddo
*
100   continue
*
      if (icent.eq.12) then
         do i=1,nder
            v(i)=-resul(i,2)-f*resul(i,3)
            if (icorps.ne.11) v(i)=v(i)+resul(i,1)
         enddo
      endif
*
      if (icent.eq.11) then
         do i=1,nder
            if (icorps.eq.12) v(i)=resul(i,2)+f*resul(i,3)
            if (icorps.ne.12) v(i)=resul(i,1)
         enddo
      endif
*
      do n=0,ider
         k=3*n
         do i=1,3
            r(i+k)=0.d0
            do j=1,3
               r(i+k)=r(i+k)+rfk5(i,j)*v(j+k)
            enddo
         enddo
      enddo
*
      ierr=0
      nul0=nul
      return
*
200   continue
      nul0=0
      return
      end
*
*
*
      double precision function OBLIMO (tdj)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul de l'obliquite moyenne.
*
*     En entree :
*     tdj         date julienne TDB (reel dp).
*
*     En sortie :
*     OBLIMO      obliquite moyenne en radian (reel dp).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      common/t2000/t2000
      common/sjul/sjul
      common/sdrad/sdrad
*
      t=(tdj-t2000)/sjul
      t2=t*t
      t3=t2*t
*
      OBLIMO=(84381.448d0-46.815d0*t-0.59d-3*t2+0.1813d-2*t3)*sdrad
*
      return
      end
*
*
*
      subroutine QPREC (ti,tf,kprc,qp)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul de quantites liees aux deplacements precessionnels de
*     l'equateur et de l'ecliptique entre deux epoques donnees.
*
*     En entree :
*     ti          date julienne TDB de l'epoque initale (reel dp).
*     tf          date julienne TDB de l'epoque finale (reel dp).
*     kprc        code precession (entier).
*                 kprc=1 precession equatoriale.
*                 kprc=2 precession ecliptique.
*
*     En sortie :
*     qp(3)       quantites liees au deplacement precessionnel des
*                 plans de reference entre l'epoque initale et
*                 l'epoque finale (reel dp).
*
*     Notations :
*     precession equatoriale          precession ecliptique
*     qp(1) = theta (radian).         qp(1) = petit pi (radian).
*     qp(2) = dzeta (radian).         qp(2) = grand pi (radian).
*     qp(3) = z     (radian).         qp(3) = p        (radian).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension qp(3)
      common/t2000/t2000
      common/sjul/sjul
      common/sdrad/sdrad
      data ti0/0.d0/,tf0/0.d0/,kp0/0/
      data q1/0.d0/,q2/0.d0/,q3/0.d0/
*
      if (ti.eq.tf.or.(kprc.ne.1.and.kprc.ne.2)) then
         qp(1)=0.d0
         qp(2)=0.d0
         qp(3)=0.d0
         return
      endif
*
      if (ti.eq.ti0.and.tf.eq.tf0.and.kprc.eq.kp0) then
         qp(1)=q1
         qp(2)=q2
         qp(3)=q3
         return
      endif
*
      if (kprc.eq.1) then
         t=(ti-t2000)/sjul
         t2=t*t
         theta1=(2004.3109d0-0.8533d0*t-0.217d-3*t2)*sdrad
         theta2=(-0.42665d0-0.217d-3*t)*sdrad
         theta3=-0.41833d-1*sdrad
         dzeta1=(2306.2181d0+1.39656d0*t-0.139d-3*t2)*sdrad
         dzeta2=(0.30188d0-0.345d-3*t)*sdrad
         dzeta3=0.17998d-1*sdrad
         za1=dzeta1
         za2=(1.09468d0+0.66d-4*t)*sdrad
         za3=0.18203d-1*sdrad
         pt=(tf-ti)/sjul
         pt2=pt*pt
         pt3=pt*pt2
         q1=theta1*pt+theta2*pt2+theta3*pt3
         q2=dzeta1*pt+dzeta2*pt2+dzeta3*pt3
         q3=za1*pt+za2*pt2+za3*pt3
         qp(1)=q1
         qp(2)=q2
         qp(3)=q3
      endif
*
      if (kprc.eq.2) then
         t=(ti-t2000)/sjul
         t2=t*t
         ppia1=(47.0029d0-0.6603d-1*t+0.598d-3*t2)*sdrad
         ppia2=(-0.3302d-1+0.598d-3*t)*sdrad
         ppia3=0.6d-4*sdrad
         gpia0=(629554.982d0+3289.4789d0*t+0.60622d0*t2)*sdrad
         gpia1=(-869.8089d0-0.50491d0*t)*sdrad
         gpia2=0.3536d-1*sdrad
         pa1=(5029.0966d0+2.22226d0*t-0.42d-4*t2)*sdrad
         pa2=(1.11113d0-0.42d-4*t)*sdrad
         pa3=-0.6d-5*sdrad
         pt=(tf-ti)/sjul
         pt2=pt*pt
         pt3=pt*pt2
         q1=ppia1*pt+ppia2*pt2+ppia3*pt3
         q2=gpia0+gpia1*pt+gpia2*pt2
         q3=pa1*pt+pa2*pt2+pa3*pt3
         qp(1)=q1
         qp(2)=q2
         qp(3)=q3
      endif
*
      ti0=ti
      tf0=tf
      kp0=kprc
*
      return
      end
*
*
*
      subroutine NUTA (tdj,knut,dpsi,deps)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul de la nutation en longitude et en obliquite.
*
*     En entree :
*     tdj         date julienne TDB (reel dp).
*     knut        code nutation (entier).
*                 knut=1 nutation en longitude uniquement.
*                 knut=2 nutation en obliquite uniquement.
*                 knut=3 nutation en longitude et en obliquite.
*
*     En sortie :
*     dpsi        nutation en longitude en radian (reel dp).
*     deps        nutation en obliquite en radian (reel dp).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      double precision l
      common/t2000/t2000
      common/sjul/sjul
      common/sdrad/sdrad
      common/dgrad/dgrad
      data t1900/2415020.0d0/
      data dpsi0/0.d0/,deps0/0.d0/
      data t/0.d0/,tj0/0.d0/,tj1/0.d0/,tj2/0.d0/
      data l/0.d0/,p/0.d0/,f/0.d0/,d/0.d0/,g/0.d0/
*
      if (knut.lt.1.or.knut.gt.3) then
         dpsi=0.d0
         deps=0.d0
         return
      endif
*
      if (tdj.ne.tj0) then
         tj=tdj-t1900
         s=tj/sjul
         s2=s*s
         s3=s2*s
         l=296.104608d0+13.0649924465d0*tj+0.9192d-2*s2+0.143d-4*s3
         l=mod(l,360.d0)*dgrad
         p=358.475833d0+0.9856002669d0*tj-0.149d-3*s2-0.331d-5*s3
         p=mod(p,360.d0)*dgrad
         f=11.250889d0+13.2293504490d0*tj-0.3211d-2*s2-0.341d-6*s3
         f=mod(f,360.d0)*dgrad
         d=350.737681d0+12.1907491914d0*tj-0.1436d-2*s2+0.19d-5*s3
         d=mod(d,360.d0)*dgrad
         g=259.183275d0-0.0529539222d0*tj+0.2078d-2*s2+0.2d-5*s3
         g=mod(g,360.d0)*dgrad
         tj0=tdj
      endif
*
*
      if ((knut.eq.1.or.knut.eq.3).and.tdj.ne.tj1) then
*
         t=(tdj-t2000)/sjul/10.d0
*
         u=(-171996-1742*t)*sin(g)
     .    +(2062+2*t)*sin(2*g)
     .    +46*sin(-2*l+2*f+g)
     .    +11*sin(2*l-2*f)
     .    -3*sin(-2*l+2*f+2*g)
     .    -3*sin(l-p-d)
     .    -2*sin(-2*p+2*f-2*d+g)
     .    +sin(2*l-2*f+g)
     .    +(-13187-16*t)*sin(2*f-2*d+2*g)
     .    +(1426-34*t)*sin(p)
*
         u=u+(-517+12*t)*sin(p+2*f-2*d+2*g)
     .    +(217-5*t)*sin(-p+2*f-2*d+2*g)
     .    +(129+t)*sin(2*f-2*d+g)
     .    +48*sin(2*l-2*d)
     .    -22*sin(2*f-2*d)
     .    +(17-t)*sin(2*p)
     .    -15*sin(p+g)
     .    +(-16+t)*sin(2*p+2*f-2*d+2*g)
     .    -12*sin(-p+g)
     .    -6*sin(-2*l+2*d+g)
*
         u=u-5*sin(-p+2*f-2*d+g)
     .    +4*sin(2*l-2*d+g)
     .    +4*sin(p+2*f-2*d+g)
     .    -4*sin(l-d)
     .    +sin(2*l+p-2*d)
     .    +sin(-2*f+2*d+g)
     .    -sin(p-2*f+2*d)
     .    +sin(p+2*g)
     .    +sin(-l+d+g)
     .    -sin(p+2*f-2*d)
*
         u=u+(-2274-2*t)*sin(2*f+2*g)
     .    +(712+t)*sin(l)
     .    +(-386-4*t)*sin(2*f+g)
     .    -301*sin(l+2*f+2*g)
     .    -158*sin(l-2*d)
     .    +123*sin(-l+2*f+2*g)
     .    +63*sin(2*d)
     .    +(63+t)*sin(l+g)
     .    +(-58-t)*sin(-l+g)
     .    -59*sin(-l+2*f+2*d+2*g)
*
         u=u-51*sin(l+2*f+g)
     .    -38*sin(2*f+2*d+2*g)
     .    +29*sin(2*l)
     .    +29*sin(l+2*f-2*d+2*g)
     .    -31*sin(2*l+2*f+2*g)
     .    +26*sin(2*f)
     .    +21*sin(-l+2*f+g)
     .    +16*sin(-l+2*d+g)
     .    -13*sin(l-2*d+g)
     .    -10*sin(-l+2*f+2*d+g)
*
         u=u-7*sin(l+p-2*d)
     .    +7*sin(p+2*f+2*g)
     .    -7*sin(-p+2*f+2*g)
     .    -8*sin(l+2*f+2*d+2*g)
     .    +6*sin(l+2*d)
     .    +6*sin(2*l+2*f-2*d+2*g)
     .    -6*sin(2*d+g)
     .    -7*sin(2*f+2*d+g)
     .    +6*sin(l+2*f-2*d+g)
     .    -5*sin(-2*d+g)
*
         u=u+5*sin(l-p)
     .    -5*sin(2*l+2*f+g)
     .    -4*sin(p-2*d)
     .    +4*sin(l-2*f)
     .    -4*sin(d)
     .    -3*sin(l+p)
     .    +3*sin(l+2*f)
     .    -3*sin(l-p+2*f+2*g)
     .    -3*sin(-l-p+2*f+2*d+2*g)
     .    -2*sin(-2*l+g)
*
         u=u-3*sin(3*l+2*f+2*g)
     .    -3*sin(-p+2*f+2*d+2*g)
     .    +2*sin(l+p+2*f+2*g)
     .    -2*sin(-l+2*f-2*d+g)
     .    +2*sin(2*l+g)
     .    -2*sin(l+2*g)
     .    +2*sin(3*l)
     .    +2*sin(2*f+d+2*g)
     .    +sin(-l+2*g)
     .    -sin(l-4*d)
*
         u=u+sin(-2*l+2*f+2*d+2*g)
     .    -2*sin(-l+2*f+4*d+2*g)
     .    -sin(2*l-4*d)
     .    +sin(l+p+2*f-2*d+2*g)
     .    -sin(l+2*f+2*d+g)
     .    -sin(-2*l+2*f+4*d+2*g)
     .    +sin(-l+4*f+2*g)
     .    +sin(l-p-2*d)
     .    +sin(2*l+2*f-2*d+g)
     .    -sin(2*l+2*f+2*d+2*g)
*
         u=u-sin(l+2*d+g)
     .    +sin(4*f-2*d+2*g)
     .    +sin(3*l+2*f-2*d+2*g)
     .    -sin(l+2*f-2*d)
     .    +sin(p+2*f+g)
     .    +sin(-l-p+2*d+g)
     .    -sin(-2*f+g)
     .    -sin(2*f-d+2*g)
     .    -sin(p+2*d)
     .    -sin(l-2*f-2*d)
*
         u=u-sin(-p+2*f+g)
     .    -sin(l+p-2*d+g)
     .    -sin(l-2*f+2*d)
     .    +sin(2*l+2*d)
     .    -sin(2*f+4*d+2*g)
     .    +sin(p+d)
*
         dpsi=u*sdrad/1.d4
         tj1=tdj
         dpsi0=dpsi
*
      else
*
         if (knut.ne.2) dpsi=dpsi0
*
      endif
*
*
      if ((knut.eq.2.or.knut.eq.3).and.tdj.ne.tj2) then
*
         t=(tdj-t2000)/sjul/10.d0
*
         v=(92025+89*t)*cos(g)
     .    +(-895+5*t)*cos(2*g)
     .    -24*cos(-2*l+2*f+g)
     .    +cos(-2*l+2*f+2*g)
     .    +cos(-2*p+2*f-2*d+g)
     .    +(5736-31*t)*cos(2*f-2*d+2*g)
     .    +(54-t)*cos(p)
     .    +(224-6*t)*cos(p+2*f-2*d+2*g)
     .    +(-95+3*t)*cos(-p+2*f-2*d+2*g)
     .    -70*cos(2*f-2*d+g)
*
         v=v+cos(2*l-2*d)
     .    +9*cos(p+g)
     .    +7*cos(2*p+2*f-2*d+2*g)
     .    +6*cos(-p+g)
     .    +3*cos(-2*l+2*d+g)
     .    +3*cos(-p+2*f-2*d+g)
     .    -2*cos(2*l-2*d+g)
     .    -2*cos(p+2*f-2*d+g)
     .    +(977-5*t)*cos(2*f+2*g)
     .    -7*cos(l)
*
         v=v+200*cos(2*f+g)
     .    +(129-t)*cos(l+2*f+2*g)
     .    -cos(l-2*d)
     .    -53*cos(-l+2*f+2*g)
     .    -2*cos(2*d)
     .    -33*cos(l+g)
     .    +32*cos(-l+g)
     .    +26*cos(-l+2*f+2*d+2*g)
     .    +27*cos(l+2*f+g)
     .    +16*cos(2*f+2*d+2*g)
*
         v=v-cos(2*l)
     .    -12*cos(l+2*f-2*d+2*g)
     .    +13*cos(2*l+2*f+2*g)
     .    -cos(2*f)
     .    -10*cos(-l+2*f+g)
     .    -8*cos(-l+2*d+g)
     .    +7*cos(l-2*d+g)
     .    +5*cos(-l+2*f+2*d+g)
     .    -3*cos(p+2*f+2*g)
     .    +3*cos(-p+2*f+2*g)
*
         v=v+3*cos(l+2*f+2*d+2*g)
     .    -3*cos(2*l+2*f-2*d+2*g)
     .    +3*cos(2*d+g)
     .    +3*cos(2*f+2*d+g)
     .    -3*cos(l+2*f-2*d+g)
     .    +3*cos(-2*d+g)
     .    +3*cos(2*l+2*f+g)
     .    +cos(l-p+2*f+2*g)
     .    +cos(-l-p+2*f+2*d+2*g)
     .    +cos(-2*l+g)
*
         v=v+cos(3*l+2*f+2*g)
     .    +cos(-p+2*f+2*d+2*g)
     .    -cos(l+p+2*f+2*g)
     .    +cos(-l+2*f-2*d+g)
     .    -cos(2*l+g)
     .    +cos(l+2*g)
     .    -cos(2*f+d+2*g)
     .    -cos(-l+2*g)
     .    -cos(-2*l+2*f+2*d+2*g)
     .    +cos(-l+2*f+4*d+2*g)
*
         v=v-cos(l+p+2*f-2*d+2*g)
     .    +cos(l+2*f+2*d+g)
     .    +cos(-2*l+2*f+4*d+2*g)
     .    -cos(2*l+2*f-2*d+g)
*
         deps=v*sdrad/1.d4
         tj2=tdj
         deps0=deps
*
      else
*
         if (knut.ne.1) deps=deps0
*
      endif
*
      return
      end
*
*
*
      subroutine REQUAM (v,qp,ipva,w)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul de coordonnees rectangulaires equatoriales moyennes d'une
*     epoque a partir de coordonnees rectangulaires equatoriales
*     moyennes d'une autre epoque.
*
*     En entree :
*     v(9)        coordonnees rectangulaires equatoriales moyennes
*                 de l'epoque initiale (reel dp).
*                 v(i),i=1,3 : positions (ua).
*                 v(i),i=4,6 : vitesses (ua/j).
*                 v(i),i=7,9 : accelerations (ua/j/j).
*     qp(3)       quantites liees au deplacement precessionnel de
*                 l'equateur entre l'epoque initiale et l'epoque
*                 finale (reel dp).
*     ipva        indice position / vitesse / acceleration (entier).
*                 1 : positions.
*                 2 : positions et vitesses.
*                 3 : positions, vitesses et accelerations.
*
*     En sortie :
*     w(9)        coordonnees rectangulaires equatoriales moyennes
*                 de l'epoque finale (reel dp).
*                 w(i),i=1,3 : positions (ua).
*                 w(i),i=4,6 : vitesses (ua/j).
*                 w(i),i=7,9 : accelerations (ua/j/j).
*
*     Notations :
*     qp(1) = theta (radian).
*     qp(2) = dzeta (radian).
*     qp(3) = z     (radian).
*
*     Remarque :
*     Veiller a la coherence d'epoques pour les quantites v et qp.
*     Les dimensions de v et w sont egales a 3*ipva (maximum 9).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension v(9),qp(3),w(9),a(3,3)
      data q1/1.d-12/,q2/1.d-12/,q3/1.d-12/
*
      if (qp(1).ne.q1.or.qp(2).ne.q2.or.qp(3).ne.q3) then
         c=cos(qp(1))
         s=sin(qp(1))
         c1=cos(qp(2))
         c2=cos(qp(3))
         s1=sin(qp(2))
         s2=sin(qp(3))
         cc=c1*c2
         cs=c1*s2
         sc=s1*c2
         ss=s1*s2
         a(1,1)=-ss+cc*c
         a(1,2)=-cs-sc*c
         a(1,3)=-s*c2
         a(2,1)=sc+cs*c
         a(2,2)=cc-ss*c
         a(2,3)=-s*s2
         a(3,1)=c1*s
         a(3,2)=-s1*s
         a(3,3)=c
         q1=qp(1)
         q2=qp(2)
         q3=qp(3)
      endif
*
      do l=1,ipva
         k=3*(l-1)
         do i=1,3
            w(i+k)=0.d0
            do j=1,3
               w(i+k)=w(i+k)+a(i,j)*v(j+k)
            enddo
         enddo
      enddo
*
      return
      end
*
*
*
      subroutine RECLIM (v,qp,ipva,w)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul de coordonnees rectangulaires ecliptiques moyennes d'une
*     epoque a partir de coordonnees rectangulaires ecliptiques
*     moyennes d'une autre epoque.
*
*     En entree :
*     v(9)        coordonnees rectangulaires ecliptiques moyennes
*                 de l'epoque initiale (reel dp).
*                 v(i),i=1,3 : positions (ua).
*                 v(i),i=4,6 : vitesses (ua/j).
*                 v(i),i=7,9 : accelerations (ua/j/j).
*     qp(3)       quantites liees au deplacement precessionnel de
*                 l'ecliptique entre l'epoque initiale et l'epoque
*                 finale (reel dp).
*     ipva        indice position / vitesse / acceleration (entier).
*                 1 : positions.
*                 2 : positions et vitesses.
*                 3 : positions, vitesses et accelerations.
*
*     En sortie :
*     w(9)        coordonnees rectangulaires ecliptiques moyennes
*                 de l'epoque finale (reel dp).
*                 w(i),i=1,3 : positions (ua).
*                 w(i),i=4,6 : vitesses (ua/j).
*                 w(i),i=7,9 : accelerations (ua/j/j).
*
*     Notations :
*     qp(1) = petit pi (radian).
*     qp(2) = grand pi (radian).
*     qp(3) = p (precession generale en longitude en radian).
*
*     Remarque :
*     Veiller a la coherence d'epoques pour les quantites v et qp.
*     Les dimensions de v et w sont egales a 3*ipva (maximum 9).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension v(9),qp(3),w(9),b(3,3)
      data q1/1.d-12/,q2/1.d-12/,q3/1.d-12/
*
      if (qp(1).ne.q1.or.qp(2).ne.q2.or.qp(3).ne.q3) then
         gl=qp(2)+qp(3)
         c1=cos(qp(2))
         c2=cos(gl)
         s1=sin(qp(2))
         s2=sin(gl)
         cc=c1*c2
         cs=c1*s2
         sc=s1*c2
         ss=s1*s2
         c=cos(qp(1))
         s=sin(qp(1))
         b(1,1)=cc+ss*c
         b(1,2)=sc-cs*c
         b(1,3)=-s2*s
         b(2,1)=cs-sc*c
         b(2,2)=ss+cc*c
         b(2,3)=c2*s
         b(3,1)=s1*s
         b(3,2)=-c1*s
         b(3,3)=c
         q1=qp(1)
         q2=qp(2)
         q3=qp(3)
      endif
*
      do l=1,ipva
         k=3*(l-1)
         do i=1,3
            w(i+k)=0.d0
            do j=1,3
               w(i+k)=w(i+k)+b(i,j)*v(j+k)
            enddo
         enddo
      enddo
*
      return
      end
*
*
*
      subroutine REQUAV (v,epsm,dpsi,deps,ipva,w)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul de coordonnees rectangulaires equatoriales vraies a partir
*     de coordonnees rectangulaires equatoriales moyennes.
*
*     En entree :
*     v(9)        coordonnees rectangulaires equatoriales moyennes
*                 (reel dp).
*                 v(i),i=1,3 : positions (ua).
*                 v(i),i=4,6 : vitesses (ua/j).
*                 v(i),i=7,9 : accelerations (ua/j/j).
*     epsm        obliquite moyenne en radian (reel dp).
*     dpsi        nutation en longitude (radian).
*     deps        nutation en obliquite (radian).
*     ipva        indice position / vitesse / acceleration (entier).
*                 1 : positions.
*                 2 : positions et vitesses.
*                 3 : positions, vitesse et accelerations.
*
*     En sortie :
*     w(9)        coordonnees rectangulaires equatoriales vraies
*                 (reel dp).
*                 w(i),i=1,3 : positions (ua).
*                 w(i),i=4,6 : vitesses (ua/j).
*                 w(i),i=7,9 : accelerations (ua/j/j).
*
*     Remarque :
*     Veiller a la coherence de dates pour les quantites v, epsm, dpsi
*     et deps.
*     Les dimensions de v et w sont egales a 3*ipva (maximum 9).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension a(3,3),v(9),w(9)
      data x/0.d0/,y/0.d0/,z/0.d0/
*
      if (epsm.ne.x.or.dpsi.ne.y.or.deps.ne.z) then
         epsv=epsm+deps
         c1=cos(epsv)
         c2=cos(dpsi)
         s1=sin(epsv)
         s2=sin(dpsi)
         cc=c1*c2
         cs=c1*s2
         sc=s1*c2
         ss=s1*s2
         s=sin(epsm)
         c=cos(epsm)
         a(1,1)=c2
         a(1,2)=-s2*c
         a(1,3)=-s2*s
         a(2,1)=cs
         a(2,2)=cc*c+s1*s
         a(2,3)=cc*s-s1*c
         a(3,1)=ss
         a(3,2)=sc*c-c1*s
         a(3,3)=sc*s+c1*c
         x=epsm
         y=dpsi
         z=deps
      endif
*
      do ii=1,ipva
         k=3*(ii-1)
         do i=1,3
            w(i+k)=0.d0
            do j=1,3
               w(i+k)=w(i+k)+a(i,j)*v(j+k)
            enddo
         enddo
      enddo
*
      return
      end
*
*
*
      subroutine RECLIV (v,dpsi,ipva,w)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Calcul de coordonnees rectangulaires ecliptiques vraies a partir
*     de coordonnees rectangulaires ecliptiques moyennes.
*
*     En entree :
*     v(9)        coordonnees rectangulaires ecliptiques moyennes
*                 (reel dp).
*                 v(i),i=1,3 : positions (ua).
*                 v(i),i=4,6 : vitesses (ua/j).
*                 v(i),i=7,9 : accelerations (ua/j/j).
*     dpsi        nutation en longitude en radian (reel dp).
*     ipva        indice position / vitesse / acceleration (entier).
*                 1 : positions.
*                 2 : positions et vitesses.
*                 3 : positions, vitesses et accelerations.
*
*     En sortie :
*     w(9)        coordonnees rectangulaires ecliptiques vraies
*                 (reel dp).
*                 w(i),i=1,3 : positions (ua).
*                 w(i),i=4,6 : vitesses (ua/j).
*                 w(i),i=7,9 : accelerations (ua/j/j).
*
*     Remarque :
*     Veiller a la coherence de dates pour les parametres v et dspi.
*     Les dimensions de v et w sont egales a 3*ipva (maximum 9).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension v(9),w(9)
      c=cos(dpsi)
      s=sin(dpsi)
*
      do j=1,ipva
         i=(j-1)*3+1
         w(i)=v(i)*c-v(i+1)*s
         w(i+1)=v(i+1)*c+v(i)*s
         w(i+2)=v(i+2)
      enddo
*
      return
      end
*
*
*
      subroutine ECLIEQ (eps,v,ipva,w)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Passage de coordonnees rectangulaires ecliptiques a des
*     coordonnees rectangulaires equatoriales.
*
*     En entree :
*     eps         obliquite en radian (reel dp).
*     v(9)        coordonnees rectangulaires ecliptiques (reel dp).
*                 v(i),i=1,3 : positions (ua).
*                 v(i),i=4,6 : vitesses (ua/j).
*                 v(i),i=7,9 : accelerations (ua/j/j).
*     ipva        indice position / vitesse / acceleration (entier).
*                 1 : positions.
*                 2 : positions et vitesses.
*                 3 : positions, vitesses et accelerations.
*
*     En sortie :
*     w(9)        coordonnees rectangulaires equatoriales (reel dp).
*                 w(i),i=1,3 : positions (ua).
*                 w(i),i=4,6 : vitesses (ua/j).
*                 w(i),i=7,9 : accelerations (ua/j/j).
*
*     Remarque :
*     Veiller a la coherence de dates pour les quantites eps et v.
*     Les dimensions de v et w sont egales a 3*ipva (maximum 9).
*     La valeur de l'obliquite eps est positive.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension v(9),w(9)
*
      c=cos(eps)
      s=sin(eps)
*
      do j=1,ipva
         i=(j-1)*3+1
         w(i)=v(i)
         w(i+1)=v(i+1)*c-v(i+2)*s
         w(i+2)=v(i+2)*c+v(i+1)*s
      enddo
*
      return
      end
*
*
*
      subroutine RECSPH (v,a1,a2,dm)
*-----------------------------------------------------------------------
*
*     ref : bdl-gf9504
*
*     Objet :
*     Passage de coordonnees rectangulaires a des coordonnees
*     spheriques.
*
*     En entree :
*     v(3)        coordonnees rectangulaires (reel dp).
*                 v(i),i=1,3 : positions x, y, z (ua).
*
*     En sortie :
*     a1          angle de longitude (rad) 0<=a1<=2pi (reel dp).
*     a2          angle de latitude (rad) -pi/2<=a2<=+pi/2 (reel dp).
*     dm          distance au centre en unites astronomiques (reel dp).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      dimension v(3)
      common/dpi/dpi
*
      dm=sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
      if (dm.eq.0.d0) then
         a1=0.d0
         a2=0.d0
         return
      endif
*
      if (v(1).eq.0.d0.and.v(2).eq.0.d0) then
         a1=0.d0
         a2=dpi/4.d0*v(3)/abs(v(3))
         return
      endif
*
      a1=atan2(v(2),v(1))
      a1=mod(a1,dpi)
      if (a1.lt.0.0) a1=a1+dpi
      a2=asin(v(3)/dm)
*
      return
      end
