ComputeMeanState with low inclination values

Hi all.
I have been trying to run a code which uses the function propagator.ComputeMeanState() to convert the osculating state to mean state.
Assuming the input state is 100% correct and that the code I’m using is also 100% correct (have been using for long time for many other simulations and it has always worked), I am facing an “orbit inside Brillouin sphere” error when it tries to convert an orbit of 5.2deg instead of the polar that I usually use.
So if I run it with a polar orbit (e.g. 97.5deg), the code works smoothly.

The orbit is: 580km Apo/Peri altitude, 5.2deg inclination (other params are irrelevant)
PV (EME2000): {2027-01-01T00:05:00.000, P(6588955.511785004, 2227138.6947943233, 202951.03949671745), V(-2432.0212411631164, 7137.77048133387, 649.4769855354031), A(-1.0059492320579517, -0.3430979974953157, 0.0026681858374665913)}.

I have tried to check the orbit PV and the orbit is correct, it really seems a problem of the ComputeMeanState function, when inclination values are low.
Is there a reason why?
Can you help me with this?

Thank you in advance,

Hi @jordan_93

Could you provide a runnable code in order to investigate your problem?

Thank you,

Hi Bryan,

thank you for the help.
Sharing with you the code would be extremely hard since it reserved material plus a very big simulator with a lot of functions.

The key however is this:

reference_state = osc_propagator.getInitialState()

mean_propagator = proplib.dsst_prop_builder(spacecraft_data=spacecraft_data, propagation_data=propagation_data)

mean_ref_state = mean_propagator.computeMeanState(
LofOffset(inertialFrame, LOFType.TNW, RotationOrder.YXY, pi/2, 0.0, 0.0),

It enters the computeMeanState with the correct state I posted before, and after some seconds, it exits with the error.

If I run it but with a different state of a normal polar orbit (haven’t tried with other orbits), it works.

import yaml
from typing import Union, Dict
from datetime import datetime, timezone
from math import pi
import libs.orekit_setup
from orekit.pyhelpers import datetime_to_absolutedate
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from import NRLMSISE00
from import MarshallSolarActivityFutureEstimation
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.propagation import SpacecraftState

Osculating Propagator

from org.orekit.propagation.numerical import NumericalPropagator

Mean Propagator

from org.orekit.propagation.semianalytical.dsst import DSSTPropagator
from org.orekit.propagation.semianalytical.dsst.forces import DSSTAtmosphericDrag, DSSTSolarRadiationPressure,
DSSTThirdBody, DSSTZonal, DSSTTesseral
from org.orekit.propagation import PropagationType
from org.orekit.orbits import KeplerianOrbit, PositionAngle, OrbitType
from org.hipparchus.geometry.euclidean.threed import RotationOrder
from org.orekit.frames import FramesFactory, LOFType
from org.orekit.bodies import CelestialBodyFactory
from org.orekit.utils import Constants, IERSConventions
from org.orekit.bodies import OneAxisEllipsoid
from org.orekit.attitudes import LofOffset


sun = CelestialBodyFactory.getSun()
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
FramesFactory.getITRF(IERSConventions.IERS_2010, True))
moon = CelestialBodyFactory.getMoon()
mu = Constants.WGS84_EARTH_MU
velocity_pointing = LofOffset(FramesFactory.getEME2000(), LOFType.TNW, RotationOrder.YXY, pi/2, 0.0, 0.0)

spacecraft_data = {‘info’: {‘platform’: ‘Ciaociao’, ‘name’: ‘Ciaociao’, ‘COSPAR_ID’: ‘2027-001G’, ‘NORAD_ID’: ‘88888’},
‘disturbances’: {‘SRP’: {‘crossSection’: 1.5, ‘reflectivityCoefficient’: 1.3}, ‘drag’: {‘crossSection’: 1.5, ‘dragCoefficient’: 2.2}, ‘geometry’: {‘spacecraft_size_x’: 1.0509, ‘spacecraft_size_y’: 1.19269, ‘spacecraft_size_z’: 1.1}},
‘platform’: {‘mass_dry_empty’: 277.486, ‘cog_dry_empty’: [-0.00363, -0.00201, 0.45725], ‘I_dry_empty’: [64.8, 0.311, -0.5, 0.311, 68.9, 0.219, -0.502, 0.219, 22.8]},
‘deployables’: {‘1’: {‘mass’: 50.4, ‘cog’: [-0.00135, -0.00115, 0.05935], ‘cgeom’: [0.25425, -0.10375, 0.87176], ‘I’: [0.0113, 0.004, 0.0113], ‘disturbances’: {‘SRP’: {‘crossSection’: 0.03, ‘reflectivityCoefficient’: 1.8}, ‘drag’: {‘crossSection’: 0.03, ‘dragCoefficient’: 2.2}}}},
‘propellant’: {‘reference_date’: datetime(2022, 11, 17, 0, 0, tzinfo=timezone.utc), ‘mass_ox’: 28.8, ‘mass_fu’: 6.9, ‘cog_ox’: [0, 0, 0.34353], ‘cog_fu’: [0, 0, 0.35067], ‘I_ox’: [4.14, 0, 0, 0, 4.14, 0, 0, 0, 0.931], ‘I_fu’: [2.14, 0, 0, 0, 0.861, 0, 0, 0, 1.31]}
propagation_data = {‘computeMeanState’: {‘convergence_threshold’: ‘1E-12’, ‘max_iteration’: ‘1E13’},
‘integrator’: {‘minStep’: 1e-07, ‘maxStep’: 2700.0, ‘absTol’: 1e-11, ‘relTol’: 1e-11},
‘gravity’: {‘harmonics_degree’: 64, ‘harmonics_order’: 64}, ‘drag’: {‘solar_activity_level’: ‘AVERAGE’}, ‘elevationDetector’: {‘maxcheck’: 60.0, ‘threshold’: 0.001},
‘attitude’: {‘mode’: ‘NONE’, ‘details’: {‘mode’: ‘NONE’, ‘details’: {‘sunPointingAxis’: ‘+X’, ‘phasingAxis’: ‘+Z’}}},
‘propagator_type’: {‘cubesat_deployment_analysis’: ‘OSCULATING’, ‘orbit_transfer_analysis’: ‘OSCULATING’}}

def _add_dsst_force_models(dsst_propagator: DSSTPropagator, spacecraft_data, propagation_data) → DSSTPropagator:
Adds semi-analytical force models to DSST Propagator. Drag, SRP and Gravity models are configured according to
spacecraft data and propagation configurations
:param dsst_propagator: DSST propagator where forces will be added
:param spacecraft_data: dictionary or path to yml containing ‘drag’ (‘dragCoefficient’ and ‘crossSection’) and ‘SRP’
(‘reflectivityCoefficient’ and ‘crossSection’) data
:param propagation_data: dictionary or path to yml containing ‘gravity’ (‘harmonics_degree’ and ‘harmonics_order’) data
:return: DSST Propagator with Drag, SRP, Zonal Harmonics, Moon and Sun Third-Body force models
propagation_conf = propagation_data
spacecraft_conf = spacecraft_data
cd = spacecraft_conf[‘disturbances’][‘drag’][‘dragCoefficient’]
cr = spacecraft_conf[‘disturbances’][‘SRP’][‘reflectivityCoefficient’]
drag_area = spacecraft_conf[‘disturbances’][‘drag’][‘crossSection’]
srp_area = spacecraft_conf[‘disturbances’][‘SRP’][‘crossSection’]
harmonics_degree = propagation_conf[‘gravity’][‘harmonics_degree’]
harmonics_order = propagation_conf[‘gravity’][‘harmonics_order’]
if propagation_conf[‘drag’][‘solar_activity_level’] == ‘AVERAGE’:
strength_level = MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE
elif propagation_conf[‘drag’][‘solar_activity_level’] == ‘STRONG’:
strength_level = MarshallSolarActivityFutureEstimation.StrengthLevel.STRONG
elif propagation_conf[‘drag’][‘solar_activity_level’] == ‘WEAK’:
strength_level = MarshallSolarActivityFutureEstimation.StrengthLevel.WEAK
print(“Unrecognized setting for solar_activity_level in propagation_conf. Setting AVERAGE.”)
strength_level = MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE

supported_names = MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES
input_data = MarshallSolarActivityFutureEstimation(supported_names, strength_level)
dsst_atmosphere = NRLMSISE00(input_data, sun, earth)
dsst_drag = DSSTAtmosphericDrag(dsst_atmosphere, cd, drag_area, mu)

dsst_srp = DSSTSolarRadiationPressure(cr, srp_area, sun, R_E, mu)

# Third Body
dsst_moon = DSSTThirdBody(moon, mu)
dsst_sun = DSSTThirdBody(sun, mu)

# Zonal and Tesseral Harmonics
provider = GravityFieldFactory.getUnnormalizedProvider(harmonics_degree, harmonics_order)
dsst_zonal = DSSTZonal(provider)

dsst_tesseral = DSSTTesseral(FramesFactory.getITRF(IERSConventions.IERS_2010, True), w_E, provider)
return dsst_propagator

minStep = 0.0000001 # [s] smallest integration step allowed
maxStep = 2700.0 # [s] longest integration step allowed
absTol = 9.999999999999999e-12 # absolute tolerance on integration error
relTol = 9.999999999999999e-12 # relative tolerance on difference between solutions
integrator = DormandPrince853Integrator(minStep, maxStep, absTol, relTol)
dsst_propagator_temp = DSSTPropagator(integrator, PropagationType.MEAN)
propagator = _add_dsst_force_models(dsst_propagator_temp, spacecraft_data, propagation_data)
orbit = KeplerianOrbit(float(6958000), float(0.001), float(0.090757), float(0.01), float(0.01), float(0.01), PositionAngle.TRUE,
FramesFactory.getEME2000(), datetime_to_absolutedate(datetime.strptime(“2027-01-01T00:05:00.000Z”,“%Y-%m-%dT%H:%M:%S.%fz”)), mu)
orbit = OrbitType.convertType(OrbitType.EQUINOCTIAL,orbit)
state = SpacecraftState(orbit)
mean_state = propagator.computeMeanState(state, velocity_pointing, propagator.getAllForceModels(),

@bcazabonne did it work?


I don’t know if this will help but in your code on my machine (orekit 11.3.3) exit with with error et

mean_state = propagator.computeMeanState(state, velocity_pointing, propagator.getAllForceModels(),

JavaError: <super: <class ‘JavaError’>, >
Java stacktrace:
org.orekit.errors.OrekitIllegalArgumentException: hyperbolic orbits cannot be handled as org.orekit.orbits.EquinoctialOrbit instances
at org.orekit.orbits.EquinoctialOrbit.(
at org.orekit.orbits.OrbitType$3.mapArrayToOrbit(
at org.orekit.orbits.OrbitType$3.mapArrayToOrbit(
at org.orekit.propagation.semianalytical.dsst.DSSTPropagator.computeOsculatingOrbit(
at org.orekit.propagation.semianalytical.dsst.DSSTPropagator.computeMeanOrbit(
at org.orekit.propagation.semianalytical.dsst.DSSTPropagator.computeMeanState(

Hi, unfortunately it doesn’t help. It still proofs that ComputeMeanState is failing somehow.
I have also tried with equinotial orbit but the reult is still the same.

Anyone able to help?

Not sure this is of any help, but you integration steps bounds seem odd to me.
DSST is a semi-analytical model, and the integration step corresponds to the mean elements which vary very slowly. When using DSST, integration steps are more often several hours (much larger than orbital period), or even one day.

I’m not sure that changing the integration step will fix the problem because it is not involved in the oscillating to mean conversion.
However, Luc’s remark is important is you perform DSST-based orbit propagation. That’s the key point of the DSST compared to a numerical propagation: because it integrates numerically the mean elements, you can use big steps! For LEO satellites, you can have a min step equal to the orbital period and a max step of about 10 orbital periods. For GEO, a fixed step of 12 hours is recommended.

Back to your problem, did you try adding the force models sequentially in order to highlight the first failing configuration. You can try to (1) use Zonal/Tesseral, (2) conf (1) + third body attraction, (3) conf (2) + drag, and (4) conf (3) + SRP.
I guess the exception will be throw at configuration (3) but it is interesting to confirm this assumption.

Also, what is your Orekit version? I’m surprised that you don’t have the same exception as @AlexS . I’m asking that because Orekit 11.3.3 fix a lot of DSST-related issues. It is a very important version for Orekit DSST. Please see the release notes. I understand that according to @AlexS there is still an exception, but it is better to fix the same exception.

Regarding your issue, I’ve a last comment.
Because computeMeanState() is a static method, it doesn’t depend on DSST instantiation. In other words, it depends on the method parameters. Therefore, you can include the call of this method into a try/catch. Like that you can catch the exception and use another osculating to mean conversion method if the first one didn’t work.
I like to say that they are hundreds of oscillating to mean conversion methods. Brouwer-Lyddane is one, Eckstein-Hechler is another one, Dst J2-only is also another one, etc. Therefore, if the first one you use doesn’t works, catch the exception and use another one. Depending the results of configurations (1), (2), (3), and (4) you can use the “last one that worked”.

To be honest, I don’t see issues in your code. Therefore, the only solution I see is the try/catch method with different oscillating to mean conversion methods.

Best regards,

1 Like

A last thing to test: could you try using HarrisPriester instead of NRLMSISE00?


The problem is solved removing only the tesseral contribution.
It is also OK with 16x16 potential and all forces…

1 Like

Indeed, changing to 32X32 potential (all forces) produces:
<SpacecraftState: SpacecraftState{orbit=equinoctial parameters: {a: 6957951.422865678; ex: -3.7341735341071933E-4; ey: 2.3961789334324766E-5; hx: 0.045368108300881727; hy: 4.4019287197056135E-4; lv: 1.7125950633725817;}, attitude=org.orekit.attitudes.Attitude@23e44287, mass=1000.0, additional={}, additionalDot={}}>