c++++++++++++++++++++++++
c
c      subroutine fsizer1(nrecl,ksize,nrfile,namfil)
c
c++++++++++++++++++++++++
c
c   version 1.0 uses the inquire statement to find out the the record length 
c   of the direct access file before opening it.  this procedure is non-standard, 
c   but seems to work for vax machines. 
c
c  the subroutine also sets the values of  nrecl, nrfile, and namfil.
c
c  *****************************************************************
c  *****************************************************************
c
c  the parameters namfil, nrecl, and nrfile are to be set by the user
c
c  *****************************************************************
c
c  namfil is the external name of the binary ephemeris file
c
c      character*80  namfil
c
c      namfil=
c
c  *****************************************************************
c
c  nrecl=1 if "recl" in the open statement is the record length in s.p. words
c  nrecl=4 if "recl" in the open statement is the record length in bytes
c  (for a vax, it is probably 1)
c
c      nrecl= 
c
c  *****************************************************************
c
c  nrfile is the internal unit number used for the ephemeris file
c
c      nrfile=12
c
c  *****************************************************************
c
c   find the record size using the inquire statement  
c
c
c     irecsz=0
c
c      inquire(file=namfil,recl=irecsz)
c
c if 'inquire' does not work, usually irecsz will be left at 0 
c
c      if(irecsz .le. 0) write(*,*)
c     . ' inquire statement probably did not work'
c
c      ksize=irecsz/nrecl
c
c      return
c
c      end






      subroutine fsizer2(nrecl,ksize,nrfile,namfil)
      implicit double precision(a-h,o-z)
      save

c  this subroutine opens the file, 'namfil', with a phony record length, reads 
c  the first record, and uses the info to compute ksize, the number of single 
c  precision words in a record.  
c
c  the subroutine also sets the values of  nrecl, nrfile, and namfil.

      character*6 ttl(14,3),cnam(400)
      character*80 namfil,word,string
      logical      ibhere

      dimension ss(3)

      integer ipt(3,13)

c..for picking which ephemeris to use
      integer          which_eph
      common /dop/     which_eph



c  *****************************************************************
c  *****************************************************************
c
c  the parameters nrecl, nrfile, and namfil are to be set by the user
c
c  *****************************************************************




c  nrecl=1 if "recl" in the open statement is the record length in s.p. words
c  nrecl=4 if "recl" in the open statement is the record length in bytes

c..use nrecl=1 for sgi boxes f77
c..use nrecl=4 for mac os x f77 or f90  or sgi f90

      nrecl=1
c      nrecl=4

c  nrfile is the internal unit number used for the ephemeris file

      nrfile=12


c  namfil is the external name of the binary ephemeris file
c..pick the ephemeis to use from the common block integer which_eph
      if (which_eph .eq. 405) then
       namfil = 'unxp1600-2200.405'
       word   = 'JPL405'

      else if (which_eph .eq. 406) then
       namfil = 'unxm3000-p3000.406'
       word   = 'JPL406'

c..open something
      else
       which_eph = 405
       namfil    = 'unxp1600-2200.405'
       word      = 'JPL405'
      end if


c..see if the file is in the local directory
      inquire(file=namfil, exist=ibhere)
      if (.not.ibhere) stop 'could not find JPLDATA file'

c  **  open the direct-access file and get the pointers in order to 
c  **  determine the size of the ephemeris record

      mrecl=nrecl*1000

        open(nrfile,
     *       file=namfil,
     *       access='direct',
     *       form='unformatted',
     *       recl=mrecl,
     *       status='old')

      read(nrfile,rec=1)ttl,cnam,ss,ncon,au,emrat,
     * ((ipt(i,j),i=1,3),j=1,12),numde,(ipt(i,13),i=1,3)

      close(nrfile)


c  find the number of ephemeris coefficients from the pointers

      kmx = 0
      khi = 0

      do i = 1,13
         if (ipt(1,i) .gt. kmx) then
            kmx = ipt(1,i)
            khi = i
         endif
      enddo

      nd = 3
      if (khi .eq. 12) nd=2

      ksize = 2*(ipt(1,khi)+nd*ipt(2,khi)*ipt(3,khi)-1)

      return
      end







      subroutine fsizer3(nrecl,ksize,nrfile,namfil)
c
c++++++++++++++++++++++++
c
c  the subroutine sets the values of  nrecl, ksize, nrfile, and namfil.
c
      save

      character*80 namfil

c  *****************************************************************
c  *****************************************************************
c
c  the parameters nrecl, nrfile, and namfil are to be set by the user
c
c  *****************************************************************
c
c  nrecl=1 if "recl" in the open statement is the record length in s.p. words
c  nrecl=4 if "recl" in the open statement is the record length in bytes
c
c      nrecl=1
      nrecl=4

c  *****************************************************************
c
c  nrfile is the internal unit number used for the ephemeris file (default: 12)

      nrfile=12
c
c  *****************************************************************
c
c  namfil is the external name of the binary ephemeris file
      namfil='unxp2000.405'


c  *****************************************************************
c
c  ksize must be set by the user according to the ephemeris to be read
c
c  for  de200, set ksize to 1652
c  for  de405, set ksize to 2036
c  for  de406, set ksize to 1456
c
      ksize = 2036

      return
      end






c++++++++++++++++++++++++++
c
      subroutine pleph ( et, ntarg, ncent, rrd )
c
c++++++++++++++++++++++++++
c  note : over the years, different versions of pleph have had a fifth argument:
c  sometimes, an error return statement number; sometimes, a logical denoting
c  whether or not the requested date is covered by the ephemeris.  we apologize
c  for this inconsistency; in this present version, we use only the four necessary 
c  arguments and do the testing outside of the subroutine.
c
c
c
c     this subroutine reads the jpl planetary ephemeris
c     and gives the position and velocity of the point 'ntarg'
c     with respect to 'ncent'.
c
c     calling sequence parameters:
c
c       et = d.p. julian ephemeris date at which interpolation
c            is wanted.
c
c       ** note the entry dpleph for a doubly-dimensioned time **
c          the reason for this option is discussed in the 
c          subroutine state
c
c     ntarg = integer number of 'target' point.
c
c     ncent = integer number of center point.
c
c            the numbering convention for 'ntarg' and 'ncent' is:
c
c                1 = mercury           8 = neptune
c                2 = venus             9 = pluto
c                3 = earth            10 = moon
c                4 = mars             11 = sun
c                5 = jupiter          12 = solar-system barycenter
c                6 = saturn           13 = earth-moon barycenter
c                7 = uranus           14 = nutations (longitude and obliq)
c                            15 = librations, if on eph file
c
c             (if nutations are wanted, set ntarg = 14. for librations,
c              set ntarg = 15. set ncent=0.)
c
c      rrd = output 6-word d.p. array containing position and velocity
c            of point 'ntarg' relative to 'ncent'. the units are au and
c            au/day. for librations the units are radians and radians
c            per day. in the case of nutations the first four words of
c            rrd will be set to nutations and rates, having units of
c            radians and radians/day.
c
c            the option is available to have the units in km and km/sec.
c            for this, set km=.true. in the stcomx common block.
c

      implicit double precision (a-h,o-z)

      dimension rrd(6),et2z(2),et2(2),pv(6,13),pnut(4)
      dimension ss(3),cval(400),pvsun(6)

      logical bsave,km,bary
      logical first
      data first/.true./

      integer list(12),ipt(39),denum


      common/ephhdr/cval,ss,au,emrat,denum,ncon,ipt

      common/stcomx/km,bary,pvsun

c     initialize et2 for 'state' and set up component count
c
      et2(1)=et
      et2(2)=0.d0
      go to 11

c     entry point 'dpleph' for doubly-dimensioned time argument 
c          (see the discussion in the subroutine state)

      entry dpleph(et2z,ntarg,ncent,rrd)

      et2(1)=et2z(1)
      et2(2)=et2z(2)

  11  do i=1,6
      rrd(i)=0.d0
      enddo

c      if(first) call state(0.d0,0,0,0)
      if(first) call state(et2,list,pv,rrd)
      first=.false.

  96  if(ntarg .eq. ncent) return

      do i=1,12
      list(i)=0
      enddo

c     check for nutation call

      if(ntarg.ne.14) go to 97
        if(ipt(35).gt.0) then
          list(11)=2
          call state(et2,list,pv,rrd)
          return
        else
          write(6,297)
  297     format(' *****  no nutations on the ephemeris file  *****')
          stop
        endif

c     check for librations

  97  if(ntarg.ne.15) go to 98
        if(ipt(38).gt.0) then
          list(12)=2
          call state(et2,list,pv,rrd)
          do i=1,6
          rrd(i)=pv(i,11)
          enddo
          return
        else
          write(6,298)
  298     format(' *****  no librations on the ephemeris file  *****')
          stop
        endif

c       force barycentric output by 'state'

  98  bsave=bary
      bary=.true.

c       set up proper entries in 'list' array for state call

      do i=1,2
      k=ntarg
      if(i .eq. 2) k=ncent
      if(k .le. 10) list(k)=2
      if(k .eq. 10) list(3)=2
      if(k .eq. 3) list(10)=2
      if(k .eq. 13) list(3)=2
      enddo

c       make call to state

      call state(et2,list,pv,rrd)

      if(ntarg .eq. 11 .or. ncent .eq. 11) then
      do i=1,6
      pv(i,11)=pvsun(i)
      enddo
      endif

      if(ntarg .eq. 12 .or. ncent .eq. 12) then
      do i=1,6
      pv(i,12)=0.d0
      enddo
      endif

      if(ntarg .eq. 13 .or. ncent .eq. 13) then
      do i=1,6
      pv(i,13)=pv(i,3)
      enddo
      endif

      if(ntarg*ncent .eq. 30 .and. ntarg+ncent .eq. 13) then
      do i=1,6
      pv(i,3)=0.d0
      enddo
      go to 99
      endif

      if(list(3) .eq. 2) then
      do i=1,6
      pv(i,3)=pv(i,3)-pv(i,10)/(1.d0+emrat)
      enddo
      endif

      if(list(10) .eq. 2) then
      do i=1,6
      pv(i,10)=pv(i,3)+pv(i,10)
      enddo
      endif

  99  do i=1,6
      rrd(i)=pv(i,ntarg)-pv(i,ncent)
      enddo

      bary=bsave

      return
      end
c+++++++++++++++++++++++++++++++++
c
      subroutine interp(buf,t,ncf,ncm,na,ifl,pv)
c
c+++++++++++++++++++++++++++++++++
c
c     this subroutine differentiates and interpolates a
c     set of chebyshev coefficients to give position and velocity
c
c     calling sequence parameters:
c
c       input:
c
c         buf   1st location of array of d.p. chebyshev coefficients of position
c
c           t   t(1) is dp fractional time in interval covered by
c               coefficients at which interpolation is wanted
c               (0 .le. t(1) .le. 1).  t(2) is dp length of whole
c               interval in input time units.
c
c         ncf   # of coefficients per component
c
c         ncm   # of components per set of coefficients
c
c          na   # of sets of coefficients in full array
c               (i.e., # of sub-intervals in full interval)
c
c          ifl  integer flag: =1 for positions only
c                             =2 for pos and vel
c
c
c       output:
c
c         pv   interpolated quantities requested.  dimension
c               expected is pv(ncm,ifl), dp.
c
c
      implicit double precision (a-h,o-z)
c
      save
c
      double precision buf(ncf,ncm,*),t(2),pv(ncm,*),pc(18),vc(18)

c
      data np/2/
      data nv/3/
      data twot/0.d0/
      data pc(1),pc(2)/1.d0,0.d0/
      data vc(2)/1.d0/
c
c       entry point. get correct sub-interval number for this set
c       of coefficients and then get normalized chebyshev time
c       within that subinterval.
c
      dna=dble(na)
      dt1=dint(t(1))
      temp=dna*t(1)
      l=idint(temp-dt1)+1

c         tc is the normalized chebyshev time (-1 .le. tc .le. 1)

      tc=2.d0*(dmod(temp,1.d0)+dt1)-1.d0

c       check to see whether chebyshev time has changed,
c       and compute new polynomial values if it has.
c       (the element pc(2) is the value of t1(tc) and hence
c       contains the value of tc on the previous call.)

      if(tc.ne.pc(2)) then
        np=2
        nv=3
        pc(2)=tc
        twot=tc+tc
      endif
c
c       be sure that at least 'ncf' polynomials have been evaluated
c       and are stored in the array 'pc'.
c
      if(np.lt.ncf) then
        do 1 i=np+1,ncf
        pc(i)=twot*pc(i-1)-pc(i-2)
    1   continue
        np=ncf
      endif
c
c       interpolate to get position for each component
c
      do 2 i=1,ncm
      pv(i,1)=0.d0
      do 3 j=ncf,1,-1
      pv(i,1)=pv(i,1)+pc(j)*buf(j,i,l)
    3 continue
    2 continue
      if(ifl.le.1) return
c
c       if velocity interpolation is wanted, be sure enough
c       derivative polynomials have been generated and stored.
c
      vfac=(dna+dna)/t(2)
      vc(3)=twot+twot
      if(nv.lt.ncf) then
        do 4 i=nv+1,ncf
        vc(i)=twot*vc(i-1)+pc(i-1)+pc(i-1)-vc(i-2)
    4   continue
        nv=ncf
      endif
c
c       interpolate to get velocity for each component
c
      do 5 i=1,ncm
      pv(i,2)=0.d0
      do 6 j=ncf,2,-1
      pv(i,2)=pv(i,2)+vc(j)*buf(j,i,l)
    6 continue
      pv(i,2)=pv(i,2)*vfac
    5 continue
c
      return
c
      end

c+++++++++++++++++++++++++
c
      subroutine split(tt,fr)
c
c+++++++++++++++++++++++++
c
c     this subroutine breaks a d.p. number into a d.p. integer
c     and a d.p. fractional part.
c
c     calling sequence parameters:
c
c       tt = d.p. input number
c
c       fr = d.p. 2-word output array.
c            fr(1) contains integer part
c            fr(2) contains fractional part
c
c            for negative input numbers, fr(1) contains the next
c            more negative integer; fr(2) contains a positive fraction.
c
c       calling sequence declarations
c
      implicit double precision (a-h,o-z)

      dimension fr(2)

c       main entry -- get integer and fractional parts

      fr(1)=dint(tt)
      fr(2)=tt-fr(1)

      if(tt.ge.0.d0 .or. fr(2).eq.0.d0) return

c       make adjustments for negative input number

      fr(1)=fr(1)-1.d0
      fr(2)=fr(2)+1.d0

      return

      end





      subroutine state(et2,list,pv,pnut)
c
c++++++++++++++++++++++++++++++++
c
c this subroutine reads and interpolates the jpl planetary ephemeris file
c
c     calling sequence parameters:
c
c     input:
c
c         et2   dp 2-word julian ephemeris epoch at which interpolation
c               is wanted.  any combination of et2(1)+et2(2) which falls
c               within the time span on the file is a permissible epoch.
c
c                a. for ease in programming, the user may put the
c                   entire epoch in et2(1) and set et2(2)=0.
c
c                b. for maximum interpolation accuracy, set et2(1) =
c                   the most recent midnight at or before interpolation
c                   epoch and set et2(2) = fractional part of a day
c                   elapsed between et2(1) and epoch.
c
c                c. as an alternative, it may prove convenient to set
c                   et2(1) = some fixed epoch, such as start of integration,
c                   and et2(2) = elapsed interval between then and epoch.
c
c        list   12-word integer array specifying what interpolation
c               is wanted for each of the bodies on the file.
c
c                         list(i)=0, no interpolation for body i
c                                =1, position only
c                                =2, position and velocity
c
c               the designation of the astronomical bodies by i is:
c
c                         i = 1: mercury
c                           = 2: venus
c                           = 3: earth-moon barycenter
c                           = 4: mars
c                           = 5: jupiter
c                           = 6: saturn
c                           = 7: uranus
c                           = 8: neptune
c                           = 9: pluto
c                           =10: geocentric moon
c                           =11: nutations in longitude and obliquity
c                           =12: lunar librations (if on file)
c
c
c     output:
c
c          pv   dp 6 x 11 array that will contain requested interpolated
c               quantities.  the body specified by list(i) will have its
c               state in the array starting at pv(1,i).  (on any given
c               call, only those words in 'pv' which are affected by the
c               first 10 'list' entries (and by list(12) if librations are
c               on the file) are set.  the rest of the 'pv' array
c               is untouched.)  the order of components starting in
c               pv(1,i) is: x,y,z,dx,dy,dz.
c
c               all output vectors are referenced to the earth mean
c               equator and equinox of j2000 if the de number is 200 or
c               greater; of b1950 if the de number is less than 200. 
c
c               the moon state is always geocentric; the other nine states 
c               are either heliocentric or solar-system barycentric, 
c               depending on the setting of common flags (see below).
c
c               lunar librations, if on file, are put into pv(k,11) if
c               list(12) is 1 or 2.
c
c         nut   dp 4-word array that will contain nutations and rates,
c               depending on the setting of list(11).  the order of
c               quantities in nut is:
c
c                        d psi  (nutation in longitude)
c                        d epsilon (nutation in obliquity)
c                        d psi dot
c                        d epsilon dot
c
c           *   statement # for error return, in case of epoch out of
c               range or i/o errors.
c
c
c     common area stcomx:
c
c          km   logical flag defining physical units of the output
c               states. km = .true., km and km/sec
c                          = .false., au and au/day
c               default value = .false.  (km determines time unit
c               for nutations and librations.  angle unit is always radians.)
c
c        bary   logical flag defining output center.
c               only the 9 planets are affected.
c                        bary = .true. =\ center is solar-system barycenter
c                             = .false. =\ center is sun
c               default value = .false.
c
c       pvsun   dp 6-word array containing the barycentric position and
c               velocity of the sun.
c
c
      implicit double precision (a-h,o-z)

      save

      dimension et2(2),pv(6,12),pnut(4),t(2),pjd(4),buf(1500),
     . ss(3),cval(400)


cxxx is the original code, below it is the repaired code
cxxx to stop the array pvsun from wandering out of the declared bounds.

cxxx      double precision pvsun(3,2)
      double precision pvsun(6),pvsun_eq(3,2)



      integer list(12),ipt(3,13)

      logical first
      data first/.true./

      character*6 ttl(14,3),cnam(400)
      character*80 namfil

      logical km,bary

      common/ephhdr/cval,ss,au,emrat,numde,ncon,ipt
      common/chrhdr/cnam,ttl
      common/stcomx/km,bary,pvsun


cxxx this is my trick for having both a pvsun(6) and a pvsun(2,3)

      equivalence(pvsun,pvsun_eq)


c
c       entry point - 1st time in, get pointer data, etc., from eph file
c
      if(first) then
        first=.false.

c ************************************************************************
c ************************************************************************

c the user must select one of the following by deleting the 'c' in column 1

c ************************************************************************

c        call fsizer1(nrecl,ksize,nrfile,namfil)
        call fsizer2(nrecl,ksize,nrfile,namfil)
c        call fsizer3(nrecl,ksize,nrfile,namfil)

      if(nrecl .eq. 0) write(*,*)'  ***** fsizer is not working *****'

c ************************************************************************
c ************************************************************************

      irecsz=nrecl*ksize


      ncoeffs=ksize/2

        open(nrfile,
     *       file=namfil,
     *       access='direct',
     *       form='unformatted',
     *       recl=irecsz,
     *       status='old')



      read(nrfile,rec=1)ttl,cnam,ss,ncon,au,emrat,
     . ((ipt(i,j),i=1,3),j=1,12),numde,(ipt(i,13),i=1,3)

      read(nrfile,rec=2)cval

      nrl=0

      endif


c       ********** main entry point **********


      if(et2(1) .eq. 0.d0) return

      s=et2(1)-.5d0
      call split(s,pjd(1))
      call split(et2(2),pjd(3))
      pjd(1)=pjd(1)+pjd(3)+.5d0
      pjd(2)=pjd(2)+pjd(4)
      call split(pjd(2),pjd(3))
      pjd(1)=pjd(1)+pjd(3)

c       error return for epoch out of range

      if(pjd(1)+pjd(4).lt.ss(1) .or. pjd(1)+pjd(4).gt.ss(2)) go to 98

c       calculate record # and relative time in interval

      nr=idint((pjd(1)-ss(1))/ss(3))+3
      if(pjd(1).eq.ss(2)) nr=nr-1
      t(1)=((pjd(1)-(dble(nr-3)*ss(3)+ss(1)))+pjd(4))/ss(3)

c       read correct record if not in core

      if(nr.ne.nrl) then
        nrl=nr
        read(nrfile,rec=nr,err=99)(buf(k),k=1,ncoeffs)
      endif

      if(km) then
      t(2)=ss(3)*86400.d0
      aufac=1.d0
      else
      t(2)=ss(3)
      aufac=1.d0/au
      endif

c   interpolate ssbary sun


cxxx      call interp(buf(ipt(1,11)),t,ipt(2,11),3,ipt(3,11),2,pvsun)
      call interp(buf(ipt(1,11)),t,ipt(2,11),3,ipt(3,11),2,pvsun_eq)

cxxx      do i=1,6
cxxx      pvsun(i,1)=pvsun(i,1)*aufac
cxxx      enddo

      do i=1,6
       pvsun(i)=pvsun(i)*aufac
      enddo

c   check and interpolate whichever bodies are requested

      do 4 i=1,10
      if(list(i).eq.0) go to 4

      call interp(buf(ipt(1,i)),t,ipt(2,i),3,ipt(3,i),
     & list(i),pv(1,i))

      do j=1,6
       if(i.le.9 .and. .not.bary) then
cxxx       pv(j,i)=pv(j,i)*aufac-pvsun(j,1)
       pv(j,i)=pv(j,i)*aufac-pvsun(j)
       else
       pv(j,i)=pv(j,i)*aufac
       endif
      enddo

   4  continue

c       do nutations if requested (and if on file)

      if(list(11).gt.0 .and. ipt(2,12).gt.0)
     * call interp(buf(ipt(1,12)),t,ipt(2,12),2,ipt(3,12),
     * list(11),pnut)

c       get librations if requested (and if on file)

      if(list(12).gt.0 .and. ipt(2,13).gt.0)
     * call interp(buf(ipt(1,13)),t,ipt(2,13),3,ipt(3,13),
     * list(12),pv(1,11))

      return

  98  write(*,198)et2(1)+et2(2),ss(1),ss(2)
 198  format(' ***  requested jed,',f12.2,
     * ' not within ephemeris limits,',2f12.2,'  ***')

      stop

   99 write(*,'(2f12.2,a80)')et2,'error return in state'

      stop

      end




      subroutine const(nam,val,sss,n)
c
c+++++++++++++++++++++++++++++
c
c     this entry obtains the constants from the ephemeris file
c
c     calling seqeunce parameters (all output):
c
c       nam = character*6 array of constant names
c
c       val = d.p. array of values of constants
c
c       sss = d.p. jd start, jd stop, step of ephemeris
c
c         n = integer number of entries in 'nam' and 'val' arrays
c
      implicit double precision (a-h,o-z)

      save

      character*6 nam(*),ttl(14,3),cnam(400)

      double precision val(*),sss(3),ss(3),cval(400)

      double precision et2(2),pv(6,12),rrd(4)
      integer          list(12)

      integer ipt(3,13),denum

      common/ephhdr/cval,ss,au,emrat,denum,ncon,ipt
      common/chrhdr/cnam,ttl

c  call state to initialize the ephemeris and read in the constants

c      call state(0.d0,0,0,0)
      call state(et2,list,pv,rrd)

      n=ncon

      do i=1,3
      sss(i)=ss(i)
      enddo

      do i=1,n
      nam(i)=cnam(i)
      val(i)=cval(i)
      enddo

      return

      end