STM propagation

Orekit Python Wrapper usage
1

Hey all,

I am working on an OD project where I write my own estimation code. I am now working on the classic batch least squares algorithm and I am running into some issues because it does not seem to converge properly. It used to work fine until I changed my propagation method. First I used solve_ivp and the equations of motion to propagate the state and state transition matrix, now I switched to the STM propagator from orekit because of a higher fidelity dynamics model and it does not seem to converge properly.

I was wondering if I use the STM propagator properly. This is my code:
the ‘stm’ input is a 1D array of the STM and time is a two element array with the initial time and final time of propagation. I use this function in the prediction step.

def propagate_STM(stm, time):
    n = 6
    state = stm[:n]
    phi_mat = stm[n:].reshape((n,n))  

    STM_hipp = MatrixUtils.createRealMatrix(6,6)
    
    for i in range(6):
        for j in range(6):
            STM_hipp.setEntry(i,j, float(phi_mat[i,j]))
    
    initial_state = state2TPV(state, time[0])
       
    t0_absolute = datetime_to_absolutedate(time[0])
    t1_absolute = datetime_to_absolutedate(time[1])
    
    initialOrbit = CartesianOrbit(initial_state, frames['eci'], EARTH_MU)
    initialSpacecraftState = SpacecraftState(initialOrbit, satellite['mass'])
    
    integrator = getDORPIIntegrator(initialOrbit)

    
    if propagatorSettings['propagation_method'] == 'Kepler':
        propagator = getKeplerPropagator(integrator, initialSpacecraftState)
    elif propagatorSettings['propagation_method'] == 'Advanced':
        propagator = getAdvancedPropagator(integrator, initialSpacecraftState)
    
    
    harvester = propagator.setupMatricesComputation("STM", STM_hipp, None)
    
    time_array = [t0_absolute, t1_absolute]
    
    state = [propagator.propagate(tt) for tt in time_array]
    propagated_state = np.array(
    [state[-1].getPVCoordinates().getPosition().x,
     state[-1].getPVCoordinates().getPosition().y,
     state[-1].getPVCoordinates().getPosition().z,
     state[-1].getPVCoordinates().getVelocity().x,
     state[-1].getPVCoordinates().getVelocity().y,
     state[-1].getPVCoordinates().getVelocity().z])

    stm = harvester.getStateTransitionMatrix(state[-1])
    stm_np = np.zeros((n,n))
    for i in range(n):
        for j in range(n):
            stm_np[i,j] = stm.getEntry(i,j)
    
    stm_out = np.hstack([propagated_state.reshape(1,-1), stm_np.reshape((1, n**2))])
    return stm_out

If anything is unclear, please let me know! Thanks for your help in advance!

Cheers,

Niek

2

Hi there,

In what coordinates do you want your STM? Cartesian?
Here you’re getting it whatever OrbitType your propagator uses. If you didn’t change it, the default is EQUINOCTIAL.

Cheers,
Romain.

3

Hi Romain,

The propagator is defined using a the cartesian orbit type:

def getAdvancedPropagator(integrator, initialSpacecraftState):
    propagator = NumericalPropagator(integrator)
    propagator.setOrbitType(OrbitType.CARTESIAN)
    propagator.setInitialState(initialSpacecraftState)
    
    isotropic_drag = IsotropicDrag(satellite['A'], satellite['CD'])
    drag = DragForce(environment['atmosphere'], isotropic_drag)
    gravityProvider = GravityFieldFactory.getNormalizedProvider(10,10)
    radiation = IsotropicRadiationClassicalConvention(satellite['A'], satellite['abs_coef'], satellite['CR'])
    srp = SolarRadiationPressure(environment['sun'], environment['earth'], radiation)
    
    #   Add forces to propagator
    propagator.addForceModel(HolmesFeatherstoneAttractionModel(frames['ecef'], gravityProvider))
    propagator.addForceModel(ThirdBodyAttraction(environment['moon']))    
    propagator.addForceModel(ThirdBodyAttraction(environment['sun']))
    propagator.addForceModel(ThirdBodyAttraction(environment['mars']))
    propagator.addForceModel(ThirdBodyAttraction(environment['jupiter']))
    propagator.addForceModel(ThirdBodyAttraction(environment['saturn']))
    propagator.addForceModel(drag)
    propagator.addForceModel(srp)
    return propagator

Do I have to specify it elsewhere?

Cheers,
Niek

4

No it’s fine

5

Great, does that mean my implementation is fine this way? Or do you have any other suggestions why it might struggle to converge?

6

Does your OD converge if you use Keplerian, analytical motion? In other words is your problem specific to the numerical propagator?