Hi!
Once again, I’m coming back with problems and questions regarding orbit determination. What I aim to do is observe a satellite with a wide-FOV survey telescope. This gives me observation arcs with lengths of about 60–120 seconds. I perform IOD on these observations. Some time later (it can be a few hours if the next pass happens to be in the FOV, or 24+ hours later), I observe the same satellite again.
I would like to use these new observations in combination with the previous ones to obtain an even better orbit estimate. If possible, I’d like to avoid storing the observation data. After IOD, I’d prefer to store only the orbit information (for example, the epoch and state vector).
The problem I have is that, with angle-only observations, I believe I get a very poor estimate of the range (and therefore eccentricity), which causes the batch LS estimator not to converge when I use the initial orbit together with the new data. Is there an iterative way to improve an orbit solution every time new measurements arrive?
Important code snippets of what I have right now:
ARCSEC_TO_RAD = math.pi / (180.0 * 3600.0)
MU = Constants.WGS84_EARTH_MU
RE = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
utc = TimeScalesFactory.getUTC()
eme2000 = FramesFactory.getEME2000()
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(
Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
itrf
)
station = TopocentricFrame(
earth,
GeodeticPoint(math.radians(lat), math.radians(lon), el),
"station"
)
gs = GroundStation(station)
sat = ObservableSatellite(0)
...
measurements = []
for dt, ra_deg, dec_deg in obs:
meas = AngularRaDec(
gs,
eme2000,
dt,
[math.radians(ra_deg), math.radians(dec_deg)],
[rms, rms],
[1.0, 1.0],
sat
)
measurements.append(meas)
iod = IodGauss(Constants.WGS84_EARTH_MU)
orbit_iod = iod.estimate(
eme2000,
measurements[0],
measurements[len(measurements) // 2],
measurements[-1],
)
integrator_builder = DormandPrince853IntegratorBuilder(
0.001,
300.0,
10.0,
)
builder = NumericalPropagatorBuilder(
orbit_iod,
integrator_builder,
PositionAngleType.MEAN,
1000.0,
)
gravity_provider = GravityFieldFactory.getNormalizedProvider(20, 20)
gravity = HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravity_provider)
builder.addForceModel(gravity)
sun = CelestialBodyFactory.getSun()
atmosphere = HarrisPriester(sun, earth)
cross_section_m2 = 5
cd = 2.2
drag = DragForce(atmosphere, IsotropicDrag(cross_section_m2, cd))
builder.addForceModel(drag)
cr = 1.2
srp_spacecraft = IsotropicRadiationSingleCoefficient(cross_section_m2, cr)
srp = SolarRadiationPressure(sun, earth, srp_spacecraft)
builder.addForceModel(srp)
optimizer = LevenbergMarquardtOptimizer()
estimator = BatchLSEstimator(optimizer, builder)
estimator.setParametersConvergenceThreshold(1.0e-3)
estimator.setMaxIterations(50)
estimator.setMaxEvaluations(200)
for m in measurements:
estimator.addMeasurement(m)
estimated_propagators = estimator.estimate()
estimated_propagator = estimated_propagators[0]
orbit_refined = estimated_propagator.getInitialState().getOrbit()
That’s how I perform IOD, later I’m trying (where new_obs are ~24hours after the IOD):
integrator_builder = DormandPrince853IntegratorBuilder(
0.001,
300.0,
1.0,
)
builder = NumericalPropagatorBuilder(
orbit,
integrator_builder,
PositionAngleType.MEAN,
1000.0,
)
builder.addForceModel(gravity)
builder.addForceModel(drag)
builder.addForceModel(srp)
...
sat = ObservableSatellite(0)
new_measurements = []
for dt, ra, dec in new_obs:
meas = AngularRaDec(
gs,
eme2000,
dt,
[math.radians(ra), math.radians(dec)],
[new_rms, new_rms],
[1.0, 1.0],
sat
)
new_measurements.append(meas)
optimizer = LevenbergMarquardtOptimizer()
estimator = BatchLSEstimator(
optimizer,
builder
)
estimator.setParametersConvergenceThreshold(1.0e-3)
estimator.setMaxIterations(25)
estimator.setMaxEvaluations(35)
for m in new_measurements:
estimator.addMeasurement(m)
estimated_propagators = estimator.estimate()
estimated_propagator = estimated_propagators[0]
updated_orbit = estimated_propagator.getInitialState().getOrbit()
Of course, this fails to converge. Changing the number of iterations or evaluations doesn’t seem to make much difference.
Is there a way to perform orbit determination for LEO satellites using such short arcs that are separated by 24+ hours?
I’ve attached the data I use for testing. The RMS of this data should be approximately 5 arcseconds in both RA and Dec.
data1_2026-03-13.csv (42.3 KB)
data2_2026-03-14.csv (32.8 KB)