Asteroid propagation not matching MPC results

Orekit usage
1

Hello fellow Orekit users!

I have been recently working on a small program, that will propagate asteroid orbit from MPCORB database data and then give telescope pointing.

At first, I wanted to compare heliocentric state vectors after propagation, if they are the same. X value is always matching MPC, then Y value deviates, and Z value is completely not matching those calculated by MPC.
I have prepared simple Junit test.

What I do is :

VECTORS: Heliocentric vectors/AU
Obj JD_TT X Y Z
00825 2458209.5000000 0.234729353 2.135914513 0.881546394

From Orkeit code:
0.2344986504264944 2.310074947284788 -0.040804051106836874

Clearly X axis value is correct, and Z totally different.
I wonder, what other coordinate frames could be used by MPC ? Is it something wrong with my code?

This is my code:

public class TestEphem {
public static File orekitData = new File("src/test/resources/orekit-data");

@Test
public void computeVectorTest() throws OrekitException {
	DirectoryCrawler crawler = new DirectoryCrawler(orekitData);
	System.out.println(orekitData.getAbsolutePath());
	DataProvidersManager.getInstance().addProvider(crawler);
	// sun body and sun-centric inertial frame

	CelestialBody sunBody = CelestialBodyFactory.getSun();

	// epoch of keplerian coordinates

	Double epoch = 2458600.5;

	int julianDay = epoch.intValue();

	double secondsSinceNoon = (epoch - julianDay) * 24 * 3600;

	AbsoluteDate date = AbsoluteDate.createJDDate(julianDay, secondsSinceNoon, TimeScalesFactory.getUTC());

	// semi-major axis [m]

	double a = 2.225774 * Constants.IAU_2012_ASTRONOMICAL_UNIT;

	// orbital eccentricity

	double e = 0.0749275;

	// inclination [rad]

	double i = FastMath.toRadians(3.39960);

	// Argument of perihelion [rad]

	double peri = FastMath.toRadians(111.78637);

	// Longitude of the ascending node [rad]

	double node = FastMath.toRadians(101.40875);

	// Mean anomaly [rad]

	double m = FastMath.toRadians(353.95770);

	Orbit orb = new KeplerianOrbit(a, e, i, peri, node, m, PositionAngle.MEAN, sunBody.getInertiallyOrientedFrame(),

			date, sunBody.getGM());

	NumericalPropagator prop = getPropagator(orb);

	AbsoluteDate dat = new AbsoluteDate(2018, 04, 01, 0, 0, 0.0, TimeScalesFactory.getUTC());

	TimeStampedPVCoordinates cord = prop.propagate(dat).getPVCoordinates();

	Vector3D pos2 = cord.getPosition();

	System.out.println(pos2.getX() / Constants.IAU_2012_ASTRONOMICAL_UNIT + " "

			+ pos2.getY() / Constants.IAU_2012_ASTRONOMICAL_UNIT + " "

			+ pos2.getZ() / Constants.IAU_2012_ASTRONOMICAL_UNIT);
}

public NumericalPropagator getPropagator(Orbit orb) throws OrekitException {
	double[][] tolerances = NumericalPropagator.tolerances(1, orb, OrbitType.CARTESIAN);
	AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 3600, tolerances[0], tolerances[1]);
	NumericalPropagator propagator = new NumericalPropagator(integrator);
	propagator.setOrbitType(OrbitType.CARTESIAN);

	propagator.addForceModel(new NewtonianAttraction(orb.getMu()));
	propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getEarth()));
	propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
	propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getJupiter()));
	propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMars()));
	propagator.setInitialState(new SpacecraftState(orb));

	return propagator;
}

}

I have used newest EOP, planetary ephem and leap secodns files.

I would be very grateful for any hints provided! I would like to learn using Orekit. Propagation was looking as quite simple task, but here I am clearly getting something wrong.

Best regards,

2

Hello @AstroUser,

You code seems fine.
There is definitely an issue with the frames that are different between MPES and Orekit.
I read MPES guide but couldn’t find any further explanation on their “heliocentric position”.
The problem arises right at the start of the propagation, on 2019-04-27.
With the Keplerian inputs you provided, Orekit gives the following position:

2019-04-27T00:00:00.000 - Sun/inertial [AU]: pos = -1.8456785136544058 -0.9074482625561042 0.11813717411260657

While MPC’s ephemeris service gives:

    Obj           JD_TT               X               Y               Z            
00825        2458600.5000000     -1.845713900    -0.879496676   -0.252544599

You should probably try to find out which frame is used in the MPC data.

One thing seems weird to me though. The inclination is 3.39960° and the aphelion is 2.393 AU. So I thought that the Z component of the position should never be bigger (in absolute value) than 2.393*sin(3.39960°) = 0.142 AU.
However MEPS gives Z = -0.25 AU on 2019-04-27 and Z = 0.88 AU on 2018-04-01; which seems way too big for the given inclination.
I’m afraid the Keplerian elements are not given in the same frame than the Cartesian elements provided by MPES.

I selected “heliocentric position/velocity vector” in MPES to initialize the Orekit orbit with Cartesian coordinates.
This gives, on 2019-04-27:

VECTORS: Heliocentric vectors/AU
    Obj           JD_TT               X               Y               Z               X'              Y'              Z'
00825        2458600.5000000     -1.845713900    -0.879496676    -0.252544599     0.005565512    -0.010111635    -0.004594933

Beware the velocities are in AU/day.
Also it seems the date in MPC are given in TT timescale and not in UTC.

If you initialize orekit with these you will get, on 2018-04-01:

2018-04-01T00:00:00.000 - Sun/inertial [AU]: pos = 0.23454281838359037 2.1356736004084103 0.8815181137312449

Which is much closer than what you were expecting in your first message.

To sum up, for now you need to use the Cartesian coordinates to initialize Orekit.
If you want to use the Keplerian elements you should find out the frame they are given in in MPC.

Cheers,
Maxime

3

Hello MaximeJ,

I am very grateful for your investigation!

I am planning to use MPCORB database (https://minorplanetcenter.net/data) where JSON files contain keplerian elements. I will contact MPC and ask about the Frame.

Very clever idea with position/velocity initialization, now at least I know that this is not my code’s issue.

Thanks again, I will try to contribute to this community as much as I can.

Best Regards,
Piotr

4

After some hours of experimentation, I found the solution.

If anyone wonders, MPC is using Sun-center ecliptic coordinate system.

This is a code to create such frame:
EclipticProvider provider = new EclipticProvider(IERSConventions.IERS_2010);
Frame eclipticFrame = new Frame(sunBody.getInertiallyOrientedFrame(), provider, “eclipticFrame”,
true);

When this frame is used, for 2018-04-01 results are as follows:
x y z
0.23447562474898115 2.1357105548615065 0.8814336294268899
velX velY velZ
-0.011012475973598617 1.6335475377973047E-4 7.641261980640823E-4

1 Like
5

Good catch @AstroUser
Thank you for the hint.
I’m sure other users will be happy to find this out at some point.

6

Actually I have to make another update.
As my above example works for propagation around Sun, transforms between other planets give bad results.
EclipticProvider might not be made for transforming frames not centered in Earth.

I took another approach - created my own provider :

public class TransformSun implements TransformProvider {
	public Transform getTransform(AbsoluteDate date) throws OrekitException {
		Vector3D vector = FramesFactory.getEcliptic(IERSConventions.IERS_2010)
				.getTransformTo(CelestialBodyFactory.getSun().getInertiallyOrientedFrame(), date).getTranslation();
		return new Transform(date, vector);
	}
	public <T extends RealFieldElement<T>> FieldTransform<T> getTransform(FieldAbsoluteDate<T> date) {
		return null;
	}
}

The frame creation looks now :

TransformSun providerSun = new TransformSun();
Frame ecliptic = new Frame(FramesFactory.getEcliptic(IERSConventions.IERS_2010), providerSun, “ecliptic”, true);

Basically, I am taking pre-defined ecliptic frame centered in Earth and using a linear translation to offset it into center of Sun.
That works not only for Sun propagation, but then transformation of coordinates form sun-centric to earth-centric works just fine.

Regards,

7

There is another update:
I have seen that “x” coordinate has quite big error (0.8%). So I went to look, how X axis is oriented in ecliptic frame. Orekit states only that it is created out of MOD, but it is not exactly said.

From MPCORB format specification document, hosted on MPC’s database website, it is said that keplerian elements are provided in J2000, and inclination is “inclination to ecliptic”.

We should have then: J2000 MOD frame, aligned with ecliptic, and centered at Sun.
Whole chain looks now like that:

EclipticProvider provider = new EclipticProvider(IERSConventions.IERS_2010);
Frame earthEcliptic = new Frame(FramesFactory.getEME2000(), provider, "eclipticEarth", true);
TransformSun providerSun = new TransformSun(earthEcliptic);
Frame ecliptic = new Frame(earthEcliptic, providerSun, "eclipticSun", true);

At first, I am aligning EME2000 frame with ecliptic. We want to have X axis aligned with J2000.
Then this earth-ecliptic frame is translated to the center of a Sun, giving final Sun-center J2000 eclitpic frame.

public class TransformSun implements TransformProvider {
private Frame origin;

public TransformSun(Frame frameOrigin) {
	this.origin = frameOrigin;
}
public Transform getTransform(AbsoluteDate date) throws OrekitException {
	Vector3D vector = origin.getTransformTo(CelestialBodyFactory.getSun().getInertiallyOrientedFrame(), date)
			.getTranslation();
	return new Transform(date, vector);
}
public <T extends RealFieldElement<T>> FieldTransform<T> getTransform(FieldAbsoluteDate<T> date) {
	return null;
}
}

The results:
Using Tanina asteroid, with Keplerian elements from MPC website, and propagating it to 2018-04-01 the relative error of GEOCENTRIC (EME2000) position is:
(x-xref)/xref = 2.1208674664370107E-5
(y-yref)/yref = 1.562721112305075E-5
(z-zref)/zref = 9.12274035126E-5

Reference position has been taken from NASA JPL Horizons, as they give more description to their output.
( https://ssd.jpl.nasa.gov/horizons.cgi#results )

The propagator is the NumericalPropagator taken from :https://www.orekit.org/site-orekit-9.2/tutorial/propagation.html
Included is NewtonianAttraction of Sun, 3rdBodyAttraction of Earth, Mars, Moon , Jupiter.

Regards,

1 Like
8

The frame was still not valid, please refer to topic:

9

Was this ever resolved or is there still a discrepancy?

11

Awesome to hear!