See below. Note that in the old Python wrapper there’s something stopping me calling getPVCoordinates on the new TLEPropagator at the end, but that only affects the printing of the PV so doesn’t affect the rest of the functionality. I’ve also changed the satellite to a GPS satellite so atmospheric drag isn’t involved.
Edit: added a .propagate() call that I missed
The output should be:
Initial TLE Epoch:
2024-12-11T18:33:12.231072Z
Starting TLE:
1 24876U 97035A 24346.77305823 .00000035 00000+0 00000+0 0 9991
2 24876 55.7332 120.1869 0086157 54.9246 305.8513 2.00561920200673
Input values from TLE in the ECI Frame
Px = -1.315473354461232E7
Py = 2.292137467376963E7
Pz = 31062.321250132172
Vx = -1879.722804250604
Vy = -1115.6288201491861
Vz = 3221.9184921375195
Estimated values based upon the Numerical Propagator in the ECI Frame
Px = -1.3282248886262149E7
Py = 2.2847572940558236E7
Pz = 216.22092157859151
Vx = -1881.2306313275792
Vy = -1126.2502392826345
Vz = 3217.3494285220254
Epoch (this is the Epoch of the input TLE!):
2024-12-11T18:33:12.231072Z
New TLE:
1 24876U 97035A 24346.77305823 .00000035 00000-0 00000-0 0 9993
2 24876 55.6132 120.5525 0086153 54.8501 305.8441 2.00564010200675
PVT of the TLE at the final time point
{2024-12-13T18:33:12.231072, P(-1.4142619061821254E7, 2.225304832353743E7, 1583440.2617032689), V(-1736.7012930178475, -1363.8053685657765, 3209.0587988536736), A(0.3058772009740569, -0.4812899163947692, -0.03424671193421054)}
PVT of the numerical propagator at the final time point
{2024-12-13T18:33:12.231072, P(-1.4140121819620019E7, 2.2255010028368253E7, 1577001.1469659884), V(-1737.449970286972, -1363.0427786061248, 3208.9899113618135), A(0.3058534479627173, -0.4813808599033754, -0.034117217119593236)}
try:
import orekit_jpype as orekit
vm = orekit.initVM()
from orekit_jpype.pyhelpers import (
setup_orekit_curdir,
absolutedate_to_datetime,
)
except ModuleNotFoundError:
import orekit
vm = orekit.initVM()
from orekit.pyhelpers import setup_orekit_curdir, absolutedate_to_datetime
from pathlib import Path
setup_orekit_curdir(str(Path.cwd() / "orekit-data"))
from org.orekit.propagation.analytical.tle import TLE, TLEPropagator
try:
from org.orekit.orbits import (
PositionAngle,
)
except ImportError:
from org.orekit.orbits import PositionAngleType
from org.orekit.propagation.conversion import (
TLEPropagatorBuilder,
FiniteDifferencePropagatorConverter,
)
from org.orekit.propagation.analytical.tle.generation import (
LeastSquaresTleGenerationAlgorithm,
)
from org.hipparchus.linear import QRDecomposer
from org.hipparchus.optim.nonlinear.vector.leastsquares import (
GaussNewtonOptimizer,
)
from org.orekit.estimation.leastsquares import BatchLSEstimator
from org.orekit.estimation.measurements import PV, ObservableSatellite
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.frames import FramesFactory
from org.orekit.utils import IERSConventions
from org.orekit.orbits import CartesianOrbit
from org.orekit.propagation.conversion import DormandPrince853IntegratorBuilder
from org.orekit.propagation.conversion import NumericalPropagatorBuilder
from org.orekit.models.earth import ReferenceEllipsoid
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from org.orekit.propagation import SpacecraftState
import numpy as np
rng = np.random.default_rng(seed=42)
rnorm = rng.normal
# Using GPS BIIR-2 (PRN 13) for purpose of demonstration
line1 = "1 24876U 97035A 24346.77305823 .00000035 00000+0 00000+0 0 9991"
line2 = "2 24876 55.7332 120.1869 0086157 54.9246 305.8513 2.00561920200673"
tle = TLE(line1, line2)
print("Initial TLE Epoch:")
print(tle.getDate())
print("Starting TLE:")
print(tle.toString())
# Getting some input data (Processing of laser data/GPS data would go here) ---
propagator = TLEPropagator.selectExtrapolator(tle)
epoch = tle.getDate()
end_date = absolutedate_to_datetime(epoch.shiftedBy(86400.0))
step = 100.0
cur_date = epoch
# Generate a set of input PV data
pvcoordinates = []
dates = []
while absolutedate_to_datetime(cur_date) < end_date:
pvcoordinates.append(propagator.getPVCoordinates(cur_date))
dates.append(cur_date)
cur_date = cur_date.shiftedBy(step)
# Now we have some data, set up the numerical propagator builder --------------
gcrf = FramesFactory.getGCRF()
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
# Selecting frames to use for OD
eci = gcrf
ecef = itrf
tleInitialState = propagator.getInitialState()
tleEpoch = tleInitialState.getDate()
tleOrbit_TEME = tleInitialState.getOrbit()
tlePV_ECI = tleOrbit_TEME.getPVCoordinates(eci)
wgs84Ellipsoid = ReferenceEllipsoid.getWgs84(ecef)
tleOrbit_ECI = CartesianOrbit(tlePV_ECI, eci, wgs84Ellipsoid.getGM())
# Propagator position error for adaptive step sizing
prop_position_error = 1e-3
# Propagator step size limits
prop_min_step = 0.1
prop_max_step = 600.0
# Getting an integrator builder
integratorBuilder = DormandPrince853IntegratorBuilder(
prop_min_step, prop_max_step, prop_position_error
)
# Getting a propagator builder
measurement_error = 100.0
propagatorBuilder = NumericalPropagatorBuilder(
tleOrbit_ECI, integratorBuilder, PositionAngleType.MEAN, measurement_error
)
propagatorBuilder.setMass(419725.0) # ISS Mass
# Adding a gravity model (more perturbations should be used)
gravityProvider = GravityFieldFactory.getConstantNormalizedProvider(
64, 64, tleEpoch
)
gravityAttractionModel = HolmesFeatherstoneAttractionModel(
ecef, gravityProvider
)
propagatorBuilder.addForceModel(gravityAttractionModel)
# Set estimation to use all the parameters (PV)
print("Input values from TLE in the ECI Frame")
for element in propagatorBuilder.getOrbitalParametersDrivers().getDrivers():
print(element)
element.setSelected(True)
# Setting up the estimator ----------------------------------------------------
convergence_threshold = 1e-3
max_iterations = 100
max_evaluations = 10
matrixDecomposer = QRDecomposer(1e-11)
optimizer = GaussNewtonOptimizer(matrixDecomposer, False)
estimator = BatchLSEstimator(optimizer, propagatorBuilder)
estimator.setParametersConvergenceThreshold(convergence_threshold)
estimator.setMaxIterations(max_iterations)
estimator.setMaxEvaluations(max_evaluations)
# Adding measurements to the estimator ----------------------------------------
observableSatellite = ObservableSatellite(0) # Propagator index = 0
for date, pv in zip(dates, pvcoordinates):
pos_err = rnorm(scale=measurement_error, size=3)
vel_err = rnorm(scale=measurement_error / 10, size=3)
orekit_position = PV(
date,
# Add noise to the measurement
pv.getPosition().add(Vector3D(*[float(x) for x in pos_err])),
pv.getVelocity().add(Vector3D(*[float(x) for x in vel_err])),
measurement_error,
measurement_error / 10, # Assume velocity error 10% of position error
1.0, # Base weight
observableSatellite,
)
estimator.addMeasurement(orekit_position)
# Now performing the estimation -----------------------------------------------
propagator = estimator.estimate()[0]
print("Estimated values based upon the Numerical Propagator in the ECI Frame")
for element in estimator.getOrbitalParametersDrivers(True).getDrivers():
print(element)
for element in estimator.getPropagatorParametersDrivers(True).getDrivers():
print(element)
# So we've fitted to the first day's worth of data, now we want to go to
# one day later, then fit a TLE to that.
final_date = epoch.shiftedBy(86400.0 * 2)
propagator.propagate(final_date) # Added this line to put the propagator where we want to be
# Now we have a propagator, we want to fit a TLE to the estimated propagator
# Final end date is the final point at which we want the TLE fitted to
tle_builder = TLEPropagatorBuilder(
tle,
PositionAngleType.MEAN,
measurement_error,
LeastSquaresTleGenerationAlgorithm(),
)
fitter = FiniteDifferencePropagatorConverter(
tle_builder, convergence_threshold, max_iterations
)
# Fitting a TLE to the numerical propagator -----------------------------------
# The span is set to a small number to make this quick - to get a good fit
# a larger number should be used.
fitted = fitter.convert(propagator, convergence_threshold, 60, "BSTAR")
try:
new_tle = fitted.getTLE()
except AttributeError:
# Older version of Orekit
new_tle = TLEPropagator.cast_(fitted).getTLE()
print("Epoch (this is the Epoch of the input TLE!):")
print(new_tle.getDate())
print("New TLE:")
print(new_tle.toString())
new_tle_propagator = TLEPropagator.selectExtrapolator(new_tle)
print("PVT of the TLE at the final time point")
# In the old orekit wrapper, this line bugs out
pv_tle = new_tle_propagator.getPVCoordinates(
final_date, propagator.getFrame()
)
print(pv_tle)
print("PVT of the numerical propagator at the final time point")
pv = propagator.getPVCoordinates(
final_date, propagator.getFrame()
)
print(pv)