Propagation Error: Keplerian propagator and BLS estimator

Hi, Orekit community!

I am using a Keplerian propagator and BLS estimator to get an ephemeris with gps data given.
And I keep getting this error.

    spacecraftState = estimatedPropagator.propagate(date_current)
orekit.JavaError: <super: <class 'JavaError'>, <JavaError object>>
    Java stacktrace:
org.orekit.errors.OrekitIllegalArgumentException: true anomaly 2,101.372 out of hyperbolic range (e = 36.013, -1.599 < v < 1.599)
	at org.orekit.orbits.KeplerianOrbit.<init>(
	at org.orekit.orbits.KeplerianOrbit.<init>(
	at org.orekit.orbits.OrbitType$4.mapArrayToOrbit(
	at org.orekit.orbits.OrbitType$4.mapArrayToOrbit(
	at org.orekit.propagation.numerical.NumericalPropagator$OsculatingMapper.mapArrayToState(
	at org.orekit.propagation.integration.StateMapper.mapArrayToState(
	at org.orekit.propagation.integration.AbstractIntegratedPropagator$ConvertedMainStateEquations.computeDerivatives(
	at org.hipparchus.ode.ExpandableODE.computeDerivatives(
	at org.hipparchus.ode.AbstractIntegrator.computeDerivatives(
	at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(

I am not sure what “org.orekit.errors.OrekitIllegalArgumentException: true anomaly 2,101.372 out of hyperbolic range (e = 36.013, -1.599 < v < 1.599)” this part of the error is saying.

If any of my code should be provided to read more context of the error, please let me know.

Thank you!

Hi @psh0078

Could you provide a runnable test to reproduce your issue?

Thank you,

Hi Bryan @bcazabonne ,

I figured it out. I do not understand how it works in detail, but I was fitting the data over like 5 days or something.
Then, I ran the script on a smaller dataset over a period of 24 hours. Unfortunately, I got stuck on the “KEPLERIAN is singluar with current orbit”, but I was able to solve the error by converting Keplerian parameters to cartesian.

Now, I got stuck on this error below, and I have no clue what this is saying.


org.orekit.errors.OrekitException: NaN appears during integration near time -197.642
	at org.orekit.errors.OrekitException.unwrap(
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(
	at org.orekit.propagation.PropagatorsParallelizer.propagate(
	at org.orekit.estimation.leastsquares.AbstractBatchLSModel.value(
	at org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresFactory$LocalLeastSquaresProblem.evaluate(
	at org.orekit.estimation.leastsquares.BatchLSEstimator$TappedLSProblem.evaluate(
	at org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.optimize(
	at org.orekit.estimation.leastsquares.BatchLSEstimator.estimate(
Caused by: org.hipparchus.exception.MathIllegalStateException: NaN appears during integration near time -197.642
	at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(
	at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(
	... 7 more

Following are the code snippets that might help:

i = math.floor(len(points_df.index) / 3)
points_for_iod = points_df.iloc[::i, :]

pos_1 = points_for_iod.iloc[0]
vector_1 = Vector3D(pos_1[['x', 'y', 'z']].to_list())
date_1 = datetime_to_absolutedate(points_for_iod.index[0])

pos_2 = points_for_iod.iloc[1]
vector_2 = Vector3D(pos_2[['x', 'y', 'z']].to_list())
date_2 = datetime_to_absolutedate(points_for_iod.index[1])

pos_3 = points_for_iod.iloc[2]
vector_3 = Vector3D(pos_3[['x', 'y', 'z']].to_list())
date_3 = datetime_to_absolutedate(points_for_iod.index[2])

from org.orekit.orbits import OrbitType

iod_gibbs = IodGibbs(orekit_constants.EIGEN5C_EARTH_MU)
orbit_first_guess = iod_gibbs.estimate(gcrf,
                                      vector_1, date_1,
                                      vector_2, date_2,
                                      vector_3, date_3)
initialCartOrbit = OrbitType.CARTESIAN.convertType(orbit_first_guess)
kepler_propagator_iod = KeplerianPropagator(initialCartOrbit)

integratorBuilder = DormandPrince853IntegratorBuilder(prop_min_step, prop_max_step, prop_position_error)
propagatorBuilder = NumericalPropagatorBuilder(initialCartOrbit, integratorBuilder, PositionAngle.TRUE, estimator_position_scale)

# Adding drag
cd = 0.01 # drag coefficient
area = 1.5 # cross sectional area of an object
cdEstimated = True

msafe = MarshallSolarActivityFutureEstimation(MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES, MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)
manager = DataContext.getDefault().getDataProvidersManager()
manager.feed(msafe.getSupportedNames(), msafe)
atmosphere = NRLMSISE00(msafe, CelestialBodyFactory.getSun(), earth)

drag = DragForce(atmosphere, IsotropicDrag(area, cd))

for driver in drag.getParametersDrivers():
    if driver.getName() == DragSensitive.DRAG_COEFFICIENT:
        if cdEstimated:


matrix_decomposer = QRDecomposer(1e-11)
optimizer = GaussNewtonOptimizer(matrix_decomposer, False)

estimator = BatchLSEstimator(optimizer, propagatorBuilder)

observableSatellite = ObservableSatellite(0) # Propagator index = 0

for timestamp, pv_gcrf in points_df.iterrows():
    date = datetime_to_absolutedate(timestamp)
    vector_pos = Vector3D(pv_gcrf[['x', 'y', 'z']].to_list())
    vector_vel = Vector3D(pv_gcrf[['vx', 'vy', 'vz']].to_list())
    orekit_pv = PV(
        1.0,  # Base weight

estimatedPropagatorArray = estimator.estimate()

points_df is the pandas data frame that contains all the GPS data. (time, lat, long, alt, positions, velocities)

Thank you for the code snippets. Please find below some questions/remarks that could help you:

  1. Why considering only the drag perturbation? I think you could also add the HolmesFeatherstoneAttractionModel to improve the consistency between your GPS data and the propagation model. Adding luni-solar point masses could also be useful (i.e. ThirdBodyAttraction class)

  2. In which frame are expressed your GPS data? Are you sure they are expressed in GCRF? I ask the question because usually GPS data are in ITRF.

  3. What are the units of vector_pos and vector_vel? Orekit uses S.I units (i.e. meters for position and meters per second for velocity)


Hi Bryan @bcazabonne,

  1. I added the gravitational field model to the propagator, and the overall performance improved significantly!

  2. I have been using GTOD for my GPS data. Given no metadata about my GPS data I have, I am not sure what frame the GPS data are in. I tried both time frames, but the difference was not noticeable.
    I know I should use the right one though. Is using ITRF more of a standard way?

gtod = FramesFactory.getGTOD(True)
GPSTime = GNSSDate(int(GNSSWeeks), float(GNSSSeconds)*1000, gps)
utcDate = GPSTime.getDate()
vector_gtod = Vector3D([xgtod, ygtod, zgtod])
                    velocity_gtod = Vector3D([vx, vy, vz])
                    pv = TimeStampedPVCoordinates(utcDate, vector_gtod, velocity_gtod)
                    pv_gcrf = gtod.getTransformTo(gcrf, utcDate).transformPVCoordinates(pv)

                    TLESpacecraftStates.add(SpacecraftState(CartesianOrbit(pv_gcrf, gcrf, utcDate, Constants.WGS84_EARTH_MU)))

Thanks for your help,