I am cross posting this from the python-wrapper issue tracker based on sage advice from Maxime…
I have been shopping around for an efficient numerical or SGP4 propagator for use with a python Gym environment I am working on. The environment will keep a true state for each space vehicle (SV) orbiting the earth as well as an Unscented Kalman Filter (UKF) propagation from the initial state. At each time-step, the an observation of one or more of the SVs will be made, updating their UKF state. I think Orekit might be the solution I am looking for, but I seem to be falling short in my implementation so far. What I am trying to accomplish is to have a method I can pass a current state (PV or Kep) to along with a time step and get a single step integration to the new step with the updated coordinates (PV or Kep - same as used for current state). My experience so far is that the propagation step is plenty fast (7e-5s per step), but the process of initializing each propagator slows the system down to more like 1.8e-1s per initialization and step. Code included below if you want to see how I am implementing Orekit currently to test the computation expense of each sigma point in the UKF.
I have started looking at Shiva-Iyer’s orbdetpy (/github.com/ut-astria/orbdetpy) based on feedback from Maxime, but at first pass their implementation is different than I would do myself. I had planned to keep the variables in memory vice writing them out to files (json), but I am still looking deeper into it as I am sure they had sound reasons for choosing that direction.
1.) Is Orekit even the right option for my propagator? I have run similar tests with PyEphem (SGP) (/rhodesmill.org/pyephem/) and the PyPi SGP4 (/pypi.org/project/sgp4/)
2.) Is there a faster / lighter wait way to initialize this?
3.) Does anyone have any good sample codes for working with Orekit and UKFs beyond orbdetpy?
4.) Any other advise?
//////////////////////// Begin Code Samples //////////////////////////
/////////////////////////////// Prep Orekit //////////////////////////////////
import orekit
orekit.initVM()
#%% Setup Orekit working directory
from orekit.pyhelpers import setup_orekit_curdir
setup_orekit_curdir()
orekit.pyhelpers.download_orekit_data_curdir()
#%% imports for orekit & datetime
from java.util import Arrays
from orekit import JArray_double
from org.orekit.bodies import OneAxisEllipsoid, GeodeticPoint
from org.orekit.frames import FramesFactory, TopocentricFrame
from org.orekit.data import DataProvidersManager, ZipJarCrawler
from org.orekit.time import TimeScalesFactory, AbsoluteDate, DateComponents, TimeComponents
from org.orekit.orbits import KeplerianOrbit, CartesianOrbit
from org.orekit.utils import Constants
from org.orekit.propagation.analytical import EcksteinHechlerPropagator, KeplerianPropagator
from org.orekit.propagation.analytical.tle import TLE, TLEPropagator
from org.orekit.propagation.conversion import FiniteDifferencePropagatorConverter
from org.orekit.propagation.conversion import TLEPropagatorBuilder
from datetime import datetime
from org.orekit.propagation import SpacecraftState
from org.orekit.orbits import OrbitType, PositionAngle
from org.orekit.propagation.numerical import NumericalPropagator
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.hipparchus.ode.nonstiff import ClassicalRungeKuttaIntegrator
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from org.orekit.utils import IERSConventions, PVCoordinates, Constants
from java.io import File
#%% imports for matplotlib, math, and numpy
import matplotlib.dates as mdates
from math import radians, pi
import numpy as np
import matplotlib.pyplot as plt
import time
#%% Some constants, ae, mu, utc
ae = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
mu = Constants.WGS84_EARTH_MU
utc = TimeScalesFactory.getUTC()
/////////////////////////////// Numerical Prop ///////////////////////////////
#%% Initial osculating orbits
numsat = 1000
states = []
# Inertial frame where the satellite is defined
inertialFrame = FramesFactory.getEME2000()
PV=True
if PV:
for i in range(0,numsat):
epochDate = AbsoluteDate(2012, 1, 26, 16, 0, 00.000, utc)
# Orbit construction as Keplerian
p=Vector3D(-1780159.5047552094 + np.random.normal(0,1000),
-1164835.539234647 + np.random.normal(0,1000),
6842941.118987356 + np.random.normal(0,1000))
v=Vector3D(-7205.19342727098 + np.random.normal(0,1),
-341.0428975779148 + np.random.normal(0,1),
-1924.0083602739878 + np.random.normal(0,1))
initialOrbit = CartesianOrbit(PVCoordinates(p, v),
inertialFrame,
epochDate, mu)
satellite_mass = 100.0 # kg
initialState = SpacecraftState(initialOrbit, satellite_mass)
states.append(initialState)
print("states defined in cartesian coords")
else:
for i in range(0,numsat):
ra = 814 * 1000 + np.random.normal(0,1000) # km Apogee!
rp = 786 * 1000 + np.random.normal(0,1000) # km Perigee!
i = radians(98.55 + np.random.normal(0,0.5)) # inclination
omega = radians(90.0 + np.random.normal(0,0.5)) # perigee argument
raan = radians(5.1917 + np.random.normal(0,0.5)) # right ascension of ascending node
lv = radians(0.0567634 + np.random.normal(0,0.5)) # True anomaly
epochDate = AbsoluteDate(2012, 1, 26, 16, 0, 00.000, utc)
a = (rp + ra + 2 * ae) / 2.0 # semi major axis in km
e = 1.0 - (rp + ae) / a
# Orbit construction as Keplerian
initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lv,
PositionAngle.TRUE,
inertialFrame, epochDate, mu)
satellite_mass = 100.0 # kg
initialState = SpacecraftState(initialOrbit, satellite_mass)
states.append(initialState)
print("states defined in Keplerian elements")
#%% Set up propagator
integrator = ClassicalRungeKuttaIntegrator(30.0)
orbitType = OrbitType.CARTESIAN
propagator_num = NumericalPropagator(integrator)
propagator_num.setOrbitType(orbitType)
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True) # International Terrestrial Reference Frame, earth fixed
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
itrf)
gravityProvider = GravityFieldFactory.getNormalizedProvider(8, 8)
propagator_num.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider))
#%% Propagate one orbit numerically numsat steps forward
## Create time vector
startDate = AbsoluteDate(2012, 1, 26, 11, 0, 00.000, utc)
step_time = 30
# Time array in orekit AbsoluteDate format
dt = startDate.shiftedBy(float(step_time))
states2 = []
start = time.time()
for initialState in states:
propagator_num.setInitialState(initialState)
state = propagator_num.propagate(dt)
pv_org = state.getPVCoordinates()
states2.append(state)
end = time.time()
print("total time: ", end-start)
print("per sat time: ", (end-start)/numsat)
# total time: 18.845316886901855
# per sat time: 0.018845316886901855
#%% Propagate numsat orbits numerically forward one step
TestState = states[0]
start = time.time()
propagator_num.setInitialState(TestState)
dt = startDate
propagator_num.setInitialState(TestState)
for i in range(1000):
dt = dt.shiftedBy(float(step_time))
state = propagator_num.propagate(dt)
pv_org = state.getPVCoordinates()
end = time.time()
print("total time: ", end-start)
print("per sat time: ", (end-start)/numsat)
# total time: 0.07098388671875
# per sat time: 7.098388671875e-05
//////////////////////////// Keplarian Prop Only /////////////////////////////
#%%
inertialFrame = FramesFactory.getEME2000()
utc = TimeScalesFactory.getUTC()
initialDate = AbsoluteDate(2004, 1, 1, 23, 30, 00.000, utc)
mu = 3.986004415e+14
numsat = 1000
keplers = []
for step in range(numsat):
ra = 814 * 1000 + np.random.normal(0,1000) # km Apogee!
rp = 786 * 1000 + np.random.normal(0,1000) # km Perigee!
i = radians(98.55 + np.random.normal(0,0.5)) # inclination
omega = radians(90.0 + np.random.normal(0,0.5)) # perigee argument
raan = radians(5.1917 + np.random.normal(0,0.5)) # right ascension of ascending node
lv = radians(0.0567634 + np.random.normal(0,0.5)) # True anomaly
epochDate = AbsoluteDate(2012, 1, 26, 16, 0, 00.000, utc)
a = (rp + ra + 2 * ae) / 2.0 # semi major axis in km
e = 1.0 - (rp + ae) / a
# Orbit construction as Keplerian
initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lv,
PositionAngle.TRUE,
inertialFrame, epochDate, mu)
kepler = KeplerianPropagator(initialOrbit)
kepler.setSlaveMode()
keplers.append(kepler)
stepT = 60.
start = time.time()
for step in range(numsat):
#kepler = keplers[step]
extrapDate = initialDate.shiftedBy(stepT*(1+step))
currentState = kepler.propagate(extrapDate)
#print("step :", step+1)
#print(" time : ", currentState.getDate())
#print(" ", currentState.getOrbit())
end = time.time()
print("total time: ", end-start)
print("per sat time: ", (end-start)/numsat)
# total time: 0.008999824523925781
# per sat time: 8.99982452392578e-06'''