*  novas fortran vers f2.0 (1 nov 98)
*  standard set of subroutines

************************************************************************
*                                                                      *
*                              n o v a s                               *
*           naval observatory vector astrometry subroutines            *
*                                                                      *
*                            g. h. kaplan                              *
*                        u.s. naval observatory                        *
*                                                                      *
************************************************************************



      subroutine apstar (tjd,n,ram,decm,pmra,pmdec,parlax,radvel,ra,dec)
c
c     this subroutine computes the apparent place of a star,
c     given its mean place, proper motion, parallax, and radial
c     velocity for j2000.0.  see kaplan, et al. (1989) astronomical
c     journal 97, 1197-1210.
c
c          tjd    = tt julian date for apparent place (in)
c          n      = body identification number for the earth (in)
c          ram    = mean right ascension j2000.0 in hours (in)
c          decm   = mean declination j2000.0 in degrees (in)
c          pmra   = proper motion in ra in time seconds per julian
c                   century (in)
c          pmdec  = proper motion in dec in arcseconds per julian
c                   century (in)
c          parlax = parallax in arcseconds (in)
c          radvel = radial velocity in kilometers per second (in)
c          ra     = apparent right ascension in hours, referred to
c                   true equator and equinox of date (out)
c          dec    = apparent declination in degrees, referred to
c                   true equator and equinox of date (out)
c
c
      double precision tjd,ujd,ram,decm,pmra,pmdec,parlax,radvel,
     .     glon,glat,ht,ra,dec,t0,t1,tlast,
     .     x,secdif,eqeq,st,gast,rm,dm,pmr,pmd,pi,rv,tlight,r,d,
     .     peb,veb,pes,ves,pog,vog,pb,vb,ps,
     .     pos1,vel1,pos2,vel2,pos3,pos4,pos5,pos6,pos7,
     .     dabs,dmod
      dimension peb(3), veb(3), pes(3), ves(3), pog(3), vog(3),
     .     pb(3), vb(3), ps(3),
     .     pos1(3), vel1(3), pos2(3), vel2(3),
     .     pos3(3), pos4(3), pos5(3), pos6(3), pos7(3)
      save
c
      data t0 / 2451545.00000000d0 /
c     t0 = tdb julian date of epoch j2000.0
      data tlast / 0.0d0 /
c
c     compute t1, the tdb julian date corresponding to tjd
      call novas_times(tjd,x,secdif)
      t1 = tjd + secdif / 86400.0d0
      if (dabs(tjd-tlast).lt.1.0d-8) go to 20
c
c     get position and velocity of the earth wrt barycenter of
c     solar system and wrt center of sun
      call solsys (t1,n,0,peb,veb,ierr)
      if (ierr.ne.0) go to 40
      call solsys (t1,n,1,pes,ves,ierr)
      if (ierr.ne.0) go to 40
      tlast = tjd
c
   20 do 22 j=1,3
      pb(j) = peb(j)
      vb(j) = veb(j)
      ps(j) = pes(j)
   22 continue
      rm = ram
      dm = decm
      pmr = pmra
      pmd = pmdec
      pi = parlax
      rv = radvel
c
c     compute apparent place
   30 call vectrs (rm,dm,pmr,pmd,pi,rv,pos1,vel1)
      call propmo (t0,pos1,vel1,t1,pos2)
      call geocen (pos2,pb,pos3,tlight)
      call sunfld (pos3,ps,pos4)
      call aberat (pos4,vb,tlight,pos5)
      call preces (t0,pos5,t1,pos6)
      call nutate (t1,pos6,pos7)
      call angles (pos7,r,d)
c
      ra = r
      dec = d
      return
c
   40 ra = 0.0d0
      dec = 0.0d0
      tlast = 0.0d0
      return
c
c
      entry tpstar (ujd,glon,glat,ht,ra,dec)
c
c     this entry computes the topocentric place of a star,
c     given the location of the observer.  this entry assumes apstar
c     was previously called, and uses data computed by apstar.
c
c          ujd    = ut1 julian date, or equivalent greenwich apparent
c                   sidereal time in hours, for topocentric place (in)
c          glon   = geodetic longitude (east +) of observer
c                   in degrees (in)
c          glat   = geodetic latitude (north +) of observer
c                   in degrees (in)
c          ht     = height of observer in meters (in)
c          ra     = topocentric right ascension in hours, referred to
c                   true equator and equinox of date (out)
c          dec    = topocentric declination in degrees, referred to
c                   true equator and equinox of date (out)
c
c
      if (tlast.eq.0.0d0) go to 40
c
c     get position and velocity of observer wrt center of earth
   50 if (ujd.gt.100.0d0) go to 52
      gast = dmod(ujd,24.0d0)
      go to 55
   52 call sidtim (ujd,0.0d0,0,st)
      call etilt (t1,x,x,eqeq,x,x)
      gast = st + eqeq/3600.0d0
   55 call terra (glon,glat,ht,gast,pos1,vel1)
      call nutate (-t1,pos1,pos2)
      call preces (t1,pos2,t0,pog)
      call nutate (-t1,vel1,vel2)
      call preces (t1,vel2,t0,vog)
c
c     compute position and velocity of observer wrt barycenter of
c     solar system and position wrt center of sun
   60 do 62 j=1,3
      pb(j) = peb(j) + pog(j)
      vb(j) = veb(j) + vog(j)
      ps(j) = pes(j) + pog(j)
   62 continue
c
c     recompute apparent place using position and velocity of observer
      go to 30
c
      end



      subroutine applan (tjd,l,n,pout,vout,ra,dec,dis,drdt)
      implicit none
      save
c
c     this subroutine computes the apparent place of a planet or other
c     solar system body.  rectangular coordinates of solar system bodies
c     are obtained from subroutine solsys.  see kaplan, et al. (1989)
c     astronomical journal 97, 1197-1210.
c
c          tjd    = tt julian date for apparent place (in)
c          l      = body identification number for desired planet (in)
c          n      = body identification number for the earth (in)
c          pout   = position vector in au, referred to
c                   true equator and equinox of date (out)
c          vout   = velocity vector in au/day, referred to
c                   true equator and equinox of date (out)
c          ra     = apparent right ascension in hours, referred to
c                   true equator and equinox of date (out)
c          dec    = apparent declination in degrees, referred to
c                   true equator and equinox of date (out)
c          dis    = true distance from earth to planet in au (out)
c          drst   = true radial velocity from earth to planet in au/day (out)
c
c
      integer          l,n,ierr,lplan,j
      double precision tjd,ujd,glon,glat,ht,pout(3),vout(3),ra,dec,dis,
     .     drdt,t0,t1,t2,t3,tlast,c,
     .     x,secdif,eqeq,st,gast,tlight,r,d,s,
     .     peb,veb,pes,ves,pog,vog,pb,vb,ps,
     .     pos1,vel1,pos2,vel2,pos3,vel3,pos4,vel4,pos5,vel5,pos6,vel6,
     .     dabs,dmod,tdum,sdot
      dimension peb(3), veb(3), pes(3), ves(3), pog(3), vog(3),
     .     pb(3), vb(3), ps(3),
     .     pos1(3), vel1(3), pos2(3), vel2(3),
     .     pos3(3), vel3(3), pos4(3), vel4(3),
     .     pos5(3), vel5(3), pos6(3), vel6(3)


c
      data c / 173.14463348d0 /
c     c = speed of light in au/day
      data t0 / 2451545.00000000d0 /
c     t0 = tdb julian date of epoch j2000.0
      data tlast / 0.0d0 /

c
      if (l.eq.n) go to 40
c
c     compute t1, the tdb julian date corresponding to tjd
      call novas_times(tjd,x,secdif)
      t1 = tjd + secdif / 86400.0d0
      if (dabs(tjd-tlast).lt.1.0d-8) go to 20
c
c     get position and velocity of the earth wrt barycenter of
c     solar system and wrt center of sun
      call solsys (t1,n,0,peb,veb,ierr)
      if (ierr.ne.0) go to 40
      call solsys (t1,n,1,pes,ves,ierr)
      if (ierr.ne.0) go to 40
      tlast = tjd
c
   20 do 22 j=1,3
      pb(j) = peb(j)
      vb(j) = veb(j)
      ps(j) = pes(j)
   22 continue
      lplan = l

c
c     get position of planet wrt barycenter of solar system
   30 call solsys (t1,lplan,0,pos1,vel1,ierr)
      if (ierr.ne.0) go to 40
      call geocen (pos1,pb,pos2,tlight)
      s = tlight * c
      call geocen (vel1,vb,vel2,tdum)
      if (abs(s) .lt. 1.0d-12) then
       sdot = 0.0d0
      else
       sdot = (pos2(1)*vel2(1) + pos2(2)*vel2(2) + pos2(3)*vel2(3))/s
      end if
      t2 = t1 - tlight
   33 call solsys (t2,lplan,0,pos1,vel1,ierr)
      if (ierr.ne.0) go to 40
      call geocen (pos1,pb,pos2,tlight)
      call geocen (vel1,vb,vel2,tdum)
      t3 = t1 - tlight
      if (dabs(t3-t2).lt.1.0d-8) go to 35
      t2 = t3
      go to 33

c
c     finish apparent place computation
 35   continue
      call sunfld (pos2,ps,pos3)
      call aberat (pos3,vb,tlight,pos4)
      call preces (t0,pos4,t1,pos5)
      call nutate (t1,pos5,pos6)

c..form the output
      do j=1,3
       pout(j) = pos6(j)
       vout(j) = vel2(j)
      end do
      call angles (pout,r,d)
      ra   = r
      dec  = d
      dis  = s
      drdt = sdot
      return
c
c..land here if we had an error
   40 continue
      do j=1,3
       pout(j) = 0.0d0
       vout(j) = 0.0d0
      enddo
      ra    = 0.0d0
      dec   = 0.0d0
      dis   = 0.0d0
      drdt  = 0.0d0
      tlast = 0.0d0
      return

c
c
      entry tpplan (ujd,glon,glat,ht,pout,vout,ra,dec,dis,drdt)
c
c     this entry computes the topocentric place of a planet,
c     given the location of the observer.  this entry assumes applan
c     was previously called, and uses data computed by applan.
c
c          ujd    = ut1 julian date, or equivalent greenwich apparent
c                   sidereal time in hours, for topocentric place (in)
c          glon   = geodetic longitude (east +) of observer
c                   in degrees (in)
c          glat   = geodetic latitude (north +) of observer
c                   in degrees (in)
c          ht     = height of observer in meters (in)
c          pout   = position vector in au, referred to
c                   true equator and equinox of date (out)
c          vout   = velocity vector in au/day, referred to
c                   true equator and equinox of date (out)
c          ra     = topocentric right ascension in hours, referred to
c                   true equator and equinox of date (out)
c          dec    = topocentric declination in degrees, referred to
c                   true equator and equinox of date (out)
c          dis    = true distance from observer to planet in au (out)
c          drst   = true radial velocity from earth to planet in au/day (out)
c
c
      if (tlast.eq.0.0d0) go to 40
c
c     get position and velocity of observer wrt center of earth
   50 if (ujd.gt.100.0d0) go to 52
      gast = dmod(ujd,24.0d0)
      go to 55
   52 call sidtim (ujd,0.0d0,0,st)
      call etilt (t1,x,x,eqeq,x,x)
      gast = st + eqeq/3600.0d0
   55 call terra (glon,glat,ht,gast,pos1,vel1)
      call nutate (-t1,pos1,pos2)
      call preces (t1,pos2,t0,pog)
      call nutate (-t1,vel1,vel2)
      call preces (t1,vel2,t0,vog)
c
c     compute position and velocity of observer wrt barycenter of
c     solar system and position wrt center of sun
   60 do 62 j=1,3
      pb(j) = peb(j) + pog(j)
      vb(j) = veb(j) + vog(j)
      ps(j) = pes(j) + pog(j)
   62 continue
c
c     recompute apparent place using position and velocity of observer
      go to 30
c
      end



      subroutine applan2 (tjd,l,n,pout,vout,ra,dec,dis,drdt)
      implicit none
      save
c
c     this subroutine computes the apparent place of a planet or other
c     solar system body.  rectangular coordinates of solar system bodies
c     are obtained from subroutine solsys.  see kaplan, et al. (1989)
c     astronomical journal 97, 1197-1210.
c
c          tjd    = tt julian date for apparent place (in)
c          l      = body identification number for desired planet (in)
c          n      = body identification number for the earth (in)
c          pout   = position vector in au, referred to
c                   true equator and equinox of date (out)
c          vout   = velocity vector in au/day, referred to
c                   true equator and equinox of date (out)
c          ra     = apparent right ascension in hours, referred to
c                   true equator and equinox of date (out)
c          dec    = apparent declination in degrees, referred to
c                   true equator and equinox of date (out)
c          dis    = true distance from earth to planet in au (out)
c          drst   = true radial velocity from earth to planet in au/day (out)
c
c
      integer          l,n,ierr,lplan,j
      double precision tjd,ujd,glon,glat,ht,pout(3),vout(3),ra,dec,dis,
     .     drdt,t0,t1,t2,t3,tlast,c,
     .     x,secdif,eqeq,st,gast,tlight,r,d,s,
     .     peb,veb,pes,ves,pog,vog,pb,vb,ps,
     .     pos1,vel1,pos2,vel2,pos3,vel3,pos4,vel4,pos5,vel5,pos6,vel6,
     .     dabs,dmod,tdum,sdot,posvel(6)
      dimension peb(3), veb(3), pes(3), ves(3), pog(3), vog(3),
     .     pb(3), vb(3), ps(3),
     .     pos1(3), vel1(3), pos2(3), vel2(3),
     .     pos3(3), vel3(3), pos4(3), vel4(3),
     .     pos5(3), vel5(3), pos6(3), vel6(3)

c
      data c / 173.14463348d0 /
c     c = speed of light in au/day
      data t0 / 2451545.00000000d0 /
c     t0 = tdb julian date of epoch j2000.0
      data tlast / 0.0d0 /

c
c      if (l.eq.n) go to 40

c     compute t1, the tdb julian date corresponding to tjd
      call novas_times(tjd,x,secdif)
      t1 = tjd + secdif / 86400.0d0
      if (dabs(tjd-tlast).lt.1.0d-8) go to 20

c     get position and velocity of the earth wrt barycenter of
c     solar system and wrt center of sun
c      call solsys (t1,n,0,peb,veb,ierr)
c      if (ierr.ne.0) go to 40
c      call solsys (t1,n,1,pes,ves,ierr)
c      if (ierr.ne.0) go to 40
c      tlast = tjd
c

c..get the position of sun with respect to the barycenter
c      call pleph(t1,11,12,posvel)
c
c   20 do 22 j=1,3
c      pb(j) = posvel(j)
c      vb(j) = posvel(j+3)
c      ps(j) = pes(j)
c   22 continue

 20   continue
      lplan = l

c
c     get position of planet wrt center of solar system
   30 call solsys (t1,lplan,1,pos1,vel1,ierr)
      if (ierr.ne.0) go to 40

c      write(6,222) pos1

      s = dsqrt(pos1(1)**2 + pos1(2)**2 + pos1(3)**2) 
      tlight = s/c

c      call geocen (pos1,pb,pos2,tlight)
c      s = tlight * c
c      call geocen (vel1,vb,vel2,tdum)

      if (abs(s) .lt. 1.0d-12) then
       sdot = 0.0d0
      else
       sdot = (pos1(1)*vel1(1) + pos1(2)*vel1(2) + pos1(3)*vel1(3))/s
      end if

      t2 = t1 - tlight
   33 call solsys (t2,lplan,1,pos1,vel1,ierr)
      if (ierr.ne.0) go to 40

      x = dsqrt(pos1(1)**2 + pos1(2)**2 + pos1(3)**2) 
      tlight = x/c

c      call geocen (pos1,pb,pos2,tlight)
c      call geocen (vel1,vb,vel2,tdum)


      t3 = t1 - tlight
      if (dabs(t3-t2).lt.1.0d-8) go to 35
      t2 = t3
      go to 33

c
c     finish apparent place computation
 35   continue
c      call sunfld (pos1,ps,pos2)
c      call sunfld (pos1,pb,pos2)
c      call aberat (pos2,vb,tlight,pos3)
c      call preces (t0,pos3,t1,pos4)
c      call nutate (t1,pos4,pos5)

c      write(6,222) pos1
c      write(6,222) pos2
c      write(6,222) pos3
c      write(6,222) pos4
c      write(6,222) pos5
 222  format(1x,1p3e14.6)


      do j=1,3
       pout(j) = pos1(j)
       vout(j) = vel1(j)
      end do
      call angles (pout,r,d)
      ra   = r
      dec  = d
      dis  = s
      drdt = sdot
      return
c
c..land here if we had an error
   40 continue
      do j=1,3
       pout(j) = 0.0d0
       vout(j) = 0.0d0
      enddo
      ra    = 0.0d0
      dec   = 0.0d0
      dis   = 0.0d0
      drdt  = 0.0d0
      tlast = 0.0d0
      return
      end




      subroutine vpstar (tjd,n,ram,decm,pmra,pmdec,parlax,radvel,ra,dec)
c
c     this subroutine computes the virtual place of a star,
c     given its mean place, proper motion, parallax, and radial
c     velocity for j2000.0.  see kaplan, et al. (1989) astronomical
c     journal 97, 1197-1210.
c
c          tjd    = tt julian date for virtual place (in)
c          n      = body identification number for the earth (in)
c          ram    = mean right ascension j2000.0 in hours (in)
c          decm   = mean declination j2000.0 in degrees (in)
c          pmra   = proper motion in ra in time seconds per julian
c                   century (in)
c          pmdec  = proper motion in dec in arcseconds per julian
c                   century (in)
c          parlax = parallax in arcseconds (in)
c          radvel = radial velocity in kilometers per second (in)
c          ra     = virtual right ascension in hours, referred to
c                   mean equator and equinox of j2000.0 (out)
c          dec    = virtual declination in degrees, referred to
c                   mean equator and equinox of j2000.0 (out)
c
c
      double precision tjd,ujd,ram,decm,pmra,pmdec,parlax,radvel,
     .     glon,glat,ht,ra,dec,t0,t1,tlast,
     .     x,secdif,eqeq,st,gast,rm,dm,pmr,pmd,pi,rv,tlight,r,d,
     .     peb,veb,pes,ves,pog,vog,pb,vb,ps,
     .     pos1,vel1,pos2,vel2,pos3,pos4,pos5,
     .     dabs,dmod
      dimension peb(3), veb(3), pes(3), ves(3), pog(3), vog(3),
     .     pb(3), vb(3), ps(3),
     .     pos1(3), vel1(3), pos2(3), vel2(3),
     .     pos3(3), pos4(3), pos5(3)
      save
c
      data t0 / 2451545.00000000d0 /
c     t0 = tdb julian date of epoch j2000.0
      data tlast/ 0.0d0 /
c
c     compute t1, the tdb julian date corresponding to tjd
      call novas_times(tjd,x,secdif)
      t1 = tjd + secdif / 86400.0d0
      if (dabs(tjd-tlast).lt.1.0d-8) go to 20
c
c     get position and velocity of the earth wrt barycenter of
c     solar system and wrt center of sun
      call solsys (t1,n,0,peb,veb,ierr)
      if (ierr.ne.0) go to 40
      call solsys (t1,n,1,pes,ves,ierr)
      if (ierr.ne.0) go to 40
      tlast = tjd
c
   20 do 22 j=1,3
      pb(j) = peb(j)
      vb(j) = veb(j)
      ps(j) = pes(j)
   22 continue
      rm = ram
      dm = decm
      pmr = pmra
      pmd = pmdec
      pi = parlax
      rv = radvel
c
c     compute virtual place
   30 call vectrs (rm,dm,pmr,pmd,pi,rv,pos1,vel1)
      call propmo (t0,pos1,vel1,t1,pos2)
      call geocen (pos2,pb,pos3,tlight)
      call sunfld (pos3,ps,pos4)
      call aberat (pos4,vb,tlight,pos5)
      call angles (pos5,r,d)
c
      ra = r
      dec = d
      return
c
   40 ra = 0.0d0
      dec = 0.0d0
      tlast = 0.0d0
      return
c
c
      entry lpstar (ujd,glon,glat,ht,ra,dec)
c
c     this entry computes the local place of a star,
c     given the location of the observer.  this entry assumes vpstar
c     was previously called, and uses data computed by vpstar.
c
c          ujd    = ut1 julian date, or equivalent greenwich apparent
c                   sidereal time in hours, for local place (in)
c          glon   = geodetic longitude (east +) of observer
c                   in degrees (in)
c          glat   = geodetic latitude (north +) of observer
c                   in degrees (in)
c          ht     = height of observer in meters (in)
c          ra     = local right ascension in hours, referred to
c                   mean equator and equinox of j2000.0 (out)
c          dec    = local declination in degrees, referred to
c                   mean equator and equinox of j2000.0 (out)
c
c
      if (tlast.eq.0.0d0) go to 40
c
c     get position and velocity of observer wrt center of earth
   50 if (ujd.gt.100.0d0) go to 52
      gast = dmod(ujd,24.0d0)
      go to 55
   52 call sidtim (ujd,0.0d0,0,st)
      call etilt (t1,x,x,eqeq,x,x)
      gast = st + eqeq/3600.0d0
   55 call terra (glon,glat,ht,gast,pos1,vel1)
      call nutate (-t1,pos1,pos2)
      call preces (t1,pos2,t0,pog)
      call nutate (-t1,vel1,vel2)
      call preces (t1,vel2,t0,vog)
c
c     compute position and velocity of observer wrt barycenter of
c     solar system and wrt center of sun
   60 do 62 j=1,3
      pb(j) = peb(j) + pog(j)
      vb(j) = veb(j) + vog(j)
      ps(j) = pes(j) + pog(j)
   62 continue
c
c     recompute virtual place using position and velocity of observer
      go to 30
c
      end



      subroutine vpplan (tjd,l,n,ra,dec,dis)
c
c     this subroutine computes the virtual place of a planet or other
c     solar system body.  rectangular coordinates of solar system bodies
c     are obtained from subroutine solsys.  see kaplan, et al. (1989)
c     astronomical journal 97, 1197-1210.
c
c          tjd    = tt julian date for virtual place (in)
c          l      = body identification number for desired planet (in)
c          n      = body identification number for the earth (in)
c          ra     = virtual right ascension in hours, referred to
c                   mean equator and equinox of j2000.0 (out)
c          dec    = virtual declination in degrees, referred to
c                   mean equator and equinox of j2000.0 (out)
c          dis    = true distance from earth to planet in au (out)
c
c
      double precision tjd,ujd,glon,glat,ht,ra,dec,dis,
     .     t0,t1,t2,t3,tlast,c,
     .     x,secdif,eqeq,st,gast,tlight,r,d,s,
     .     peb,veb,pes,ves,pog,vog,pb,vb,ps,
     .     pos1,vel1,pos2,vel2,pos3,pos4,
     .     dabs,dmod
      dimension peb(3), veb(3), pes(3), ves(3), pog(3), vog(3),
     .     pb(3), vb(3), ps(3),
     .     pos1(3), vel1(3), pos2(3), vel2(3),
     .     pos3(3), pos4(3)
      save
c
      data c / 173.14463348d0 /
c     c = speed of light in au/day
      data t0 / 2451545.00000000d0 /
c     t0 = tdb julian date of epoch j2000.0
      data tlast / 0.0d0 /
c
      if (l.eq.n) go to 40
c
c     compute t1, the tdb julian date corresponding to tjd
      call novas_times(tjd,x,secdif)
      t1 = tjd + secdif / 86400.0d0
      if (dabs(tjd-tlast).lt.1.0d-8) go to 20
c
c     get position and velocity of the earth wrt barycenter of
c     solar system and wrt center of sun
      call solsys (t1,n,0,peb,veb,ierr)
      if (ierr.ne.0) go to 40
      call solsys (t1,n,1,pes,ves,ierr)
      if (ierr.ne.0) go to 40
      tlast = tjd
c
   20 do 22 j=1,3
      pb(j) = peb(j)
      vb(j) = veb(j)
      ps(j) = pes(j)
   22 continue
      lplan = l
c
c     get position of planet wrt barycenter of solar system
   30 call solsys (t1,lplan,0,pos1,vel1,ierr)
      if (ierr.ne.0) go to 40
      call geocen (pos1,pb,pos2,tlight)
      s = tlight * c
      t2 = t1 - tlight
   33 call solsys (t2,lplan,0,pos1,vel1,ierr)
      if (ierr.ne.0) go to 40
      call geocen (pos1,pb,pos2,tlight)
      t3 = t1 - tlight
      if (dabs(t3-t2).lt.1.0d-8) go to 35
      t2 = t3
      go to 33
c
c     finish virtual place computation
   35 continue
      call sunfld (pos2,ps,pos3)
      call aberat (pos3,vb,tlight,pos4)
      call angles (pos4,r,d)
      ra = r
      dec = d
      dis = s
      return
c
   40 ra = 0.0d0
      dec = 0.0d0
      dis = 0.0d0
      tlast = 0.0d0
      return
c
c
      entry lpplan (ujd,glon,glat,ht,ra,dec,dis)
c
c     this entry computes the local place of a planet, given
c     the location of the observer.  this entry assumes vpplan was
c     previously called, and uses data computed by vpplan.
c
c          ujd    = ut1 julian date, or equivalent greenwich apparent
c                   sidereal time in hours, for local place (in)
c          glon   = geodetic longitude (east +) of observer
c                   in degrees (in)
c          glat   = geodetic latitude (north +) of observer
c                   in degrees (in)
c          ht     = height of observer in meters (in)
c          ra     = local right ascension in hours, referred to
c                   mean equator and equinox of j2000.0 (out)
c          dec    = local declination in degrees, referred to
c                   mean equator and equinox of j2000.0 (out)
c          dis    = true distance from observer to planet in au (out)
c
c
      if (tlast.eq.0.0d0) go to 40
c
c     get position and velocity of observer wrt center of earth
   50 if (ujd.gt.100.0d0) go to 52
      gast = dmod(ujd,24.0d0)
      go to 55
   52 call sidtim (ujd,0.0d0,0,st)
      call etilt (t1,x,x,eqeq,x,x)
      gast = st + eqeq/3600.0d0
   55 call terra (glon,glat,ht,gast,pos1,vel1)
      call nutate (-t1,pos1,pos2)
      call preces (t1,pos2,t0,pog)
      call nutate (-t1,vel1,vel2)
      call preces (t1,vel2,t0,vog)
c
c     compute position and velocity of observer wrt barycenter of
c     solar system and wrt center of sun
   60 do 62 j=1,3
      pb(j) = peb(j) + pog(j)
      vb(j) = veb(j) + vog(j)
      ps(j) = pes(j) + pog(j)
   62 continue
c
c     recompute virtual place using position and velocity of observer
      go to 30
c
      end



      subroutine asstar (tjd,n,ram,decm,pmra,pmdec,parlax,radvel,ra,dec)
c
c     this subroutine computes the astrometric place of a star,
c     given its mean place, proper motion, parallax, and radial
c     velocity for j2000.0.   see kaplan, et al. (1989) astronomical
c     journal 97, 1197-1210.
c
c          tjd    = tt julian date for astrometric place (in)
c          n      = body identification number for the earth (in)
c          ram    = mean right ascension j2000.0 in hours (in)
c          decm   = mean declination j2000.0 in degrees (in)
c          pmra   = proper motion in ra in time seconds per julian
c                   century (in)
c          pmdec  = proper motion in dec in arcseconds per julian
c                   century (in)
c          parlax = parallax in arcseconds (in)
c          radvel = radial velocity in kilometers per second (in)
c          ra     = astrometric right ascension in hours, referred to
c                   mean equator and equinox of j2000.0 (out)
c          dec    = astrometric declination in degrees, referred to
c                   mean equator and equinox of j2000.0 (out)
c
c
      double precision tjd,ram,decm,pmra,pmdec,parlax,radvel,ra,dec,
     .     t0,t1,tlast,x,secdif,rm,dm,pmr,pmd,pi,rv,tlight,r,d,
     .     peb,veb,pb,pos1,vel1,pos2,pos3,dabs
      dimension peb(3), veb(3), pb(3),
     .     pos1(3), vel1(3), pos2(3), pos3(3)
      save
c
      data t0 / 2451545.00000000d0 /
c     t0 = tdb julian date of epoch j2000.0
      data tlast / 0.0d0 /
c
c     compute t1, the tdb julian date corresponding to tjd
      call novas_times(tjd,x,secdif)
      t1 = tjd + secdif / 86400.0d0
      if (dabs(tjd-tlast).lt.1.0d-8) go to 20
c
c     get position and velocity of the earth wrt barycenter of
c     solar system
      call solsys (t1,n,0,peb,veb,ierr)
      if (ierr.ne.0) go to 40
      tlast = tjd
c
   20 do 22 j=1,3
      pb(j) = peb(j)
   22 continue
      rm = ram
      dm = decm
      pmr = pmra
      pmd = pmdec
      pi = parlax
      rv = radvel
c
c     compute astrometric place
   30 call vectrs (rm,dm,pmr,pmd,pi,rv,pos1,vel1)
      call propmo (t0,pos1,vel1,t1,pos2)
      call geocen (pos2,pb,pos3,tlight)
      call angles (pos3,r,d)
      ra = r
      dec = d
      return
c
   40 ra = 0.0d0
      dec = 0.0d0
      tlast = 0.0d0
      return
c
      end



      subroutine asplan (tjd,l,n,ra,dec,dis)
c
c     this subroutine computes the astrometric place of a planet or
c     other solar system body.   rectangular coordinates of solar system
c     bodies are obtained from subroutine solsys.  see kaplan, et al.
c     astronomical journal 97, 1197-1210.
c
c          tjd    = tt julian date for astrometric place (in)
c          l      = body identification number for desired planet (in)
c          n      = body identification number for the earth (in)
c          ra     = astrometric right ascension in hours, referred to
c                   mean equator and equinox of j2000.0 (out)
c          dec    = astrometric declination in degrees, referred to
c                   mean equator and equinox of j2000.0 (out)
c          dis    = true distance from earth to planet in au (out)
c
c
      double precision tjd,ra,dec,dis,t1,t2,t3,tlast,x,secdif,
     .     c,tlight,r,d,s,
     .     peb,veb,pb,pos1,vel1,pos2,dabs
      dimension peb(3), veb(3), pb(3),
     .     pos1(3), vel1(3), pos2(3)
      save
c
      data c / 173.14463348d0 /
c     c = speed of light in au/day
      data tlast / 0.0d0 /
c
      if (l.eq.n) go to 40
c
c     compute t1, the tdb julian date corresponding to tjd
      call novas_times(tjd,x,secdif)
      t1 = tjd + secdif / 86400.0d0
      if (dabs(tjd-tlast).lt.1.0d-8) go to 20
c
c     get position and velocity of the earth wrt barycenter of
c     solar system
      call solsys (t1,n,0,peb,veb,ierr)
      if (ierr.ne.0) go to 40
      tlast = tjd
c
   20 do 22 j=1,3
      pb(j) = peb(j)
   22 continue
      lplan = l
c
c     get position of planet wrt barycenter of solar system
   30 call solsys (t1,lplan,0,pos1,vel1,ierr)
      if (ierr.ne.0) go to 40
      call geocen (pos1,pb,pos2,tlight)
      s = tlight * c
      t2 = t1 - tlight
   33 call solsys (t2,lplan,0,pos1,vel1,ierr)
      if (ierr.ne.0) go to 40
      call geocen (pos1,pb,pos2,tlight)
      t3 = t1 - tlight
      if (dabs(t3-t2).lt.1.0d-8) go to 35
      t2 = t3
      go to 33
c
c     finish astrometric place computation
   35 continue
      call angles (pos2,r,d)
      ra = r
      dec = d
      dis = s
      return
c
   40 ra = 0.0d0
      dec = 0.0d0
      dis = 0.0d0
      tlast = 0.0d0
      return
c
      end



      subroutine mpstar (tjd,n,ra,dec,ram,decm)
c
c     this subroutine computes the mean place of a star for j2000.0,
c     given its apparent place at date tjd.  proper motion, parallax,
c     and radial velocity are assumed to be zero.
c
c          tjd    = tt julian date of apparent place (in)
c          n      = body identification number for the earth (in)
c          ra     = apparent right ascension in hours, referred to
c                   true equator and equinox of date (in)
c          dec    = apparent declination in degrees, referred to
c                   true equator and equinox of date (in)
c          ram    = mean right ascension j2000.0 in hours (out)
c          decm   = mean declination j2000.0 in degrees (out)
c
c
      double precision tjd,ra,dec,ram,decm,t1,ramnew,dcmnew,
     .     ramold,dcmold,r,d,delra,deldec,dabs,dmod
c
      t1 = tjd
      ramnew = dmod(ra,24.0d0)
      if (ramnew.lt.0.0d0) ramnew = ramnew + 24.0d0
      dcmnew = dec
      iter = 0
c
   20 iter = iter + 1
      ramold = ramnew
      dcmold = dcmnew
      r = ramold
      d = dcmold
      call apstar (t1,n,r,d,0.0d0,0.0d0,0.0d0,0.0d0,r,d)
      delra = r - ramold
      deldec = d - dcmold
      if (delra.lt.-12.0d0) delra = delra + 24.0d0
      if (delra.gt.+12.0d0) delra = delra - 24.0d0
      ramnew = ra - delra
      dcmnew = dec - deldec
      if (iter.gt.20) go to 40
      if (dabs(ramnew-ramold).gt.1.0d-10) go to 20
      if (dabs(dcmnew-dcmold).gt.1.0d-09) go to 20
c
      ram = ramnew
      decm = dcmnew
      if (ram.lt. 0.0d0) ram = ram + 24.0d0
      if (ram.ge.24.0d0) ram = ram - 24.0d0
      return
c
   40 ram = 0.0d0
      decm = 0.0d0
      return
c
      end



      subroutine sidtim (tjdh,tjdl,k,gst)
c
c     this subroutine computes the greenwich sidereal time
c     (either mean or apparent) at julian date tjdh + tjdl.
c     see aoki, et al. (1982) astronomy and astropysics 105, 359-361.
c
c          tjdh   = julian date, high-order part (in)
c          tjdl   = julian date, low-order part (in)
c                   julian date may be split at any point, but
c                   for highest precision, set tjdh to be the integral
c                   part of the julian date, and set tjdl to be the
c                   fractional part
c          k      = time selection code (in)
c                   set k=0 for greenwich mean sidereal time
c                   set k=1 for greenwich apparent sidereal time
c          gst    = greenwich (mean or apparent) sidereal time
c                   in hours (out)
c
c     note:  for most applications, basis for input julian date should
c     be ut1, which results in ordinary sidereal time output in gst.
c     use of input julian date based on tdb results in 'dynamical
c     sidereal time'.
c
c
      double precision tjdh,tjdl,tjd,th,tl,t0,t,t2,t3,gst,
     .     x,eqeq,st,dmod
c
      data t0 / 2451545.00000000d0 /
c     t0 = tdb julian date of epoch j2000.0
c
      tjd = tjdh + tjdl
      th = (tjdh - t0) / 36525.0d0
      tl =  tjdl       / 36525.0d0
      t = th + tl
      t2 = t * t
      t3 = t2 * t
c
c     for apparent sidereal time, obtain equation of the equinoxes
      eqeq = 0.0d0
      if (k.eq.1) call etilt (tjd,x,x,eqeq,x,x)
c
      st = eqeq - 6.2d-6*t3 + 0.093104d0*t2 + 67310.54841d0
     .     + 8640184.812866d0 *tl
     .     + 3155760000.0d0   *tl
     .     + 8640184.812866d0 *th
     .     + 3155760000.0d0   *th
c
      gst = dmod (st / 3600.0d0, 24.0d0)
      if (gst.lt.0.0d0) gst = gst + 24.0d0
      return
c
      end



      subroutine pnsw (tjd,gast,x,y,vece,vecs)
c
c     transforms a vector from earth-fixed system to space-fixed system
c     by applying rotations for wobble, spin, nutation, and precession.
c     (combined rotation is symbolized  p n s w .)   specifically,
c     it transforms a vector from earth-fixed geographic system to
c     space-fixed system based on mean equator and equinox of j2000.0.
c
c          tjd    = tt julian date (in)
c          gast   = greenwich apparent sidereal time, in hours (in)
c          x      = conventionally-defined x coordinate of celestial
c                   ephemeris pole with respect to iers reference
c                   pole, in arcseconds (in)
c          y      = conventionally-defined y coordinate of celestial
c                   ephemeris pole with respect to iers reference
c                   pole, in arcseconds (in)
c          vece   = vector in geocentric rectangular
c                   earth-fixed system, referred to geographic
c                   equator and greenwich meridian (in)
c          vecs   = vector in geocentric rectangular
c                   space-fixed system, referred to mean equator
c                   and equinox of j2000.0 (out)
c
c     note:  tjd=0.d0 means no precession/nutation transformation,
c     gast=0.d0 means no spin transformation, x=y=0.d0 means no
c     wobble transformation.
c
c
      double precision tjd,gast,x,y,vece,vecs,t0,t1,z,secdif,
     .     v1,v2,v3
      dimension vece(3), vecs(3), v1(3), v2(3), v3(3)
c
      data t0 / 2451545.00000000d0 /
c     t0 = tdb julian date of epoch j2000.0
c
c     compute t1, the tdb julian date corresponding to tjd
      if (tjd.eq.0.0d0) go to 20
      call novas_times(tjd,z,secdif)
      t1 = tjd + secdif / 86400.0d0
c
   20 if (x.eq.0.0d0 .and. y.eq.0.0d0) go to 25
      call wobble (x,y,vece,v1)
      go to 30
   25 do 28 j=1,3
   28 v1(j) = vece(j)
c
   30 if (gast.eq.0.0d0) go to 35
      call spin (gast,v1,v2)
      go to 40
   35 do 38 j=1,3
   38 v2(j) = v1(j)
c
   40 if (tjd.eq.0.0d0) go to 45
      call nutate (-t1,v2,v3)
      call preces (t1,v3,t0,vecs)
      go to 50
   45 do 48 j=1,3
   48 vecs(j) = v2(j)
c
   50 return
c
      end



      subroutine gethip (rah,dech,pmrah,pmdech,parxh,rvh,
     .                   ra2,dec2,pmra2,pmdec2,parx2,rv2)
c
c     this subroutine converts hipparcos data at epoch j1991.25
c     to epoch j2000.0 and fk5-style units.  to be used only for
c     hipparcos or tycho stars with linear space motion.
c
c          rah    = hipparcos right ascension in degrees (in)
c          dech   = hipparcos declination in degrees (in)
c          pmrah  = hipparcos proper motion in ra * cos(dech)
c                   in milliarcseconds per year (in)
c          pmdech = hipparcos proper motion in dec
c                   in milliarcseconds per year (in)
c          parxh  = hipparcos parallax in milliarcseconds (in)
c          rvh    = radial velocity at hipparcos epoch
c                   in kilometers per second (in)
c          ra2    = right ascension at j2000.0 in hours (out)
c          dec2   = declination at j2000.0 in degrees (out)
c          pmra2  = proper motion in ra at j2000.0
c                   in time seconds per julian century (out)
c          pmdec2 = proper motion in dec at j2000.0
c                   in arcseconds per julian century (out)
c          parx2  = parallax at j2000.0 in arcseconds (out)
c          rv2    = radial velocity at j2000.0 in kilometers
c                   per second (out)
c
c
      double precision rah,dech,pmrah,pmdech,parxh,rvh,
     .     ra2,dec2,pmra2,pmdec2,parx2,rv2,
     .     radcon,epoch1,epoch2,
     .     ra1,dec1,pmra1,pmdec1,parx1,rv1,dcos
c
      data radcon / 0.0174532925199433d0 /,
     .     epoch1, epoch2 / 2448349.0625d0, 2451545.0000d0 /
c
      ra1 = rah / 15.d0
      dec1 = dech
      pmra1 = pmrah / ( 150.d0 * dcos ( dec1 * radcon ) )
      pmdec1 = pmdech / 10.d0
      parx1 = parxh / 1000.d0
      rv1 = rvh
c
      call catran ( 1,epoch1,ra1,dec1,pmra1,pmdec1,parx1,rv1,
     .                epoch2,ra2,dec2,pmra2,pmdec2,parx2,rv2 )
c
      return
c
      end



      subroutine catran (it,date1,ra1,dec1,pmra1,pmdec1,parx1,rv1,
     .                      date2,ra2,dec2,pmra2,pmdec2,parx2,rv2)
c
c     this subroutine transforms a star's catalog quantities for
c     a change of epoch and/or equator and equinox.
c
c          it     = transformation option (in)
c                   set it=1 to change epoch (same equator and equinox)
c                   set it=2 to change equator and equinox (same epoch)
c                   set it=3 to change equator and equinox and epoch
c          date1  = tt julian date, or year, of original catalog
c                   data (the following six arguments) (in)
c          ra1    = original mean right ascension in hours (in)
c          dec1   = original mean declination in degrees (in)
c          pmra1  = original proper motion in ra
c                   in time seconds per julian century (in)
c          pmdec1 = original proper motion in dec
c                   in arcseconds per julian century (in)
c          parx1  = original parallax in arcseconds (in)
c          rv1    = original radial velocity in kilometers
c                   per second (out)
c          date2  = tt julian date, or year, for transformed
c                   output data (the following six arguments) (in)
c          ra2    = transformed mean right ascension in hours (out)
c          dec2   = transformed mean declination in degrees (out)
c          pmra2  = transformed proper motion in ra
c                   in time seconds per julian century (out)
c          pmdec2 = transformed proper motion in dec
c                   in arcseconds per julian century (out)
c          parx2  = transformed parallax in arcseconds (out)
c          rv2    = transformed radial velocity in kilometers
c                   per second (out)
c
c     note 1:  date1 and date2 may be specified either as a julian
c     date (e.g., 2433282.5d0) or a julian year and fraction
c     (e.g., 1950.0d0).  values less than 10000 are assumed to
c     be years.
c
c     note 2:  it=1 updates the star's data to account for
c     the star's space motion between the first and second dates,
c     within a fixed reference frame.  it=2 applies a rotation
c     of the reference frame corresponding to precession between
c     the first and second dates, but leaves the star fixed in space.
c     it=3 provides both transformations.
c
c     note 3:  this subroutine cannot be properly used to bring data
c     from old (pre-fk5) star catalogs into the modern system, because
c     old catalogs were compiled using a set of constants that are
c     incompatible with the iau (1976) system.
c
c
      double precision date1,ra1,dec1,pmra1,pmdec1,parx1,rv1,
     .     date2,ra2,dec2,pmra2,pmdec2,parx2,rv2,
     .     seccon,kmau,tjd1,pos1,vel1,tjd2,pos2,vel2,
     .     paralx,dist,r,d,cra,sra,cdc,sdc,pmr,pmd,rvl,
     .     xyproj,dcos,dsin,datan2
      integer it,j
c
      dimension pos1(3), vel1(3), pos2(3), vel2(3)
c
      data seccon / 206264.8062470964d0 /,   kmau / 1.49597870d8 /
c
c --- if necessary, compute julian dates -----------------------------
c
c     subroutine uses tdb julian dates internally, but no
c     distinction between tdb and tt is necessary
c
      if ( date1 .lt. 10000.d0 ) then
           tjd1 = 2451545.0d0 + ( date1 - 2000.d0 ) * 365.25d0
      else
           tjd1 = date1
      end if
      if ( date2 .lt. 10000.d0 ) then
           tjd2 = 2451545.0d0 + ( date2 - 2000.d0 ) * 365.25d0
      else
           tjd2 = date2
      end if
c
c --- convert input angular components to vectors --------------------
c
c     if parallax is unknown, undetermined, or zero, set it to 1e-7
c     arcsecond, corresponding to a distance of 10 megaparsecs
      paralx = parx1
      if ( paralx .le. 0.d0 ) paralx = 1.d-7
c
c     convert right ascension, declination, and parallax to position
c     vector in equatorial system with units of au
      dist = seccon / paralx
      r = ra1 * 54000.d0 / seccon
      d = dec1 * 3600.d0 / seccon
      cra = dcos(r)
      sra = dsin(r)
      cdc = dcos(d)
      sdc = dsin(d)
      pos1(1) = dist * cdc * cra
      pos1(2) = dist * cdc * sra
      pos1(3) = dist * sdc
c
c     convert proper motion and radial velocity to orthogonal
c     components of motion, in spherical polar system at star's
c     original position, with units of au/day
      pmr = pmra1 * 15.d0 * cdc / ( paralx * 36525.d0 )
      pmd = pmdec1              / ( paralx * 36525.d0 )
      rvl = rv1 * 86400.d0 / kmau
c
c     transform motion vector to equatorial system
      vel1(1) = - pmr * sra - pmd * sdc * cra + rvl * cdc * cra
      vel1(2) =   pmr * cra - pmd * sdc * sra + rvl * cdc * sra
      vel1(3) =               pmd * cdc       + rvl * sdc
c
c --- update star's position vector for space motion -----------------
c     (only if it=1 or it=3)
c
      if ( it .eq. 1 .or. it .eq. 3 ) then
          do 22 j=1,3
              pos2(j) = pos1(j) + vel1(j) * ( tjd2 - tjd1 )
              vel2(j) = vel1(j)
   22     continue
      else
          do 24 j=1,3
              pos2(j) = pos1(j)
              vel2(j) = vel1(j)
   24     continue
      end if
c
c --- precess position and velocity vectors --------------------------
c     (only if it=2 or it=3)
c
      if ( it .eq. 2 .or. it .eq. 3 ) then
          do 32 j=1,3
              pos1(j) = pos2(j)
              vel1(j) = vel2(j)
   32     continue
          call preces ( tjd1, pos1, tjd2,    pos2 )
          call preces ( tjd1, vel1, tjd2,    vel2 )
      end if
c
c --- convert vectors back to angular components for output ----------
c
c     from updated position vector, obtain star's new position
c     expressed as angular quantities
      xyproj = dsqrt ( pos2(1)**2 + pos2(2)**2 )
      r = 0.d0
      if ( xyproj .gt. 0.d0 ) r = datan2 ( pos2(2), pos2(1) )
      ra2 = r * seccon / 54000.d0
      if ( ra2 .lt.  0.d0 ) ra2 = ra2 + 24.d0
      if ( ra2 .ge. 24.d0 ) ra2 = ra2 - 24.d0
      d = datan2 ( pos2(3), xyproj  )
      dec2 = d * seccon / 3600.d0
      dist = dsqrt ( pos2(1)**2 + pos2(2)**2 + pos2(3)**2 )
      paralx = seccon / dist
      parx2 = paralx
c
c     transform motion vector back to spherical polar system at star's
c     new position
      cra = dcos(r)
      sra = dsin(r)
      cdc = dcos(d)
      sdc = dsin(d)
      pmr = - vel2(1) * sra       + vel2(2) * cra
      pmd = - vel2(1) * cra * sdc - vel2(2) * sra * sdc + vel2(3) * cdc
      rvl =   vel2(1) * cra * cdc + vel2(2) * sra * cdc + vel2(3) * sdc
c
c     convert components of motion to from au/day to normal
c     catalog units
      pmra2  = pmr * paralx * 36525.d0 / ( 15.d0 * cdc )
      pmdec2 = pmd * paralx * 36525.d0
      rv2    = rvl * kmau / 86400.d0
c
c     take care of zero-parallax case
      if ( parx2 .le. 1.01d-7 ) then
          parx2 = 0.d0
          rv2 = rv1
      end if
c
      return
c
      end



      subroutine zdaz (ujd,x,y,glon,glat,ht,ra,dec,irefr,
     .                 zd,az,rar,decr)
c
c     this subroutine transforms topocentric right ascension and
c     declination to zenith distance and azimuth.  this routine uses
c     a method that properly accounts for polar motion, which is
c     significant at the sub-arcsecond level.  this subroutine
c     can also adjust coordinates for atmospheric refraction.
c
c          ujd    = ut1 julian date, or equivalent greenwich apparent
c                   sidereal time in hours (in)
c          x      = conventionally-defined x coordinate of celestial
c                   ephemeris pole with respect to iers reference
c                   pole, in arcseconds (in)
c          y      = conventionally-defined y coordinate of celestial
c                   ephemeris pole with respect to iers reference
c                   pole, in arcseconds (in)
c          glon   = geodetic longitude (east +) of observer
c                   in degrees (in)
c          glat   = geodetic latitude (north +) of observer
c                   in degrees (in)
c          ht     = height of observer in meters (in)
c          ra     = topocentric right ascension of object of interest,
c                   in hours, referred to true equator and equinox
c                   of date (in)
c          dec    = topocentric declination of object of interest,
c                   in degrees, referred to true equator and equinox
c                   of date (in)
c          irefr  = atmospheric refraction option (in):
c                   set irefr=0  for no refraction
c                   set irefr=1  to include refraction
c          zd     = topocentric zenith distance in degrees,
c                   affected by refraction if irefr=1 (out)
c          az     = topocentric azimuth (measured east from north)
c                   in degrees (out)
c          rar    = topocentric right ascension of object of interest,
c                   in hours, referred to true equator and equinox
c                   of date, affected by refraction if irefr=1 (out)
c          decr   = topocentric declination of object of interest,
c                   in degrees, referred to true equator and equinox
c                   of date, affected by refraction if irefr=1 (out)
c
c     note 1:  ujd may be specified either as a ut1 julian date
c     (e.g., 2451251.823d0) or an hour and fraction of greenwich
c     apparent sidereal time (e.g., 19.1846d0).  x and y can be
c     set to zero if sub-arcsecond accuracy is not needed.
c     ht is used only for refraction, if irefr=1.  ra and dec can
c     be obtained from tpstar or tpplan.
c
c     note 2:  the directons zd=0 (zenith) and az=0 (north) are
c     here considered fixed in the terrestrial frame.  specifically,
c     the zenith is along the geodetic normal, and north is toward
c     the iers reference pole.
c
c     note 3:  if irefr=0, then rar=ra and decr=dec.
c
c
      double precision ujd,x,y,glon,glat,ht,ra,dec,zd,az,rar,decr,
     .     pi,degrad,raddeg,gast,
     .     sinlat,coslat,sinlon,coslon,sindc,cosdc,sinra,cosra,
     .     uze,une,uwe,uz,un,uw,p,pr,pz,pn,pw,proj,
     .     zd0,zd1,refr,cosr,prlen,rlen,
     .     dsin,dcos,dsqrt,datan2
      dimension uze(3), une(3), uwe(3), uz(3), un(3), uw(3),
     .     p(3), pr(3)
c
      data pi / 3.141592653589793d0 /
c
      degrad = pi / 180.d0
      raddeg = 180.d0 / pi
c
      if ( ujd .gt. 100.d0 ) then
          call sidtim ( ujd, 0.0d0, 1,   gast )
      else
          gast = dmod ( ujd, 24.0d0 )
      end if
c
      rar    = ra
      decr   = dec
      sinlat = dsin ( glat * degrad )
      coslat = dcos ( glat * degrad )
      sinlon = dsin ( glon * degrad )
      coslon = dcos ( glon * degrad )
      sindc  = dsin ( dec * degrad )
      cosdc  = dcos ( dec * degrad )
      sinra  = dsin ( ra * 15.0d0 * degrad )
      cosra  = dcos ( ra * 15.0d0 * degrad )
c
c --- set up orthonormal basis vectors in local earth-fixed system ----
c
c     define vector toward local zenith in earth-fixed system (z axis)
      uze(1) =  coslat * coslon
      uze(2) =  coslat * sinlon
      uze(3) =  sinlat
c
c     define vector toward local north in earth-fixed system (x axis)
      une(1) = -sinlat * coslon
      une(2) = -sinlat * sinlon
      une(3) =  coslat
c
c     define vector toward local west in earth-fixed system (y axis)
      uwe(1) =  sinlon
      uwe(2) = -coslon
      uwe(3) =  0.d0
c
c --- obtain vectors in celestial system ------------------------------
c
c     rotate earth-fixed orthonormal basis vectors to celestial system
c     (wrt equator and equinox of date)
      call pnsw ( 0.d0, gast, x, y, uze,   uz )
      call pnsw ( 0.d0, gast, x, y, une,   un )
      call pnsw ( 0.d0, gast, x, y, uwe,   uw )
c
c     define unit vector p toward object in celestial system
c     (wrt equator and equinox of date)
      p(1) = cosdc * cosra
      p(2) = cosdc * sinra
      p(3) = sindc
c
c --- compute coordinates of object wrt orthonormal basis -------------
c
c     compute components of p -- projections of p onto rotated
c     earth-fixed basis vectors
      pz = p(1) * uz(1) + p(2) * uz(2) + p(3) * uz(3)
      pn = p(1) * un(1) + p(2) * un(2) + p(3) * un(3)
      pw = p(1) * uw(1) + p(2) * uw(2) + p(3) * uw(3)
c
c     compute azimuth and zenith distance
      proj = dsqrt ( pn**2 + pw**2 )
      az = 0.d0
      if ( proj .gt. 0.d0 ) az = -datan2 ( pw, pn ) * raddeg
      if ( az .lt.   0.d0 ) az = az + 360.d0
      if ( az .ge. 360.d0 ) az = az - 360.d0
      zd = datan2 ( proj, pz ) * raddeg
c
c --- apply atmospheric refraction if requested -----------------------
c
      if ( irefr .eq. 1 ) then
c
c         get refraction in zenith distance
c         iterative process required because refraction algorithms are
c         always a function of observed (not computed) zenith distance
          zd0 = zd
  40      zd1 = zd
          call refrac ( ht, zd,   refr )
          zd = zd0 - refr
c         require convergence to 0.2 arcsec (actual accuracy less)
          if ( dabs ( zd - zd1 ) .gt. 5.d-5 ) go to 40
c
c         apply refraction to celestial coordinates of object
          if ( refr .gt. 0.d0 .and. zd .gt. 0.01d0 ) then
c
c             shift position vector of object in celestial system
c             to account for for refraction (see usno/aa technical
c             note 9)
              cosr  = dcos ( refr * degrad )
              prlen = dsin ( zd0 * degrad ) / dsin ( zd * degrad )
              rlen  = dsqrt ( 1.d0 + prlen**2 - 2.d0 * prlen * cosr )
c             add small refraction displacement vector to p
              do 50 j = 1, 3
  50          pr(j) = ( p(j) + rlen * uz(j) ) / prlen
c
c             compute refracted right ascension and declination
              proj = dsqrt ( pr(1)**2 + pr(2)**2 )
              rar = 0.d0
              if ( proj .gt. 0.d0 ) rar = datan2 ( pr(2), pr(1) )
     .                                    * raddeg / 15.d0
              if ( rar .lt.  0.d0 ) rar = rar + 24.d0
              if ( rar .ge. 24.d0 ) rar = rar - 24.d0
              decr = datan2 ( pr(3), proj ) * raddeg
c
          end if
c
      end if
c
c ---------------------------------------------------------------------
c
      return
c
      end



      subroutine vectrs (ra,dec,pmra,pmdec,parllx,rv,pos,vel)
c
c     this subroutine converts angular quantities to vectors.
c
c          ra     = right ascension in hours (in)
c          dec    = declination in degrees (in)
c          pmra   = proper motion in ra in time seconds per
c                   julian century (in)
c          pmdec  = proper motion in dec in arcseconds per
c                   julian century (in)
c          parllx = parallax in arcseconds (in)
c          rv     = radial velocity in kilometers per second (in)
c          pos    = position vector, equatorial rectangular coordinates,
c                   components in au (out)
c          vel    = velocity vector, equatorial rectangular coordinates,
c                   components in au/day (out)
c
c
      double precision ra,dec,pmra,pmdec,parllx,rv,pos,vel,
     .     seccon,kmau,paralx,dist,r,d,cra,sra,cdc,sdc,pmr,pmd,rvl,
     .     dcos,dsin
      dimension pos(3), vel(3)
c
      data seccon / 206264.8062470964d0 /,     kmau / 1.49597870d8 /
c
c     if parallax is unknown, undetermined, or zero, set it to 1e-7
c     arcsecond, corresponding to a distance of 10 megaparsecs
      paralx = parllx
      if (paralx.le.0.0d0) paralx = 1.0d-7
c
c     convert right ascension, declination, and parallax to position
c     vector in equatorial system with units of au
      dist = seccon / paralx
      r = ra * 54000.0d0 / seccon
      d = dec * 3600.0d0 / seccon
      cra = dcos(r)
      sra = dsin(r)
      cdc = dcos(d)
      sdc = dsin(d)
      pos(1) = dist * cdc * cra
      pos(2) = dist * cdc * sra
      pos(3) = dist * sdc
c
c     convert proper motion and radial velocity to orthogonal components
c     of motion with units of au/day
      pmr = pmra * 15.0d0 * cdc / (paralx * 36525.0d0)
      pmd = pmdec / (paralx * 36525.0d0)
      rvl = rv * 86400.0d0 / kmau
c
c     transform motion vector to equatorial system
      vel(1) = - pmr * sra   - pmd * sdc * cra   + rvl * cdc * cra
      vel(2) =   pmr * cra   - pmd * sdc * sra   + rvl * cdc * sra
      vel(3) =                 pmd * cdc         + rvl * sdc
c
      return
c
      end



      subroutine angles (pos,ra,dec)
c
c     this subroutine converts a vector to angular quantities.
c
c          pos = position vector, equatorial rectangular
c                coordinates (in)
c          ra  = right ascension in hours (out)
c          dec = declination in degrees (out)
c
c
      double precision pos,ra,dec,seccon,xyproj,r,d,dsqrt,datan2
      dimension pos(3)
c
      data seccon / 206264.8062470964d0 /
c
      xyproj = dsqrt(pos(1)**2 + pos(2)**2)
      r = 0.d0
      if (xyproj.gt.0.d0) r = datan2(pos(2),pos(1))
      ra = r * seccon / 54000.0d0
      if (ra.lt. 0.0d0) ra = ra + 24.0d0
      if (ra.ge.24.0d0) ra = ra - 24.0d0
      d = datan2(pos(3),xyproj)
      dec = d * seccon / 3600.0d0
      return
c
      end



      subroutine propmo (tjd1,pos1,vel1,tjd2,pos2)
c
c     this subroutine applies proper motion, including foreshortening
c     effects, to a star's position.
c
c          tjd1 = tdb julian date of first epoch (in)
c          pos1 = position vector at first epoch (in)
c          vel1 = velocity vector at first epoch (in)
c          tjd2 = tdb julian date of second epoch (in)
c          pos2 = position vector at second epoch (out)
c
c
      double precision tjd1,pos1,vel1,tjd2,pos2
      dimension pos1(3), vel1(3), pos2(3)
c
      do 20 j=1,3
   20 pos2(j) = pos1(j) + vel1(j) * (tjd2 - tjd1)
      return
c
      end



      subroutine geocen (pos1,pe,pos2,tlight)
c
c     this subroutine moves the origin of coordinates from the
c     barycenter of the solar system to the center of mass of the
c     earth, i.e., this subroutine corrects for parallax.
c
c          pos1   = position vector, referred to origin at solar system
c                   barycenter, components in au (in)
c          pe     = position vector of center of mass of the earth,
c                   referred to origin at solar system barycenter,
c                   components in au (in)
c          pos2   = position vector, referred to origin at center of
c                   mass of the earth, components in au (out)
c          tlight = light time from body to earth in days (out)
c
c
      double precision pos1,pe,pos2,tlight,c,dsqrt
      dimension pos1(3), pe(3), pos2(3)
c
      data c / 173.14463348d0 /
c     c = speed of light in au/day
c
      do 20 j=1,3
   20 pos2(j) = pos1(j) - pe(j)
      tlight = dsqrt(pos2(1)**2 + pos2(2)**2 + pos2(3)**2) / c
      return
c
      end



      subroutine sunfld (pos1,pe,pos2)
c
c     subroutine sunfld version 1.
c     this subroutine corrects position vector for the deflection
c     of light in the gravitational field of the sun.  see misner,
c     thorne, and wheeler (1973), gravitation, pp. 184-185.  this
c     subroutine valid for bodies within the solar system as well as
c     for stars.
c
c          pos1 = position vector, referred to origin at center of mass
c                 of the earth, components in au (in)
c          pe   = position vector of center of mass of the earth,
c                 referred to origin at center of mass of
c                 the sun, components in au (in)
c          pos2 = position vector, referred to origin at center of mass
c                 of the earth, corrected for gravitational deflec-
c                 tion, components in au (out)
c
c
      double precision pos1,pe,pos2,p1hat,pehat,p1mag,pemag,mau,gs,c,f,
     .     cosd,sind,b,bm,pqmag,zfinl,zinit,xifinl,xiinit,
     .     delphi,delphp,delp,dabs,dsqrt
      dimension pos1(3), pe(3), pos2(3), p1hat(3), pehat(3)
c
      data mau / 1.49597870d11 /
c     mau = number of meters per au
      data gs / 1.32712438d20 /
c     gs = heliocentric gravitational constant
      data c / 299792458.0d0 /
c     c = speed of light
c
      f = 0.0d0
c
c     compute vector magnitudes and unit vectors
      p1mag = dsqrt (pos1(1)**2 + pos1(2)**2 + pos1(3)**2)
      pemag = dsqrt (  pe(1)**2 +   pe(2)**2 +   pe(3)**2)
      do 20 j=1,3
      p1hat(j) = pos1(j) / p1mag
   20 pehat(j) =   pe(j) / pemag
c
c     compute geometrical quantities
c     cosd and sind are cosine and sine of d, the angular separation
c     of the body from the sun as viewed from the earth
      cosd = - pehat(1)*p1hat(1) - pehat(2)*p1hat(2) - pehat(3)*p1hat(3)
      if (dabs(cosd).gt.0.9999999999d0) go to 40
      sind = dsqrt (1.0d0 - cosd**2)
c     b is the impact parameter for the ray
      b = pemag * sind
      bm = b * mau
c     pqmag is the distance of the body from the sun
      pqmag = dsqrt (p1mag**2 + pemag**2 - 2.0d0 * p1mag * pemag * cosd)
c
c     compute delphi, the angle of deflection of the ray
      zfinl = pemag * cosd
      zinit = -p1mag + zfinl
      xifinl = zfinl / b
      xiinit = zinit / b
      delphi = 2.0d0*gs/(bm*c*c) * (xifinl / dsqrt (1.0d0 + xifinl**2)
     .                            - xiinit / dsqrt (1.0d0 + xiinit**2))
c
c     compute delphp, the change in angle as seen at the earth
      delphp = delphi / (1.0d0 + (pemag / pqmag))
c
c     fix up position vector
c     pos2 is pos1 rotated through angle delphp in plane defined
c     by pos1 and pe
      f = delphp * p1mag / sind
   40 do 50 j=1,3
      delp = f * (cosd * p1hat(j) + pehat(j))
   50 pos2(j) = pos1(j) + delp
c
      return
c
      end



      subroutine aberat (pos1,ve,tlight,pos2)
c
c     this subroutine corrects position vector for aberration of light.
c     algorithm includes relativistic terms.  see murray (1981)
c     mon. notices royal ast. society 195, 639-648.
c
c          pos1   = position vector, referred to origin at center of
c                   mass of the earth, components in au (in)
c          ve     = velocity vector of center of mass of the earth,
c                   referred to origin at solar system barycenter,
c                   components in au/day (in)
c          tlight = light time from body to earth in days (in)
c                   if tlight = 0.0d0, this subroutine will compute
c          pos2   = position vector, referred to origin at center of
c                   mass of the earth, corrected for aberration,
c                   components in au (out)
c
c
      double precision pos1,ve,tlight,pos2,c,tl,p1mag,vemag,
     .     beta,dot,cosd,gammai,p,q,r,dsqrt
      dimension pos1(3), ve(3), pos2(3)
c
      data c / 173.14463348d0 /
c     c = speed of light in au/day
c
      tl = tlight
      p1mag = tl * c
      if (tl.ne.0.0d0) go to 20
      p1mag = dsqrt(pos1(1)**2 + pos1(2)**2 + pos1(3)**2)
      tl = p1mag / c
   20 vemag = dsqrt(ve(1)**2 + ve(2)**2 + ve(3)**2)
      beta = vemag / c
      dot = pos1(1)*ve(1) + pos1(2)*ve(2) + pos1(3)*ve(3)
      cosd = dot / (p1mag * vemag)
      gammai = dsqrt(1.0d0 - beta**2)
      p = beta * cosd
      q = (1.0d0 + p / (1.0d0 + gammai)) * tl
      r = 1.0d0 + p
c
      do 30 j=1,3
   30 pos2(j) = (gammai * pos1(j) + q * ve(j)) / r
      return
c
      end



      subroutine preces (tjd1,pos1,tjd2,pos2)
c
c     this subroutine precesses equatorial rectangular coordinates from
c     one epoch to another.  the coordinates are referred to the mean
c     equator and equinox of the two respective epochs.  see
c     explanatory supplement to the astronomical almanac, pp. 103-104,
c     lieske, et al. (1977) astronomy and astrophysics 58, 1-16, and
c     lieske (1979) astronomy and astrophysics 73, 282-284.
c
c          tjd1 = tdb julian date of first epoch (in)
c          pos1 = position vector, geocentric equatorial rectangular
c                 coordinates, referred to mean equator and equinox of
c                 first epoch (in)
c          tjd2 = tdb julian date of second epoch (in)
c          pos2 = position vector, geocentric equatorial rectangular
c                 coordinates, referred to mean equator and equinox of
c                 second epoch (out)
c
c
      double precision tjd1,tjd2,t0,t,t02,t2,t3,pos1,pos2,seccon,
     .     zeta0,zee,theta,czeta0,szeta0,czee,szee,ctheta,stheta,
     .     xx,yx,zx,xy,yy,zy,xz,yz,zz,t1last,t2last,dabs,dcos,dsin
      dimension pos1(3), pos2(3)
      save
c
      data seccon / 206264.8062470964d0 /
      data t1last,t2last / 0.0d0,0.0d0 /
c
      if (dabs(tjd1-t1last).lt.1.0d-8.and.dabs(tjd2-t2last).lt.1.0d-8)
     .     go to 20
      if (dabs(tjd1-t2last).lt.1.0d-8.and.dabs(tjd2-t1last).lt.1.0d-8)
     .     go to 30
c
c     t0 and t below correspond to lieske's big t and little t
c     time scale is assumed to be tdb
      t0 = (tjd1 - 2451545.00000000d0) / 36525.0d0
      t = (tjd2 - tjd1) / 36525.0d0
      t02 = t0 * t0
      t2 = t * t
      t3 = t2 * t
c     zeta0, zee, and theta below correspond to lieske's zeta-sub-a,
c     z-sub-a, and theta-sub-a
      zeta0 = (2306.2181d0 + 1.39656d0*t0 - 0.000139d0*t02) * t
     .      + (0.30188d0 - 0.000344d0*t0) * t2
     .      +  0.017998d0 * t3
      zee   = (2306.2181d0 + 1.39656d0*t0 - 0.000139d0*t02) * t
     .      + (1.09468d0 + 0.000066d0*t0) * t2
     .      +  0.018203d0 * t3
      theta = (2004.3109d0 - 0.85330d0*t0 - 0.000217d0*t02) * t
     .      + (-0.42665d0 - 0.000217d0*t0) * t2
     .      -  0.041833d0 * t3
      zeta0 = zeta0 / seccon
      zee = zee / seccon
      theta = theta / seccon
      czeta0 = dcos(zeta0)
      szeta0 = dsin(zeta0)
      czee = dcos(zee)
      szee = dsin(zee)
      ctheta = dcos(theta)
      stheta = dsin(theta)
c
c     precession rotation matrix follows
      xx = czeta0*ctheta*czee - szeta0*szee
      yx = -szeta0*ctheta*czee - czeta0*szee
      zx = -stheta*czee
      xy = czeta0*ctheta*szee + szeta0*czee
      yy = -szeta0*ctheta*szee + czeta0*czee
      zy = -stheta*szee
      xz = czeta0*stheta
      yz = -szeta0*stheta
      zz = ctheta
      t1last = tjd1
      t2last = tjd2
c
c     perform rotation
   20 pos2(1) = xx*pos1(1) + yx*pos1(2) + zx*pos1(3)
      pos2(2) = xy*pos1(1) + yy*pos1(2) + zy*pos1(3)
      pos2(3) = xz*pos1(1) + yz*pos1(2) + zz*pos1(3)
      go to 50
c
c     perform inverse rotation
   30 pos2(1) = xx*pos1(1) + xy*pos1(2) + xz*pos1(3)
      pos2(2) = yx*pos1(1) + yy*pos1(2) + yz*pos1(3)
      pos2(3) = zx*pos1(1) + zy*pos1(2) + zz*pos1(3)
c
   50 return
c
      end



      subroutine nutate (tjd,pos1,pos2)
c
c     this subroutine nutates equatorial rectangular coordinates from
c     mean equator and equinox of epoch to true equator and equinox of
c     epoch.  see explanatory supplement to the astronomical almanac,
c     pp. 114-115.
c
c          tjd    = tdb julian date of epoch (in)
c          pos1   = position vector, geocentric equatorial rectangular
c                   coordinates, referred to mean equator and equinox
c                   of epoch (in)
c          pos2   = position vector, geocentric equatorial rectangular
c                   coordinates, referred to true equator and equinox
c                   of epoch (out)
c
c     note:  if tjd is negative, inverse nutation (true to mean)
c     is applied.
c
c
      double precision tjd,pos1,pos2,tjd1,seccon,oblm,oblt,eqeq,
     .     dpsi,deps,cobm,sobm,cobt,sobt,cpsi,spsi,
     .     xx,yx,zx,xy,yy,zy,xz,yz,zz,dabs,dcos,dsin
      dimension pos1(3), pos2(3)
c
      data seccon / 206264.8062470964d0 /
c
      tjd1 = dabs(tjd)
c
      call etilt (tjd1,oblm,oblt,eqeq,dpsi,deps)
      oblm = oblm * 3600.0d0 / seccon
      oblt = oblt * 3600.0d0 / seccon
      dpsi = dpsi / seccon
      deps = deps / seccon
      cobm = dcos(oblm)
      sobm = dsin(oblm)
      cobt = dcos(oblt)
      sobt = dsin(oblt)
      cpsi = dcos(dpsi)
      spsi = dsin(dpsi)
c
c     nutation rotation matrix follows
      xx = cpsi
      yx = -spsi*cobm
      zx = -spsi*sobm
      xy = spsi*cobt
      yy = cpsi*cobm*cobt + sobm*sobt
      zy = cpsi*sobm*cobt - cobm*sobt
      xz = spsi*sobt
      yz = cpsi*cobm*sobt - sobm*cobt
      zz = cpsi*sobm*sobt + cobm*cobt
   10 if (tjd.lt.0.0d0) go to 30
c
c     perform rotation
   20 pos2(1) = xx*pos1(1) + yx*pos1(2) + zx*pos1(3)
      pos2(2) = xy*pos1(1) + yy*pos1(2) + zy*pos1(3)
      pos2(3) = xz*pos1(1) + yz*pos1(2) + zz*pos1(3)
      go to 50
c
c     perform inverse rotation
   30 pos2(1) = xx*pos1(1) + xy*pos1(2) + xz*pos1(3)
      pos2(2) = yx*pos1(1) + yy*pos1(2) + yz*pos1(3)
      pos2(3) = zx*pos1(1) + zy*pos1(2) + zz*pos1(3)
c
   50 return
c
      end



      subroutine spin (st,pos1,pos2)
c
c     this subroutine transforms geocentric rectangular coordinates
c     from rotating system based on rotational equator and orthogonal
c     reference meridian to non-rotating system based on true equator
c     and equinox of date.
c
c          st     = local apparent sidereal time at reference meridian
c                   in hours (in)
c          pos1   = vector in geocentric rectangular
c                   rotating system, referred to rotational equator
c                   and orthogonal reference meridian (in)
c          pos2   = vector in geocentric rectangular
c                   non-rotating system, referred to true equator
c                   and equinox of date (out)
c
c
      double precision st,pos1,pos2,seccon,tlast,str,cosst,sinst,
     .     xx,yx,zx,xy,yy,zy,xz,yz,zz,dabs,dcos,dsin
      dimension pos1(3), pos2(3)
      save
c
      data seccon / 206264.8062470964d0 /
      data tlast / -999.0d0 /
c
      if (dabs(st-tlast).lt.1.0d-12) go to 10
c
      str   = st * 15.0d0 * 3600.0d0 / seccon
      cosst = dcos(str)
      sinst = dsin(str)
c
c     sidereal time rotation matrix follows
      xx =  cosst
      yx = -sinst
      zx =  0.0d0
      xy =  sinst
      yy =  cosst
      zy =  0.0d0
      xz =  0.0d0
      yz =  0.0d0
      zz =  1.0d0
      tlast = st
   10 continue
c
c     perform rotation
   20 pos2(1) = xx*pos1(1) + yx*pos1(2) + zx*pos1(3)
      pos2(2) = xy*pos1(1) + yy*pos1(2) + zy*pos1(3)
      pos2(3) = xz*pos1(1) + yz*pos1(2) + zz*pos1(3)
c
   50 return
c
      end



      subroutine wobble (x,y,pos1,pos2)
c
c     this subroutine corrects earth-fixed geocentric rectangular
c     coordinates for polar motion.  it transforms a vector from
c     earth-fixed geographic system to rotating system based on
c     rotational equator and orthogonal greenwich meridian through
c     axis of rotation.
c
c          x      = conventionally-defined x coordinate of celestial
c                   ephemeris pole with respect to iers reference
c                   pole, in arcseconds (in)
c          y      = conventionally-defined y coordinate of celestial
c                   ephemeris pole with respect to iers reference
c                   pole, in arcseconds (in)
c          pos1   = vector in geocentric rectangular
c                   earth-fixed system, referred to geographic
c                   equator and greenwich meridian (in)
c          pos2   = vector in geocentric rectangular
c                   rotating system, referred to rotational equator
c                   and orthogonal greenwich meridian (out)
c
c
      double precision x,y,pos1,pos2,seccon,xpole,ypole,
     .     xx,yx,zx,xy,yy,zy,xz,yz,zz
      dimension pos1(3), pos2(3)
c
      data seccon / 206264.8062470964d0 /
c
      xpole = x / seccon
      ypole = y / seccon
c
c     wobble rotation matrix follows
      xx =  1.0d0
      yx =  0.0d0
      zx = -xpole
      xy =  0.0d0
      yy =  1.0d0
      zy =  ypole
      xz =  xpole
      yz = -ypole
      zz =  1.0d0
   10 continue
c
c     perform rotation
   20 pos2(1) = xx*pos1(1) + yx*pos1(2) + zx*pos1(3)
      pos2(2) = xy*pos1(1) + yy*pos1(2) + zy*pos1(3)
      pos2(3) = xz*pos1(1) + yz*pos1(2) + zz*pos1(3)
c
   50 return
c
      end



      subroutine terra (glon,glat,ht,st,pos,vel)
c
c     this subroutine computes the position and velocity vectors of
c     a terrestrial observer with respect to the center of the earth.
c
c          glon   = longitude of observer with respect to reference
c                   meridian (east +) in degrees (in)
c          glat   = geodetic latitude (north +) of observer
c                   in degrees (in)
c          ht     = height of observer in meters (in)
c          st     = local apparent sidereal time at reference meridian
c                   in hours (in)
c          pos    = position vector of observer with respect to center
c                   of earth, equatorial rectangular coordinates,
c                   referred to true equator and equinox of date,
c                   components in au (out)
c          vel    = velocity vector of observer with respect to center
c                   of earth, equatorial rectangular coordinates,
c                   referred to true equator and equinox of date,
c                   components in au/day (out)
c
c     note:  if reference meridian is greenwich and st=0.d0, pos
c     is effectively referred to equator and greenwich.
c
c
      double precision glon,glat,ht,st,pos,vel,seccon,erad,f,omega,
     .     kmau,df2,phi,sinphi,cosphi,c,s,ach,ash,stlocl,sinst,cosst,
     .     dsqrt,dcos,dsin
      dimension pos(3), vel(3)
c
      data seccon / 206264.8062470964d0 /
c
      data erad / 6378.140d0 /,   f / 0.00335281d0 /
c     erad = radius of earth in km, f = earth ellipsoid flattening
      data omega / 7.292115d-5 /
c     omega = rotational angular velocity of earth in radians/sec
      data kmau / 1.49597870d8 /
c     kmau = kilometers per astronomical unit
c
c     compute parameters relating to geodetic to geocentric conversion
      df2 = (1.0d0 - f)**2
      phi = glat * 3600.0d0 / seccon
      sinphi = dsin(phi)
      cosphi = dcos(phi)
      c = 1.0d0 / dsqrt ( cosphi**2 + df2 * sinphi**2 )
      s = df2 * c
      ach = erad * c + ht/1000.0d0
      ash = erad * s + ht/1000.0d0
c
c     compute local sidereal time factors
      stlocl = (st * 54000.0d0 + glon * 3600.0d0) / seccon
      sinst = dsin(stlocl)
      cosst = dcos(stlocl)
c
c     compute position vector components in km
      pos(1) = ach * cosphi * cosst
      pos(2) = ach * cosphi * sinst
      pos(3) = ash * sinphi
c
c     compute velocity vector components in km/sec
      vel(1) = -omega * ach * cosphi * sinst
      vel(2) =  omega * ach * cosphi * cosst
      vel(3) =  0.0d0
c
c     convert position and velocity components to au and au/day
      do 20 j=1,3
      pos(j) = pos(j) / kmau
      vel(j) = vel(j) / kmau * 86400.0d0
   20 continue
c
      return
c
      end



      subroutine novas_times(tdbjd,ttjd,secdif)
c
c     this subroutine computes the terrestrial time (tt) julian date
c     corresponding to a barycentric dynamical time (tdb) julian date.
c     expressions used in this version are approximations resulting
c     in accuracies of about 20 microseconds.  see explanatory
c     supplement to the astronomical almanac, pp. 42-44 and 316.
c
c          tdbjd  = tdb julian date (in)
c          ttjd   = tt julian date (out)
c          secdif = difference tdbjd-ttjd, in seconds (out)
c
c
      double precision tdbjd,ttjd,secdif,seccon,rev,t0,ecc,
     .     tdays,m,l,lj,e,dsin
c
      data seccon / 206264.8062470964d0 /,   rev / 1296000.d0 /
      data t0 / 2451545.00000000d0 /
c     t0 = tdb julian date of epoch j2000.0
      data ecc / 0.01671022d0 /
c     ecc = eccentricity of earth-moon barycenter orbit
c
      tdays = tdbjd - t0
      m  = ( 357.51716d0 + 0.985599987d0 * tdays ) * 3600.d0
      l  = ( 280.46435d0 + 0.985609100d0 * tdays ) * 3600.d0
      lj = (  34.40438d0 + 0.083086762d0 * tdays ) * 3600.d0
      m  = dmod (  m, rev ) / seccon
      l  = dmod (  l, rev ) / seccon
      lj = dmod ( lj, rev ) / seccon
      e  = m + ecc * dsin ( m ) + 0.5d0 * ecc**2 * dsin ( 2.d0 * m )
      secdif =   1.658d-3 * dsin ( e )
     .         + 20.73d-6 * dsin ( l - lj )
      ttjd = tdbjd - secdif / 86400.d0
c
      return
c
      end



      subroutine etilt (tjd,oblm,oblt,eqeq,dpsi,deps)
c
c     this subroutine computes quantities related to the orientation
c     of the earth's rotation axis at julian date tjd.
c     implements equation of the equinoxes definition as per
c     iau resolution c7 of 1994.
c
c          tjd    = tdb julian date for orientation parameters (in)
c          oblm   = mean obliquity of the ecliptic in degrees at
c                   date tjd (out)
c          oblt   = true obliquity of the ecliptic in degrees at
c                   date tjd (out)
c          eqeq   = equation of the equinoxes in time seconds at
c                   date tjd (out)
c          dpsi   = nutation in longitude in arcseconds at
c                   date tjd (out)
c          deps   = nutation in obliquity in arcseconds at
c                   date tjd (out)
c
c
      double precision tjd,t0,t,t2,t3,tlast,oblm,oblt,eqeq,dpsi,deps,
     .     seccon,obm,obt,ee,delpsi,deleps,psi,eps,x,omega,
     .     ddpsi,ddeps,psicor,epscor,dabs,dcos
      save
c
      data t0 / 2451545.00000000d0 /
c     t0 = tdb julian date of epoch j2000.0
      data seccon / 206264.8062470964d0 /
      data tlast / 0.0d0 /,   psicor,epscor / 0.0d0,0.0d0 /
c
      t = (tjd - t0) / 36525.0d0
      t2 = t * t
      t3 = t2 * t
c
c     obtain nutation parameters in arcseconds
      if (dabs(tjd-tlast).lt.1.0d-8) go to 20
      call funarg (t,x,x,x,x,omega)
      call nod (t,delpsi,deleps)
      tlast = tjd
   20 psi = delpsi + psicor
      eps = deleps + epscor
c
c     compute mean obliquity of the ecliptic in arcseconds
      obm = 84381.4480d0 - 46.8150d0*t - 0.00059d0*t2
     .    + 0.001813d0*t3
c
c     compute true obliquity of the ecliptic in arcseconds
      obt = obm + eps
c
c     compute equation of the equinoxes in arcseconds, time seconds
c     (iau 1994 and iers 1996 definition)
      ee = psi * dcos (obm/seccon)
     .     + 0.00264d0 * dsin(omega) + 0.000063d0 * dsin(2.d0*omega)
      ee = ee / 15.d0
c
c     convert obliquity values to degrees
      obm = obm / 3600.0d0
      obt = obt / 3600.0d0
c
      oblm = obm
      oblt = obt
      eqeq = ee
      dpsi = psi
      deps = eps
c
      return
c
c
      entry celpol (ddpsi,ddeps)
c
c     this entry allows for the specification of celestial pole
c     offsets for high-precision applications.  these are added
c     to the nutation parameters delta psi and delta epsilon.
c     daily values of the offsets are published, for example,
c     in iers bulletins a and b.  this entry, if used, should
c     be called before any other routines for a given date.
c     values of the pole offsets specified via a call to this
c     entry will be used until explicitly changed.
c
c          ddpsi  = value of offset in delta psi (dpsi)
c                   in arcseconds (in)
c          ddeps  = value of offset in delta epsilon (deps)
c                   in arcseconds (in)
c
c
      psicor = ddpsi
      epscor = ddeps
      return
c
      end



      subroutine funarg (t,el,elprim,f,d,omega)
c
c     this subroutine computes fundamental arguments (mean elements)
c     of the sun and moon.  see seidelmann (1982) celestial
c     mechanics 27, 79-106 (1980 iau theory of nutation).
c
c          t      = tdb time in julian centuries since j2000.0 (in)
c          el     = mean anomaly of the moon in radians
c                   at date tjd (out)
c          elprim = mean anomaly of the sun in radians
c                   at date tjd (out)
c          f      = mean longitude of the moon minus mean longitude
c                   of the moon's ascending node in radians
c                   at date tjd (out)
c          d      = mean elongation of the moon from the sun in
c                   radians at date tjd (out)
c          omega  = mean longitude of the moon's ascending node
c                   in radians at date tjd (out)
c
c
      double precision t,tlast,el,elprim,f,d,omega,arg,seccon,rev,dmod
      dimension arg(5)
      save
c
      data seccon / 206264.8062470964d0 /,   rev / 1296000.d0 /
      data tlast / 0.0d0 /
c
      if (dabs(t-tlast).lt.1.0d-12) go to 40
c
c     compute fundamental arguments in arcseconds
c
      arg(1) = ((+0.064d0 * t + 31.310d0) * t + 715922.633d0) * t
     .     + 485866.733d0 + dmod(1325.0d0*t,1.0d0) * rev
      arg(1) = dmod(arg(1),rev)
c
      arg(2) = ((-0.012d0 * t - 0.577d0) * t + 1292581.224d0) * t
     .     + 1287099.804d00 + dmod(99.0d0*t,1.0d0) * rev
      arg(2) = dmod(arg(2),rev)
c
      arg(3) = ((+0.011d0 * t - 13.257d0) * t + 295263.137d0) * t
     .     + 335778.877d0 + dmod(1342.0d0*t,1.0d0) * rev
      arg(3) = dmod(arg(3),rev)
c
      arg(4) = ((+0.019d0 * t - 6.891d0) * t + 1105601.328d0) * t
     .     + 1072261.307d0 + dmod(1236.0d0*t,1.0d0) * rev
      arg(4) = dmod(arg(4),rev)
c
      arg(5) = ((0.008d0 * t + 7.455d0) * t - 482890.539d0) * t
     .     + 450160.280d0  - dmod(5.0d0*t,1.0d0) * rev
      arg(5) = dmod(arg(5),rev)
c
c     convert arguments to radians
      do 30 i=1,5
      arg(i) = dmod(arg(i),rev)
      if (arg(i).lt.0.d0) arg(i) = arg(i) + rev
      arg(i) = arg(i) / seccon
   30 continue
      tlast = t
c
   40 el     = arg(1)
      elprim = arg(2)
      f      = arg(3)
      d      = arg(4)
      omega  = arg(5)
c
      return
c
      end



      subroutine refrac (height,zdobs,refr)
c
c     this subroutine computes atmospheric refraction in zenith
c     distance.  this version computes approximate refraction for
c     optical wavelengths.  it can be used for planning observations
c     or telescope pointing, but should not be used for the reduction
c     of precise observations.  basic algorithm is described in the
c     explanatory supplement to the astronomical almanac, p. 144,
c     and is an adaptation of a formula in bennett (1982), journal
c     of navigation (royal institute) 35, 255-259.
c
c          height = height of observer in meters (in)
c          zdobs  = observed zenith distance in degrees (in)
c          refr   = atmospheric refraction in degrees (out)
c
c     note:  height is not used if entry refdat has been called
c     to specify atmospheric pressure.
c
c
      double precision height,zdobs,refr,pi,s,degrad,
     .     pobs,tobs,dobs,wlobs,obsp,obst,obsd,obswl,p,t,d,wl,h,r,
     .     dexp,dtan
      save
c
      data pi / 3.141592653589793d0 /
      data pobs, tobs, dobs, wlobs / 4 * -999.d0 /
      data s / 9.1d3 /
c     s is approximate scale height of atmosphere in meters
c
      degrad = pi / 180.d0
c
c     compute refraction only for zenith distances
c     between 0.1 and 91 degrees
      if ( zdobs .lt. 0.1d0 .or. zdobs .gt. 91.d0 ) then
          refr = 0.d0
          go to 77
      end if
c
c     if observed weather data are available, use them
c     otherwise, use crude estimates of average conditions
      if ( pobs .ge. 1.d0 .and. tobs .gt. -100.d0 ) then
          p  = pobs
          t  = tobs
          d  = dobs
          wl = wlobs
      else
          p  = 1010.d0 * dexp ( -height / s )
          t  = 10.d0
          d  =  0.d0
          wl =  0.5d0
      end if
c     d and wl not used in this version
c
      h = 90.d0 - zdobs
      r = 0.016667d0 / dtan ( ( h +  7.31d0 / ( h + 4.4d0 ) ) * degrad )
      refr = r * ( 0.28d0 * p / ( t + 273.d0 ) )
c
  77  return
c
c
      entry refdat (obsp,obst,obsd,obswl)
c
c     this entry allows for the specification of weather observations
c     and other data to be used in the atmospheric refraction
c     calculation.  this entry, if used, should be called before
c     subroutine refrac or zdaz for a given date/time.  data specified
c     via a call to this entry will be used until explicitly changed.
c
c          obsp   = observed atmospheric pressure in millibars (in)
c          obst   = observed temperature in degrees celsius (in)
c          obsd   = observed dew point in degrees celsius (in)
c          obswl  = observing wavelength in microns (in)
c
c     note:  obsd and obswl are not used in this version's refraction
c     algorithm, and can be set to any value.
c
c
      pobs  = obsp
      tobs  = obst
      dobs  = obsd
      wlobs = obswl
      return
c
      end



      subroutine juldat (i,m,k,h,tjd)
c
c     this subroutine computes julian date, given calendar date and
c     time.  input calendar date must be gregorian.  input time value
c     can be in any ut-like time scale (utc, ut1, tt, etc.) - output
c     julian date will have same basis.  algorithm by fliegel and
c     van flandern.
c
c          i      = year (in)
c          m      = month number (in)
c          k      = day of month (in)
c          h      = ut hours (in)
c          tjd    = julian date (out)
c
c
      double precision h,tjd
c
c     jd=julian day no for day beginning at greenwich noon on given date
      jd = k-32075+1461*(i+4800+(m-14)/12)/4+367*(m-2-(m-14)/12*12)/12
     .     -3*((i+4900+(m-14)/12)/100)/4
      tjd = jd - 0.5d0 + h/24.d0
c
      return
      end



      subroutine caldat (tjd,i,m,k,h)
c
c     this subroutine computes calendar date and time, given julian
c     date.  input julian date can be based on any ut-like time scale
c     (utc, ut1, tt, etc.) - output time value will have same basis.
c     output calendar date will be gregorian.  algorithm by fliegel and
c     van flandern.
c
c          tjd    = julian date (in)
c          i      = year (out)
c          m      = month number (out)
c          k      = day of month (out)
c          h      = ut hours (out)
c
c
      double precision tjd,h,djd,dmod
c
      djd = tjd + 0.5d0
      jd = djd
      h = dmod (djd,1.d0) * 24.d0
c     jd=julian day no for day beginning at greenwich noon on given date
      l = jd + 68569
      n = 4*l/146097
      l = l - (146097*n+3)/4
c     i=year, m=month, k=day
      i = 4000*(l+1)/1461001
      l = l - 1461*i/4 + 31
      m = 80*l/2447
      k = l - 2447*m/80
      l = m / 11
      m = m + 2 - 12*l
      i = 100*(n-49) + i + l
c
      return
      end



      subroutine nod (t,dpsi,deps)
c
c     subroutine nod version 1.
c     this subroutine evaluates the nutation series and returns the
c     values for nutation in longitude and nutation in obliquity.
c     wahr nutation series for axis b for gilbert & dziewonski earth
c     model 1066a.  see seidelmann (1982) celestial mechanics 27,
c     79-106.  1980 iau theory of nutation.
c
c          t    = tdb time in julian centuries since j2000.0 (in)
c          dpsi = nutation in longitude in arcseconds (out)
c          deps = nutation in obliquity in arcseconds (out)
c
c
      double precision t,dpsi,deps,l,lp,f,d,om,arg,dble,dsin,dcos
      dimension x(9,106),x1(90),x2(90),x3(90),x4(90),x5(90),x6(90),
     .     x7(90),x8(90),x9(90),xa(90),xb(54)
      equivalence(x(1,  1),x1(1))
      equivalence(x(1, 11),x2(1))
      equivalence(x(1, 21),x3(1))
      equivalence(x(1, 31),x4(1))
      equivalence(x(1, 41),x5(1))
      equivalence(x(1, 51),x6(1))
      equivalence(x(1, 61),x7(1))
      equivalence(x(1, 71),x8(1))
      equivalence(x(1, 81),x9(1))
      equivalence(x(1, 91),xa(1))
      equivalence(x(1,101),xb(1))
c
c***********************************************************************
c
c
c     table of multiples of arguments and coefficients
c
c                   multiple of            longitude        obliquity
c              l    l'   f    d  omega   coeff. of sin    coeff. of cos
      data x1/ 0.,  0.,  0.,  0.,  1., -171996., -174.2,  92025.,  8.9,
     .         0.,  0.,  2., -2.,  2.,  -13187.,   -1.6,   5736., -3.1,
     .         0.,  0.,  2.,  0.,  2.,   -2274.,   -0.2,    977., -0.5,
     .         0.,  0.,  0.,  0.,  2.,    2062.,    0.2,   -895.,  0.5,
     .         0.,  1.,  0.,  0.,  0.,    1426.,   -3.4,     54., -0.1,
     .         1.,  0.,  0.,  0.,  0.,     712.,    0.1,     -7.,  0.0,
     .         0.,  1.,  2., -2.,  2.,    -517.,    1.2,    224., -0.6,
     .         0.,  0.,  2.,  0.,  1.,    -386.,   -0.4,    200.,  0.0,
     .         1.,  0.,  2.,  0.,  2.,    -301.,    0.0,    129., -0.1,
     .         0., -1.,  2., -2.,  2.,     217.,   -0.5,    -95.,  0.3/
      data x2/ 1.,  0.,  0., -2.,  0.,    -158.,    0.0,     -1.,  0.0,
     .         0.,  0.,  2., -2.,  1.,     129.,    0.1,    -70.,  0.0,
     .        -1.,  0.,  2.,  0.,  2.,     123.,    0.0,    -53.,  0.0,
     .         1.,  0.,  0.,  0.,  1.,      63.,    0.1,    -33.,  0.0,
     .         0.,  0.,  0.,  2.,  0.,      63.,    0.0,     -2.,  0.0,
     .        -1.,  0.,  2.,  2.,  2.,     -59.,    0.0,     26.,  0.0,
     .        -1.,  0.,  0.,  0.,  1.,     -58.,   -0.1,     32.,  0.0,
     .         1.,  0.,  2.,  0.,  1.,     -51.,    0.0,     27.,  0.0,
     .         2.,  0.,  0., -2.,  0.,      48.,    0.0,      1.,  0.0,
     .        -2.,  0.,  2.,  0.,  1.,      46.,    0.0,    -24.,  0.0/
      data x3/ 0.,  0.,  2.,  2.,  2.,     -38.,    0.0,     16.,  0.0,
     .         2.,  0.,  2.,  0.,  2.,     -31.,    0.0,     13.,  0.0,
     .         2.,  0.,  0.,  0.,  0.,      29.,    0.0,     -1.,  0.0,
     .         1.,  0.,  2., -2.,  2.,      29.,    0.0,    -12.,  0.0,
     .         0.,  0.,  2.,  0.,  0.,      26.,    0.0,     -1.,  0.0,
     .         0.,  0.,  2., -2.,  0.,     -22.,    0.0,      0.,  0.0,
     .        -1.,  0.,  2.,  0.,  1.,      21.,    0.0,    -10.,  0.0,
     .         0.,  2.,  0.,  0.,  0.,      17.,   -0.1,      0.,  0.0,
     .         0.,  2.,  2., -2.,  2.,     -16.,    0.1,      7.,  0.0,
     .        -1.,  0.,  0.,  2.,  1.,      16.,    0.0,     -8.,  0.0/
      data x4/ 0.,  1.,  0.,  0.,  1.,     -15.,    0.0,      9.,  0.0,
     .         1.,  0.,  0., -2.,  1.,     -13.,    0.0,      7.,  0.0,
     .         0., -1.,  0.,  0.,  1.,     -12.,    0.0,      6.,  0.0,
     .         2.,  0., -2.,  0.,  0.,      11.,    0.0,      0.,  0.0,
     .        -1.,  0.,  2.,  2.,  1.,     -10.,    0.0,      5.,  0.0,
     .         1.,  0.,  2.,  2.,  2.,      -8.,    0.0,      3.,  0.0,
     .         0., -1.,  2.,  0.,  2.,      -7.,    0.0,      3.,  0.0,
     .         0.,  0.,  2.,  2.,  1.,      -7.,    0.0,      3.,  0.0,
     .         1.,  1.,  0., -2.,  0.,      -7.,    0.0,      0.,  0.0,
     .         0.,  1.,  2.,  0.,  2.,       7.,    0.0,     -3.,  0.0/
      data x5/-2.,  0.,  0.,  2.,  1.,      -6.,    0.0,      3.,  0.0,
     .         0.,  0.,  0.,  2.,  1.,      -6.,    0.0,      3.,  0.0,
     .         2.,  0.,  2., -2.,  2.,       6.,    0.0,     -3.,  0.0,
     .         1.,  0.,  0.,  2.,  0.,       6.,    0.0,      0.,  0.0,
     .         1.,  0.,  2., -2.,  1.,       6.,    0.0,     -3.,  0.0,
     .         0.,  0.,  0., -2.,  1.,      -5.,    0.0,      3.,  0.0,
     .         0., -1.,  2., -2.,  1.,      -5.,    0.0,      3.,  0.0,
     .         2.,  0.,  2.,  0.,  1.,      -5.,    0.0,      3.,  0.0,
     .         1., -1.,  0.,  0.,  0.,       5.,    0.0,      0.,  0.0,
     .         1.,  0.,  0., -1.,  0.,      -4.,    0.0,      0.,  0.0/
      data x6/ 0.,  0.,  0.,  1.,  0.,      -4.,    0.0,      0.,  0.0,
     .         0.,  1.,  0., -2.,  0.,      -4.,    0.0,      0.,  0.0,
     .         1.,  0., -2.,  0.,  0.,       4.,    0.0,      0.,  0.0,
     .         2.,  0.,  0., -2.,  1.,       4.,    0.0,     -2.,  0.0,
     .         0.,  1.,  2., -2.,  1.,       4.,    0.0,     -2.,  0.0,
     .         1.,  1.,  0.,  0.,  0.,      -3.,    0.0,      0.,  0.0,
     .         1., -1.,  0., -1.,  0.,      -3.,    0.0,      0.,  0.0,
     .        -1., -1.,  2.,  2.,  2.,      -3.,    0.0,      1.,  0.0,
     .         0., -1.,  2.,  2.,  2.,      -3.,    0.0,      1.,  0.0,
     .         1., -1.,  2.,  0.,  2.,      -3.,    0.0,      1.,  0.0/
      data x7/ 3.,  0.,  2.,  0.,  2.,      -3.,    0.0,      1.,  0.0,
     .        -2.,  0.,  2.,  0.,  2.,      -3.,    0.0,      1.,  0.0,
     .         1.,  0.,  2.,  0.,  0.,       3.,    0.0,      0.,  0.0,
     .        -1.,  0.,  2.,  4.,  2.,      -2.,    0.0,      1.,  0.0,
     .         1.,  0.,  0.,  0.,  2.,      -2.,    0.0,      1.,  0.0,
     .        -1.,  0.,  2., -2.,  1.,      -2.,    0.0,      1.,  0.0,
     .         0., -2.,  2., -2.,  1.,      -2.,    0.0,      1.,  0.0,
     .        -2.,  0.,  0.,  0.,  1.,      -2.,    0.0,      1.,  0.0,
     .         2.,  0.,  0.,  0.,  1.,       2.,    0.0,     -1.,  0.0,
     .         3.,  0.,  0.,  0.,  0.,       2.,    0.0,      0.,  0.0/
      data x8/ 1.,  1.,  2.,  0.,  2.,       2.,    0.0,     -1.,  0.0,
     .         0.,  0.,  2.,  1.,  2.,       2.,    0.0,     -1.,  0.0,
     .         1.,  0.,  0.,  2.,  1.,      -1.,    0.0,      0.,  0.0,
     .         1.,  0.,  2.,  2.,  1.,      -1.,    0.0,      1.,  0.0,
     .         1.,  1.,  0., -2.,  1.,      -1.,    0.0,      0.,  0.0,
     .         0.,  1.,  0.,  2.,  0.,      -1.,    0.0,      0.,  0.0,
     .         0.,  1.,  2., -2.,  0.,      -1.,    0.0,      0.,  0.0,
     .         0.,  1., -2.,  2.,  0.,      -1.,    0.0,      0.,  0.0,
     .         1.,  0., -2.,  2.,  0.,      -1.,    0.0,      0.,  0.0,
     .         1.,  0., -2., -2.,  0.,      -1.,    0.0,      0.,  0.0/
      data x9/ 1.,  0.,  2., -2.,  0.,      -1.,    0.0,      0.,  0.0,
     .         1.,  0.,  0., -4.,  0.,      -1.,    0.0,      0.,  0.0,
     .         2.,  0.,  0., -4.,  0.,      -1.,    0.0,      0.,  0.0,
     .         0.,  0.,  2.,  4.,  2.,      -1.,    0.0,      0.,  0.0,
     .         0.,  0.,  2., -1.,  2.,      -1.,    0.0,      0.,  0.0,
     .        -2.,  0.,  2.,  4.,  2.,      -1.,    0.0,      1.,  0.0,
     .         2.,  0.,  2.,  2.,  2.,      -1.,    0.0,      0.,  0.0,
     .         0., -1.,  2.,  0.,  1.,      -1.,    0.0,      0.,  0.0,
     .         0.,  0., -2.,  0.,  1.,      -1.,    0.0,      0.,  0.0,
     .         0.,  0.,  4., -2.,  2.,       1.,    0.0,      0.,  0.0/
      data xa/ 0.,  1.,  0.,  0.,  2.,       1.,    0.0,      0.,  0.0,
     .         1.,  1.,  2., -2.,  2.,       1.,    0.0,     -1.,  0.0,
     .         3.,  0.,  2., -2.,  2.,       1.,    0.0,      0.,  0.0,
     .        -2.,  0.,  2.,  2.,  2.,       1.,    0.0,     -1.,  0.0,
     .        -1.,  0.,  0.,  0.,  2.,       1.,    0.0,     -1.,  0.0,
     .         0.,  0., -2.,  2.,  1.,       1.,    0.0,      0.,  0.0,
     .         0.,  1.,  2.,  0.,  1.,       1.,    0.0,      0.,  0.0,
     .        -1.,  0.,  4.,  0.,  2.,       1.,    0.0,      0.,  0.0,
     .         2.,  1.,  0., -2.,  0.,       1.,    0.0,      0.,  0.0,
     .         2.,  0.,  0.,  2.,  0.,       1.,    0.0,      0.,  0.0/
      data xb/ 2.,  0.,  2., -2.,  1.,       1.,    0.0,     -1.,  0.0,
     .         2.,  0., -2.,  0.,  1.,       1.,    0.0,      0.,  0.0,
     .         1., -1.,  0., -2.,  0.,       1.,    0.0,      0.,  0.0,
     .        -1.,  0.,  0.,  1.,  1.,       1.,    0.0,      0.,  0.0,
     .        -1., -1.,  0.,  2.,  1.,       1.,    0.0,      0.,  0.0,
     .         0.,  1.,  0.,  1.,  0.,       1.,    0.0,      0.,  0.0/
c
c***********************************************************************
c
c
c     get fundamental arguments
c
      call funarg (t,l,lp,f,d,om)
c
c***********************************************************************
c
c
c     sum nutation series terms, from smallest to largest
c
      dpsi = 0.d0
      deps = 0.d0
c
      do 10 j=1,106
      i = 107 - j
c     formation of multiples of arguments
      arg = dble(x(1,i)) * l
     .    + dble(x(2,i)) * lp
     .    + dble(x(3,i)) * f
     .    + dble(x(4,i)) * d
     .    + dble(x(5,i)) * om
c     evaluate nutation
      dpsi = (dble(x(6,i)) + dble(x(7,i))*t) * dsin(arg) + dpsi
      deps = (dble(x(8,i)) + dble(x(9,i))*t) * dcos(arg) + deps
   10 continue
c
c
      dpsi = dpsi * 1.0d-4
      deps = deps * 1.0d-4
c
c***********************************************************************
c
c
      return
c
      end





      subroutine solsys (tjd,body,origin, pos,vel,ierr)
c
c     subroutine solsys version 2.
c
c-----------------------------------------------------------------------
c
c---purpose: this is solsys version 2.  it is intended to provide
c            an interface between the jpl binary direct-access solar
c            system ephemerides and the 'novas' astrometric subroutine
c            library.
c
c---reference:  documentation by e. m. standish on cd-rom
c            'jpl planetary and lunar ephemerides' (c)1997 by jpl
c            and published by willmann-bell, inc., richmond, va, usa.
c
c---input arguments:     tjd = julian date of the desired time, on
c                              the tdb time scale (double precision).
c                       body = body identification number for the
c                              solar system object of interest;
c                              mercury= 1,...,pluto= 9, sun= 10,
c                              moon= 11 (integer).
c                     origin = origin code; solar system barycenter= 0,
c                              center of mass of the sun= 1 (integer).
c
c---output arguments:    pos = position vector of 'body' at tjd;
c                              equatorial rectangular coordinates in
c                              au referred to the mean equator and
c                              equinox of j2000.0 (double precision).
c                        vel = velocity vector of 'body' at tjd;
c                              equatorial rectangular system referred
c                              to the mean equator and equinox of
c                              j2000.0, in au/day (double precision).
c                       ierr = 0 ... everything ok
c                            = 1 ... 'tjd' before first ephemeris date
c                            = 2 ... 'tjd' after last ephemeris date
c                            = 3 ... invalid value of 'body' or
c                                    'origin' (integer).
c
c---common blocks: none.
c
c---subroutines called: subroutine const   (jpl)
c                       subroutine pleph   (jpl)
c                       subroutine auxpos  (supplied)
c
c---version/date/programmer: v1/02-90/jab
c                            v2/07-91/ghk
c                            v3/05-98/ghk
c
c---notes: 1. subroutine pleph is a jpl-supplied routine that calls
c             a variety of other jpl subroutines.
c          2. this routine is for use with 1997 version of jpl
c             ephemeris software, e.g., as distributed on cd-rom
c             'jpl planetary and lunar ephemerides' (c)1997 by jpl.
c             it may not work properly with previous versions of the
c             jpl software.
c          3. for body identification numbers outside the range 1-11,
c             subroutine 'auxpos' will be called to supply positions
c             and velocities from sources external to the jpl
c             ephemerides.  a dummy version of this routine is
c             provided, which can be replaced by the user.
c
c-----------------------------------------------------------------------
c
      implicit none
      save

      integer          body,origin,ierr,targ,cent,i,n,jerr
      double precision tjd,pos(3),vel(3),posvel(6)

      include 'ephem.dek'



c---initialize output arguments
      pos(1) =  0.d0
      pos(2) =  0.d0
      pos(3) = 99.d0
      vel(1) =  0.d0
      vel(2) =  0.d0
      vel(3) =  0.d0
      ierr = 0


c---perform sanity checks on the input body and origin.
      if ( ( origin .lt. 0 ) .or. ( origin .gt. 1 ) ) then
          ierr = 3
          go to 99


c..call auxpos for auxiliary bodies (if any)
      else if ( ( body .lt. 1 ) .or. ( body .gt. nplanets ) ) then


c..for the martian moons
        if (body .ge. nplanets+1  .and. 
     1      body .le. nsum1) then
         call marsmoonpos ( tjd, body, origin,   pos, vel, jerr )


c..for the jovian moons
        else if (body .ge. nsum1+1  .and.  
     1           body .le. nsum2) then
         call jupmoonpos ( tjd, body, origin,   pos, vel, jerr )


c..for the saturn moons
        else if (body .ge. nsum2+1  .and.  
     1           body .le. nsum3) then
         call satmoonpos ( tjd, body, origin,   pos, vel, jerr )


c..for the uranus moons
        else if (body .ge. nsum3+1  .and.  
     1           body .le. nsum4) then
         call uramoonpos ( tjd, body, origin,   pos, vel, jerr )


c..for comets and asteroids
        else 
         call auxpos ( tjd, body, origin,   pos, vel, jerr )
        end if
        ierr = jerr
        go to 99
      endif


c---check that requested julian date is within range of ephemeris.
      if ( tjd .lt. begjd ) then
          ierr = 1
          go to 99
      else if ( tjd .gt. endjd ) then
          ierr = 2
          go to 99
      endif


c---select 'targ' according to value of 'body'.
      if ( body .eq. 10 ) then
          targ = 11
      else if ( body .eq. 11 ) then
          targ = 10
      else
          targ = body
      endif


c---select 'cent' according to the value of 'origin'.
      if ( origin .eq. 0 ) then
          cent = 12
      else if ( origin .eq. 1 ) then
          cent = 11
      endif


c---call jpl routine 'pleph' to obtain position and velocity array 'posvel'

      call pleph ( tjd, targ, cent,   posvel )

      do 10 i = 1, 3
       pos(i) = posvel(i)
       vel(i) = posvel(i+3)
   10 continue


   99 continue
      return
      end




      subroutine marsmoonpos (tjd,m,k, pos,vel,jerr)
      implicit none
      save
c
c          tjd  = tdb julian date of desired epoch (in)
c          m    = body identification number (in)
c          k    = origin selection code (in)
c                 set k=0 for origin at solar system barycenter
c                 set k=1 for origin at center of sun
c          pos  = position vector, equatorial rectangular
c                 coordinates, referred to mean equator and equinox
c                 of j2000.0, components in au (out)
c          vel  = velocity vector, equatorial rectangular
c                 coordinates, referred to mean equator and equinox
c                 of j2000.0, components in au/day (out)
c          jerr = error indicator (out)
c                 jerr=0 means everything ok
c                 jerr=1 means tjd before first valid date
c                 jerr=2 means tjd after last valid date
c                 jerr=3 means invalid value of m or k
c
c     generally, it is necessary for this routine
c     to provide only barycentric positions for the input jd.
c     apparent place routines such as applan need velocities and
c     heliocentric positions only for the earth (which will be
c     obtained from the jpl ephemeris).  also, do not use fortran
c     i/o unit number 12, which is used by the jpl routines.
c


c..declare the pass
      integer          m,k,jerr
      double precision tjd,pos(3),vel(3)


c..for communicating the orbital plane solution
      double precision rorb(6)
      common /oplanec/ rorb

c..local variables
      integer          i,ksat
      double precision xp(3),vp(3),posvel(6),kmau,kmaui
      parameter        (kmau  = 1.49597870d8,
     1                  kmaui = 1.0d0/kmau)


      include 'ephem.dek'


c..martain satellite number
      ksat = m - nplanets


c..get the satellite state vector
      if (ksat .eq. 1) then
       call ephob(tjd,xp,vp,pos,vel)   
      else if (ksat .eq. 2) then
       call edeim(tjd,xp,vp,pos,vel) 
      endif

c..convert from km to au
      do i=1,3
       xp(i)  = xp(i) * kmaui
       vp(i)  = vp(i) * kmaui
       pos(i) = pos(i) * kmaui
       vel(i) = vel(i) * kmaui
      enddo

c..store the orbital plane solution
      do i=1,3
       rorb(i)   = xp(i)
       rorb(i+3) = vp(i)
      enddo


c..get the position of mars with respect to the barycenter
      call pleph(tjd,4,12,posvel)


c..the heliocentric equatorial state vector 
      do i=1,3
       pos(i) = posvel(i) + pos(i)      
       vel(i) = posvel(i+3) + vel(i)
      enddo

      jerr   = 0
      return
      end





      subroutine jupmoonpos (tjd,m,k, pos,vel,jerr)
      implicit none
      save
c
c
c          tjd  = tdb julian date of desired epoch (in)
c          m    = body identification number (in)
c          k    = origin selection code (in)
c                 set k=0 for origin at solar system barycenter
c                 set k=1 for origin at center of sun
c          pos  = position vector, equatorial rectangular
c                 coordinates, referred to mean equator and equinox
c                 of j2000.0, components in au (out)
c          vel  = velocity vector, equatorial rectangular
c                 coordinates, referred to mean equator and equinox
c                 of j2000.0, components in au/day (out)
c          jerr = error indicator (out)
c                 jerr=0 means everything ok
c                 jerr=1 means tjd before first valid date
c                 jerr=2 means tjd after last valid date
c                 jerr=3 means invalid value of m or k
c
c     generally, it is necessary for this routine
c     to provide only barycentric positions for the input jd.
c     apparent place routines such as applan need velocities and
c     heliocentric positions only for the earth (which will be
c     obtained from the jpl ephemeris).  also, do not use fortran
c     i/o unit number 12, which is used by the jpl routines.
c

c..declare the pass
      integer          m,k,jerr
      double precision tjd,pos(3),vel(3)

c..local variables
      integer          i,ksat
      double precision r(6),posvel(6)


c..for communicating the orbital plane solution
      double precision rorb(6)
      common /oplanec/ rorb

      include 'ephem.dek'



c..get the satellite state vector
      ksat = m - nsum1
      call galsat(r,rorb,tjd,ksat,2)


c..get the position of jupiter with respect to the barycenter
      call pleph(tjd,5,12,posvel)


c..the heliocentric equatorial state vector 
      do i=1,3
       pos(i) = posvel(i) + r(i)      
       vel(i) = posvel(i+3) + r(i+3)
      enddo

      jerr   = 0
      return
      end






      subroutine satmoonpos (tjd,m,k, pos,vel,jerr)
      implicit none
      save
c
c
c          tjd  = tdb julian date of desired epoch (in)
c          m    = body identification number (in)
c          k    = origin selection code (in)
c                 set k=0 for origin at solar system barycenter
c                 set k=1 for origin at center of sun
c          pos  = position vector, equatorial rectangular
c                 coordinates, referred to mean equator and equinox
c                 of j2000.0, components in au (out)
c          vel  = velocity vector, equatorial rectangular
c                 coordinates, referred to mean equator and equinox
c                 of j2000.0, components in au/day (out)
c          jerr = error indicator (out)
c                 jerr=0 means everything ok
c                 jerr=1 means tjd before first valid date
c                 jerr=2 means tjd after last valid date
c                 jerr=3 means invalid value of m or k
c
c     generally, it is necessary for this routine
c     to provide only barycentric positions for the input jd.
c     apparent place routines such as applan need velocities and
c     heliocentric positions only for the earth (which will be
c     obtained from the jpl ephemeris).  also, do not use fortran
c     i/o unit number 12, which is used by the jpl routines.
c


c..declare the pass
      integer          m,k,jerr
      double precision tjd,pos(3),vel(3),kmau,kmaui,dayi,tjd2000,
     1                 oblm,oblt,eqeq,dpsi,deps
      parameter        (kmau    = 1.49597870d8,
     1                  kmaui   = 1.0d0/kmau,
     2                  dayi    = 1.0d0/365.25d0,
     3                  tjd2000 = 2451545.00000000d0)


c..for communicating the orbital plane solution
      double precision rorb(6)
      common /oplanec/ rorb

c..local variables
      integer          i,ksat
      double precision posvel(6)

      include 'ephem.dek'



c..saturn satellite number and ecliptic state vector
      ksat = m - nsum2
      call posired(tjd,ksat,pos,vel)


c..convert from au/year to au/day
      do i=1,3
       vel(i)  = vel(i) *dayi
      enddo


c..tass17 returns the j2000 mean ecliptic solution
c..so we have to rotate to the j2000 mean equator
      call etilt (tjd2000,oblm,oblt,eqeq,dpsi,deps)
      call rotax(pos,-oblm,pos)
      call rotax(vel,-oblm,vel)


c..store the orbital plane solution
      do i=1,3
       rorb(i)   = pos(i)  
       rorb(i+3) = vel(i)
      enddo


c..get the position of saturn with respect to the barycenter
      call pleph(tjd,6,12,posvel)


c..the heliocentric equatorial state vector 
      do i=1,3
       pos(i) = posvel(i) + pos(i)      
       vel(i) = posvel(i+3) + vel(i)
      enddo

      jerr   = 0
      return
      end




      subroutine uramoonpos (tjd,m,k, pos,vel,jerr)
      implicit none
      save
c
c
c          tjd  = tdb julian date of desired epoch (in)
c          m    = body identification number (in)
c          k    = origin selection code (in)
c                 set k=0 for origin at solar system barycenter
c                 set k=1 for origin at center of sun
c          pos  = position vector, equatorial rectangular
c                 coordinates, referred to mean equator and equinox
c                 of j2000.0, components in au (out)
c          vel  = velocity vector, equatorial rectangular
c                 coordinates, referred to mean equator and equinox
c                 of j2000.0, components in au/day (out)
c          jerr = error indicator (out)
c                 jerr=0 means everything ok
c                 jerr=1 means tjd before first valid date
c                 jerr=2 means tjd after last valid date
c                 jerr=3 means invalid value of m or k
c
c     generally, it is necessary for this routine
c     to provide only barycentric positions for the input jd.
c     apparent place routines such as applan need velocities and
c     heliocentric positions only for the earth (which will be
c     obtained from the jpl ephemeris).  also, do not use fortran
c     i/o unit number 12, which is used by the jpl routines.
c


c..declare the pass
      integer          m,k,jerr
      double precision tjd,pos(3),vel(3),kmau,kmaui
      parameter        (kmau  = 1.49597870d8,
     1                  kmaui = 1.0d0/kmau)


c..for communicating the orbital plane solution
      double precision rorb(6)
      common /oplanec/ rorb

c..local variables
      integer          i,ksat
      double precision posvel(6)

      include 'ephem.dek'



c..saturn satellite number and state vector
      ksat = m - nsum3

      call gust86(tjd,ksat,pos,vel)

c..convert from km and km/s to au and au/day
      do i=1,3
       pos(i) = pos(i) * kmaui
       vel(i) = vel(i) * kmaui * 86400.0d0
      enddo


c..store the orbital plane solution
      do i=1,3
       rorb(i)   = pos(i)
       rorb(i+3) = vel(i)
      enddo


c..get the position of uranus with respect to the barycenter
      call pleph(tjd,7,12,posvel)


c..the heliocentric j2000 equatorial state vector 
      do i=1,3
       pos(i) = posvel(i) + pos(i)      
       vel(i) = posvel(i+3) + vel(i)
      enddo

      jerr   = 0
      return
      end








      subroutine auxpos (tjd,m,k, pos,vel,jerr)
      implicit none
      save
c
c          tjd  = tdb julian date of desired epoch (in)
c          m    = body identification number (in)
c          k    = origin selection code (in)
c                 set k=0 for origin at solar system barycenter
c                 set k=1 for origin at center of sun
c          pos  = position vector, equatorial rectangular
c                 coordinates, referred to mean equator and equinox
c                 of j2000.0, components in au (out)
c          vel  = velocity vector, equatorial rectangular
c                 coordinates, referred to mean equator and equinox
c                 of j2000.0, components in au/day (out)
c          jerr = error indicator (out)
c                 jerr=0 means everything ok
c                 jerr=1 means tjd before first valid date
c                 jerr=2 means tjd after last valid date
c                 jerr=3 means invalid value of m or k
c
c     generally, it is necessary for this routine
c     to provide only barycentric positions for the input jd.
c     apparent place routines such as applan need velocities and
c     heliocentric positions only for the earth (which will be
c     obtained from the jpl ephemeris).  also, do not use fortran
c     i/o unit number 12, which is used by the jpl routines.



c..declare the pass
      integer          m,k,jerr
      double precision tjd,pos(3),vel(3)


c..math constants
      double precision pi,a2rad,rad2a,twopi
      parameter        (pi     = 3.1415926535897932384d0, 
     1                  a2rad  = pi/180.0d0,  
     2                  rad2a  = 180.0d0/pi,
     3                  twopi  = 2.0d0*pi)


c..physical constants
      double precision kgauss,kearth,gm
      parameter        (kgauss = 0.01720209895d0,
     1                  kearth = 0.07436680d0,
     2                  gm     = kgauss*kgauss)

      include 'ephem.dek'



c..for integrating the orbits
      external          derv,jaco,bsstep
      integer           xdim,ydim,kount,nok,nbad,iprint,i
      parameter         (xdim=10, ydim=6)
      double precision  xrk(xdim),yrk(ydim,xdim),bc(ydim),dydx,stptry,
     1                  stpmin,tol,start,sstop,dxsav,odescal
      parameter         (tol = 1.0d-8, odescal=1.0d0,
     1                   dxsav=0.0d0, iprint=0)


c..communication between routine derv and this routine for the origin
      integer           icenter
      common /auxint1/  icenter



c..for picking how to do the bodies
c..if imethod = 1, then static orbital elements are used
c..                this is fast, but not very accurate far from
c..                the epoch for which the elements are given
c..
c..if imethod = 2, then the orbit is determined by numerical integration
c..                this is slower, but more accurate far from
c..                the epoch for which the elements are given

      integer          imethod
      parameter        (imethod = 1)




c..local variables
      integer          j,kk,kk1,iyear,imonth,iday,ibeg
      double precision xrate,pecl(3),vecl(3),x,y,z,
     1                 oblm,oblt,eqeq,dpsi,deps,
     2                 rhour,tjdp,tjd2000,tjd_start,qper,
     3                 manom,semi_maj,ecen,asend,xinclin,argper,tjdm,
     4                 manom_int,pdum(3),vdum(3),
     5                 tjde,secdif
      parameter        (tjd2000 = 2451545.0d0)



c..popular format statements
 11     format(a)
 222    format(1x,1p6e16.8)




c..load the arguments: perihelion distance, semimajor axis, 
c..longitude of ascending node, inclination, eccentricity, 
c..the argument of perhelion and either the perihelion date
c..or the mean anomaly. the given elements are slightly different for
c..each class of minor objects, we'll handle their differences here.

c..for comets
      if (m .ge. nsum4 + 1  .and. 
     1    m .le. nsum4 + ncomet) then
       j        = m - nsum4 
       qper     = comet_q(j)
       asend    = comet_node(j)
       xinclin  = comet_i(j)
       ecen     = comet_e(j)
       argper   = comet_w(j)
       iyear    = comet_tp_year(j)
       imonth   = comet_tp_month(j)
       iday     = comet_tp_day(j)
       rhour    = comet_tp_hour(j)
       call juldat2(iyear,imonth,iday,rhour,tjdp)

c..start the orbit determination from the perihelion date
c..where the mean anomaly is zero.
       tjd_start = tjdp
       manom_int = 0.0d0



c..for the numbered asteroids
      else if (m .ge. nsum4 + ncomet + 1 .and. 
     1         m .le. nsum4 + ncomet + nnaster) then
       j        = m - nsum4 - ncomet
       semi_maj = naster_a(j)
       asend    = naster_node(j)
       xinclin  = naster_i(j)
       ecen     = naster_e(j)
       argper   = naster_w(j)
       qper     = semi_maj * abs(1.0d0 - ecen)
       tjdp     = 0.0d0

c..start from the given epoch and mean anomaly
       tjd_start = naster_epoch(j) + 2400000.5d0
       manom_int = naster_m(j)



c..else its a blown body number
      else
       goto 1000
      end if



c..put all times on the tdb scale
       call deltat(tjd_start,tjde,secdif)
       tjd_start = tjde
       call novas_times(tjd_start,x,secdif)
       tjd_start = tjd_start + secdif / 86400.0d0

       call deltat(tjdp,tjde,secdif)
       tjdp = tjde
       call novas_times(tjdp,x,secdif)
       tjdp = tjdp + secdif / 86400.0d0




c..use the static keplerian elements for a fast evaluation

      if (imethod .eq. 1) then

c..get the initial ecliptic state vector from the orbital elements
       call kep_to_state(gm,
     1   tjdp,tjd_start,manom_int,qper,ecen,asend,xinclin,argper,
     2   pdum,vdum)


c..and get the exact solution to the two body problem
       call twobody(gm,
     1       tjd_start,pdum(1),pdum(2),pdum(3),vdum(1),vdum(2),vdum(3),
     2       tjd,pecl(1),pecl(2),pecl(3),vecl(1),vecl(2),vecl(3))


c..convert the mean j2000 ecliptic state vector to a 
c..mean j2000 equatorial state vector
       call etilt (tjd2000,oblm,oblt,eqeq,dpsi,deps)

       call recl_to_requ(pecl(1),pecl(2),pecl(3),oblm,
     1                   pos(1),pos(2),pos(3))

       call recl_to_requ(vecl(1),vecl(2),vecl(3),oblm,
     1                   vel(1),vel(2),vel(3))

       jerr = 0
       return




c..do an integration of the objects orbit 
      else if (imethod .eq. 2) then


c..get the initial ecliptic state vector from the orbital elements
       call kep_to_state(gm,
     1      tjdp,tjd_start,manom_int,qper,ecen,asend,xinclin,argper,
     2      pecl,vecl)


c..load the state vector as the initial conditions
       bc(1) = pecl(1)
       bc(2) = pecl(2)
       bc(3) = pecl(3)
       bc(4) = vecl(1)
       bc(5) = vecl(2)
       bc(6) = vecl(3)


c..set the start, step and stop times
       start   = tjd_start
       stptry  = 1.0d-3 * sign(1.0d0,(tjd - tjd_start))
       stpmin  = 1.0d-6 
       sstop   = tjd

c..set if this pass is the solar system barycenter or center of the sun
       if (k .eq. 0) then
        icenter = 12
       else if (k .eq. 1) then
        icenter = 11
       else
        write(6,*)
        write(6,*) 'impossible center for integration in routine auxpos'
        stop 'bad center in routine auxpos'
       end if
        

c..and integrate
       call cometint(start,stptry,stpmin,sstop,bc,
     1               tol,dxsav,xdim,
     2               xrk,yrk,xdim,ydim,xdim,ydim,
     3               nok,nbad,kount,odescal,iprint,
     4               derv,jaco,bsstep)


c..convert the ecliptic state vector to an equatorial state vector
       call etilt (tjd,oblm,oblt,eqeq,dpsi,deps)

       call recl_to_requ(bc(1),bc(2),bc(3),oblt,
     1                   pos(1),pos(2),pos(3))

       call recl_to_requ(bc(4),bc(5),bc(6),oblt,
     1                   vel(1),vel(2),vel(3))

       jerr = 0
       return
      end if



c..land here if we encountered an error
 1000 continue
      write(6,*)
      write(6,*) 'bad body number in routine auxpos',m
      write(6,*) 'setting zero positions and velocities and returning'
      write(6,*) 
      pos(1) =  0.d0
      pos(2) =  0.d0
      pos(3) = 99.d0
      vel(1) =  0.d0
      vel(2) =  0.d0
      vel(3) =  0.d0
      jerr   = 3
      return

      end