Thank you very much @MaximeJ !
All of your suggestions worked! I have also added some force models, but TBH I didn’t notce any significant changes in estimated orbits (I’m working with LEO satellites with altitudes <1000km so I feel like I did something wrong).
Right now I’m dealing with another problem that I don’t know if it is easily resolvable. I’m working mainly with pretty short (but dense) arcs (for example I have ~100 measurments for 20 sec arc) and I don’t know if this is enough to estimate orbit with enough accuracy.
In my data I’ve found one satellite that I have observations of in 2 consecutive nights. So I performed OD on measurments from first night and propagated it to the next night. Unfortunately results from my propagation are way off compared to measurments:
Red dots are measurments used for orbit determination (blue lines). Green dots are measurments from next night, purple dots are calculated positions for next night.
Code that im working on:
from org.orekit.frames import FramesFactory
from org.orekit.utils import IERSConventions
gcrf = FramesFactory.getGCRF()
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
from org.orekit.time import TimeScalesFactory
utc = TimeScalesFactory.getUTC()
from org.orekit.models.earth import ReferenceEllipsoid
wgs84Ellipsoid = ReferenceEllipsoid.getWgs84(itrf)
from org.orekit.bodies import CelestialBodyFactory
moon = CelestialBodyFactory.getMoon()
sun = CelestialBodyFactory.getSun()
Transforming data from itrf to gcrf (I feel like there is more efficient and pythonic way):
from org.orekit.estimation.measurements import Position
import math
i = math.ceil(len(dane.index) / 2)
points_for_iod = dane.iloc[[0, i, -1]]
from org.hipparchus.geometry.euclidean.threed import Vector3D
from orekit.pyhelpers import datetime_to_absolutedate
from org.orekit.utils import PVCoordinates
pos_1 = points_for_iod.iloc[0]
p_vector_1 = Vector3D(pos_1[['x_', 'y_', 'z_']].to_list())
v_vector_1 = Vector3D(pos_1[["vx_", "vy_", "vz_"]].to_list())
date_1 = datetime_to_absolutedate(points_for_iod.index[0])
tmp_pv1 = PVCoordinates(p_vector_1, v_vector_1)
pv1 = itrf.getTransformTo(gcrf, date_1).transformPVCoordinates(tmp_pv1)
pos_2 = points_for_iod.iloc[1]
p_vector_2 = Vector3D(pos_2[['x_', 'y_', 'z_']].to_list())
v_vector_2 = Vector3D(pos_2[["vx_", "vy_", "vz_"]].to_list())
date_2 = datetime_to_absolutedate(points_for_iod.index[1])
tmp_pv2 = PVCoordinates(p_vector_2, v_vector_2)
pv2 = itrf.getTransformTo(gcrf, date_2).transformPVCoordinates(tmp_pv2)
pos_3 = points_for_iod.iloc[2]
p_vector_3 = Vector3D(pos_3[['x_', 'y_', 'z_']].to_list())
v_vector_3 = Vector3D(pos_3[["vx_", "vy_", "vz_"]].to_list())
date_3 = datetime_to_absolutedate(points_for_iod.index[2])
tmp_pv3 = PVCoordinates(p_vector_3, v_vector_3)
pv3 = itrf.getTransformTo(gcrf, date_3).transformPVCoordinates(tmp_pv3)
Performing IOD:
from org.orekit.estimation.iod import IodGibbs
from org.orekit.utils import Constants as orekit_constants
iod_gibbs = IodGibbs(orekit_constants.EIGEN5C_EARTH_MU)
initialOrbit = iod_gibbs.estimate(gcrf,
pv1.getPosition(), date_1,
pv2.getPosition(), date_2,
pv3.getPosition(), date_3)
display(initialOrbit)
#<KeplerianOrbit: Keplerian parameters: {a: 5665948.126048018; e: 0.2763193273172972; i: 81.17959707499644; pa: -44.55326269968; raan: 119.82447080790052; v: 178.94313046625413;}>
Setting up estimator and propagator:
# Estimator parameters
estimator_position_scale = 1.0 # m
estimator_convergence_thres = 1e-2
estimator_max_iterations = 25
estimator_max_evaluations = 35
# Orbit propagator parameters
prop_min_step = 0.001 # s
prop_max_step = 300.0 # s
prop_position_error = 10.0 # m
from org.orekit.propagation.conversion import DormandPrince853IntegratorBuilder
integratorBuilder = DormandPrince853IntegratorBuilder(prop_min_step, prop_max_step, prop_position_error)
from org.orekit.propagation.conversion import NumericalPropagatorBuilder
from org.orekit.orbits import PositionAngle, OrbitType
initialCartOrbit = OrbitType.CARTESIAN.convertType(initialOrbit)
propagatorBuilder = NumericalPropagatorBuilder(initialCartOrbit,
integratorBuilder,
PositionAngle.TRUE,
estimator_position_scale)
from org.orekit.forces.gravity.potential import GravityFieldFactory
gravityProvider = GravityFieldFactory.getConstantNormalizedProvider(64, 64)
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
gravityAttractionModel = HolmesFeatherstoneAttractionModel(itrf, gravityProvider)
propagatorBuilder.addForceModel(gravityAttractionModel)
from org.orekit.forces.gravity import ThirdBodyAttraction
moon_3dbodyattraction = ThirdBodyAttraction(moon)
propagatorBuilder.addForceModel(moon_3dbodyattraction)
sun_3dbodyattraction = ThirdBodyAttraction(sun)
propagatorBuilder.addForceModel(sun_3dbodyattraction)
from org.orekit.forces.radiation import IsotropicRadiationSingleCoefficient
isotropicRadiationSingleCoeff = IsotropicRadiationSingleCoefficient(100.0, 1.0);
from org.orekit.forces.radiation import SolarRadiationPressure
solarRadiationPressure = SolarRadiationPressure(sun, wgs84Ellipsoid.getEquatorialRadius(),
isotropicRadiationSingleCoeff)
propagatorBuilder.addForceModel(solarRadiationPressure)
from org.orekit.models.earth.atmosphere.data import MarshallSolarActivityFutureEstimation
from orekit import JArray
msafe = MarshallSolarActivityFutureEstimation(
'(?:jan|feb|mar|apr|may|jun|jul|aug|sep|oct|nov|dec)\p{Digit}\p{Digit}\p{Digit}\p{Digit}f10\.(?:txt|TXT)',
MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)
DM.feed(msafe.getSupportedNames(), msafe)
from org.orekit.models.earth.atmosphere import NRLMSISE00
atmosphere = NRLMSISE00(msafe, sun, wgs84Ellipsoid)
from org.orekit.forces.drag import IsotropicDrag
isotropicDrag = IsotropicDrag(100.0, 1.0)
from org.orekit.forces.drag import DragForce
dragForce = DragForce(atmosphere, isotropicDrag)
propagatorBuilder.addForceModel(dragForce)
from org.hipparchus.linear import QRDecomposer
matrix_decomposer = QRDecomposer(1e-11)
from org.hipparchus.optim.nonlinear.vector.leastsquares import GaussNewtonOptimizer
optimizer = GaussNewtonOptimizer(matrix_decomposer, False)
from org.orekit.estimation.leastsquares import BatchLSEstimator
estimator = BatchLSEstimator(optimizer, propagatorBuilder)
estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
estimator.setMaxIterations(estimator_max_iterations)
estimator.setMaxEvaluations(estimator_max_evaluations)
Adding measurments to the estimator and estimating:
from org.orekit.estimation.measurements import Position, ObservableSatellite
observableSatellite = ObservableSatellite(0) # Propagator index = 0
for timestamp, pv_gcrf in dane.iterrows():
tmp_p = Vector3D(pv_gcrf[['x_', 'y_', 'z_']].to_list())
tmp_v = Vector3D(pv_gcrf[['vx_', 'vy_', 'vz_']].to_list())
tmp_pvc = PVCoordinates(tmp_p, tmp_v)
pvc = itrf.getTransformTo(gcrf, datetime_to_absolutedate(timestamp)).transformPVCoordinates(tmp_pvc)
orekit_position = Position(
datetime_to_absolutedate(timestamp),
pvc.getPosition(),
1000.0,
1.0, # Base weight
observableSatellite
)
estimator.addMeasurement(orekit_position)
estimatedPropagatorArray = estimator.estimate()
Getting estimated orbit:
from org.orekit.orbits import KeplerianOrbit
estimatedPropagator = estimatedPropagatorArray[0]
estimatedInitialState = estimatedPropagator.getInitialState()
tmp_orbit = OrbitType.KEPLERIAN.convertType(estimatedInitialState.getOrbit())
orbit = KeplerianOrbit.cast_(tmp_orbit)
display(orbit)
#<KeplerianOrbit: Keplerian parameters: {a: 7129764.816986412; e: 0.01604576521836157; i: 81.19281051652376; pa: -18.340653585441842; raan: 119.8104073176335; v: 152.7333235809643;}>
Setting up bounded propagator:
dt = 60.0
date_start = datetime_to_absolutedate(dane.index[0])
date_end = datetime_to_absolutedate(dane.index[-1]).shiftedBy(86400.0)
ephmerisGenerator = estimatedPropagator.getEphemerisGenerator()
estimatedPropagator.propagate(date_start, date_end)
boundedPropagator = ephmerisGenerator.getGeneratedEphemeris()
Propagating to the next night, saving results to “pos_df”:
from orekit.pyhelpers import absolutedate_to_datetime
pos_df = pd.DataFrame(columns=['x', 'y', 'z', 'vx', 'vy', 'vz'])
date_current = date_start
while date_current.compareTo(date_end) <= 0:
datetime_current = absolutedate_to_datetime(date_current)
spacecraftState = boundedPropagator.propagate(date_current)
pv = spacecraftState.getPVCoordinates(itrf)
pos = np.array(pv.getPosition().toArray())
pos_df.loc[datetime_current] = np.concatenate(
(pos,
np.array(pv.getVelocity().toArray()))
)
date_current = date_current.shiftedBy(dt)
Calculating next night position in roughy same times as my measurments (“dane2” is dataframe with measurments from the next night):
pos_df_next = pd.DataFrame(columns=['x', 'y', 'z', 'vx', 'vy', 'vz'])
date_start = datetime_to_absolutedate(dane2.index[0])
date_end = datetime_to_absolutedate(dane2.index[-1])
dt = 0.2
display(date_start, date_end)
date_current = date_start
while date_current.compareTo(date_end) <= 0:
datetime_current = absolutedate_to_datetime(date_current)
spacecraftState = boundedPropagator.propagate(date_current)
pv = spacecraftState.getPVCoordinates(itrf)
pos = np.array(pv.getPosition().toArray())
pos_df_next.loc[datetime_current] = np.concatenate(
(pos,
np.array(pv.getVelocity().toArray()))
)
date_current = date_current.shiftedBy(dt)
Plotting (plot uploaded above code):
import plotly.graph_objects as go
fig_data = data=[
go.Scatter3d(x=pos_df['x'], y=pos_df['y'], z=pos_df['z'],
mode='lines',
name='Batch least squares solution'),
go.Scatter3d(x=dane['x_'], y=dane['y_'], z=dane['z_'],
mode='markers',
name='Measurements'),
go.Scatter3d(x=dane2['x_'], y=dane2['y_'], z=dane2['z_'],
mode='markers',
name='Measurements (next day)'),
go.Scatter3d(x=pos_df_next['x'], y=pos_df_next['y'], z=pos_df_next['z'],
mode='markers',
name='Batch least squares solution (next day)'),
]
scene=dict(
aspectmode='data',
)
layout = dict(
scene=scene
)
fig = go.Figure(
data=fig_data,
layout=layout)
fig.show()
Is there something really wrong with this code or it’s just not possible to determine accurate enough orbit with such short arcs?