Position of the Moon (topocentric) using Astrometry.jl (SOFA/ERFA) routines

tamasgal · February 13, 2026, 10:25am

OK I finally ask for help here because I am losing my head

I want to calculate the Moon coordinates in an equatorial coordinate system for an observer on Earth at a given time using Astrometry.jl which uses SOFA/ERFA routines..

What I got so far is calculating the Moon coordinates in a geocentric coordinate system using this tiny function

function moon(dt::DateTime)
    ht = 0.0  # sitting on Earht's surface
    jd = (datetime2julian(dt), 0.0)
    tt = SOFA.taitt(SOFA.utctai(jd...)...)
    pm, vm = SOFA.moon98(tt...)
    rc, dc = SOFA.c2s(pm)  # cartesian to spherical
    ra, dec = mod2pi(rc), dc
    return (ra=ra, dec=dec)
end

and it matches AstroPy’s get_body(...) when called without a location.

Fine, the problem I have now is to calculate the coordinates for a given location, which means I need to calculate it in topocentric system.

Using AstroPy (btw. I get the same results using SLALIB’s RDPLAN routine) I obtain the following values for a given location and time:

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

location = EarthLocation(
    lat=0.7470155743891863*u.rad,
    lon=0.10510862459055528*u.rad,
    height=0.0*u.m
)

d = Time(["2028-11-08T23:34:54.00"], scale="utc")

print(get_body("moon", d, location))
WARNING: Tried to get polar motions for times after IERS data is valid. Defaulting to polar motion from the 50-yr mean for those. This may affect precision at the arcsec level. Please check your astropy.utils.iers.conf.iers_auto_url and point it to a newer version if necessary. [astropy.coordinates.builtin_frames.utils]

<SkyCoord (GCRS: obstime=['2028-11-08T23:34:54.000'], obsgeoloc=[(3144039.81787578, 3484619.75485466, 4300693.72117526)] m, obsgeovel=[(-254.1002625, 228.37771983, 0.71891401)] m / s): (ra, dec, distance) in (deg, deg, AU)
    [(127.74528269, 16.71982861, 0.00250404)]>

So we get

RA = 127.74528269
DEC = 16.71982861

Note that our moon function above yields coordinates roughly 1deg off due to the parallax of the Moon (it’s in the geocentric system and not taking account any location on Earth)

julia> dt = DateTime("2028-11-08T23:34:54.00")
2028-11-08T23:34:54

julia> pos = moon(dt)
(ra = 2.2167352184921754, dec = 0.30209193883117474, azimuth = 1.4716666159254528, elevation = 0.34136871844980865)

julia> rad2deg.([pos.ra, pos.dec])
2-element Vector{Float64}:
 127.00957231761205
  17.30859311995054

Now comes my function (see below) which tries to do everything right, but doesn’t, and I can’t figure out where the problem is. I tried to comment everything so that it’s maybe a bit easier to spot the error.

This is what I get

julia> pos = moon_topocentric(lon, lat, dt)
(ra = 2.21935735847888, dec = 0.299922564561618, distance_au = 0.0025181965734329967)

julia> rad2deg.([pos.ra, pos.dec])
2-element Vector{Float64}:
 127.15980987214272
  17.18429713012066

Which is still off by around half a degree or so!

Does anybody know what’s wrong here? I got great help two years ago from @DaveMacMahon here Calculate positions of astronomical objects from a terrestrial observatory where I calculated the Sun coordinates but that did not involve any parallax corrections.

I currently think that the nutation/precession correction is the last missing piece but cannot get it right.

Here is the function I am sitting at since yesterday

using Astrometry.SOFA
using Dates

function moon_topocentric(lon, lat, dt::DateTime; height=0.0)
    # UTC as JD
    utc1, utc2 = datetime2julian(dt), 0.0
    
    # Earth orientation (use IERS values for high precision)
    dut1 = 0.0
    xp = 0.0
    yp = 0.0
    
    # Disable refraction, we are observing neutrinos so we don't care right now ;)
    phpa = 0.0   # pressure (hPa)
    tc   = 0.0   # temperature (C)
    rh   = 0.0   # relative humidity
    wl   = 0.0   # wavelength (micron)
    
    # Build astrometry context
    astrom, eo = SOFA.apco13(
        utc1, utc2,
        dut1,
        lon, lat, height,
        xp, yp,
        phpa, tc, rh, wl
    )
    
    # TT for Moon ephemeris
    tai1, tai2 = SOFA.utctai(utc1, utc2)
    tt1, tt2   = SOFA.taitt(tai1, tai2)
    
    # Moon geocentric GCRS position (AU and AU/day)
    pm, vm = SOFA.moon98(tt1, tt2)
    
    # Calculate distance in AU
    distance_au = sqrt(pm[1]^2 + pm[2]^2 + pm[3]^2)
    
    # Parallax in arcseconds: π = 1 AU / distance
    parallax_arcsec = 1.0 / distance_au
    
    # Cartesian -> spherical (GCRS)
    ra_gc, dec_gc = SOFA.c2s(pm)
    
    # Calculate radial velocity (AU/day)
    rv_au_day = (pm[1]*vm[1] + pm[2]*vm[2] + pm[3]*vm[3]) / distance_au
    # Convert to km/s (1 AU/day ≈ 1731.46 km/s)
    rv_km_s = rv_au_day * 1731.46
    
    # GCRS -> topocentric apparent place with parallax
    # Note: For Moon, proper motion during the observation is negligible
    # compared to parallax effect, so we can use 0.0 for pmra/pmdec
    ra_topo, dec_topo = SOFA.atciq(
        ra_gc, dec_gc,
        0.0, 0.0,           # proper motion (negligible for Moon)
        parallax_arcsec, 
        rv_km_s,            # radial velocity (km/s)
        astrom
    )
    
    return (
        ra  = mod2pi(ra_topo),
        dec = dec_topo,
        distance_au = distance_au
    )
end

tamasgal · February 13, 2026, 1:29pm

I added the alt/az calculation and compared to AstroPy by using AstroPy to generate a dataset of 24 hours.

What you can see is the moon function (geocentric) which ignores the Earth location and the slightly extended moon_topocentric function (topocentric). The differences are periodic and it seems to me that I am not at the right “location” on Earth.

The azimuth and right ascension are correct twice a day

Here is the current code

function moon_topocentric2(lon, lat, dt::DateTime; height=0.0)
    # UTC as JD
    utc1, utc2 = datetime2julian(dt), 0.0
    
    # Earth orientation (use IERS values for high precision)
    dut1 = 0.0
    xp = 0.0
    yp = 0.0
    
    # Disable refraction
    phpa = 0.0   # pressure (hPa)
    tc   = 0.0   # temperature (C)
    rh   = 0.0   # relative humidity
    wl   = 0.0   # wavelength (micron)
    
    # Build astrometry context
    astrom, eo = SOFA.apco13(
        utc1, utc2,
        dut1,
        lon, lat, height,
        xp, yp,
        phpa, tc, rh, wl
    )
    
    # TT for Moon ephemeris
    tai1, tai2 = SOFA.utctai(utc1, utc2)
    tt1, tt2   = SOFA.taitt(tai1, tai2)
    
    # Moon geocentric GCRS position (AU and AU/day)
    pm, vm = SOFA.moon98(tt1, tt2)
    
    # Calculate distance in AU
    distance_au = sqrt(pm[1]^2 + pm[2]^2 + pm[3]^2)
    
    # Parallax in arcseconds: π = 1 AU / distance
    parallax_arcsec = 1.0 / distance_au
    
    # Cartesian → spherical (GCRS)
    ra_gc, dec_gc = SOFA.c2s(pm)
    
    # Calculate radial velocity (AU/day)
    rv_au_day = (pm[1]*vm[1] + pm[2]*vm[2] + pm[3]*vm[3]) / distance_au
    # Convert to km/s (1 AU/day ≈ 1731.46 km/s)
    rv_km_s = rv_au_day * 1731.46
    
    # GCRS → topocentric apparent place with parallax
    # Note: For Moon, proper motion during the observation is negligible
    # compared to parallax effect, so we can use 0.0 for pmra/pmdec
    ra_topo, dec_topo = SOFA.atciq(
        ra_gc, dec_gc,
        0.0, 0.0,           # proper motion (negligible for Moon)
        parallax_arcsec,    # parallax (arcsec) - THIS IS THE KEY!
        rv_km_s,            # radial velocity (km/s)
        astrom
    )
    
    moon_dist_au = sqrt(sum(pm.^2))
    # 180 × 3600 / π, for reference ;)
    parallax = 1.0 / (moon_dist_au * 206264.806247)  # parallax in arcseconds
    moonobs = SOFA.atco13(ra_topo, dec_topo, 0.0, 0.0, parallax, 0.0, utc1, utc2, 0.0, lon, lat, height, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

    az, el = moonobs.azi, π/2 - moonobs.zen

    return (ra=mod2pi(ra_topo), dec=dec_topo, distance_au = distance_au, azimuth=az, elevation=el)
end

cormullion · February 13, 2026, 2:02pm

Can NASA software help?

mihalybaci · February 13, 2026, 2:25pm

Piggybacking off this, I would say that Python’s skyfield package is much better designed for ephemeris calculations than even Astropy is. It might be good to take inspiration from that package, it says it is “pure python” so the calculations may be exposed in a way you could translate to Julia.

tamasgal · February 13, 2026, 4:51pm

Thanks @mihalybaci I will have a look but as a side note; we also have a Slalib based package which agrees with the dataset generated by AstroPy within 0.0001 degree. The fact that I cannot get the Moon right within 0.7deg is very likely a user mistake on my side

tamasgal · February 13, 2026, 4:58pm

I will check this out as well, as a fourth reference

mihalybaci · February 13, 2026, 5:29pm

Yeah, it’s tricky and I am certainly not well-versed in planetary motion. If I had to guess, and I don’t, I would take a long look at the lunar parallax again as you already mentioned. In your code, the parallax is simply a function of the distance, but the parallax will be a function of time and location on the Earth as well. Your code may already take that into account and I missed it.

Edit: It looks like SPICE.jl may be the current best way to implement this kind of calculation, unless you do want/need to roll your own.

Topic		Replies	Views
Calculate positions of astronomical objects from a terrestrial observatory Astro/Space question	8	983	December 20, 2023
Sun elevation and moon phase Astro/Space question	6	1377	January 16, 2024
[ANN] EarthOrientation.jl Astro/Space announcement	8	1863	August 28, 2020
Using Julia to calc distance to moon with UNCERTAINTIES Offtopic	2	984	September 6, 2019
Sun's trajectory around earth General Usage	19	1438	October 20, 2025

Position of the Moon (topocentric) using Astrometry.jl (SOFA/ERFA) routines

Related topics