Possible issue with Moon positions

@mkbrewer

@mkbrewer, @dhomeier - I'm moving our discussion in #11753 to a new issue, since I convinced myself that at least it has nothing to do specifically with the use of moon98. But there does seem something weird with the Moon. Somewhat similarly to @dhomeier's worry, I had just been checking what happens if I adjust the time a bit, and I find the match is significantly better for a time offset of 0.65 s. It is independent of whether I use the new moon98 or our own implementation, and does not seem to depend on time span I try to cover. Note that there is no such offset for body='earth'.

This is my test programme:

import numpy as np
from astropy.time import Time
from astropy import units as u
from astropy.coordinates import get_body_barycentric_posvel, get_body_barycentric

t = Time('1990-03-20T12:30:00') + np.linspace(0, 20*u.yr, 800)
time_diffs = np.linspace(-1*u.s, 1.*u.s, 21)
body = 'moon'
dp, dv = get_body_barycentric_posvel(body, t, ephemeris='de432s')
bp = get_body_barycentric(body, t, ephemeris='builtin')
for time_diff in time_diffs:
    pos_diff = (bp + time_diff * dv - dp).norm().to(u.km)
    print(f'Δt={time_diff:5.2f}, rms={np.sqrt((pos_diff**2).mean()):.2f}')

(Note I use the velocity from JPL so I can also check current master - the velocities are close enough it does not matter much.)

Output

Δt=-1.00 s, rms=53.25 km
Δt=-0.90 s, rms=50.52 km
Δt=-0.80 s, rms=47.82 km
Δt=-0.70 s, rms=45.16 km
Δt=-0.60 s, rms=42.54 km
Δt=-0.50 s, rms=39.96 km
Δt=-0.40 s, rms=37.45 km
Δt=-0.30 s, rms=35.02 km
Δt=-0.20 s, rms=32.67 km
Δt=-0.10 s, rms=30.43 km
Δt= 0.00 s, rms=28.33 km
Δt= 0.10 s, rms=26.41 km
Δt= 0.20 s, rms=24.69 km
Δt= 0.30 s, rms=23.22 km
Δt= 0.40 s, rms=22.07 km
Δt= 0.50 s, rms=21.27 km
Δt= 0.60 s, rms=20.87 km
Δt= 0.70 s, rms=20.90 km
Δt= 0.80 s, rms=21.34 km
Δt= 0.90 s, rms=22.18 km
Δt= 1.00 s, rms=23.37 km

From the dates in #11753 (comment) it looks quite evident that the distance offsets are highest at first quarter, close to 0 around full/new Moon, and probably smallest (negative) at last quarter (the first date is only a few day before FM).

Thinking on this further, in example 47.a in his book, Meeus uses his algorithm to calculate the apparent position of the Moon directly, so his algorithm must include the effect of light time delay on the Moon's orbit around the Sun. This makes me wonder why SOFA would leave out the effect of light time delay on the Moon's longitude. It would be interesting to look at the distance error as a function of the Moon's elongation from the Sun. When the Moon is leading the Earth in its orbit around the Sun, it appears closer to the Earth than when it is trailing due to the motion of the Earth during the light travel time. If, as I suspect, this motion is taken into account in the algorithm, then adding in the light time delay in _get_apparent_body_position() will exacerbate the situation by effectively adding in the Earth's motion twice.

distance error as a function of the Moon's elongation from the Sun. When the Moon is leading the Earth in its orbit around the Sun, it appears closer to the Earth than when it is trailing due to the motion of the Earth during the light travel time. If, as I

Yes, that would match the positive distance offsets I mentioned above, largest around 2017-03-19 (first quarter) and smaller towards 2016-03-20 (some days before FM).

Similarly, adding in the light time delay twice, should lead to approximately +/-15" error in RA of the full and new Moon.

@mkbrewer

@mkbrewer - yes, you are right, adding the light travel delay causes an increase in RA, Dec error. Amusingly, it also decreases the error in distance, to really quite minimal. New test routine:

import numpy as np
from astropy.time import Time
from astropy import units as u
from astropy.coordinates import get_body


t = Time('2016-03-20T12:30:00') + np.arange(80.) * 10 * u.day
body = 'moon'
builtin = get_body(body, t, ephemeris='builtin')
de432 = get_body(body, t, ephemeris='de432s')

print(f'rms_ra={(builtin.ra - de432.ra).std().to(u.arcsec)}\n'
      f'rms_dec={(builtin.dec - de432.dec).std().to(u.arcsec)}\n'
      f'rms_distance={(builtin.distance - de432.distance).std()}')

With current master (or simple replacement of calc_moon with erfa.moon84):

rms_ra=3.0923007201476302 arcsec
rms_dec=1.5814318611792149 arcsec
rms_distance=28.42713318597424 km

With the replacement including light travel time correction:

rms_ra=14.78075407708036 arcsec
rms_dec=5.091270306592033 arcsec
rms_distance=3.7637321170709894 km

EDIT: For completeness, with the straight moon98 replacement (i.e., assuming light travel time is included already), I do get slightly better agreement than with our own implementation on master:

rms_ra=3.0794952046457307 arcsec
rms_dec=1.4965303053714698 arcsec
rms_distance=28.42473375418744 km

Well, that's curious. I wonder if the time delay is included in the calculation of longitude and latitude, but not distance. That would seem rather inconsistent.

I think in order to get to the bottom of this, we'd need plots of the difference in distance and right ascension as a function of the elongation of the Moon (difference between ecliptic longitude of the Sun and Moon) both with and without light time delay added.

I'd suggest doing this in one hour intervals over the period of one lunar month.

Weird indeed. I can get very good agreement by just removing earth motion from the distance, by adjusting calc_moon from #11753 as follows (and then not applying any further light travel time corrections):

    pv = erfa.moon98(*get_jd12(t, 'tt'))
    pos = CartesianRepresentation(pv['p'] << u.AU, unit=u.km, xyz_axis=-1)
    vel = CartesianRepresentation(pv['v'] << u.AU/u.day, xyz_axis=-1)
    # Adjust distance for light travel time effect due to Earth moving,
    # which appears to have been missed in the Meeus algorithm.
    light_travel_time = pos.norm() / speed_of_light
    _, earth_vel = _get_body_barycentric_posvel('earth', t, ephemeris='builtin')
    pos_normal = pos / pos.norm()
    pos -= light_travel_time * earth_vel.dot(pos_normal) * pos_normal
    # Create GCRS sky coordinate.
    gcrs_rep = pos.with_differentials(
        {'s': CartesianDifferential.from_cartesian(vel)})
    return SkyCoord(GCRS(gcrs_rep, obstime=t))

With that, the rms differences I get are:

rms_lon=3.031427731940901 arcsec
rms_lat=1.359440972999706 arcsec
rms_distance=3.7794434832812143 km

p.s. Yes, plots would probably be a better idea.

OK. I made some plots that reveal the problem. First I computed the ICRS position of the Moon using the JPL DE405 ephemeris (the reference for ELP/MPP02) without any time delay and subtracted the position of the Earth. This is merely the position given by ephemeris and can be considered as the astrometric position wrt the geocenter at the time that the light was emitted. Then I compared this to the output of Moon98() again without any time delay. Here are the results for the difference in RA/Dec and distance over one lunar month:

Then I added the light travel time and aberration to yield the GCRS position for the JPL result and compared this to the light time delayed output of Moon98(). Here are those results:

While the difference in distance is much larger due to the change in the position of the Earth during the light travel time, the difference in RA/Dec is the same to within 10 mas.

This indicates that annual aberration exactly cancels the component of planetary aberration due to the motion of the Earth. What's left is merely the component of planetary aberration due to the change in the position of the Moon. Since this is taken into account in both the time delayed JPL and Moon98() positions, it also cancels in the difference.

The upshot of this is that in treating the output of Moon98() as a GCRS postion and then transforming to ICRS, you are removing annual aberration, but not the component of planetary aberration from the motion of the Earth during the light travel time which is already built in to Moon98(). Then when you add in the light time delay and transform back to GCRS, you are adding in this planetary aberration twice while only adding back annual aberration once.

The solution is to treat the output of Moon98() as an astrometric position wrt the geocenter, then just add in the position of the Earth wrt the SSB at the time of observation to get the equivalent ICRS ephemeris position. Then adding in the time delay will add planetary aberration, adjust for apparent distance, and transforming the result to GCRS will add in annual aberration canceling the planetary aberration due to the motion of the Earth. This is why, in their paper, The Lunar Ephemeris ELP 2000, referenced by Meeus, Chapront-Touzé and Chapront do a direct comparison of their ephemeris with JPL LE200.

Btw, if you want to reproduce these plots, I used the SOFA function iauFad03() for the elongation of the Moon used as the X axis.

@mkbrewer

@mkbrewer - thanks so much, that makes a lot of sense! And your suggestion for the implementation will be straightforward to use.

Just for completeness in case someone stumbles on this in the future: #11753 now uses moon98 correctly and will thus fix this issue.

@mhvk

@mhvk I believe that this issue can now be closed.

Navigation Menu

Repository navigation

Possible issue with Moon positions #11755

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Possible issue with Moon positions #11755

Description

Activity

dhomeier commented on May 18, 2021

mkbrewer commented on May 18, 2021

dhomeier commented on May 18, 2021

mkbrewer commented on May 18, 2021

mhvk commented on May 18, 2021

mkbrewer commented on May 18, 2021

mkbrewer commented on May 18, 2021

mkbrewer commented on May 18, 2021

mhvk commented on May 18, 2021

mhvk commented on May 18, 2021

mkbrewer commented on May 19, 2021

mkbrewer commented on May 20, 2021

mhvk commented on May 20, 2021

mhvk commented on May 20, 2021

mkbrewer commented on May 21, 2021

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions