*     PLANETAP
*-------------------------------------------------------------------------------
*
*     Object :
*     Example program for subroutine PLANETAP (Fortran 77, DOS 5.0).
*     Read comments in subroutine PLANETAP.
*
*-------------------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
      character*14 corps(8)
      character*7 escscr
      dimension xyz(6)
      data corps/'MERCURY','VENUS','E-M BARYCENTER','MARS',
     .           'JUPITER','SATURN','URANUS','NEPTUNE'/
      escscr=char(27)//char(91)//'2J'//char(27)//char(91)//'H'
*
*     Compute the longitudes, latitudes, radius vectors and rectangular
*     coordinates with velocities of the eight major planets by step of
*     -150000 days for 4 dates from dj0=2451545.0d0
*
      dj0=2451545.0d0
      pas=-150000.d0
      itmax=3
*
      do 110 ipla=1,8
         do 100 it=0,itmax
            dj=dj0+it*pas
C           write (*,'(2x,a)') escscr
            write (*,1000) corps(ipla)
            call PLANETAP (dj,ipla,dlong,dlat,rvec,xyz,ierr)
            if (ierr.eq.1) stop ' Wrong planet index'
            if (ierr.eq.2) stop ' Error in Kepler equation'
            write (*,1001) dj,dlong,dlat,rvec,(xyz(k),k=1,6)
            if (ierr.eq.3) write (*,1002)
C            pause ' Enter a blank line to continue'
100      continue
110   continue
*
      stop1248
*
1000  format (2x,'PROGRAM PLANETAP    PLANET : ',a///
     .        2x,'Longitude (L), Latitude (B), Radius vector (R)'/
     .        2x,'Rectangular heliocentric coordinates X  Y  Z'/
     .        2x,'Rectangular heliocentric velocities  X'' Y'' Z'''/
     .        2x,'Equinox and ecliptic J2000'///)
1001  format (2x,'JD= ',f14.5//
     .        2x,'L = ',f12.7,' rad ',3x,'B = ',f12.7,' rad ',
     .        3x,'R = ',f12.7,' au'//
     .        2x,'X = ',f12.7,' au  ',3x,'Y = ',f12.7,' au  ',
     .        3x,'Z = ',f12.7,' au'//
     .        2x,'X''= ',f12.8,' au/d',3x,'Y''= ',f12.8,' au/d',
     .        3x,'Z''= ',f12.8,' au/d'//)
1002  format (2x,'You are outside of the interval 1000-3000;'/
     .        2x,'the precision decreases outside of that interval.')
*
      end
*
*
*
      subroutine PLANETAP (dj,ipla,dlong,dlat,rvec,xyz,ierr)
*-----------------------------------------------------------------------
*
*     Object :
*     Computation of approximative ephemerides of the 8 major planets.
*     - longitude, latitude, radius vector.
*     - rectangular heliocentric coordinates and velocities.
*     Reference frame : Equinox and ecliptic J2000 (JD2451545.0).
*
*     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 :
*     dj        Julian date tdb (double real).
*     ipla      planet index (integer).
*               ipla=1 Mercury                  ipla=5 Jupiter
*               ipla=2 Venus                    ipla=6 Saturn
*               ipla=3 Earth-Moon Barycenter    ipla=7 Uranus
*               ipla=4 Mars                     ipla=8 Neptune
*
*     Output parameters :
*     dlong     longitude in radians (double real).
*     dlat      latitude  in radians (double real).
*     rvec      radius vector in au  (double real).
*     xyz(6)    rectangular heliocentric coordinates (double real).
*               xyz(1) : X in au     xyz(4) : X' in au/day.
*               xyz(2) : Y in au     xyz(5) : Y' in au/day.
*               xyz(3) : Z in au     xyz(6) : Z' in au/day.
*     ierr      error index (integer)
*               ierr=0 no error.
*               ierr=1 error in planet index ipla.
*               ierr=2 error in Kepler equation.
*               ierr=3 date outside of the interval (1000-3000).
*
*     Data tables :
*
*     The tables  a, dlm, e, pi, dinc, omega  give the mean elements of
*     classical Keplerian variables of planets limited to t**2 terms.
*     - a:      semi-major axis,
*     - dlm :   mean longitude (degree and arcsecond),
*     - e :     eccentricity,
*     - pi :    longitude of the perihelion (degree and arcsecond),
*     - dinc :  inclination (degree and arcsecond),
*     - omega : longitude of the ascending node (degree and arcsecond).
*
*     The tables  kp, ca, sa  give the trigonometric terms to be added
*     to the mean elements of the semi-major axes.
*
*     The tables  kq, cl, sl  give the trigonometric terms to be added
*     to the mean elements of the mean longitudes.
*
*     Precision :
*     Over the interval 1800-2200, the precision is :
*     for the longitude of Mercury          4"
*     for the longitude of Venus            5"
*     for the longitude of E-M barycenter   6"
*     for the longitude of Mars            17"
*     for the longitude of Jupiter         71"
*     for the longitude of Saturn          81"
*     for the longitude of Uranus          86"
*     for the longitude of Neptune         11"
*     Over the interval 1000-3000, the precision is better than 1.5
*     times the precision over 1800-2200.
*     Outside of the interval 1000-3000, the precision decreases.
*
*     Subroutine called : ELLIPX.
*
*-----------------------------------------------------------------------
*
*     Declarations.
*
      implicit double precision (a-h,o-z)
      dimension xyz(6),xxp(6),amas(8)
      dimension a(3,8),dlm(3,8),e(3,8),pi(3,8),dinc(3,8),omega(3,8)
      dimension kp(9,8),ca(9,8),sa(9,8),kq(10,8),cl(10,8),sl(10,8)
      data amas/6023600.d0,408523.5d0,328900.5d0,3098710.d0,
     .          1047.355d0,  3498.5d0,  22869.d0,  19314.d0/
      data sdrad/0.4848136811095360d-5/
      data dpi/6.283185307179586d0/
      ierr=0
*
*     Tables for mean elements.
*
      data a/
     .  0.3870983098d0,            0.d0,      0.d0,
     .  0.7233298200d0,            0.d0,      0.d0,
     .  1.0000010178d0,            0.d0,      0.d0,
     .  1.5236793419d0,          3.d-10,      0.d0,
     .  5.2026032092d0,      19132.d-10,  -39.d-10,
     .  9.5549091915d0, -0.0000213896d0,  444.d-10,
     . 19.2184460618d0 ,     -3716.d-10,  979.d-10,
     . 30.1103868694d0 ,    -16635.d-10,  686.d-10/
*
      data dlm/
     . 252.25090552d0, 5381016286.88982d0,  -1.92789d0,
     . 181.97980085d0, 2106641364.33548d0,   0.59381d0,
     . 100.46645683d0, 1295977422.83429d0,  -2.04411d0,
     . 355.43299958d0,  689050774.93988d0,   0.94264d0,
     .  34.35151874d0,  109256603.77991d0, -30.60378d0,
     .  50.07744430d0,   43996098.55732d0,  75.61614d0,
     . 314.05500511d0,   15424811.93933d0,  -1.75083d0,
     . 304.34866548d0,    7865503.20744d0,   0.21103d0/
*
      data e/
     . 0.2056317526d0,  0.0002040653d0,     -28349.d-10,
     . 0.0067719164d0, -0.0004776521d0,      98127.d-10,
     . 0.0167086342d0, -0.0004203654d0, -0.0000126734d0,
     . 0.0934006477d0,  0.0009048438d0,     -80641.d-10,
     . 0.0484979255d0,  0.0016322542d0, -0.0000471366d0,
     . 0.0555481426d0, -0.0034664062d0, -0.0000643639d0,
     . 0.0463812221d0, -0.0002729293d0,  0.0000078913d0,
     . 0.0094557470d0,  0.0000603263d0,            0.d0/
*
      data pi/
     .  77.45611904d0,  5719.11590d0,   -4.83016d0,
     . 131.56370300d0,   175.48640d0, -498.48184d0,
     . 102.93734808d0, 11612.35290d0,   53.27577d0,
     . 336.06023395d0, 15980.45908d0,  -62.32800d0,
     .  14.33120687d0,  7758.75163d0,  259.95938d0,
     .  93.05723748d0, 20395.49439d0,  190.25952d0,
     . 173.00529106d0,  3215.56238d0,  -34.09288d0,
     .  48.12027554d0,  1050.71912d0,   27.39717d0/
*
      data dinc/
     . 7.00498625d0, -214.25629d0,   0.28977d0,
     . 3.39466189d0,  -30.84437d0, -11.67836d0,
     .         0.d0,  469.97289d0,  -3.35053d0,
     . 1.84972648d0, -293.31722d0,  -8.11830d0,
     . 1.30326698d0,  -71.55890d0,  11.95297d0,
     . 2.48887878d0,   91.85195d0, -17.66225d0,
     . 0.77319689d0,  -60.72723d0,   1.25759d0,
     . 1.76995259d0,    8.12333d0,   0.08135d0/
*
      data omega/
     .  48.33089304d0,  -4515.21727d0,  -31.79892d0,
     .  76.67992019d0, -10008.48154d0,  -51.32614d0,
     . 174.87317577d0,  -8679.27034d0,   15.34191d0,
     .  49.55809321d0, -10620.90088d0, -230.57416d0,
     . 100.46440702d0,   6362.03561d0,  326.52178d0,
     . 113.66550252d0,  -9240.19942d0,  -66.23743d0,
     .  74.00595701d0,   2669.15033d0,  145.93964d0,
     . 131.78405702d0,   -221.94322d0 ,  -0.78728d0/
*
*     Tables for trigonometric terms
*
      data  kp/
     . 69613, 75645, 88306, 59899, 15746, 71087, 142173,  3086,    0,
     . 21863, 32794, 26934, 10931, 26250, 43725,  53867, 28939,    0,
     . 16002, 21863, 32004, 10931, 14529, 16368,  15318, 32794,    0,
     . 6345,   7818, 15636,  7077,  8184, 14163,   1107,  4872,    0,
     . 1760,   1454,  1167,   880,   287,  2640,     19,  2047, 1454,
     .  574,      0,   880,   287,    19,  1760,   1167,   306,  574,
     .  204,      0,   177,  1265,     4,   385,    200,   208,  204,
     .    0,    102,   106,     4,    98,  1367,    487,   204,    0/
*
      data  ca/
     .       4,    -13,    11,    -9,    -9,    -3,    -1,     4,    0,
     .    -156,     59,   -42,     6,    19,   -20,   -10,   -12,    0,
     .      64,   -152,    62,    -8,    32,   -41,    19,   -11,    0,
     .     124,    621,  -145,   208,    54,   -57,    30,    15,    0,
     .  -23437,  -2634,  6601,  6259, -1507, -1821,  2620, -2115,-1489,
     .   62911,-119919, 79336, 17814,-24241, 12068,  8306, -4893, 8902,
     .  389061,-262125,-44088,  8387,-22976, -2093,  -615, -9720, 6633,
     . -412235,-157046,-31430, 37817, -9740,   -13, -7449,  9644,    0/
*
      data  sa/
     .     -29,     -1,     9,     6,    -6,     5,     4,     0,    0,
     .     -48,   -125,   -26,   -37,    18,   -13,   -20,    -2,    0,
     .    -150,    -46,    68,    54,    14,    24,   -28,    22,    0,
     .    -621,    532,  -694,   -20,   192,   -94,    71,   -73,    0,
     .  -14614, -19828, -5869,  1881, -4372, -2255,   782,   930,  913,
     .  139737,      0, 24667, 51123, -5102,  7429, -4095, -1976,-9566,
     . -138081,      0, 37205,-49039,-41901,-33872,-27037,-12474,18797,
     .       0,  28492,133236, 69654, 52322,-49577,-26430, -3593,    0/
*
      data  kq/
     . 3086,  15746, 69613, 59899, 75645, 88306,  12661,  2658,  0,   0,
     .21863,  32794, 10931,    73,  4387, 26934,   1473,  2157,  0,   0,
     .   10,  16002, 21863, 10931,  1473, 32004,   4387,    73,  0,   0,
     .   10,   6345,  7818,  1107, 15636,  7077,   8184,   532, 10,   0,
     .   19,   1760,  1454,   287,  1167,   880,    574,  2640, 19,1454,
     .   19,    574,   287,   306,  1760,    12,     31,    38, 19, 574,
     .    4,    204,   177,     8,    31,   200,   1265,   102,  4, 204,
     .    4,    102,   106,     8,    98,  1367,    487,   204,  4, 102/
*
      data  cl/
     .     21,   -95, -157,   41,   -5,   42,   23,   30,     0,    0,
     .   -160,  -313, -235,   60,  -74,  -76,  -27,   34,     0,    0,
     .   -325,  -322,  -79,  232,  -52,   97,   55,  -41,     0,    0,
     .   2268,  -979,  802,  602, -668,  -33,  345,  201,   -55,    0,
     .   7610, -4997,-7689,-5841,-2617, 1115, -748, -607,  6074,  354,
     . -18549, 30125,20012, -730,  824,   23, 1289, -352,-14767,-2062,
     .-135245,-14594, 4197,-4030,-5630,-2898, 2540, -306,  2939, 1986,
     .  89948,  2103, 8963, 2695, 3682, 1648,  866, -154, -1963, -283/
*
      data  sl/
     .   -342,   136,  -23,   62,   66,  -52,  -33,   17,     0,    0,
     .    524,  -149,  -35,  117,  151,  122,  -71,  -62,     0,    0,
     .   -105,  -137,  258,   35, -116,  -88, -112,  -80,     0,    0,
     .    854,  -205, -936, -240,  140, -341,  -97, -232,   536,    0,
     . -56980,  8016, 1012, 1448,-3024,-3710,  318,  503,  3767,  577,
     . 138606,-13478,-4964, 1441,-1319,-1482,  427, 1236, -9167,-1918,
     .  71234,-41116, 5334,-4935,-1848,   66,  434,-1748,  3780, -701,
     . -47645, 11647, 2166, 3194,  679,    0, -244, -419, -2531,   48/
*
*     Wrong planet index.
*
      if (ipla.lt.1.or.ipla.gt.8) then
         ierr=1
      else
*
*     Time parameters.
*
         t=(dj-2451545.d0)/365250.d0
         t2=t**2
*
*     Computation of the mean elements.
*
         da=a(1,ipla)+t*a(2,ipla)+t2*a(3,ipla)
         dl=(3600.d0*dlm(1,ipla)+t*dlm(2,ipla)+t2*dlm(3,ipla))*sdrad
         de=e(1,ipla)+t*e(2,ipla)+t2*e(3,ipla)
         dp=(3600.d0*pi(1,ipla)+t*pi(2,ipla)+t2*pi(3,ipla))*sdrad
         dp=mod(dp,dpi)
         if (dp.lt.0.d0) dp=dp+dpi
         di=(3600.d0*dinc(1,ipla)+t*dinc(2,ipla)+t2*dinc(3,ipla))*sdrad
         do=(3600.d0*omega(1,ipla)+t*omega(2,ipla)
     .      +t2*omega(3,ipla))*sdrad
         do=mod(do,dpi)
         if (do.lt.0.d0) do=do+dpi
*
*     Computation of the trigonometric terms.
*
         dmu=0.35953620*t
         do 100 j=1,8
            arga=kp(j,ipla)*dmu
            argl=kq(j,ipla)*dmu
            da=da+(ca(j,ipla)*cos(arga)+sa(j,ipla)*sin(arga))*1.d-7
            dl=dl+(cl(j,ipla)*cos(argl)+sl(j,ipla)*sin(argl))*1.d-7
100      continue
         arga=kp(9,ipla)*dmu
         da=da+t*(ca(9,ipla)*cos(arga)+sa(9,ipla)*sin(arga))*1.d-7
         do 200 j=9,10
            argl=kq(j,ipla)*dmu
            dl=dl+t*(cl(j,ipla)*cos(argl)+sl(j,ipla)*sin(argl))*1.d-7
200      continue
         dl=mod(dl,dpi)
         if (dl.lt.0.d0) dl=dl+dpi
*
*     Computation of the rectangular coordinates and velocities.
*
         dmas=1.d0/amas(ipla)
         call ELLIPX (da,dl,de,dp,di,do,xxp,dmas,ierr)
         if (ierr.eq.0) then
            do 300 i=1,6
               xyz(i)=xxp(i)
300         continue
*
*     Date outside of the interval (1000-3000).
*
            if (dj.gt.2817151.5d0.or.dj.lt.2086307.5d0) ierr=3
*
*     Computation of longitude, latitude, radius vector.
*
            dlong=datan2(xxp(2),xxp(1))
            rvec=sqrt(xxp(1)**2+xxp(2)**2+xxp(3)**2)
            dlat=dasin(xxp(3)/rvec)
*
         else
            ierr=2
            dlong=0.d0
            dlat=0.d0
            rvec=0.d0
            xyz(1)=0.d0
            xyz(2)=0.d0
            xyz(3)=0.d0
         endif
      endif
*
      return
      end
*
*
*
      subroutine ELLIPX (a,dlm,e,p,dia,omega,xxp,dmas,ierr)
*-----------------------------------------------------------------------
*
*     Object :
*     This subroutine computes the variables x, y, z, x',y',z'
*     from the elliptic variables.
*
*     Error index :
*     ierr=0 : no error.
*     ierr=1 : error in Kepler equation.
*
*     Intrinsic complex functions :
*     dcmplex : complex number (*16) with two real numbers (*8).
*     dimag :   imaginary part (*8) of complex number (*16).
*     dreal :   real part (*8) of complex number (*16).
*
*     Subroutine called : KEPLER.
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-y)
      complex*16 z,zto,dcmplx
      dimension xxp(6)
      ierr=0
*
      xa=a
      xl=dlm
      xk=e*cos(p)
      xh=e*sin(p)
      xq=sin(dia/2)*cos(omega)
      xp=sin(dia/2)*sin(omega)
*
      asr=648000/3.141592653589793d0
      gk=3548.1876069651d0
      xfi=sqrt(1.d0-xk**2-xh**2)
      xki=sqrt(1.d0-xq**2-xp**2)
      u=1.d0/(1.d0+xfi)
      z=dcmplx(xk,xh)
      call KEPLER (xl,z,zto,r,u,ierr)
      if (ierr.ne.0) return
      xcw=dreal(zto)
      xsw=dimag(zto)
      xm=xp*xcw-xq*xsw
      xxx=xa*r
      xxp(1)=xxx*(xcw-2*xp*xm)
      xxp(2)=xxx*(xsw+2*xq*xm)
      xxp(3)=-2*xxx*xki*xm
      xms=xa*(xh+xsw)/xfi
      xmc=xa*(xk+xcw)/xfi
      xn=gk*((1.d0+dmas)/xa**3)**(0.5d0)/asr
      xxp(4)=xn*((-1.d0+2*xp**2)*xms+2*xp*xq*xmc)
      xxp(5)=xn*((1.d0-2*xq**2)*xmc-2*xp*xq*xms)
      xxp(6)=2*xn*xki*(xp*xms+xq*xmc)
*
      return
      end
*
*
*
      subroutine KEPLER (al,z,zto,r,u,ierr)
*-----------------------------------------------------------------------
*
*     Object :
*     Kepler equation.
*
*     Error index :
*     ierr=0 : no error.
*     ierr=1 : error in algorithm.
*
*     Intrinsic complex functions :
*     dcmplex : complex number (*16) with two real numbers (*8).
*     dimag :   imaginary part (*8) of complex number (*16).
*     dreal :   real part (*8) of complex number (*16).
*     dconjg :  conjugate (*16) of complex number (*16).
*     cdabs :   absolute value (*8) of complex number (*16).
*     cdexp :   exponent (*16) of complex number (*16).
*
*-----------------------------------------------------------------------
*
      integer ierr,k,kmax
      complex*16 dcmplx,dconjg,cdexp
      double precision dreal,dimag
      double precision ex,al,e,dl,u,r
      complex*16 z1,z2,z3,zteta,z,zto
      data kmax/10/
      k=0
      ierr=0
*
      al=mod(al,2*3.141592653589793d0)
      ex=cdabs(z)
      z1=dconjg(z)
      e=al+(ex-0.125*ex**3)*sin(al)+0.5*ex**2*sin(2*al)
     .  +0.375*ex**3*sin(al*3)
*
1     continue
      k=k+1
      z2=dcmplx(0.d0,e)
      zteta=cdexp(z2)
      z3=z1*zteta
      dl=al-e+dimag(z3)
      r=1-dreal(z3)
      if (abs(dl).lt.1.d-12) then
         z1=u*z*dimag(z3)
         z2=dcmplx(dimag(z1),-dreal(z1))
         zto=(-z+zteta+z2)/r
         return
      else
         if (k.gt.kmax) then
            ierr=1
            return
         endif
      endif
      e=e+dl/r
      goto 1
*
      end
