Add usage example for converting from UTC to HJD and BJD

[icweaver]

Not sure where this should live, so just sharing a modified snippet from our slack convo here for now. Will hopefully be useful for future docs, happy to make any edits.

Background

Given a time in UTC, we want to get its resulting light travel time lags and corresponding times in HJD and BJD. AstroLib.jl can get us part of the way there with HJD, but it would be great to have a comprehensive example using more modern tools like Astrometry.jl. Below is some comparison code showing how we might accomplish this (+ AstroTime.jl to help make dealing with timescales a bit less painful):

Setup

using Dates: DateTime
using AstroTime: TDBEpoch, TDB, to_utc, from_utc, days, julian, value

# Observation params
t_utc = DateTime(2026, 02, 21, 23, 0, 0)
ra, dec = 122.39896, 44.47156 # Degrees

AstroLib.jl

using AstroLib: helio_jd

function t_corrected_astrolib(t, ra, dec)
    rjd = value(julian(from_utc(t; scale = TDB))) - 2_400_000.0  # Reduced JD
    rhjd = helio_jd(rjd, ra, dec) # Reduced HJD
    hjd = rhjd + 2_400_000.0 # HJD
    ep = TDBEpoch(hjd * days; origin = :julian)
    return (; hjd, ep, utc = to_utc(DateTime, ep))
end

t_helio_astrolib = t_corrected_astrolib(t_utc, ra, dec)

(
    hjd = 2.461093463238254e6,
    ep = 2026-02-21T23:07:03.785 TDB,
    utc = 2026-02-21T23:05:54.599,
)

Astrometry.jl

using Astrometry: SOFA
using LinearAlgebra: ⋅

function t_corrected_astrom(t, ra, dec; kind = :heliocentric)
    # Astrometry
    jdtdb = from_utc(t; scale = TDB) |> julian |> value
        # Can also use Astrometry.SOFA built-ins
        # jdutc = datetime2julian(t)
        # jdtai = SOFA.utctai(jdutc, 0.0)
        # jdtt = SOFA.taitt(jdtai...)
        # dtr = SOFA.dtdb(jdtt..., 0.0, 0.0, 0.0, 0.0)
        # jdtdb = SOFA.tttdb(jdtt..., dtr)
    astrom = SOFA.apcg13(jdtdb, 0.0)
    gcra, gcdec = SOFA.atciq(deg2rad(ra), deg2rad(dec), 0.0, 0.0, 0.0, 0.0, astrom)

    # Heliocentric light travel time delay
    r⃗ = if kind == :heliocentric
        astrom.em * astrom.eh # Earth -- Sun vec, AU
    else
        astrom.eb # Earth -- Solar System barycenter vec, AU
    end
    n̂ = SOFA.s2c(gcra, gcdec) # Earth -- source unit vec
    c = SOFA.LIGHTSPEED / SOFA.AU * SOFA.SECPERDAY # AU/day
    Δt = r⃗ ⋅ n̂ / c # days

    # Store results
    bjd = jdtdb + Δt
    ep = TDBEpoch(bjd * days; origin = :julian)
    utc = to_utc(DateTime, ep)

    return if kind == :heliocentric
        (; hjd = bjd, ep, utc)
    else
        (; bjd, ep, utc)
    end
end

t_helio_astrom = t_corrected_astrom(t_utc, ra, dec)
t_bary_astrom = t_corrected_astrom(t_utc, ra, dec; kind = :barycentric)

# helio
(
    hjd = 2.4610934632390006e6,
    ep = 2026-02-21T23:07:03.849 TDB,
    utc = 2026-02-21T23:05:54.664,
)

# bary
(
    bjd = 2.4610934632187826e6,
    ep = 2026-02-21T23:07:02.102 TDB,
    utc = 2026-02-21T23:05:52.917,
)

Comparison

• The Astrometry.jl helio method is within 65 ms of AstroLib.jl
• The Astrometry.jl bary method is within 23 ms of the OSU calculator

# Value from https://astroutils.astronomy.osu.edu/time/
# for t_utc = 2026-02-21T23:00:00
bjd_OSU = 2461093.463218510 # BJD_TDB

utc_OSU = let
    ep = TDBEpoch(bjd_OSU * days; origin = :julian)
    to_utc(DateTime, ep)
end

# 2026-02-21T23:05:52.894

t_helio_astrom.hjd - t_helio_astrolib.hjd
# 7.46455043554306e-7

t_helio_astrom.utc - t_helio_astrolib.utc
# 65 milliseconds

t_bary_astrom.bjd - bjd_OSU
# 2.7241185307502747e-7

t_bary_astrom.utc - utc_OSU
# 23 milliseconds