STM propagation with Orekit (python wrapped version)

Orekit Python Wrapper usage
1

Hi Everyone! First time posting on here, but long time user of Orekit :slight_smile:

I am trying to propagate an STM using orekit (python wrapped version) and I’m running into an error (which I’ve spent hours trying to debug/resolve with no avail). The error occurs when I call setInitialJacobians(). Below is the code; everything works up until the aforementioned line…

Thanks in advance for any help!
-Shaylah

gravity_field = GravityFieldFactory.getNormalizedProvider(int(70), int(70))
earth = OneAxisEllipsoid(gravity_field.getAe(),
Constants.WGS84_EARTH_FLATTENING,
FramesFactory.getITRF(IERSConventions.IERS_2010, True))

Set up the initial date

initial_date = AbsoluteDate(2023, 4, 21, 0, 0, 0.0, TimeScalesFactory.getUTC())

Set up the inertial frame

inertial_frame = FramesFactory.getEME2000()

Create a Cartesian orbit

orbitPV = np.array([1813.337619, 7203.761092, -90.135766,
-0.708160485, 0.786128638, 7.181618172])*1e3
position = Vector3D(float(orbitPV[0]),float(orbitPV[1]),float(orbitPV[2]))
velocity = Vector3D(float(orbitPV[3]),float(orbitPV[4]),float(orbitPV[5]))
pvCoordinates = PVCoordinates(position, velocity)
orbit = CartesianOrbit(pvCoordinates, inertial_frame, initial_date, Constants.WGS84_EARTH_MU)

Set up the integrator

min_step = 10.0 # seconds
max_step = 300.0 # seconds
absolute_tolerance = 1.0e-4 # Absolute tolerance
relative_tolerance = 1.0e-4 # Relative tolerance
integrator = DormandPrince853Integrator(min_step, max_step, absolute_tolerance, relative_tolerance)

Set up the numerical propagator

propagator = NumericalPropagator(integrator)
propagator.setOrbitType(OrbitType.CARTESIAN)

propagator.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravity_field))

partialDerivatives = PartialDerivativesEquations(“STM”, propagator)
state_dimension = 6 # state vector dimension (position and velocity)
initialDyDy0 = Array2DRowRealMatrix(state_dimension, state_dimension)

initial_state_with_derivatives = partialDerivatives.setInitialJacobians(initial_state, initialDyDy0) #initialDyDP = None b/c no param estimated

Error for the last line:
InvalidArgsError: (<class ‘org.orekit.propagation.numerical.PartialDerivativesEquations’>, ‘setInitialJacobians’, (<SpacecraftState: SpacecraftState{orbit=Cartesian parameters: {P(1813337.619, 7203761.092, -90135.766), V(-708.160485, 786.128638, 7181.618172)}, attitude=org.orekit.attitudes.Attitude@59942b48, mass=10.0, additional={}, additionalDot={}}>, <Array2DRowRealMatrix: Array2DRowRealMatrix{{0.0,0.0,0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0,0.0,0.0}}>))

2

I use this code snippet to initialize covariance propagation, so I’m not 100% sure that this one snippet alone will solve your problem. But here goes nothing:

# Create NumericalPropagator object
backProp = numericPropBuilder. \
    buildPropagator(numericPropBuilder.getSelectedNormalizedParameters())
  
# Add calculation of STM to propagator
# Note: use of None keyword for input 2 will initialize STM as
# an identity matrix.  Still not 100% sure about input 3.
harvester = backProp.setupMatricesComputation("stm", None, None)

# Figured I might as well leave the covariance code in

# Create covariance
p_m_m = # 6x6 RealMatrix object

# Append covariance to propagator
covInit = StateCovariance(p_m_m, startDate, frame, OrbitType.CARTESIAN, PositionAngleType.TRUE)
provider = StateCovarianceMatrixProvider("covariance", "stm", harvester, covInit)
backProp.addAdditionalStateProvider(provider)

Link to method description

1 Like
3

Thank you! This was immensely helpful!

Here’s what ultimately ended up working (based on your snippet of code and extended):

partialDerivatives = PartialDerivativesEquations(“STM”, propagator)

Set initial state

initial_state = SpacecraftState(orbit, sat_mass)

initial_state_with_derivatives = partialDerivatives.setInitialJacobians(initial_state)
propagator.setInitialState(initial_state_with_derivatives)

Add calculation of STM to propagator

harvester = propagator.setupMatricesComputation(“stm”, None, None)

Create covariance

p_m_m = Array2DRowRealMatrix(6, 6)# 6x6 RealMatrix object

Append covariance to propagator

covInit = StateCovariance(p_m_m, initial_date, eciFrame, OrbitType.CARTESIAN, PositionAngle.TRUE)
provider = StateCovarianceMatrixProvider(“covariance”, “stm”, harvester, OrbitType.CARTESIAN, PositionAngle.TRUE, covInit)
propagator.addAdditionalStateProvider(provider)

current_state = propagator.propagate(target_date)
stm = current_state.getAdditionalState(“stm”)