Hi,
I’m currently using org.orekit.propagation.numerical.cr3bp to write a function for CR3BP propagation in Python.
I’ve seen many examples written in Java, but I’m working in Python and have been struggling to find any Python-specific examples on using the CR3BP class.
My main question is whether I’m setting up the propagation correctly or possibly misusing the class because I’m getting trajectory results that are way off and clearly incorrect. I based my code on this thread and the Java code shared there: Propagation CR3BP
My current code :
def CR3BPropagator(initialState,
cr3bpSystem: CR3BPSystem,
attitudeProvider = None,
stateCovariance = None,
integration = Integrator()):
# Validate input state type
if not isinstance(initialState, SpacecraftState):
raise TypeError("initialState must be a SpacecraftState object")
# Verify state is in barycenter rotating frame
if not initialState.getFrame().equals(cr3bpSystem.getRotatingFrame()):
raise ValueError("Initial state must be in the barycenter rotating frame")
normalizedState = initialState
# Create the numerical propagator with non-Keplerian settings
propagator = NumericalPropagator(integration)
# 1. Set orbit type to None for non-Keplerian system
propagator.setOrbitType(None) # Non-Keplerian system
propagator.setIgnoreCentralAttraction(True) # CR3BP doesn't use central attraction
# 2. Add force models
propagator.addForceModel(CR3BPForceModel(cr3bpSystem))
# 3. Set up matrices computation if covariance is needed
scmProvider = None
if stateCovariance is not None:
if not isinstance(stateCovariance, StateCovariance):
raise TypeError("stateCovariance must be a StateCovariance object")
# Check whether the state and covariance frames match
if not stateCovariance.getFrame().equals(initialState.getFrame()):
raise ValueError("State and Covariance must be in the same frame")
# Check whether the state and covariance are defined at the same epoch
if not stateCovariance.getDate().equals(initialState.getDate()):
warnings.warn("State and Covariance are not defined at the same Epoch.\nTime-shifting covariance to state epoch using linearized model.", stacklevel=3)
stateCovariance = stateCovariance.shiftedBy(initialState.getDate())
# Set up matrices computation before setting initial state
stm0 = MatrixUtils.createRealIdentityMatrix(6)
scmHarvester = propagator.setupMatricesComputation("stm", stm0, None)
if scmHarvester is None:
raise RuntimeError("STM setup failed — scmHarvester is None")
# Now set initial state (with STM embedded)
propagator.setInitialState(normalizedState)
# Register covariance propagation
scmProvider = StateCovarianceMatrixProvider("covariance", "stm", scmHarvester, stateCovariance)
else:
# No covariance: just set the initial state directly
propagator.setInitialState(normalizedState)
# 5. Set attitude provider if provided (optional)
if attitudeProvider is not None:
if not isinstance(attitudeProvider, AttitudeProvider):
raise TypeError("attitudeProvider must be an AttitudeProvider object")
propagator.setAttitudeProvider(attitudeProvider)
return (propagator, scmProvider)
A side note, I’m also having trouble setting up the STM when I try to use it. I keep running into a NullPointerException, which might be a separate issue, but I wanted to mention it as well if anyone had any insight if my set up for STM is incorrect.
Here is how I am calling my function:
utc = TimeScalesFactory.getUTC()
earth = CelestialBodyFactory.getEarth()
moon = CelestialBodyFactory.getMoon()
a = 384400e3 # meters
cr3bpSystem = CR3BPSystem(earth, moon, a)
initialEpoch = AbsoluteDate(2025, 5, 20, 12, 0, 0.0, utc)
mass = 1000.0 # kg
# JPL Halo Orbit around L1 (normalized units)
x_norm = 8.2339081983651485E-1
y_norm = -1.9017764504099543E-28
z_norm = 9.8941366235910004E-4
vx_norm = -2.3545391932685812E-15
vy_norm = 1.2634272983881797E-1
vz_norm = 2.2367029429442455E-16
normalizedPV = AbsolutePVCoordinates(
cr3bpSystem.getRotatingFrame(), # CR3BP rotating frame
initialEpoch,
Vector3D(x_norm, y_norm, z_norm),
Vector3D(vx_norm, vy_norm, vz_norm))
# Convert to physical units using Orekit
physicalPV = cr3bpSystem.getRealAPV(
normalizedPV,
initialEpoch,
cr3bpSystem.getRotatingFrame())
initialState = SpacecraftState(physicalPV, mass)
minStep = .5
maxStep = 1000.0
positionTolerance = 1.0e-3
integrator = DormandPrince853Integrator(minStep, maxStep, [positionTolerance] * 3, [positionTolerance] * 3)
# State Covariance
covMat = DiagonalMatrix([1e4, 1e4, 1e4, 0.01, 0.01, 0.01]) # m^2 and m^2/s^2
covCR3BP = StateCovariance(covMat, initialEpoch, cr3bpSystem.getRotatingFrame(), OrbitType.CARTESIAN, PositionAngleType.TRUE)
cr3bProp,_ = CR3BPropagator(initialState=initialState,
cr3bpSystem=cr3bpSystem,
stateCovariance=covCR3BP)
# Propagate for one orbit
T = cr3bpSystem.getTdim()
time = [initialEpoch.shiftedBy(float(dt)) for dt in np.arange(0.0, T, 600)] # 10-minute steps
cr3bpStates = [cr3bProp.propagate(t) for t in time]
But my plot is clearly incorrect:
I know this is a long post, but I’ve been working through this for a while and would really appreciate any insight. If anyone has experience using the CR3BP class in Python or has run into similar issues, any help would be greatly appreciated.Thank you!!
(Using orekit 13.0.1 and python 3.12.2)