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/mailinglistarchives/orekitusers/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 2darray 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:

define a numerical propagator (twobody 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

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,