Quick snippet to use astropy to convert from heliocentric times to barycentric and vice-versa

helio_bary.py

from astropy.coordinates import SkyCoord, EarthLocation
from astropy import units as u
from astropy.time import Time


def helio_to_bary(coords, hjd, obs_name):
    helio = Time(hjd, scale='utc', format='jd')
    obs = EarthLocation.of_site(obs_name)
    star = SkyCoord(coords, unit=(u.hour, u.deg)) 
    ltt = helio.light_travel_time(star, 'heliocentric', location=obs)
    guess = helio - ltt
    # if we assume guess is correct - how far is heliocentric time away from true value?
    delta = (guess + guess.light_travel_time(star, 'heliocentric', obs)).jd  - helio.jd
    # apply this correction
    guess -= delta * u.d

    ltt = guess.light_travel_time(star, 'barycentric', obs)
    return guess.tdb + ltt


def bary_to_helio(coords, bjd, obs_name):
    bary = Time(bjd, scale='tdb', format='jd')
    obs = EarthLocation.of_site(obs_name)
    star = SkyCoord(coords, unit=(u.hour, u.deg))
    ltt = bary.light_travel_time(star, 'barycentric', location=obs) 
    guess = bary - ltt
    delta = (guess + guess.light_travel_time(star, 'barycentric', obs)).jd  - bary.jd
    guess -= delta * u.d
    
    ltt = guess.light_travel_time(star, 'heliocentric', obs)
    return guess.utc + ltt

Hi Ian,

What you've done is slightly incorrect. There are two issues; one is that the observer-provided t0 is on the TDB timescale, and you should initialise as such:

t = Time(dt, scale='tdb',  location=(str(INT_Lon) + 'd', str(INT_Lat) + 'd'))

The second is that the call to light_travel_time will not quite get the correct time, as it is meant to be used the other way around. This is what the second iteration in my functions is meant to correct for; guess the light_travel_time, use the guess to correct to UTC, and see if the corrected UTC converts to the given TDB. A version of your code using these fixes would look like:

def convert_tdb_to_utc(dt, targ):
    """
    Convert TDB to UTC
    'targ' is a SkyCoord object of the target star
    """
    t = Time(dt, scale='tdb',  location=(str(INT_Lon) + 'd', str(INT_Lat) + 'd'))
    td = t.light_travel_time(targ)
    guess = t - td
    # if guess is correct, what's the difference between the inferred BJD and the given BJD?
    delta = (guess + guess.light_travel_time(targ)).jd  - t.jd
    # apply this correction
    guess -= delta * u.d
    # return on UTC scale as date time
    return guess.utc.datetime

However, the correction involves a second call to light_travel_time and only provides a minor change in precision - around 1/100th of a second. So a cruder but faster version (x2) would just be

def convert_tdb_to_utc(dt, targ):
    """
    Convert TDB to UTC
    'targ' is a SkyCoord object of the target star
    """
    t = Time(dt, scale='tdb',  location=(str(INT_Lon) + 'd', str(INT_Lat) + 'd'))
    td = t.light_travel_time(targ)
    correct = t - td
    return correct.utc.datetime

Thanks Stuart, I like the idea of checking the guess and correcting it slightly. We don’t have to run this process too many times in the code, only in setup in fact, so I might put this in for completeness. Thanks for your help, it’s much appreciated. Ian