Unexpectedly, they do not produce the same results. If I give the epoch time as propto, I expect that the result will be the same as the initial state. With the first form, it is, but with the second, it is different by hundreds of meters (for one LEO tested), and the difference between the two persists around that magnitude for the propagation over one hour. Since the first form gives the correct answer for zero propagation time, I trust it, but I am wondering what the reason is for the discrepancy.
Those two are not exactly the same thing. The propagator can only propagate in one way: with its propagate method (i.e. the second of the 2 you wrote). The other one is not actually propagating anything with the propagator. The getEphemerisGenerator method returns you the object is responsible for the storing of the computation of the integrator (i.e. the integration steps), which you are then retrieving with getGeneratedEphemeris. You now have a BoundedPropagator object which allows you to get the state at any point within the bounds of the arc it spans (the propagate method in this case does exactly that, it interpolates the orbit so that you can get the state at any timestamp within the given window). This means that, in general, the state you get from the Ephemeris is an interpolated state (unless you happen to query the timestamps of the integration steps), while the one you obtain from the NumericalPropagator is an integrated one.
Now, to address your concerns. I find it very strange that the first one is performing better than the second, since again, it’s an interpolation. If you’re querying the initial epoch of the integration, I would expect the two approaches to return the same results, since no interpolation should happen, but in general, the integrated state should always* be closer to the true state.
A couple of side notes:
The line you provided for the EphemerisGenerator should not work. That code is:
creating the propagator
extracting the ephemeris generator
extracting the ephemeris
interpolating
Now, since no propagation has actually happened, no ephemeris should be extracted and no interpolation should be allowed. The correct way to use this would be to:
create the propagator
create the ephemeris generator
propagate
extract the ephemeris from the generator
interpolate
Note that the ephemeris generator only stores the data for the latest propagation.
I ran a quick test and for me the second approach is returning the expected result (i.e. the initial state)
*This is not necessarily/strictly true, both integration and interpolation are approximate solutions. However, the second is an approximation based on the first one, while the first one is an approximation of the actual physical model, therefore, it should be more accurate and closer to the truth. The interpolation is just math, while the integration considers math (integration scheme) and physics (dynamics). However, it’s theoretically possible that the error introduced by the interpolation compensates the error introduced by the integration, resulting in the interpolation being closer to the true path. However, for a properly setup integration scheme, I think this discussion becomes academic/philosophical. If your integrator is properly capturing the dynamics of the system that means that the interpolation will also be “close enough” to the real trajectory and you won’t really see big differences between interpolation and integration error
Thanks for the thorough explanation. For your side note 1, note that I posted an oversimplification of my code to get to the point; in fact I do a propagation before I extract the ephemeris, as you indicate. What you say makes sense, I’ll look at this later with this in mind to see if I can come up with a better understanding of the result.
In fact, in the numerical propagation case, the interpolation used in the ephemeris is really based on the full stored integrator states. We use the continuous output feature of the integrator. This feature is what allows getting intermediate states even during the initial propagation (when using a step handler that uses a step interpolator instead of using a fixed step handler). What is stored in the ephemeris are really the step interpolators. This implies that ephemeris and integration should always return the same results.
Do you recreate both the integrator, the force models and the propagator between your two attempts or do you reconfigure already built objects in the second run? You should recreate objects because from step to step some values are updated and kept, so when reusing objects, there is a risk you do not restart exactly with the same configuration. What value did you use for tolerance (if you used an adaptive step integrator like Dormand-Prince 8(5,3) ? Note that for a long time, we did use extremely crude values (like 10m on position) and this may still appear in some tests and examples, this was a very bad advice, by now we strongly recommend using much smaller tolerances (like 1mm or even 0.1mm).
Just an additional remark to what’s already been said.
If your NumericalPropagator has true for the resetAtEnd flag, then it will keep its last integration state as the new initial condition. This means that integrating back and forth in time will give you slightly different results due to integration noise.
I have been successful in reproducing my results in a standalone block of code
import orekit_jpype as orekit
orekit.initVM()
from orekit_jpype.pyhelpers import setup_orekit_data
setup_orekit_data()
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.orbits import CartesianOrbit, OrbitType
from org.orekit.frames import FramesFactory
from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.numerical import NumericalPropagator
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.orekit.utils import PVCoordinates, TimeStampedPVCoordinates, IERSConventions
from org.orekit.propagation import Propagator, BoundedPropagator
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from orekit_jpype.pyhelpers import datetime_to_absolutedate
epoch = AbsoluteDate(2025, 1, 1, 0, 0, 00.000, TimeScalesFactory.getUTC())
vecp = Vector3D([5740132.68349, 3314067.15, 0.])
vecv = Vector3D([-2750.82684, 4764.57184, 5501.65367])
initstate = TimeStampedPVCoordinates(epoch, vecp, vecv)
gp = GravityFieldFactory.getNormalizedProvider(0, 0)
initstate = CartesianOrbit(initstate, FramesFactory.getEME2000(), gp.getMu())
tolerances = NumericalPropagator.tolerances(1.0, initstate, initstate.getType())
integrator = DormandPrince853Integrator(0.001, 1000.0, tolerances[0], tolerances[1])
integrator.setInitialStepSize(60.0)
initialState = SpacecraftState(initstate, 100.0)
propagator = NumericalPropagator(integrator)
propagator.setOrbitType(OrbitType.CARTESIAN)
propagator.setInitialState(initialState)
propagator.addForceModel(HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True), gp))
ephgen = propagator.getEphemerisGenerator()
propagator.propagate(initstate.getDate(), initstate.getDate().shiftedBy(86400.0))
gge = ephgen.getGeneratedEphemeris()
gge.propagate(initstate.getDate()) # This gets the answer I expect, same as initstate
# 5740132.68349, 3314067.15, 0.0
propagator.propagate(initstate.getDate()) # This is several hundred meters off
# 5740386.498610164, 3313626.921144008, -508.15693273954093
I found that the discrepancy shrinks to around 1m if the initial propagation is one hour instead of one day. I’m not sure how to interpret these results in light of the previous comments; with the exact code now explicit, what is happening here?
So the fact that the two results are not the same comes from what I explained: you propagate first with the reset flag at true (by default) so when you call propagate on your original initial date, it does backward propagation there. To avoid this do: propagator.setResetAtEnd(false)
Moreover, the difference you’re getting is rather big I think because of high integration noise in the call to tolerances. Use 1e-3 or 1e-4 as the first argument.
Thanks. I first reduced the tolerance to 1e-3 and that reduced the discrepancy to around 1m. Then I set the .setResetAtEnd and that gave agreement exactly (in all digits).
I had the impression that propagation from the generated ephemeris (which I presumed was some kind of interpolator) was preferable to re-integrating. This should be faster, correct?
Yes, it is correct. As usual, we are trading performances against memory consumption because storing all integrator internal states can be huge when tolerance is small and hence steps are small.