*     PRCBDL94
*-----------------------------------------------------------------------
*
*     Object :
*     Example program for subroutine PRCBDL94 (Fortran 77, DOS 5.0).
*     Read comments in subroutine PRCBDL94.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      character*7 escscr
      dimension qp(12)
      data sdrad/0.4848136811095360d-5/
      escscr=char(27)//char(91)//'2J'//char(27)//char(91)//'H'
*
*     Fixed epoch : 4 dates since J2000 with -200000 days step.
*     Date epoch :  4 dates since fixed epoch with +100000 days step.
*
      do 300 ief=1,4
         ef=2451545.0d0-2.d5*(ief-1)
         do 200 ied=1,4
            ed=ef+1.d5*(ied-1)
            call PRCBDL94 (ef,ed,qp)
            do 100 i=1,12
               qp(i)=qp(i)/sdrad
100         continue
C           write (*,'(2x,a)') escscr
            write (*,1000) ed,ef,qp
C           pause ' Enter a blank line to continue'
200      continue
300   continue
      stop1248
*
1000  format (2x,'PRECESSION QUANTITIES (BDL94)'//
     .        2x,'Date : JD',f14.5,
     .        20x,'Fixed epoch : JD',f14.5//
     .        2x,'sin(pi.a)*sin(PI.a) :',f12.5,' "'/
     .        2x,'sin(pi.a)*cos(PI.a) :',f12.5,' "'//
     .        2x,'pi.a :               ',f12.5,' "'/
     .        2x,'PI.a :               ',f12.5,' "'/
     .        2x,'p.a :                ',f12.5,' "'//
     .        2x,'theta.a :            ',f12.5,' "'/
     .        2x,'zeta.a :             ',f12.5,' "'/
     .        2x,'z.a :                ',f12.5,' "'//
     .        2x,'epsilon.a :          ',f12.5,' "'//
     .        2x,'omega.a :            ',f12.5,' "'/
     .        2x,'psi.a :              ',f12.5,' "'/
     .        2x,'chi.a :              ',f12.5,' "'/)
      end
*
*
*
      subroutine PRCBDL94 (epocf,epocd,qp)
*-----------------------------------------------------------------------
*
*     Object :
*     Computation of precession quantities with formulae BDL94.
*     Precessional displacements from a fixed epoch "epocf" to the epoch
*     of the date "epocd".
*
*     Reference :
*     Astron. Astrophys. 282, 663 (1994).
*
*     Authors :
*     J.L. Simon, P. Bretagnon, J. Chapront, M. Chapront-Touze,
*     G. Francou, J. Laskar (Bureau des Longitudes, Paris, France).
*
*     Input parameters :
*     epocf       Julian date fixed epoch (real double).
*     epocd       Julian date epoch of the date (real double).
*
*     Output parameters :
*     qp(12)      precessional quantities in radians (real double).
*
*     Notations :
*     For each array hereunder, the indices correspond respectively to
*     coefficients of power of time :
*               T = (fixed epoch - J2000 epoch)
*          and  t = (epoch of the date - fixed epoch)
*          in thousands of Julian years of 365250 days.
*     These coefficients are double real expressed in 1/10000 ".
*     ss    (0:5,1:5) :   sin(pi.a)*sin(PI.a)     ---->   qp(1)
*     sc    (0:5,1:5) :   sin(pi.a)*cos(PI.a)     ---->   qp(2)
*     spi   (0:5,1:4) :   pi.a                    ---->   qp(3)
*     cpi   (0:5,0:5) :   PI.a                    ---->   qp(4)
*     p     (0:5,1:6) :   p.a (precession)        ---->   qp(5)
*     theta (0:5,1:6) :   theta.a                 ---->   qp(6)
*     zeta  (0:5,1:6) :   zeta.a                  ---->   qp(7)
*     z     (0:5,1:6) :   z.a                     ---->   qp(8)
*     epsil (0:5,0:5) :   epsilon.a (obliquity)   ---->   qp(9)
*     omega (0:5,0:6) :   omega.a                 ---->   qp(10)
*     psi   (0:5,1:6) :   psi.a                   ---->   qp(11)
*     chi   (0:5,1:6) :   chi.a                   ---->   qp(12)
*
*-----------------------------------------------------------------------
*
*     Declarations.
*
      implicit double precision (a-h,o-z)
      dimension qp(12),qpl(12)
      dimension ct(0:5),st(0:6),lct(0:6,12),lst(2,12)
      dimension coef(0:5,0:6,12)
      dimension ss(0:5,5),sc(0:5,5)
      dimension spi(0:5,4),cpi(0:5,0:5),p(0:5,6)
      dimension theta(0:5,6),zeta(0:5,6),z(0:5,6)
      dimension epsil(0:5,0:5),omega(0:5,0:6),psi(0:5,6),chi(0:5,6)
      equivalence (coef(0,1,1),ss(0,1)),(coef(0,1,2),sc(0,1))
      equivalence (coef(0,1,3),spi(0,1)),(coef(0,0,4),cpi(0,0))
      equivalence (coef(0,1,5),p(0,1)),(coef(0,1,6),theta(0,1))
      equivalence (coef(0,1,7),zeta(0,1)),(coef(0,1,8),z(0,1))
      equivalence (coef(0,0,9),epsil(0,0)),(coef(0,0,10),omega(0,0))
      equivalence (coef(0,1,11),psi(0,1)),(coef(0,1,12),chi(0,1))
      save qpl,epocf0,epocd0
*
*     Initializations.
*
      data epocf0/0.d0/,epocd0/0.d0/
      data sdrad/0.4848136811095360d-5/
      data t2000/2451545.d0/
      data tunit/365250.d0/
      data lst/1,5,1,5,1,4,0,5,1,6,1,6,1,6,1,6,0,5,0,6,1,6,1,6/
      data lct/0,5,4,3,2,0,0,0,5,4,3,2,0,0,
     .         0,4,4,3,2,0,0,5,4,3,2,1,0,0,
     .         0,5,4,3,2,1,0,0,5,4,3,2,1,0,
     .         0,5,4,3,2,1,0,0,5,4,3,2,1,0,
     .         5,5,4,3,2,1,0,5,0,4,3,2,1,0,
     .         0,5,4,3,2,1,0,0,5,4,3,2,1,0/
*
*     Data for numerical expressions of precession with following constants :
*     - general precession :                   50288.200  "/1000y
*     - obliquity in J2000 :                   84381.412  "
*     - Ratio  Sun mass / Mercury mass :     6023600
*     - Ratio  Sun mass / Venus   mass :      408523.5
*     - Ratio  Sun mass / E-M Baryc. mass :   328900.5
*     - Ratio  Sun mass / Mars    mass :     3098710
*     - Ratio  Sun mass / Jupiter mass :        1047.355
*     - Ratio  Sun mass / Saturn  mass :        3498.5
*     - Ratio  Sun mass / Uranus  mass :       22869
*     - Ratio  Sun mass / Neptune mass :       19314
*
      data ss/
     .    419971.d0,  -753286.d0,   3179.d0,  3178.d0,    7.d0,   -4.d0, T0-5 t1
     .    193971.d0,     5740.d0,  -2541.d0,    -5.d0,    6.d0,    0.d0, T0-4 t2
     .     -2235.d0,      859.d0,     33.d0,    -3.d0,    0.d0,    0.d0, T0-3 t3
     .      -104.d0,       -4.d0,      2.d0,     0.d0,    0.d0,    0.d0, T0-2 t4
     .         2.d0,        0.d0,      0.d0,     0.d0,    0.d0,    0.d0/ T0   t5
      data sc/
     .  -4680927.d0,     -305.d0,  59967.d0,  -205.d0, -125.d0,   -2.d0, T0-5 t1
     .     51043.d0,   -31633.d0,   -326.d0,   138.d0,   -2.d0,    0.d0, T0-4 t2
     .      5223.d0,      318.d0,    -66.d0,    -4.d0,    0.d0,    0.d0, T0-3 t3
     .       -57.d0,       19.d0,     -1.d0,     0.d0,    0.d0,    0.d0, T0-2 t4
     .        -1.d0,        0.d0,      0.d0,     0.d0,    0.d0,    0.d0/ T0   t5
      data spi/
     .   4699729.d0,   -67011.d0,    448.d0,   -19.d0,   -1.d0,    0.d0, T0-4 t1
     .    -33505.d0,      448.d0,    -28.d0,    -2.d0,    1.d0,    0.d0, T0-4 t2
     .     -1237.d0,       -4.d0,     -2.d0,     1.d0,    0.d0,    0.d0, T0-3 t3
     .         3.d0,       -1.d0,      1.d0,     0.d0,    0.d0,    0.d0/ T0-2 t4
      data cpi/
     .6295434330.d0,329296590.d0, 953520.d0,   -50.d0,-4590.d0, -100.d0, T0-5 t0
     . -86792700.d0,  -158510.d0,  -1130.d0, -4480.d0, -190.d0,    0.d0, T0-4 t1
     .    153420.d0,     -190.d0,  -4320.d0,  -230.d0,    0.d0,    0.d0, T0-3 t2
     .        50.d0,    -2080.d0,   -150.d0,     0.d0,    0.d0,    0.d0, T0-2 t3
     .      -370.d0,      -50.d0,      0.d0,     0.d0,    0.d0,    0.d0, T0-1 t4
     .       -10.d0,        0.d0,      0.d0,     0.d0,    0.d0,    0.d0/ T0   t5
      data p/
     . 502882000.d0,  2224045.d0,   2095.d0, -9408.d0,  -90.d0,   10.d0, T0-5 t1
     .   1112022.d0,     2095.d0, -14111.d0,  -180.d0,   26.d0,    0.d0, T0-4 t2
     .       773.d0,    -9410.d0,   -180.d0,    35.d0,    0.d0,    0.d0, T0-3 t3
     .     -2353.d0,      -90.d0,     26.d0,     0.d0,    0.d0,    0.d0, T0-2 t4
     .       -18.d0,       10.d0,      0.d0,     0.d0,    0.d0,    0.d0, T0-1 t5
     .         2.d0,        0.d0,      0.d0,     0.d0,    0.d0,    0.d0/ T0   t6
      data theta/
     . 200420207.d0,  -853131.d0,  -2111.d0,  3642.d0,    8.d0,   -5.d0, T0-5 t1
     .   -426566.d0,    -2111.d0,   5463.d0,    17.d0,  -12.d0,    0.d0, T0-4 t2
     .   -418238.d0,      359.d0,     27.d0,    -1.d0,    0.d0,    0.d0, T0-3 t3
     .      -731.d0,       19.d0,      9.d0,     0.d0,    0.d0,    0.d0, T0-2 t4
     .      -127.d0,       11.d0,      0.d0,     0.d0,    0.d0,    0.d0, T0-1 t5
     .         4.d0,        0.d0,      0.d0,     0.d0,    0.d0,    0.d0/ T0   t6
      data zeta/
     . 230609097.d0,  1397495.d0,    -38.d0, -5918.d0,  -37.d0,    7.d0, T0-5 t1
     .    302226.d0,    -2523.d0,  -3840.d0,   -14.d0,    7.d0,    0.d0, T0-4 t2
     .    180183.d0,    -1326.d0,      6.d0,     5.d0,    0.d0,    0.d0, T0-3 t3
     .      -583.d0,       -1.d0,      7.d0,     0.d0,    0.d0,    0.d0, T0-2 t4
     .      -285.d0,        0.d0,      0.d0,     0.d0,    0.d0,    0.d0, T0   t5
     .        -2.d0,        0.d0,      0.d0,     0.d0,    0.d0,    0.d0/ T0   t6
      data z/
     . 230609097.d0,  1397495.d0,    -38.d0, -5918.d0,  -37.d0,    7.d0, T0-5 t1
     .   1095270.d0,     2446.d0, -13913.d0,  -134.d0,   26.d0,    0.d0, T0-4 t2
     .    182667.d0,   -11400.d0,   -173.d0,    44.d0,    0.d0,    0.d0, T0-3 t3
     .     -2821.d0,      -93.d0,     32.d0,     0.d0,    0.d0,    0.d0, T0-2 t4
     .      -301.d0,        6.d0,      0.d0,     0.d0,    0.d0,    0.d0, T0-1 t5
     .        -1.d0,        0.d0,      0.d0,     0.d0,    0.d0,    0.d0/ T0   t6
      data epsil/
     . 843814120.d0, -4680927.d0,   -152.d0, 19989.d0,  -51.d0,  -25.d0, T0-5 t0
     .  -4680927.d0,     -305.d0,  59967.d0,  -205.d0, -125.d0,   -2.d0, T0-5 t1
     .      -152.d0,    59967.d0,   -308.d0,  -250.d0,   -6.d0,    0.d0, T0-4 t2
     .     19989.d0,     -205.d0,   -250.d0,    -8.d0,    0.d0,    0.d0, T0-3 t3
     .       -51.d0,     -125.d0,     -6.d0,     0.d0,    0.d0,    0.d0, T0-2 t4
     .       -25.d0,       -2.d0,      0.d0,     0.d0,    0.d0,    0.d0/ T0-1 t5
      data omega/
     . 843814120.d0, -4680927.d0,   -152.d0, 19989.d0,  -51.d0,  -25.d0, T0-5 t0
     .         0.d0,        0.d0,      0.d0,     0.d0,    0.d0,    0.d0,      t1
     .     51294.d0,   -91954.d0,    298.d0,   389.d0,    2.d0,    0.d0, T0-4 t2
     .    -77276.d0,      235.d0,    987.d0,    -1.d0,    0.d0,    0.d0, T0-3 t3
     .       -48.d0,      954.d0,     -7.d0,     0.d0,    0.d0,    0.d0, T0-2 t4
     .       333.d0,       -9.d0,      0.d0,     0.d0,    0.d0,    0.d0, T0-1 t5
     .        -3.d0,        0.d0,      0.d0,     0.d0,    0.d0,    0.d0/ T0   t6
      data psi/
     . 503850672.d0,   492595.d0,  -1344.d0, -2115.d0,   17.d0,    3.d0, T0-5 t1
     .  -1072374.d0,   -10919.d0,  13673.d0,   137.d0,  -28.d0,    0.d0, T0-4 t2
     .    -11424.d0,    26425.d0,     87.d0,  -111.d0,    0.d0,    0.d0, T0-3 t3
     .     13279.d0,     -110.d0,   -170.d0,     0.d0,    0.d0,    0.d0, T0-2 t4
     .       -94.d0,     -123.d0,      0.d0,     0.d0,    0.d0,    0.d0, T0-1 t5
     .       -35.d0,        0.d0,      0.d0,     0.d0,    0.d0,    0.d0/ T0   t6
      data chi/
     .   1055794.d0, -1888214.d0,  -1888.d0, +7950.d0,  101.d0,   -9.d0, T0-5 t1
     .  -2381379.d0,   -10910.d0,  30291.d0,  +290.d0,  -59.d0,    0.d0, T0-4 t2
     .    -12117.d0,    39055.d0,    229.d0,  -159.d0,    0.d0,    0.d0, T0-3 t3
     .     17024.d0,      -38.d0,   -214.d0,     0.d0,    0.d0,    0.d0, T0-2 t4
     .       -77.d0,     -145.d0,      0.d0,     0.d0,    0.d0,    0.d0, T0-1 t5
     .       -40.d0,        0.d0,      0.d0,     0.d0,    0.d0,    0.d0/ T0   t6
*
      if (epocf.ne.epocf0.or.epocd.ne.epocd0) then
*
*        Time parameters
*
         ct(0)=1.d0
         st(0)=1.d0
         ct(1)=(epocf-t2000)/tunit
         st(1)=(epocd-epocf)/tunit
         do 100 i=2,5
            ct(i)=ct(1)*ct(i-1)
100      continue
         do 110 i=2,6
            st(i)=st(1)*st(i-1)
110      continue
*
*        Computation of precession quantities.
*
         do 220 k=1,12
            qpl(k)=0.d0
            jmin=lst(1,k)
            jmax=lst(2,k)
            if (epocd.eq.epocf) jmax=0
            do 210 j=jmin,jmax
               q=0.d0
               imax=lct(j,k)
               if (epocf.eq.t2000) imax=0
               do 200 i=0,imax
                  q=q+coef(i,j,k)*ct(i)
200            continue
               qpl(k)=qpl(k)+q*st(j)
210         continue
            qpl(k)=qpl(k)*1.d-4*sdrad
220      continue
         epocf0=epocf
         epocd0=epocd
*
      endif
*
*     Results.
*
      do 300 i=1,12
         qp(i)=qpl(i)
300   continue
*
      return
      end
