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.
Thanks!
org.orekit.errors.OrekitException: NaN appears during integration near time -197.642
at org.orekit.errors.OrekitException.unwrap(OrekitException.java:154)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:511)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:408)
at org.orekit.propagation.PropagatorsParallelizer.propagate(PropagatorsParallelizer.java:140)
at org.orekit.estimation.leastsquares.AbstractBatchLSModel.value(AbstractBatchLSModel.java:319)
at org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresFactory$LocalLeastSquaresProblem.evaluate(LeastSquaresFactory.java:440)
at org.orekit.estimation.leastsquares.BatchLSEstimator$TappedLSProblem.evaluate(BatchLSEstimator.java:615)
at org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.optimize(GaussNewtonOptimizer.java:163)
at org.orekit.estimation.leastsquares.BatchLSEstimator.estimate(BatchLSEstimator.java:435)
Caused by: org.hipparchus.exception.MathIllegalStateException: NaN appears during integration near time -197.642
at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:269)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:477)
... 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:
driver.setSelected(True)
propagatorBuilder.addForceModel(drag)
matrix_decomposer = QRDecomposer(1e-11)
optimizer = GaussNewtonOptimizer(matrix_decomposer, False)
estimator = BatchLSEstimator(optimizer, propagatorBuilder)
estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
estimator.setMaxIterations(estimator_max_iterations)
estimator.setMaxEvaluations(estimator_max_evaluations)
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(
date,
vector_pos,
vector_vel,
sigma_position,
sigma_velocity,
1.0, # Base weight
observableSatellite
)
estimator.addMeasurement(orekit_pv)
estimatedPropagatorArray = estimator.estimate()
points_df is the pandas data frame that contains all the GPS data. (time, lat, long, alt, positions, velocities)