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
