Implementing a state transition matrix (STM) in python

Hi Orekit users,

I would like to implement a state transition matrix (STM) for the propagation.

Now I have my own function for computing the STM.

However, I’d like to use classes, such as ‘PartialDerivativesEquations’, ‘JacobiansMapper’, etc., provided by Orekit instead of using my function.

So far, my trials were not successful and some examples (Java) seemed not to work for me (I’m using python wrapper.)

I wonder if I can find better examples of how to get the STM in python.

Thank you for any suggestions!.

:: Inkwan

Hi,

Could you describe in what way it did not work, was it the classes that was not visible in the Python environment or did you get some error in the input/output? Could you provide a pointer to the java example you are refering to?

Not sure if it is related but if you need to make java matrixes in python there are ways to do that, but requires a little bit of caution. See the example of numerical propagation at
orekit gitlab or in the pyhelpers file.

Regards

Hi, thank you for the response.

I think now I can get a state transition matrix (STM) from orekit, but am not sure if it’s a usual way.

First of all, here’s a link for the java examples.

https://www.orekit.org/mailing-list-archives/orekit-users/msg00215.html

https://www.orekit.org/forge/projects/orekit/repository/revisions/master/entry/src/main/java/org/orekit/propagation/numerical/JacobiansMapper.java#L140

One of the biggest problems was how to define 2d-array in python. And I solved it by using ‘JArray’, which would be in the second link you suggested.

The other issue was to get the state transition matrix.

I tried to use the way suggested in the first java example.

Array2DRowRealMatrix dYdY0 = new Array2DRowRealMatrix(6, 6);
PDE.getMapper().getStateJacobian(finalState, dYdY0.getDataRef());

But I got ‘Invalid input argument error.’ Even though I’ve gotten the STM, I could not use ‘getMapper()’ yet.
(Sorry I don’t have an old version, which gave me the error message.)

Here’s what I’ve done in Python:

  1. define a numerical propagator (two-body dynamics):
    def get_numerical_propagator(initial_state):
    # Propagator min step (s), max step (s), position error (m) and normalization scale (m) [0.001, 300, 10.0, 1000.0]
    min_step = float(0.001)
    max_step = float(60.0)
    position_error = float(10.0)

     # numerical propagator
     t = NumericalPropagator.tolerances(position_error, initial_state, initial_state.getType())
     # integrator
     tols = [orekit.pyhelpers.JArray('double').cast_(x) for x in t]
     dpIntegrator = DormandPrince853Integrator(min_step, max_step, tols[0], tols[1])
    
     num_propagator = NumericalPropagator(dpIntegrator)
     ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
     earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, ecef)
     # for gfile in files:
     #     GravityFieldFactory.addPotentialCoefficientsReader(ICGEMFormatReader(gfile, True))
     gravity_field = GravityFieldFactory.getNormalizedProvider(1, 0)
     attraction_model = HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravity_field)
     num_propagator.addForceModel(attraction_model)
     num_propagator.setOrbitType(initial_state.getType())
     return num_propagator
    
  2. get the STM (part of the function):`
    while (dt_sign*(extrap_date.compareTo(fin_date)) <= 0.0):
    # Numerical propagation
    tmp_fin_state = propagator.propagate(extrap_date)
    tmp_stm = np.array([[y for y in JArray(‘double’).cast_(x)] for x in (tmp_fin_state.additionalStates).values()])

What I tried was to extract additional states, i.e., integrated jacobian matrix, for getting the STM.

But still, I’m wondering how to use JacobianMapper class properly in Python.

Thank you again for your concern.

Best,

Hi,

I haven’t used Jacobians but will try to make an example.

I think the form of the setInitialJacobian as used in the example is depriciated (initialstate, 6, 6) so if this is where the problem is this is likely the cause.

Note that the python gateway to orekit can be more restrictive on the input parameters than java in order to find the right method call, if the method call states it shall be a float it needs to be a float and not an int, and if it needs an matrix one needs to provide such (by casting). I think if one uses orekit from inside java this is more automatic. Usually not a problem but when one gets the invalid args this is sometimes the case.

There are some examples of usage in java in the test case:
https://www.orekit.org/static/apidocs/org/orekit/propagation/numerical/PartialDerivativesEquations.html#setInitialJacobians-org.orekit.propagation.SpacecraftState-

Hi,

My apologies that my response is too late. Thank you very much for your response.

Yes. I totally agree with you. Sometimes I spent pretty long time to figure out how to define inputs in python.

Using java directly would be even easier to use orekit.

So far, the stm seems to work as I expected, but I’ll have more tests how to use jacobian as well.

Thank you again for your concern. I appreciate it.

Have a great day!

Best,

:: Inkwan