Earth position does not match that of JPL Horizons

[simas]

Hello, I am trying to extract Earth’s ICRF position at a specific Julian Date on the TDB time scale. For validation, I am comparing my result with JPL Horizons. The problem is that I am unable to match the values exactly. From DE441 ephemerides Horizons produces the following state:

$$SOE
2460157.500000000 = A.D. 2023-Aug-01 00:00:00.0000 TDB [del_T= 69.183312 s]
XYZ : 9.262412495155303E+07 -1.097336742893849E+08 -4.753366571502377E+07 2.291818900912849E+01 1.678010135905066E+01 7.273858345769444E+00
sigmas: n.a. n.a. n.a. n.a. n.a. n.a.
$$EOE

I am using the default orekit-data package. It contains DE-441-ephemerides. To make sure that correct ephemerides are used, I create the earth CelestialBody in the following way:

JPLEphemeridesLoader earthLoader = new JPLEphemeridesLoader("lnxp1990.441", JPLEphemeridesLoader.EphemerisType.EARTH);
CelestialBody earth = earthLoader.loadCelestialBody("Earth");

I then create my AbsoluteDate in the following way:

TDBScale TDBscale = TimeScalesFactory.getTDB();
AbsoluteDate date = AbsoluteDate.createJDDate(2460157, Constants.JULIAN_DAY/2.0d, TDBscale);

I also create the ICRF frame like this:

CelestialBody centralBody = CelestialBodyFactory.getSolarSystemBarycenter();
Frame inertialFrame = centralBody.getInertiallyOrientedFrame();

Finally, I extract the PVCoordinates of earth at JDTDB 2460157.5 like follows:

earth.getPVCoordinates(date,inertialFrame)

I get the following result:
{2023-07-31T23:58:50.8167034179512,
P(9.262412495125844E10, -1.0973367428959982E11, -4.75336657151169E10),
V(22918.18900917888, 16780.101358989006, 7273.858345741602),
A(-0.003542073372692932, 0.0041202768855387335, 0.001782459856161173)}

There are a few discrepancies that I’m not able to resolve. First, the discrepancy in position coordinates, for example, the x coordinates differ by around 30 cm:
horizons: 92624124951.55303 m
orekit: 92624124951.25844 m

Secondly, there is a discrepancy in the AbsoluteDate, which, I suspect, could be causing the discrepancy in the position. If I extract the TDB date using

date.toString(TDBscale)

I get 2023-07-31T23:59:59.99998730725312, which is not exactly 2023-08-01T00:00:00.000 as I would expect and as also seen in the Horizons output. However, even if some error is introduced during conversion of JDTDB to the date time format, it seems strange that it would cause an error in the ephemeride readings, since the ephemerides directly require the JDTDB date as input.

In conclusion, I cannot figure out why the discrepancy is there, perhaps someone has an idea? Thanks in advance.

[MaximeJ]

Hi @simas,

Welcome to the forum!

@simas

Secondly, there is a discrepancy in the AbsoluteDate, which, I suspect, could be causing the discrepancy in the position.

I think you’re right.
The discrepancy in date is 1.26927E-05 seconds, which, multiplied by Earth velocity from JPL Horizons (29321.0539427297 m/s), gives a discrepancy in position of 0.372164717 m.
This is consistent with the norm of the position error that you observe (0.376430341 m).

Could you try building the AbsoluteDate with

AbsoluteDate date = new AbsoluteDate("2023-08-01T00:00:00.000", TDBscale);

And tell us if you get a better match between Orekit and JPL Horizons ?

As for why the conversion from Julian date to AbsoluteDate in TDB introduces an error, I have no clue for now.
It does seem to come from the value of the secondsSinceNoon parameter in AbsoluteDate.createJDDate method. If it’s zero, we have a perfect match between calendar and JD date, but with secondsSinceNoon growing the error between the 2 dates grows. From 1.06E-06s at 1PM to 2.43E-05s at 11AM the next day.
Note that it happens for time scales with a bias that is not constant with TAI (i.e. TDB, TCB, TT…).
I would say it’s a bug but I’d like to wait for date experts (@luc and @evan.ward for example) to give their opinion to be sure.

Edit: the date difference is the difference of offset between TAI and TDB taken at noon (2023-07-31T12:00:00.000 TDB) and midnight (2023-08-01T00:00:00.000 TDB).

It comes from the construction of the AbsoluteDate from Julian date in Orekit.

public static AbsoluteDate createJDDate(final int jd, final double secondsSinceNoon, final TimeScale timeScale) {
        return new AbsoluteDate(new DateComponents(DateComponents.JULIAN_EPOCH, jd),
                                TimeComponents.H12, timeScale).shiftedBy(secondsSinceNoon);
    }

Here the offset from TAI is taken at noon (TimeComponents.H12), then the date is shifted by 12 * 3600 seconds.
While doing new AbsoluteDate("2023-08-01T00:00:00.000", TDBscale); takes the offset at the proper date (so exactly 2023-08-01T00:00:00.000 in TDB).

Maxime

[simas]

Hi @MaximeJ , thank you for your answer.

Creating the date with

new AbsoluteDate("2023-08-01T00:00:00.000", TDBscale)

does indeed produce a much better match. For Earth’s ICRF position I got:
2023-07-31T23:58:50.81671611069808,
P(9.262412495154932E10, -1.0973367428938684E11, -4.753366571502458E10)

However, it still differs from Horizons by about 5 mm. This is probably good enough for most cases, but I still wonder why I do not get an exact match. In my case, my Orekit app and another application I’m building need to agree exactly on the situation at a certain epoch and that’s why I am trying to get a perfect match.
When reading the ephemerides with the other application, I am able to get the exact same result as Horizons simply by inputting the exact same Julian Date. I wish there was a way to do that with my Orekit app, since it needs to receive a Julian Date and start a propagation from that epoch.
It seems that now my best bet is to convert the Julian Date to a calendar date format and then create the AbsoluteDate from that. However, I would still like to ask about the remaining difference in position. Perhaps it’s possible to get an exact match?

Thank you for your time.

[MaximeJ]

Hi @simas,

@simas

When reading the ephemerides with the other application, I am able to get the exact same result as Horizons simply by inputting the exact same Julian Date. I wish there was a way to do that with my Orekit app, since it needs to receive a Julian Date and start a propagation from that epoch.

In my opinion, the issue you experienced could be considered a bug.
Could you please open an issue on the Orekit forge about it?

I think I found a way to get a date that is close using directly the MJD.
Could you try using this function in your code?

public static AbsoluteDate createJDDate(final int jd, final double secondsSinceNoon,
                                        final TimeScale timeScale,
                                        final TimeScale pivotTimeScale) {
        final AbsoluteDate dateInPivotTimeScale =  new AbsoluteDate(new DateComponents(DateComponents.JULIAN_EPOCH, jd),
                                                                    TimeComponents.H12,
                                                                    pivotTimeScale).shiftedBy(secondsSinceNoon);
        final double offsetFromTAI = timeScale.offsetFromTAI(dateInPivotTimeScale) -  pivotTimeScale.offsetFromTAI(dateInPivotTimeScale);
        return dateInPivotTimeScale.shiftedBy(-offsetFromTAI);
    }

The idea is to introduce a pivot time scale that is close to the target timescale but has a constant bias with TAI.
In your case, you would use the TT time scale as a pivot time scale to get the date in TDB (since TT and TDB are less than 1 ms apart).
That way the difference between AbsoluteDate from MJD and calendar date should be around 1e-13s.
I’d like to know what is the difference between Orekit and Horizons using this function.

@simas

However, I would still like to ask about the remaining difference in position. Perhaps it’s possible to get an exact match?

That I wouldn’t know. If this is still a date issue, a 5mm bias in position means we still have a date bias of around 0.15 µs.
What is the other application you’re using?

Maxime

[simas]

Hi @MaximeJ .

Sure, I will open an issue today.

I tried your method and it works. The date I created with it prints as:
2023-07-31T23:59:59.99999999999979
so, indeed it’s only 2e-13 seconds off.

The Earth position I get with it is identical to the one we get when using the exact date, created with

AbsoluteDate date = new AbsoluteDate("2023-08-01T00:00:00.000", TDBscale);

That position in meters is
92624124951.54932, -109733674289.38684, -47533665715.02458
Position from Horizons:
92624124951.55303, -109733674289.3849, -47533665715.02377

So, the 5mm discrepancy remains. The other application I’m using can be found here. I’m using the C# binary Development Ephemeris reader found in there. That produces:
92624124951.553, -109733674289.385, -47533665715.0238
I have no idea where the remaining difference comes from, but I don’t think it’s a date issue anymore since both the exact date and the workaround date produce this result and the workaround can be seen to be 2e-13 sec off, which is much less than 0.15 µs.

Thanks for your help!

[MaximeJ]

Hi @simas,

@simas

Sure, I will open an issue today.

Thanks.

@simas

I tried your method and it works.

Great!

@simas

So, the 5mm discrepancy remains. The other application I’m using can be found here. I’m using the C# binary Development Ephemeris reader found in there. That produces:
92624124951.553, -109733674289.385, -47533665715.0238
I have no idea where the remaining difference comes from, but I don’t think it’s a date issue anymore since both the exact date and the workaround date produce this result and the workaround can be seen to be 2e-13 sec off, which is much less than 0.15 µs.

Thanks for the link to the other application you mentioned. There’s a bridge to Java so maybe we could try to compare the results with Orekit and see where the 5mm discrepancy comes from.
Could you please open a second issue related to this discrepancy between Orekit and JPL Horizons?

@simas

Thanks for your help!

You’re welcome.

Maxime