CR3BP Propagation - Python

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)

can u post your code as an attachment or in a code block? it will be convenient for others to find the problem.

yes, sorry, just did now. thank you for pointing that out!

for orekit’s cr3bp implementation, the normalized PV should be used as the initial input. the real PV is scaled back after the integration.

i do not fully setup the python code and run it. just a quick check, try it.

1 Like

FYI, https://forum.orekit.org/t/the-normalized-units-in-orekit-bodies-cr3bpsystem/4372

i found the orekit may be is wrong with the scale units in its cr3bp, but had no replies from the developers or any others. if u r interested, please check it out, and shared any suggestions.

1 Like

ok thank you! my issue was i was getting confused if normalized inputs were needed, thank you for clarifying that. I also realized my time wasnt normalized, once i changed that it got working!