*
*     ----------------
*     PROGRAM SERIES96   BDL-GF9612
*     ----------------
*
*     This program is an example of use of the subroutine SERIES
*     and the files of planetary series (96) in the computation of the
*     astrometric coordinates of a planet (equinox and equator J2000).
*
*     If the user wants to improve the rapidity of computation, it is
*     recommended to transform the files in direct access and to read
*     the series in memory once for all in the subroutine SERIES.
*
*     The subroutine SERIES computes the heliocentric coordinates of
*     planets (equinox and equator J2000).
*     The subroutine EARTH is necessary for geocentric coordinates.
*
*     ------------
*     Declarations
*     ------------
*
      implicit double precision (a-h,o-z)
*
      character*3 ext(9)
      character*8 pla(9)
*
      character*40 fich
*
      dimension v1(6),v2(6),w(3)
*
      data tlight/0.577551831d-2/
      data dpi/6.283185307179586d0/
*
      data ext/'mer','ven','*','mar','jup','sat','ura','nep','plu'/
      data pla/'Mercury','Venus','*','Mars','Jupiter','Saturn',
     .         'Uranus','Neptune','Pluto'/
*
*     ------------
*     Dates limits
*     ------------
*
      data t1/2415020.5d0/
      data t2/2487980.5d0/
*
*     -------------------
*     Initial julian date
*     -------------------
*
*     From : JD2415020.5d0 (1 Jan 1900 0h)
*     To   : JD2487980.5d0 (4 Oct 2099 0h)
*
      t0=2415020.5d0
*
*     ---------------
*     Number of dates
*     ---------------

      ndat=5
*
*     ------------------
*     Tab interval (day)
*     ------------------
*
      step=14500.25d0
*
*     ----------
*     Test dates
*     ----------
*
      tlim=t0+(ndat-1)*step
      if (t0.lt.t1.or.t0.gt.t2) stop '  Error : Date'
*
*     ------------
*     Planet index
*     ------------
*
*     ipla=1   Mercury
*     ipla=2   Venus
*     ipla=4   Mars
*     ipla=5   Jupiter
*     ipla=6   Saturn
*     ipla=7   Uranus
*     ipla=8   Neptune
*     ipla=9   Pluto
*
      ipla=1
*
*     -----------------
*     Test planet index
*     -----------------
*
      if (ipla.lt.1.or.ipla.eq.3.or.ipla.gt.9) stop '  Error : Planet'
*
*     ----------------
*     Open planet file
*     ----------------
*
      lupla=ipla+10
      fich='series96.'//ext(ipla)
      open (lupla,file=fich,status='old',iostat=nerr)
      if (nerr.ne.0) stop '  Error : planet File'
*
*     -------------
*     Open EMB file
*     -------------
*
      luemb=13
      fich='series96.emb'
      open (luemb,file=fich,status='old',iostat=nerr)
      if (nerr.ne.0) stop '  Error : EMB File'
*
*     ----------
*     Dates loop
*     ----------
*
      do n=1,ndat
*
         if (mod(n,10).eq.1) write (*,1001)
*
         tjd=t0+(n-1)*step
*
*        -------------------------------
*        Compute astrometric coordinates
*        -------------------------------
*
         call SERIES (tjd,3,luemb,v1,ierr)
         if (ierr.ne.0) stop '  Error : EMB series'
         call EARTH (tjd,v2)
         do i=1,3
            v1(i)=v1(i)-v2(i)
         enddo
*
         call SERIES (tjd,ipla,lupla,v2,ierr)
         if (ierr.ne.0) stop '  Error : planetary series'
*
         do k=1,2
            do i=1,3
               w(i)=v2(i)-v1(i)
            enddo
            dist=dsqrt(w(1)*w(1)+w(2)*w(2)+w(3)*w(3))
            t=tjd-dist*tlight
            call SERIES (t,ipla,lupla,v2,ierr)
            if (ierr.ne.0) stop '  Error : planetary series'
         enddo
*
         do i=1,3
            w(i)=v2(i)-v1(i)
         enddo
*
*        ---------------
*        Display results
*        ---------------
*
         alpha=datan2(w(2),w(1))
         if (alpha.lt.0.0) alpha=alpha+dpi
         p=dsqrt(w(1)*w(1)+w(2)*w(2))
         delta=datan(w(3)/p)
*
         write (*,1002) pla(ipla),tjd,alpha,delta
*
         if (mod(ndat,10).eq.1) then
            write (*,1003)
            read (*,*)
         endif
*
      enddo
*
*     -----------------
*     Close file SERIES
*     -----------------
*
      close (lupla)
      close (luemb)
*
*     ---
*     End
*     ---
*
      write (*,1004)
      read (*,*)
      stop
*
*     -------
*     Formats
*     -------
*
1001  format (////
     .        2x,'--------------------------------------------'/
     .        2x,'PLANEPH : ASTROMETRIC COORDINATES OF PLANETS'/
     .        2x,'--------------------------------------------'/)
1002  format (2x,a8,2x,'JD',f13.5,
     .        3x,'alpha: ',f12.9,' rad',
     .        3x,'delta: ',f12.9,' rad')
1003  format (//2x,'Continue : Press Enter'//)
1004  format (//2x,'Program Ended : Press Enter'//)
*
      end
*
*
*
      subroutine SERIES (tjd,ipla,lu,r,ierr)
*     ======================================
*
*
*     Ref : Bureau des Longitudes - 96.12
*           J. Chapront, G. Francou (BDL)
*
*
*     Object
*     ------
*
*     Compute planetary ephemerides with series built by a method of
*     representation using frequency analysis.
*     Ref : J. Chapront, 1995, Astron. Atrophys. Sup. Ser., 109, 181.
*
*     Planetary Ephemerides : rectangular heliocentric coordinates.
*     Source : DE403 (Jet Propulsion Laboratory).
*     Frame : Dynamical Mean Equinox and Equator J2000 (DE403).
*     Time  : Barycentric Dynamical Time (TDB).
*
*     Input
*     -----
*
*     tjd :       Julian date TDB (double real).
*
*     ipla :      Planet index (integer).
*                 1 : Mercury.
*                 2 : Venus.
*                 3 : Earth-Moon Barycenter.
*                 4 : Mars.
*                 5 : Jupiter.
*                 6 : Saturn.
*                 7 : Uranus.
*                 8 : Neptune.
*                 9 : Pluto.
*
*     lu :        Logical unit of planet file (integer).
*
*     Output
*     ------
*
*     r(6) :      Table of rectangular coordinates (double real).
*                 r(1) : X  position (au).
*                 r(2) : Y  position (au).
*                 r(3) : Z  position (au).
*                 r(4) : X' velocity (au/day).
*                 r(5) : Y' velovity (au/day).
*                 r(6) : Z' velocity (au/day).
*
*     ierr :      Error index (integer).
*                 0 : no error.
*                 1 : date error.
*                 2 : planet index error.
*
*     Files
*     -----
*
*     According to the planet the series are in a file named SERIES.xxx
*     where xxx is : mer, ven, emb, mar, jup, sat, ura, nep, plu.
*     The file is opened in the monitor (logical unit lu).
*
*
*     Declarations
*     ------------
*
      implicit double precision (a-h,o-z)
      character*3 nom(9)
      character*12 cname
*
      parameter (maxor=2,maxfq=250)
*
      dimension r(6)
      dimension v(3),w(3),ws(3)
      dimension nf(0:maxor),fq(maxfq,0:maxor)
      dimension sec(0:3,3)
      dimension ct(maxfq,0:maxor,3)
      dimension st(maxfq,0:maxor,3)
*
      data nom/'mer','ven','emb','mar','jup','sat','ura','nep','plu'/
*
*
*     Clear results table
*     -------------------
*
      do i=1,6
         r(i)=0.d0
      enddo
*
*
*     Read file parameters
*     ---------------------
*
      rewind lu
      read (lu,fmt='(a12,f10.1,f8.0,3i3)') cname,tzero,dt,mx,imax,iblock
*
*
*     Check up date
*     -------------
*
      ierr=1
      tmax=tzero+dt*iblock
      if (tjd.lt.tzero-0.5d0.or.tjd.gt.tmax+0.5d0) return
*
      nb=(tjd-tzero)/dt+1
      if (tjd.le.tzero) nb=1
      if (tjd.ge.tmax) nb=iblock
*
*
*     Check up planet
*     ---------------
*
      ierr=2
      if (ipla.lt.1.or.ipla.gt.9) return
      if (cname(10:12).ne.nom(ipla)) return
*
*
*     Read frequencies
*     ----------------
*
      do m=0,mx
         read (lu,fmt='(i4)') nf(m)
         do i=1,nf(m)
            read (lu,fmt='(f20.16)') fq(i,m)
         enddo
      enddo
*
*
*     Read separator
*     --------------
*
      read (lu,*)
*
*
*     Read series
*     -----------
*
      do k=1,nb
      do iv=1,3
*
*
*        Read block and variable numbers
*        -------------------------------
*
         read (lu,*)
*
*
*        Read secular terms
*        ------------------
*
         do i=0,imax,2
            read (lu,fmt='(2f14.0)') sec(i,iv),sec(i+1,iv)
         enddo
*
*
*        Read Poisson terms
*        ------------------
*
         do m=0,mx
            ip=mod(m,2)
            do i=1,nf(m)
               read (lu,fmt='(2f14.0)') a,b
               if (ip.eq.0) then
                  ct(i,m,iv)=a
                  st(i,m,iv)=b
               else
                  ct(i,m,iv)=b
                  st(i,m,iv)=a
               endif
            enddo
         enddo
*
      enddo
      enddo
*
*
*     Change variable
*     ---------------
*
      tdeb=tzero+(nb-1)*dt
      x=2.d0*(tjd-tdeb)/dt-1.d0
      fx=x*dt/2.d0
*
*     Compute positions (secular terms)
*     ---------------------------------
*
      do iv=1,3
         v(iv)=0.d0
         wx=1.d0
         max=2*imax-1
         do i=0,max
            v(iv)=v(iv)+sec(i,iv)*wx
            wx=wx*x
         enddo
      enddo
*
*     Compute positions (Poisson terms)
*     ---------------------------------
*
      wx=1.d0
      do m=0,mx
         nw=nf(m)
         do iv=1,3
            ws(iv)=0.d0
         enddo
         do i=1,nw
            f=fq(i,m)*fx
            cf=cos(f)
            sf=sin(f)
            do iv=1,3
               ws(iv)=ws(iv)+ct(i,m,iv)*cf+st(i,m,iv)*sf
            enddo
         enddo
         do iv=1,3
            v(iv)=v(iv)+ws(iv)*wx
         enddo
         wx=wx*x
      enddo
*
*
*     Compute velocities (secular terms)
*     ----------------------------------
*
      wt=2.d0/dt
      do iv=1,3
         w(iv)=0.d0
         wx=1.d0
         max=2*imax-1
         do i=1,max
            w(iv)=w(iv)+i*sec(i,iv)*wx
            wx=wx*x
         enddo
         w(iv)=wt*w(iv)
      enddo
*
*
*     Compute velocities (Poisson terms)
*     ----------------------------------
*
      wx=1.d0
      do m=0,mx
         nw=nf(m)
         do i=1,nw
            fw=fq(i,m)
            f=fw*fx
            cf=cos(f)
            sf=sin(f)
            do iv=1,3
               stw=st(i,m,iv)
               ctw=ct(i,m,iv)
               w(iv)=w(iv)+fw*(stw*cf-ctw*sf)*wx
               if (m.gt.0) w(iv)=w(iv)+m*wt*(ctw*cf+stw*sf)*wy
            enddo
         enddo
         wy=wx
         wx=wx*x
      enddo
*
*     Stock results
*     -------------
*
      do iv=1,3
         r(iv)=v(iv)/1.d10
         r(iv+3)=w(iv)/1.d10
      enddo
*
*     Exit
*     ----
*
      ierr=0
      return
      end
*
*
*
      subroutine EARTH (tjd,r)
*     ========================
*
*
*     Ref : Bureau des Longitudes - 96.12
*           J. Chapront, G. Francou (BDL)
*
*
*     Object
*     ------
*
*     Rectangular coordinates of geocentric Earth-Moon barycenter
*     (equinox and equateur J2000).
*
*     Input
*     -----
*
*     tjd :       Julian date TDB (double real).
*
*
*     Output
*     ------
*
*     r(3) :      Table of rectangular coordinates (double real).
*                 r(1) : X  position (au).
*                 r(2) : Y  position (au).
*                 r(3) : Z  position (au).
*
*-----------------------------------------------------------------------
*
      implicit double precision (a-h,o-z)
*
      dimension r(3),v(3),n1(3),n2(3)
      dimension c(138),s(138),f(138)
*
      data n1/01,44,93/
      data n2/43,92,138/
*
      data c/
     .-244075.,  -2965.,   8528.,   2345.,  -2486.,   1426.,    527.,
     .    -43.,   -393.,    394.,   -218.,     73.,     91.,   -173.,
     .     25.,    -20.,     75.,     72.,      6.,     72.,    -40.,
     .    -58.,     56.,    -53.,     46.,    -44.,     -5.,      0.,
     .     -1.,    -12.,     21.,     -4.,     -4.,     -2.,      9.,
     .      8.,     10.,      2.,    -10.,    -12.,    -11.,     10.,
     .    -10.,-176962., -23344., -11109.,   -922.,  -4118.,    714.,
     .  -1135.,   -601.,    299.,    564.,   -311.,    261.,   -251.,
     .    254.,    229.,    213.,   -179.,     57.,    -19.,   -125.,
     .   -113.,    -87.,    -52.,     75.,     16.,     -5.,    -42.,
     .     -4.,    -10.,      8.,    -19.,      9.,     40.,     29.,
     .    -25.,     11.,     19.,    -12.,     -5.,     18.,    -15.,
     .    -16.,     13.,     13.,      9.,     -1.,     -3.,     -5.,
     .     -1., -76714.,  25611., -10120.,   -400.,   1387.,  -1785.,
     .    310.,    580.,   -492.,   -527.,     44.,    130.,    244.,
     .   -135.,    113.,    -38.,    110.,     92.,    -78.,     25.,
     .    -54.,    -26.,    -49.,    -38.,     27.,    -23.,     32.,
     .      2.,     -2.,      1.,    -18.,     -2.,    -23.,     -4.,
     .      3.,     -8.,      4.,    -11.,     17.,      3.,    -11.,
     .     13.,      8.,    -11.,     -7.,      4./
*
      data s/
     . 192874.,  25444.,   1005.,   4489.,   -778.,   1238.,   -326.,
     .   -614.,    339.,   -285.,   -276.,   -232.,    195.,    -63.,
     .    136.,    124.,     95.,     57.,    -82.,      5.,     45.,
     .      4.,     11.,     -9.,     20.,    -10.,    -44.,    -32.,
     .     27.,    -20.,      6.,    -20.,     17.,     17.,    -14.,
     .    -14.,      4.,    -11.,    -10.,      3.,      6.,     -7.,
     .      4.,-223938.,  -2720.,    635.,   7824.,   2151.,  -2281.,
     .   1309.,   -675.,    483.,    -40.,   -360.,    362.,    327.,
     .   -200.,    204.,     67.,     84.,   -159.,   -145.,     23.,
     .    -18.,     69.,     66.,      5.,     65.,     66.,    -36.,
     .    -53.,     51.,    -48.,     42.,    -40.,     -5.,      0.,
     .     -1.,    -19.,    -11.,     17.,     20.,     -4.,     -4.,
     .     -2.,      9.,      7.,     -9.,    -13.,    -11.,    -10.,
     .     11., -97079.,  -1464.,  -1179.,   3392.,   1557.,    933.,
     .   -989.,   -754.,    567.,   -470.,    334.,    210.,    -17.,
     .   -156.,    157.,   -151.,    -87.,     29.,     36.,    -69.,
     .     10.,     43.,     -8.,     30.,    -39.,     29.,      2.,
     .     29.,     29.,    -25.,    -16.,    -23.,      3.,     22.,
     .    -21.,     18.,    -17.,     13.,     -2.,     14.,     -8.,
     .      0.,     -9.,      0.,     -8.,      9./
*
      data f/
     . 0.2299708345453799, 0.0019436907548255, 0.4579979783362081,
     . 0.0324605575663244,-0.1955665862245038, 0.4274811115247091,
     .-0.2318206046403833, 0.6555082553155372, 0.2127688645204655,
     . 0.2471728045705681, 0.6860251221267625, 0.2604877013568789,
     . 0.0496625275912389,-0.1783646161993155, 0.0172021241604381,
     . 0.0191456607800137,-0.2260834530360027, 0.4102791414995209,
     .-0.0152582792703628, 0.4407960083110198, 0.8835353991060917,
     . 0.1994539677341547, 0.2662248529612594, 0.4751999483613963,
     . 0.4427395449305955,-0.0037934608495551, 0.6383062852903491,
     . 0.4637351299405887, 0.1937168161297741, 0.0152585875411362,
     .-0.2127685562496920, 0.0000001541352498,-0.4598477484309377,
     . 0.9140522659175908, 0.2452292679512662, 0.6249913885040383,
     . 0.2300338378809035,-0.2299078312098563, 0.4446830815498973,
     . 0.8530185322948665, 0.4885148451477070, 0.0381977091707050,
     . 0.2147124011397673, 0.2299708345453799, 0.0019436907548255,
     . 0.2308957195928816, 0.4579979783362081, 0.0324605575663244,
     .-0.1955665862245038, 0.4274811115247091,-0.0028685758023272,
     .-0.2318206046403833, 0.6555082553155372, 0.2127688645204655,
     . 0.2471728045705681, 0.1946417011770021, 0.6860251221267625,
     . 0.4589228633837098, 0.2604877013568789, 0.0496625275912389,
     .-0.1783646161993155,-0.0333854426135524, 0.0172021241604381,
     . 0.0191456607800137,-0.2260834530360027, 0.4102791414995209,
     .-0.0152582792703628, 0.4284059965722108, 0.4407960083110198,
     . 0.8835353991060917, 0.1994539677341547, 0.2662248529612594,
     . 0.4751999483613963, 0.4427395449305955,-0.0037934608495551,
     . 0.6383062852903491, 0.4637351299405887, 0.1937168161297741,
     . 0.6564331403627652, 0.0152585875411362, 0.1774397311520876,
     .-0.2127685562496920, 0.0000001541352498,-0.4598477484309377,
     . 0.9140522659175908, 0.2452292679512662, 0.6249913885040383,
     . 0.4446830815498973, 0.6869500071742641, 0.8530185322948665,
     . 0.4885148451477070, 0.2251585679885010, 0.2299708345453799,
     . 0.2308957195928816, 0.0019436907548255, 0.4579979783362081,
     .-0.0028685758023272, 0.0324605575663244,-0.1955665862245038,
     . 0.1946417011770021, 0.4274811115247091, 0.4589228633837098,
     .-0.0333854426135524,-0.2318206046403833, 0.6555082553155372,
     . 0.2127688645204655, 0.2471728045705681, 0.4284059965722108,
     . 0.6860251221267625, 0.2604877013568789, 0.0496625275912389,
     .-0.1783646161993155, 0.0172021241604381, 0.6564331403627652,
     . 0.0191456607800137,-0.2260834530360027, 0.1774397311520876,
     . 0.4102791414995209,-0.0152582792703628, 0.6869500071742641,
     . 0.4407960083110198, 0.2251585679885010, 0.8835353991060917,
     . 0.1994539677341547, 0.4226688449678302, 0.2662248529612594,
     . 0.4751999483613963, 0.4427395449305955,-0.0037934608495551,
     . 0.2118436712021903, 0.6383062852903491,-0.0505874126387406,
     .-0.2614125864043805, 0.4637351299405887, 0.0200705458272416,
     . 0.1937168161297741, 0.0143333942228611,-0.0181270092079398/
*
      t=tjd-2451545.0d0
*
      do iv=1,3
         v(iv)=0.d0
         do n=n1(iv),n2(iv)
            x=f(n)*t
            cx=cos(x)
            sx=sin(x)
            v(iv)=v(iv)+c(n)*cx+s(n)*sx
         enddo
         v(iv)=v(iv)/1.d10
         r(iv)=v(iv)
      enddo
*
      return
      end
