Generating a JPL SPICE kernel using orekit ephemeris

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
  }

I’ve reproduced the ephemeris using JPL/Horizons and check with my ephemeris generation.
The are different (see below).

The interesting thing is the definition of “Ecliptic at the standard reference epoch
at the end of the file. I suppose that this is the frame that should be used in orekit
to reproduce the ephemeris.


JPL/HORIZONS 470599 (2008 OG19) 2022-Feb-06 10:14:49
Rec #: 470599 (+COV) Soln.date: 2022-Feb-05_14:51:14 # obs: 412 (2008-2021)

IAU76/J2000 helio. ecliptic osc. elements (au, days, deg., period=Julian yrs):

EPOCH= 2457813.5 ! 2017-Mar-01.00 (TDB) Residual RMS= .096413
EC= .4208246209085354 QR= 38.57289818195139 TP= 2456526.0765123828
OM= 164.0323635412571 W= 140.4370264774148 IN= 13.16355471822164
A= 66.59968564696165 MA= 2.334625544334328 ADIST= 94.6264731119719
PER= 543.52119 N= .001813409 ANGMOM= .127348201
DAN= 81.12388 DDN= 41.38049 L= 305.2165494
B= 8.3400651 MOID= 37.56719971 TP= 2013-Aug-21.5765123828

Asteroid physical parameters (km, seconds, rotational period in hours):
GM= n.a. RAD= n.a. ROTPER= 8.727
H= 4.94 G= .150 B-V= n.a.
ALBEDO= n.a. STYP= n.a.

ASTEROID comments:
1: soln ref.= JPL#16, OCC=2
2: source=ORB



Ephemeris / WWW_USER Sun Feb 6 10:14:49 2022 Pasadena, USA / Horizons


Target body name: 470599 (2008 OG19) {source: JPL#16}
Center body name: Earth (399) {source: DE441}
Center-site name: Sierra Nevada Observatory


Start time : A.D. 2014-Jan-01 00:00:00.0000 TDB
Stop time : A.D. 2014-Dec-31 18:00:00.0000 TDB
Step-size : 360 minutes


Center geodetic : 356.615300,37.0641713,2.9288151 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 356.615300,5098.04112,3824.8436 {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}
Small perturbers: Yes {source: SB441-N16}
Output units : KM-S
Output type : GEOMETRIC cartesian states
Output format : 2 (position and velocity)
EOP file : eop.220203.p220429
EOP coverage : DATA-BASED 1962-JAN-20 TO 2022-FEB-03. PREDICTS-> 2022-APR-28
Reference frame : Ecliptic of J2000.0


Initial IAU76/J2000 heliocentric ecliptic osculating elements (au, days, deg.):
EPOCH= 2457813.5 ! 2017-Mar-01.00 (TDB) Residual RMS= .096413
EC= .4208246209085354 QR= 38.57289818195139 TP= 2456526.0765123828
OM= 164.0323635412571 W= 140.4370264774148 IN= 13.16355471822164
Equivalent ICRF heliocentric cartesian coordinates (au, au/d):
X= 2.539017738013506E+01 Y=-2.827653221397832E+01 Z=-6.999439963213364E+00
VX= 2.546312651901348E-03 VY= 2.081893209906159E-03 VZ= 2.331976176974970E-04
Asteroid physical parameters (km, seconds, rotational period in hours):
GM= n.a. RAD= n.a. ROTPER= 8.727
H= 4.94 G= .150 B-V= n.a.
ALBEDO= n.a. STYP= n.a.


        JDTDB, Calendar Date (TDB),                      X,                      Y,                      Z,                     VX,                     VY,                     VZ,

JPL/HORIZONS. 2014-Jan-01 00:00:00.0000, 3.372977465257917E+09, -4.773393508083921E+09, 8.254475881520498E+08, 3.481616918924401E+01, 8.639364048871855E+00, -1.027099467026608E+00,
MY_EPHEMERID 2014 Jan 01 00:00:00.000 4.447192243699149E9 -3.7266698562568884E9 5.522928898642708E8 3.9748941075621156 3.9165281154133824 -1.1306718725573761

TIME

Barycentric Dynamical Time (“TDB” or T_eph) output was requested. This
continuous relativistic coordinate time is equivalent to the relativistic
proper time of a clock at rest in a reference frame comoving with the
solar system barycenter but outside the system’s gravity well. It is the
independent variable in the solar system relativistic equations of motion.

TDB runs at a uniform rate of one SI second per second and is independent
of irregularities in Earth’s rotation.

Calendar dates prior to 1582-Oct-15 are in the Julian calendar system.
Later calendar dates are in the Gregorian system.

REFERENCE FRAME AND COORDINATES

** Ecliptic at the standard reference epoch**

Reference epoch: J2000.0
X-Y plane: adopted Earth orbital plane at the reference epoch
           Note: IAU76 obliquity of 84381.448 arcseconds wrt ICRF X-Y plane
X-axis   : ICRF
Z-axis   : perpendicular to the X-Y plane in the directional (+ or -) sense
           of Earth's north pole at the reference epoch.

Symbol meaning:

JDTDB    Julian Day Number, Barycentric Dynamical Time
  X      X-component of position vector (km)
  Y      Y-component of position vector (km)
  Z      Z-component of position vector (km)
  VX     X-component of velocity vector (km/sec)                           
  VY     Y-component of velocity vector (km/sec)                           
  VZ     Z-component of velocity vector (km/sec)                           

ABERRATIONS AND CORRECTIONS

Geometric state vectors have NO corrections or aberrations applied.

Hi @rmorales
It is difficult to analyze the problem. There are most probably frames errors, but I don’t think they are only a missing rotation. Looking at your first post, and computing the angular separation between position and velocity given by MPC and doing the same thing from your code leads to different angles (89.95° for MPC, 85.75° for your code). The modules of the vectors are also different. The dates are consistent and correspond to the initial state, so this cannot be explained by a difference in propagation force models. I am puzzled.

One thing that seems strange to me is that if I compute the cross product of the position and velocity from MPC results to get the orbital momentum and then look at its Z coordinate, it does not correspond to a 13.11° inclination at all (but rather to 11.34°), so I am not sure about the relationship between the reference Cartesian coordinates and the orbital elements for initializing the orbit, even prior to converting any frame.

Hi Luc,

You are right, I am mixing different coordinate systems and frames.
I need further study on how to set them up in orekit, MPC and JPL SPICE mkspk.

Thanks for your time in reviewing my post.
Best regards.

Hello everybody,

After some studies, I have corrected some bugs and simplified the code to generate an asteroid ephemeris.
Now JPL and MPC produce almost the same results. I guess the differences are related to the orbital perturbations considered and the different time scale (JTDB in JPL and TT in MPC).
The configuration used in JPL and MPC are in the attached images.

My code using Orekit almost generates similar results in X and Y axis but not in Z axis. I have tried different time scales (see below): TDB, TT and UTC with almost the same results.
The problem with the Z axis is not new and was originally discussed at.

a) Asteroid propagation not matching MPC results
b) Ecliptic frame, J2000 sun-centric - #3 by evan.ward

In fact, I’m using the framework defined in the last post, but the Z axis still have a different value.
I guess now the problem is in the frame definition. This is the information I have about the frame definition in JPL, MPC and my code.

Any help/clue would be appreciated.

#-------------------------------
Frame definition

#-------------------------------
JPL
#-------------------------------
International Celestial Reference Frame (ICRF)
IAU76/J2000 heliocentric ecliptic  osc. elements (au, days, deg., period=Julian yrs)

#-------------------------------
MPC
#-------------------------------
Heliocentric vectors/AU

#-------------------------------
MyCode (extracted from the post b))
#-------------------------------
    Heliocentric Aries (J2000) ecliptic frame
    This system has its Z axis perpendicular to the plane of the Earth’s orbit
    around the Sun (positive North) and its X axis towards the First Point of
    Aries (the direction in space defined by the intersection between the Earth’s
    equatorial plane and the plane of its orbit around the Sun (the plane of the
    ecliptic). This system is (to first order) fixed with respect to the distant
    stars.

public Transform getTransform(final AbsoluteDate date) throws OrekitException {
      // translation to Sun
      Vector3D sunCoords = PARENT_FRAME
            .getTransformTo(CelestialBodyFactory.getSun().getInertiallyOrientedFrame(), date).getTranslation();
      Transform translation = new Transform(date, sunCoords);
      // Ecliptic frame is based on MOD, so it has to be fixed on J2000 epoch where
      // EME2000 = MOD
      Transform rotation = PARENT_FRAME.getTransformTo(FramesFactory.getEcliptic(IERSConventions.IERS_2010),
            AbsoluteDate.J2000_EPOCH);

      return new Transform(date, translation, rotation);
   }



#------------------------------------------------------------------------------
#JPL (JTDB)
#------------------------------------------------------------------------------

        JDTDB,            Calendar Date (TDB), 

2456658.500000000, A.D. 2014-Jan-01 00:00:00.0000,

X, Y, Z,


2.237136669896858E+01, -3.058225474419807E+01, -7.245011392344233E+00,

VX, VY, VZ,


2.679191885366865E-03, 1.910546575274903E-03, 1.927081955270723E-04,

#------------------------------------------------------------------------------
#MPC (TT)
#------------------------------------------------------------------------------
JD_TT
2456658.5000000

X Y Z
22.370745536 -30.582674516 -7.245049448

X’ Y’ Z’
0.002679209 0.001910519 0.000192702

#------------------------------------------------------------------------------
#My code (TDB)
#------------------------------------------------------------------------------
TDB date
2013-12-31T23:58:52.81607656493335

X Y Z
22.372979502879865 -30.946344505777397 5.528654433467755

X’ Y’ Z’
3.1004372671743857E-8 2.1211561117427533E-8 -6.790065020044835E-9
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
#My code (TT)
#------------------------------------------------------------------------------
TT date
2013-12-31T23:58:52.816

X Y Z
22.372979502862012 -30.946344505789487 5.528654433471626

X’ Y’ Z’
3.1004372671753214E-8 2.121156111741524E-8 -6.790065020042673E-9

#------------------------------------------------------------------------------
#My code (UTC)
#------------------------------------------------------------------------------
UTC date
2014-01-01T00:00:00.000

X Y Z
22.372979409281665 -30.946344548084404 5.528654447746767

X’ Y’ Z’
3.100437282349128E-8 2.1211560955444972E-8 -6.790065064655598E-9
#------------------------------------------------------------------------------Preformatted text

//---------------------------------------------------------------------------
  @throws[OrekitException]
  def getPropagator(orb: Orbit) = {
    val tolerances = NumericalPropagator.tolerances(1, orb, 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(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))
    propagator
  }
  //---------------------------------------------------------------------------
  //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)
      translation
    }
    override def getTransform[T <: CalculusFieldElement[T]](date: FieldAbsoluteDate[T]): FieldTransform[T] = { // method not implemented
      null
    }
  }

def OG19_Ephemeris(dbHub: DatabaseHub): 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) 2

val asteroid = dbHub.mpcDB.getOrbitalElement(470599).get
val au2m =  Constants.IAU_2012_ASTRONOMICAL_UNIT            // astronomical unit to meters
val a = asteroid.a * au2m                                   // (1.0113461830157723E13) semi-major axis [m]
val e = asteroid.e                                          // (0.42897099256515503) orbital eccentricity
val i    = FastMath.toRadians(asteroid.i)                   // (0.22889486739824355) inclination [rad]
val peri = FastMath.toRadians(asteroid.w)                   // (2.463201368056165) argument of perihelion [rad]
val node = FastMath.toRadians(asteroid.N)                   // (2.8649983406708177) longitude of the ascending node [rad]
val m    = FastMath.toRadians(asteroid.M)                   // (0.08991639943050779) mean anomaly [rad]
val ep   = Time.fromJulian(asteroid.epoch)                  //2022-01-21T00:00 epoch
//-----------------------------------------------------------------
val timeScale   = TimeScalesFactory.getUTC
val sunBody     = CelestialBodyFactory.getSun                      // heliocentric
val mu          = sunBody.getGM()                                  // central attraction coefficient (m³/s²)
val anomalyType = PositionAngle.MEAN
val epoch = new AbsoluteDate(ep.getYear, ep.getMonthValue, ep.getDayOfMonth, ep.getHour, ep.getMinute, ep.getSecond, timeScale)
//-----------------------------------------------------------------
val initialDate = new AbsoluteDate("2014-01-01T00:00:00", timeScale)
val endDate     = new AbsoluteDate("2014-12-31T00:00:00", timeScale)
val dateStep    = Time.SECONDS_IN_ONE_DAY
val ephemerisOutputFile = "output/og19_ephemeris_data.txt"
//-----------------------------------------------------------------
//sun-center ecliptic coordinate system (inertial frame)
val inertialFrame = new Frame(FramesFactory.getEME2000(), new LocalTransformProvider(), "Heliocentric Aries (J2000) ecliptic frame", true)
//-----------------------------------------------------------------
//orbit
val orbit = new KeplerianOrbit(a, e, i, peri, node, m, anomalyType, inertialFrame, epoch, mu);
//-----------------------------------------------------------------
//propagator
val ephemerisPropagator = getPropagator(orbit)
//-----------------------------------------------------------------
//file generation
val pw = new PrintWriter(new File(ephemerisOutputFile))
var currentDate = initialDate
while (currentDate.isBefore(endDate)) {
  val state = ephemerisPropagator.propagate(currentDate).getPVCoordinates
  val p = state.getPosition //m
  val v = state.getVelocity //m/s

  //write position in au, speed in au/s
  val line = currentDate.toString.replace("Z","") +
    s" ${p.getX / au2m} ${p.getY / au2m} ${p.getZ / au2m} " +
    s" ${v.getX / au2m} ${v.getY / au2m} ${v.getZ / au2m}"

  println(line)
  pw.write(line + "\n")
  currentDate = currentDate.shiftedBy(dateStep)
}
pw.close()

}

#------------------------------------------------------------------------------Preformatted text

Dear all,

Finally I’ve fix the problem with Z axis and the frame.
Now I’m using “sunBody.getInertiallyOrientedFrame()” instead the “LocalTransformProvider
and applying a rotation on X axis using the obliquity of the ecliptic at requested time.

I’m sure that the rotation can be done using the current “LocalTransformProvider”, I but I did not how to indicated it.

Now, the last problem is the speed that is completely wrong.

#-------------------------------------------------------------------------
My Code
2014-01-01T00:00 

X, Y, Z,
22.379877914416127 -30.58660629859005 -7.246486864727551  

VX, VY, VZ,
3.09832025474121E-8 2.1235393890194217E-8 -6.761285175308039E-9
#-------------------------------------------------------------------------

#-------------------------------------------------------------------------
JPL Horizons
2456658.500000000, A.D. 2014-Jan-01 00:00:00.0000,

X, Y, Z,
2.237136669896858E+01, -3.058225474419807E+01, -7.245011392344233E+00,

VX, VY, VZ,
2.679191885366865E-03, 1.910546575274903E-03, 1.927081955270723E-04,

#-------------------------------------------------------------------------

 //---------------------------------------------------------------------------
  @throws[OrekitException]
  def getPropagator(orb: Orbit) = {
    val tolerances = NumericalPropagator.tolerances(1, orb, 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(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))
    propagator
  }

  /---------------------------------------------------------------------------
  //inclination of Earth's equator (projected on celestial sphere) with respect to the ecliptic
  //https://en.wikipedia.org/wiki/Ecliptic
  def getObliquityOfTheEclipticAngle(timestamp: LocalDateTime) = OBLIQUITY_OF_THE_ECLIPTIC - 3.563E-7 * Time.getTimeScale(timestamp)

  //---------------------------------------------------------------------------
  private final val FRACTION_DAY_HOUR =   1d / HOURS_IN_A_DAY
  private final val FRACTION_DAY_MINUTE = 1d / MINUTES_IN_A_ONE_DAY
  private final val FRACTION_DAY_SECOND = 1d / SECONDS_IN_A_ONE_DAY
  //-------------------------------------------------------------------------
  def getFractionDay(h: Int, m: Int, s:Double) =
    (h * FRACTION_DAY_HOUR) + (m * FRACTION_DAY_MINUTE) + (s * FRACTION_DAY_SECOND)
  /-------------------------------------------------------------------------
  //https://www.stjarnhimlen.se/comp/ppcomp.html#3
  //http://www.stjarnhimlen.se/comp/tutorial.html
  //The time scale in these formulae are counted in days. Hours, minutes, seconds are expressed as fractions of a day.
  // Day 0.0 occurs at 2000 Jan 0.0 UT (or 1999 Dec 31, 0:00 UT)
  def getTimeScale(date: LocalDateTime) = {
    val y = date.getYear
    val m = date.getMonthValue
    val D = date.getDayOfMonth
    val d = 367 * y - 7 * ( y + (m+9)/12 ) / 4 - 3 * ( ( y + (m-9)/7 ) / 100 + 1 ) / 4 + 275*m/9 + D - 730515
    val fractionalDay = getFractionDay(date.getHour,date.getMinute,date.getSecond)
    d + fractionalDay
  }
//---------------------------------------------------------------------------
  def OG19_Ephemeris(dbHub: DatabaseHub): 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) 2
    val asteroid = dbHub.mpcDB.getOrbitalElement(470599).get
    val au2m =  Constants.IAU_2012_ASTRONOMICAL_UNIT            // astronomical unit to meters
    val a = asteroid.a * au2m                                   // (1.0113461830157723E13) semi-major axis [m]
    val e = asteroid.e                                          // (0.42897099256515503) orbital eccentricity
    val i    = FastMath.toRadians(asteroid.i)                   // (0.22889486739824355) inclination [rad]
    val peri = FastMath.toRadians(asteroid.w)                   // (2.463201368056165) argument of perihelion [rad]
    val node = FastMath.toRadians(asteroid.N)                   // (2.8649983406708177) longitude of the ascending node [rad]
    val m    = FastMath.toRadians(asteroid.M)                   // (0.08991639943050779) mean anomaly [rad]
    val ep   = Time.fromJulian(asteroid.epoch)                  //2022-01-21T00:00 epoch
    //-----------------------------------------------------------------
    val timeScale   = TimeScalesFactory.getUTC
    val sunBody     = CelestialBodyFactory.getSun                      // heliocentric
    val mu          = sunBody.getGM()                                  // central attraction coefficient (m³/s²)
    val anomalyType = PositionAngle.MEAN
    val epoch = new AbsoluteDate(ep.getYear, ep.getMonthValue, ep.getDayOfMonth, ep.getHour, ep.getMinute, ep.getSecond, timeScale)
    //-----------------------------------------------------------------
    val initialDate = new AbsoluteDate("2014-01-01T00:00:00", timeScale)
    val endDate     = new AbsoluteDate("2014-12-31T00:00:00", timeScale)
    val dateStep    = Time.SECONDS_IN_ONE_DAY
    val ephemerisOutputFile = "output/og19_ephemeris_data.txt"
    //-----------------------------------------------------------------
    //sun-center ecliptic coordinate system (inertial frame)
    val inertialFrame =  sunBody.getInertiallyOrientedFrame()
    //-----------------------------------------------------------------
    //orbit
    val orbit = new KeplerianOrbit(a, e, i, peri, node, m, anomalyType, inertialFrame, epoch, mu);
    //-----------------------------------------------------------------
    //propagator
    val ephemerisPropagator = getPropagator(orbit)
    //---------------------------------------------------------------------------
    //(x,y,z) all in au
    def getHeliocentricEquatorialCoordinates(timestamp: LocalDateTime, x: Double ,y: Double ,z: Double) = {
      //it is just a rotation on Z axis (X is fixed, Y and Z axis are rotated) of the obliquity of the ecliptic applied on heliocentric ecliptic coordinates
      val obEc = toRadians(getObliquityOfTheEclipticAngle(timestamp))
      val sinObEc = sin(obEc)
      val cosObEc = cos(obEc)

      //Celestial mechanics. Jeremy Tatum. Eq. 10.7.10
      val y2 = (y * cosObEc) + (-sinObEc * z)
      val z2 = (y * sinObEc) + (cosObEc * z)
      (x,y2,z2)
    }
    //-----------------------------------------------------------------
    //file generation
    val pw = new PrintWriter(new File(ephemerisOutputFile))
    var currentDate = initialDate
    while (currentDate.isBefore(endDate)) {
      val state = ephemerisPropagator.propagate(currentDate).getPVCoordinates
      val p = state.getPosition //m
      val v = state.getVelocity //m/s
      val timeStamp = LocalDateTime.parse(currentDate.toString.replace("Z",""))

      val (x,y,z) = getHeliocentricEquatorialCoordinates(
         timeStamp
        , p.getX / au2m
        , p.getY / au2m
        , p.getZ / au2m)

      val (vx,vy,vz) =(
          v.getX / au2m
        , v.getY / au2m
        , v.getZ / au2m)


      //write position in au, speed in au/s
      val line = timeStamp +
        s" $x $y $z " +
        s" $vx $vy $vz"


      println(line)
      pw.write(line + "\n")
      currentDate = currentDate.shiftedBy(dateStep)
    }
    pw.close()

    //now build the binary spk:
    //./mkspk -setup /home/rafa/proyecto/henosis/output/og19_setup.txt -input/home/rafa/proyecto/henosis/output/og19_ephemeris_data.txt -output /home/rafa/proyecto/m2/input/spice/kernels/spk/objects/2470599_og19.bsp
  }
  //---------------------------------------------------------------------------

It seems there is a unit error, your velocity is probably per seconds and JPL velocity is per day.
After correcting this issue, there is still a 13° error in velocity direction, despite position direction is good. Perhaps some wrong velocity composition?

Hi Luc,

You are right, it’s a problem with the units. Now, all is in KM and KM/S.
Also it is required a rotation on X axis also in speed. After the changes,
the result are close (0.01 deg). I’m assuming that the differences are related to the perturbers used, and maybe in the time scale: JPL in TBC and my code in UTC.

If some people is interested, the code can be downloaded from:

https://gitlab.com/rmorales_iaa/henosis/-/blob/main/src/main/scala/com/henosis/SandBox.scala

 val (x,y,z) = getHeliocentricEquatorialCoordinates(
         timeStamp
        , p.getX / 1000
        , p.getY / 1000
        , p.getZ / 1000)


  val (vx,vy,vz) = getHeliocentricEquatorialCoordinates(
        timeStamp
        , v.getX  / 1000
        , v.getY  / 1000
        , v.getZ  / 1000)

#--------------------------------------------------------------------------
#JPL (KM and KM/S)
 2014-Jan-01 00:00:00.0000,  
 X,                      Y,                      Z,
 3.346708822814586E+09, -4.575040190937004E+09, -1.083838277491939E+09,  

 VX,                     VY,                     VZ,
 4.638905107032425E+00,  3.308028929795170E+00,  3.336659226538110E-01,

#-------------------------------------------------------------------------- 
#my code (KM and KM/S)
2014-01-01T00:00 
X,                 Y,                    Z,
3.34798209635394E9 -4.575691168908265E9 -1.0840590044320982E9

X,               VY,               VZ,
4.63502111123374 3.316981125108725 0.3355317874377496
#--------------------------------------------------------------------------