Re: [Orekit Developers] [Orekit Users] Planetary Ephemeris

[MAISONOBE Luc <luc.maisonobe@c-s.fr>]

Hi,

[note that I have added the dev list here, if we continue
this thread it should probably be done on the dev lists
rather than the users list]

Evan Ward <evan.ward@nrl.navy.mil> a écrit :

> Hi Hank,
>
> My guess is a reference frame mismatch since the errors increase the
> further you get from Earth. Have you tried getting the positions in the
> GCRF frame in Orekit? I think Orekit assumes that the JPL ECI frame is
> GCRF. I think the slight bias between GCRF and J2000 as implemented in
> Orekit would be about the same order of magnitude for the errors you're
> seeing.

I agree this may be a problem. Take a look for example at slide 12
in the JPL NAIF tutorial about frames here:

  
<http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/Tutorials/pdf/individual_docs/17_frames_and_coordinate_systems.pdf>

This slides reads:

------------------
 The realization of ICRF was made to coincide almost
   exactly with the J2000 frame.

 –  The difference is very small–a rotation of less than 0.1 arc second
 –  These two frames are considered the same in SPICE
      »  In reality, any SPICE data said to be referenced to the J2000
         frame are actually referenced to the ICRF frame
      »  For historical and backwards compatibility reasons, only the
         name “J2000” is recognized by SPICE software as a frame
         name–not “ICRF”

•  No transformation is required to convert SPICE state
   vectors orientation data from the J2000 frame to
   the ICRF.
------------------

In Orekit, however, we do take into account this bias. As up to now
we thought JPL published everything in "ICRF aligned" frames centered
at various points and it appears it may really be in "EME2000 eligned"
frames at various points, we may need to change JPLEphemerisLoader.

There is however another issue we identified only yesterday, while
working on the IAU pole models. Up to now, Orekit did consider that
the reference dates of these models may be in some specific time scale
(most of the time TDB), the durations from these reference dates was
always computed in an Eart-based scale (i.e. with a time flow rate
equal to that of TT, TAI or UTC). This was wrong. For solar system
related models, the flow rate of time should be something less Earth-centric,
it is really defined in TDB. On the long term, there is no drift as
TDB is a small harmonic variation with respect to TT (due to Earth
velocity I guess and relativity I guess).

This was fixed yesterday for IAU pole (see commit 2fa02be). It should
also be fixed in JPL ephemeris. I have set up a first patch (not complete
because I also have to update unit tests, 13 of which now failing due
to the change). This path is attached. Could you try it and tell me if
this improves things? It will probably not fix everything as the
EME2000/ICRF bias issue is till there.

We are currently working on some new tests with reference data computed
using NAIF libraries, so I think we will eventually be able to fix
this completely.

best regards,
Luc

> Best Regards,
> Evan
>
> On 10/18/2014 03:26 PM, Hank Grabowski wrote:
> > I'm querying for planetary ephemeris and to make sure I am getting
> > expected values I'm comparing my results against the JPL HORIZONS
> > website.   Unfortunately I'm not seeing the accuracies that I was
> > expecting.   I'm querying for January 1st at midnight for the year
> > 2000 and 2015.  I'm doing this for the Sun, Moon and Jupiter.  The Sun
> > position is around 5100 meters off, the moon position is about 40
> > meters off and Jupiter's position is between 137 km and 179 km off.
> >
> > I downloaded and have Orekit using the 421 ephemerides at this point.
> > These get me closest.  I've tried 406, 423 and 430.  The 430
> > ephemerides didn't even load (I believe they changed the format
> > slightly) and 406 was worse by at least 20%.  My settings on the
> > website are:
> >
> > Ephemeris Type: Vectors
> > Target Body: Jupiter [599], Sun [10] or Moon [301]
> > Coordinate Origin: Geocentric [500] (which is Earth[399] BODY CENTER)
> > Reference plane: Earth mean equator and equinox of reference
> > Reference System: ICRF/J2000
> > Type: Geometric states (no aberration; instantaneous dynamic states)
> >
> > Times entered for vector format is always in Coordinate Time, which
> > matches Orekit's TDB timescale.  I'm generating my "actual" values
> > using code like this:
> >
> > PVCoordinatesProvider pvSun = CelestialBodyFactory.getSun();
> > Frame j2000 = FramesFactory.getEME2000();
> >
> > AbsoluteDate epoch = new AbsoluteDate(2000, 1, 1, 0, 0, 0,
> >                                               TimeScalesFactory.getTDB());
> >
> > Vector3D actualPosition = pvSun.getPVCoordinates(epoch, j2000)
> >             .getPosition();
> >
> >
> > It's not clear to me what I am doing wrong.  I've tried messing around
> > with the coordinate frame, timescale and changing the type generated
> > by HORIZONS to include light time delay and stellar aberration.  The
> > closest I've gotten to matching HORIZONS is with the above
> > configuration of the HORIZONS site and this code.
> >
> > If anything jumps out at someone, let me know.
> >
> > Thanks!
> >
> > Hank


diff --git a/src/main/java/org/orekit/bodies/JPLEphemeridesLoader.java b/src/main/java/org/orekit/bodies/JPLEphemeridesLoader.java
index 7859358..5174883 100644
--- a/src/main/java/org/orekit/bodies/JPLEphemeridesLoader.java
+++ b/src/main/java/org/orekit/bodies/JPLEphemeridesLoader.java
@@ -995,7 +995,7 @@ public class JPLEphemeridesLoader implements CelestialBodyLoader {
                 }
 
                 // build the position-velocity model for current chunk
-                entries.add(new PosVelChebyshev(chunkStart, duration, xCoeffs, yCoeffs, zCoeffs));
+                entries.add(new PosVelChebyshev(chunkStart, timeScale, duration, xCoeffs, yCoeffs, zCoeffs));
 
             }
 
diff --git a/src/main/java/org/orekit/bodies/PosVelChebyshev.java b/src/main/java/org/orekit/bodies/PosVelChebyshev.java
index 8565e31..12a79cc 100644
--- a/src/main/java/org/orekit/bodies/PosVelChebyshev.java
+++ b/src/main/java/org/orekit/bodies/PosVelChebyshev.java
@@ -21,6 +21,7 @@ import java.io.Serializable;
 import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
 import org.apache.commons.math3.util.FastMath;
 import org.orekit.time.AbsoluteDate;
+import org.orekit.time.TimeScale;
 import org.orekit.time.TimeStamped;
 import org.orekit.utils.PVCoordinates;
 
@@ -34,7 +35,10 @@ import org.orekit.utils.PVCoordinates;
 class PosVelChebyshev implements TimeStamped, Serializable {
 
     /** Serializable UID. */
-    private static final long serialVersionUID = -2220448511466595393L;
+    private static final long serialVersionUID = 20151023L;
+
+    /** Time scale in which the ephemeris is defined. */
+    private final TimeScale timeScale;
 
     /** Start of the validity range of the instance. */
     private final AbsoluteDate start;
@@ -53,6 +57,7 @@ class PosVelChebyshev implements TimeStamped, Serializable {
 
     /** Simple constructor.
      * @param start start of the validity range of the instance
+     * @param timeScale time scale in which the ephemeris is defined
      * @param duration duration of the validity range of the instance
      * @param xCoeffs Chebyshev polynomials coefficients for the X component
      * (a reference to the array will be stored in the instance)
@@ -61,14 +66,14 @@ class PosVelChebyshev implements TimeStamped, Serializable {
      * @param zCoeffs Chebyshev polynomials coefficients for the Z component
      * (a reference to the array will be stored in the instance)
      */
-    PosVelChebyshev(final AbsoluteDate start, final double duration,
-                           final double[] xCoeffs, final double[] yCoeffs,
-                           final double[] zCoeffs) {
-        this.start    = start;
-        this.duration = duration;
-        this.xCoeffs  = xCoeffs;
-        this.yCoeffs  = yCoeffs;
-        this.zCoeffs  = zCoeffs;
+    PosVelChebyshev(final AbsoluteDate start, final TimeScale timeScale, final double duration,
+                    final double[] xCoeffs, final double[] yCoeffs, final double[] zCoeffs) {
+        this.start     = start;
+        this.timeScale = timeScale;
+        this.duration  = duration;
+        this.xCoeffs   = xCoeffs;
+        this.yCoeffs   = yCoeffs;
+        this.zCoeffs   = zCoeffs;
     }
 
     /** {@inheritDoc} */
@@ -91,7 +96,7 @@ class PosVelChebyshev implements TimeStamped, Serializable {
      * @return true if the instance is the successor of the predecessor model
      */
     public boolean isSuccessorOf(final PosVelChebyshev predecessor) {
-        final double gap = start.durationFrom(predecessor.start) - predecessor.duration;
+        final double gap = start.offsetFrom(predecessor.start, timeScale) - predecessor.duration;
         return FastMath.abs(gap) < 0.001;
     }
 
@@ -100,7 +105,7 @@ class PosVelChebyshev implements TimeStamped, Serializable {
      * @return true if date is in validity range
      */
     public boolean inRange(final AbsoluteDate date) {
-        final double dt = date.durationFrom(start);
+        final double dt = date.offsetFrom(start, timeScale);
         return (dt >= -0.001) && (dt <= duration + 0.001);
     }
 
@@ -111,7 +116,7 @@ class PosVelChebyshev implements TimeStamped, Serializable {
     public PVCoordinates getPositionVelocityAcceleration(final AbsoluteDate date) {
 
         // normalize date
-        final double t = (2 * date.durationFrom(start) - duration) / duration;
+        final double t = (2 * date.offsetFrom(start, timeScale) - duration) / duration;
         final double twoT = 2 * t;
 
         // initialize Chebyshev polynomials recursion