[AstroPy] Heliocentric/Barycentric JD corrections

[Erik Tollerud]

For those who aren't monitoring astropy-dev and github, the plan is
that astropy.time will probably support these eventually, but it
requires some additional work to get it to operate smoothly in the
existing architecture.

So this is a good interim solution - thanks, Leo!

On Wed, Apr 16, 2014 at 11:11 AM, Leo Huckvale
<leo.huckvale at postgrad.manchester.ac.uk> wrote:
> Hi all,
>
>   Following recent discussion about implementing HJD on astropy-dev and
> github, I thought I might offer up one solution to correct timestamps to
> HJD/BJD. It seemed more appropriate here than on the dev mailing list.
>
>   The code below defines the function "jd_corr" which takes Modified
> Julian Dates (on UTC scale), source RA and Dec and a correction type
> "bjd" or "hjd". This returns "new_jd", which is the corrected set of
> time stamps in an astropy.time.Time object. The delta-t corrections are
> obtained from barycentric ephemerides via jplephem, as recommended in
> the discussion on github, and a corresponding data package (in this
> case, de423).
>
> Hopefully this might help someone. I'd appreciate any feedback -
> especially if I've overlooked something vital, which is entirely
> possible! I have relied in a large part on astropy.time to work its
> magic on time-standard conversions.
>
> Best regards,
> Leo
>
> ----
>
> import numpy as np
>
> import astropy.time as time
> import astropy.coordinates as coords
> import astropy.units as u
> import astropy.constants as const
>
> import jplephem
> import de423
>
> VISTA_LAT = -24.6157    # Geographic coordinates of observatory
> VISTA_LON = -70.3976
>
> def jd_corr(mjd, ra, dec, jd_type='bjd'):
>      """Return BJD or HJD for input MJD(UTC)."""
>
>      # Initialise ephemeris from jplephem
>      eph = jplephem.Ephemeris(de423)
>
>      # Source unit-vector
>      ## Assume coordinates in ICRS
>      ## Set distance to unit (kilometers)
>      src_vec = coords.ICRS(ra=ra, dec=dec,
>                            unit=(u.degree, u.degree),
>                            distance=coords.Distance(1, u.km))
>
>      # Convert epochs to astropy.time.Time
>      ## Assume MJD(UTC)
>      t = time.Time(mjd, scale='utc', format='mjd', lat=VISTA_LAT,
> lon=VISTA_LON)
>
>      # Get Earth-Moon barycenter position
>      ## NB: jplephem uses Barycentric Dynamical Time, e.g. JD(TDB)
>      ## and gives positions relative to solar system barycenter
>      barycenter_earthmoon = eph.position('earthmoon', t.tdb.jd)
>
>      # Get Moon position vectors
>      moonvector = eph.position('moon', t.tdb.jd)
>
>      # Compute Earth position vectors
>      pos_earth = (barycenter_earthmoon - moonvector * eph.earth_share)*u.km
>
>      if jd_type == 'bjd':
>          # Compute BJD correction
>          ## Assume source vectors parallel at Earth and Solar System
> Barycenter
>          ## i.e. source is at infinity
>          corr = np.dot(pos_earth.T, src_vec.cartesian.value)/const.c
>      elif jd_type == 'hjd':
>          # Compute HJD correction via Sun ephemeris
>          pos_sun = eph.position('sun', t.tdb.jd)*u.km
>          sun_earth_vec = pos_earth - pos_sun
>          corr = np.dot(sun_earth_vec.T, src_vec.cartesian.value)/const.c
>
>      # TDB is the appropriate time scale for these ephemerides
>      dt = time.TimeDelta(corr, scale='tdb', format='jd')
>
>      # Compute and return HJD/BJD as astropy.time.Time
>      new_jd = t + dt
>
>      return new_jd
>
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy



-- 
Erik T