Hi all!
I’m trying to generate a JPL SPICE kernel (SPK) using the orbit parameters (MPC or JPL Small database).
I use orekit to propagate the orbit during a time period obtaining a set of ephemeris that are saved to a file.
This file generated is finally passed to JPL “mkspk” tool to create the desired SPK file in its binary forn (BSP).
The problem is that the ephemeris that I obtain are different that I get when using MPC ephemereis service.
MPC
VECTORS: Heliocentric vectors/AU
Obj JD_TT X Y Z X’ Y’ Z’
l0599 2456658.5000000 22.370745536 -30.582674516 -7.245049448 0.002679209 0.001910519 0.000192702
MyCode X Y Z X’ Y’ Z’
2014 Jan 01 00:00:00.000 4.44719224369915E9 -3.72666985625689E9 5.52292889864271E8 3.9748941075621174 3.9165281154133837 -1.1306718725573766
Reading the posts, I suppose that the problem is in the heliocentric frame definition.
Here is my code (it is in Scala, but can easily translated to Java) to process the ephemeris of asteroid “470599 (2008 OG19)” along 2014.
Any comment will be appreciate.
Thanks in advance.
//---------------------------------------------------------------------------
//https://forum.orekit.org/t/ecliptic-frame-j2000-sun-centric/373/3
val PARENT_FRAME = FramesFactory.getEME2000()
private class LocalTransformProvider extends TransformProvider {
@throws[OrekitException]
override def getTransform(date: AbsoluteDate): Transform = { // translation to Sun
val sunCoords = PARENT_FRAME.getTransformTo(CelestialBodyFactory.getSun.getInertiallyOrientedFrame, date).getTranslation
val translation = new Transform(date, sunCoords)
// Ecliptic frame is based on MOD, so it has to be fixed on J2000 epoch where
// EME2000 = MOD
val rotation = PARENT_FRAME.getTransformTo(FramesFactory.getEcliptic(IERSConventions.IERS_2010), AbsoluteDate.J2000_EPOCH)
new Transform(date, translation, rotation)
}
override def getTransform[T <: CalculusFieldElement[T]](date: FieldAbsoluteDate[T]): FieldTransform[T] = { // method not implemented
null
}
}
//---------------------------------------------------------------------------
#-----------------------------------------------------------------------------
def oreKitEphemeris(): Unit = {
//=========================================================================
//https://www.orekit.org/site-orekit-tutorials-10.3/tutorials/propagation.html
val orekitData = new File("/home/rafa/proyecto/henosis/input/orekit-data-master/")
val manager = DataContext.getDefault.getDataProvidersManager
manager.addProvider(new DirectoryCrawler(orekitData))
//-----------------------------------------------------------------
//470599 (2008 OG19)
//https://ssd.jpl.nasa.gov/tools/sbdb_lookup.html#/?sstr=2470599
val a = 67.60470176087289 * Constants.IAU_2012_ASTRONOMICAL_UNIT // semi-major axis [m]
val e = 0.4289740851645675 // orbital eccentricity
val i = FastMath.toRadians(13.1147010952624) // inclination [rad]
val peri = FastMath.toRadians(141.1304206312195) // Argument of perihelion [rad]
val node = FastMath.toRadians(164.1523398804327) // Longitude of the ascending node [rad]
val m = FastMath.toRadians(5.152000573779065) // Mean anomaly [rad]
//-----------------------------------------------------------------
val timeScale = TimeScalesFactory.getUTC
val sunBody = CelestialBodyFactory.getSun // heliocentric
val mu = sunBody.getGM() // central attraction coefficient (m³/s²)
val anomalyType = PositionAngle.MEAN
//-----------------------------------------------------------------
val initialDate = new AbsoluteDate(2014, 1, 1, 0, 0, 00.000, timeScale) //jd 2456658.5000000
val integratorStepSize = 100
val integratorYearsToGenerate = 1
val ephemerisStepInHours = 6
val ephemerisOutputFile = "output/og19_ephemeris_data.txt"
//-----------------------------------------------------------------
//sun-center ecliptic coordinate system (inertial frame)
val inertialFrame = new Frame(PARENT_FRAME, new LocalTransformProvider(), "Heliocentric Aries (J2000) ecliptic frame", true);
//-----------------------------------------------------------------
//orbit
val orbit = new KeplerianOrbit(a, e, i, peri, node, m, anomalyType, inertialFrame, initialDate, mu);
//-----------------------------------------------------------------
//propagator
val tolerances = NumericalPropagator.tolerances(1, orbit, OrbitType.CARTESIAN)
val integrator = new DormandPrince853Integrator(0.001, 3600, tolerances(0), tolerances(1))
val propagator = new NumericalPropagator(integrator)
propagator.setOrbitType(OrbitType.CARTESIAN)
propagator.addForceModel(new NewtonianAttraction(orbit.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(orbit))
val generator = propagator.getEphemerisGenerator()
propagator.propagate(initialDate.shiftedBy(Time.daysToSeconds(Time.DAYS_PER_YEAR) * integratorYearsToGenerate))
val ephemerisPropagator = generator.getGeneratedEphemeris
//-----------------------------------------------------------------
println("Ephemeris defined from " + ephemerisPropagator.getMinDate() + " to " + ephemerisPropagator.getMaxDate())
val startTime = LocalDateTime.parse(ephemerisPropagator.getMinDate.toString().dropRight(1),DateTimeFormatter.ISO_LOCAL_DATE_TIME)
val stopTime = LocalDateTime.parse(ephemerisPropagator.getMaxDate.toString().dropRight(1),DateTimeFormatter.ISO_LOCAL_DATE_TIME)
var currentTime = startTime
//-----------------------------------------------------------------
def getAbsoluteDate(t: LocalDateTime) =
new AbsoluteDate(t.getYear
, t.getMonthValue, t.getDayOfMonth, t.getHour, t.getMinute, t.getSecond + Time.nanosToSeconds(t.getNano), timeScale)
//-----------------------------------------------------------------
//file generation
val ephemerisStepSeconds = Time.SECONDS_IN_ONE_HOUR * ephemerisStepInHours
val pw = new PrintWriter(new File(ephemerisOutputFile))
while (currentTime.isBefore(stopTime)) {
val state = ephemerisPropagator.propagate(getAbsoluteDate(currentTime)).getPVCoordinates()
val p = state.getPosition //m
val v = state.getVelocity //m/s
//write position in KM, speed in KM/s
val line = currentTime.format(GaiaToMongo.ephemerisTimeFormatter) +
s" ${p.getX / 1000} ${p.getY / 1000} ${p.getZ / 1000} " +
s" ${v.getX / 1000} ${v.getY / 1000} ${v.getZ / 1000}"
println(p.getX / Constants.IAU_2012_ASTRONOMICAL_UNIT)
println(line)
pw.write(line + "\n")
currentTime = currentTime.plusSeconds(ephemerisStepSeconds)
}
pw.close()
//now build the binary spk:
//./mkspk -setup /home/rafa/proyecto/henosis/output/og19_ephemeris_setup.txt -input /home/rafa/proyecto/henosis/output/og19_ephemeris_data.txt -output /home/rafa/proyecto/henosis/output/og19.bsp
}