Moon position comparison: VSOP87 / Horizons

[StacyMader]

I have been looking at the Moon position from PyEphem and comparing to Horizons. Unless I am doing something really dumb here, the differences in the Astrometric position relative to an observer are quite large.

Attached are difference plots of the RA and DEC Astrometric positions for the Moon from 2022-04-01, 00:00:00 to 05:38:00 UTC. The PyEphem plots are a lot more than 1 arc-second accuracy. AstroPy is bang on, though I am using DE432s.

I also add my PyEphem script below and the Horizons output. In both cases, I have disabled refraction calculations by setting pressure to zero.

Any help appreciated.

Thanks,
Stacy.

#!/usr/bin/python3

from astropy.time import Time
import astropy.units as u
import math
import ephem
from matplotlib import pylab as plt

home = ephem.Observer()
home.lon =  ephem.degrees('148.2635192')
home.lat = ephem.degrees('-32.9983897')
home.elevation = 414.666
home.horizon = ephem.degrees('30.25')
home.pressure = 0

tmin = Time("2022-04-01 00:00:00")

moon = ephem.Moon()
for hour in range(0,24):
   h = tmin + hour*u.hour
   for m in range(0,60):
      t = h + m*u.min
      home.date = ephem.Date(t.datetime)
      moon.compute(home)
      if math.degrees(moon.alt) > 30.25:
         print("{} {:.3f} {:.3f} {} {}".format(home.date, math.degrees(moon.a_ra), math.degrees(moon.a_dec), math.degrees(float(moon.az)), math.degrees(float(moon.alt))))

[brandon-rhodes]

Thanks for providing a sample script, but I'm not sure it cut-and-pasted cleanly? There's no code to call the function. And, it tries to use a name pks in moon.compute(pks) that's not defined anywhere in the program.

[StacyMader]

Apologies, updated now to run as python3 script.py. Also running v4.1.3 on macOS 12.3 (installed with pip3.)

[brandon-rhodes]

Thanks for the working script! What are the units in your x-y plot, so that I know how big a difference is being shown?

[StacyMader]

Oops, the y-axis is the decimal difference between the RA and DEC coordinates for the cases Horizons-PyEphem, Horizons-AstroPy. The x-axis has a point for everyminute from 00:00:00 to 05:38:00 UT- when the Moon is visible above 30.25 degrees. I'll draw up better plots now.

[brandon-rhodes]

Oops, the y-axis is the decimal difference between the RA and DEC coordinates…

With the two quantities expressed in degrees, or is RA in hours? (Edit: I'm guessing hours, from the script, now that I re-read it.)

[StacyMader]

[brandon-rhodes]

@StacyMader — Thanks for the updated plots!

Per the quick reference, PyEphem's a_ra and a_dec are measured from the center of the Earth — they are geocentric RA and dec. Is it possible that Horizons and your AstroPy-powered script are instead computing RA and dec from the home location instead of the center of the Earth?

[StacyMader]

@brandon-rhodes , Horizons header below. I take this to mean the coordinates are for the home location, not Earth centre?

*******************************************************************************
Ephemeris / WWW_USER Sun Apr  3 06:11:26 2022 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: Moon (301)                      {source: DE441}
Center body name: Earth (399)                     {source: DE441}
Center-site name: (user defined site below)
*******************************************************************************
Start time      : A.D. 2022-Apr-01 00:00:00.0000 UT      
Stop  time      : A.D. 2022-Apr-02 00:00:00.0000 UT      
Step-size       : 1 minutes
*******************************************************************************
Target pole/equ : MOON_ME                         {East-longitude positive}
Target radii    : 1737.4 x 1737.4 x 1737.4 km     {Equator, meridian, pole}    
Center geodetic : 148.263519,-32.998390,0.4146660 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 148.263519,5354.91983,-3454.035 {E-lon(deg),Dxy(km),Dz(km)}
Center pole/equ : ITRF93                          {East-longitude positive}
Center radii    : 6378.1 x 6378.1 x 6356.8 km     {Equator, meridian, pole}    
Target primary  : Earth
Vis. interferer : MOON (R_eq= 1737.400) km        {source: DE441}
Rel. light bend : Sun, EARTH                      {source: DE441}
Rel. lght bnd GM: 1.3271E+11, 3.9860E+05 km^3/s^2                              
Atmos refraction: NO (AIRLESS)
RA format       : DEG
Time format     : CAL 
RTS-only print  : NO       
EOP file        : eop.220331.p220624                                           
EOP coverage    : DATA-BASED 1962-JAN-20 TO 2022-MAR-31. PREDICTS-> 2022-JUN-23
Units conversion: 1 au= 149597870.700 km, c= 299792.458 km/s, 1 day= 86400.0 s 
Table cut-offs 1: Elevation ( 30.2deg=YES),Airmass (>38.000=NO), Daylight (NO )
Table cut-offs 2: Solar elongation (  0.0,180.0=NO ),Local Hour Angle( 0.0=NO )
Table cut-offs 3: RA/DEC angular rate (     0.0=NO )                           

I read here about a_ra, a_dec:

PyEphem now knows the “star-atlas position” of the body, and checks for which star-atlas epoch you want coordinates expressed in ...

So I took this to say a_ra and a_dec are J2000, and given the closeness of J2000 to ICRF (minus the small rotation in axes), thought they should be the same.

Instead of a_ra/a_dec, I plotted ra and dec (the Apparent Topocentric Position), with still an appreciative offset for PyEphem. However, I note the Az/Alt for both Horizons and PyEphem are pretty much the same (to ~ 2 decimal places).

[StacyMader]

I attach my Python script for AstroPy, see below. The home location in EarthLocation and Horizons is Geodetic, not geocentric, though I do note AstroPy stores the location internally as Geocentric.

#!/usr/bin/python3

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

home = EarthLocation(lon=148.2635192 * u.deg, lat=-32.9983897 * u.deg, height=414.666 * u.m)
tmin = Time("2022-04-01 00:00:00")

with solar_system_ephemeris.set('de432s'):
    for hour in range(0,24):
        h = tmin + hour*u.hour
        for m in range(0,60):
            t = h + m*u.min
            moon = get_body('moon', t, home)
            altAzFrame = AltAz(location=home, obstime=t, pressure=0)
            altAz = moon.transform_to(altAzFrame).to_string('decimal').split()
            if float(altAz[1]) > 30.25:
                print("{} {:.3f} {:.3f} {} {}".format(t, moon.ra, moon.dec, altAz[0], altAz[1]))

[brandon-rhodes]

I want to say:

I may not have been clear in my earlier comment: I think that a_ra and a_dec are computing where in the sky the Moon is if viewed from the center of the Earth. Whereas ra and dec instead return the different position that the Moon has when the Moon is viewed from the location on Earth that you specify with lat/lon.

Since the Moon is only about 60 Earth radii away, a given lat/lon on Earth's surface will, as the Earth rotates, see the Moon first from one side and then the other side of the Earth's center. This moves the Moon back and forth a bit against the stars, creating the difference in RA and dec that you are noticing.

But instead I'm going to ask:

What happened to the sine-wave like shape of the dec difference in the earlier diagrams? Now it looks like a flat line, which doesn't agree with the theory I just typed up — actually, I guess they should both be a sine wave. Hmm. I might need to really take time to look at this the next time I have the chance to do a deeper dive than I have time for this evening.

[brandon-rhodes]

@StacyMader — It looks like I might have a little time this month for investigating old PyEphem issues, so I wanted to check back about where things stand with this issue. Did you ever, in the meantime, resolve the question about the difference you were seeing between the output of PyEphem and the output of your other tools?

[StacyMader]

Sorry @brandon-rhodes I've not had a chance to look at this again, but I did not resolve the difference.