# State transition matrix does not resemble I + F dt

I’m trying to produce a state transition matrix for a propagation. I’ve looked through most of the other posts on STMs and Jacobians and can’t seem to identify my problem.

``````class MyHandler(PythonOrekitFixedStepHandler):
mapper         = None

def init(self, x0, t, step):
self.x = []
self.t = []
self.jPhi = Array2DRowRealMatrix(6,6)

def handleStep(self, x0, is_last):
pv = x0.getPVCoordinates()
t  = pv.getDate().durationFrom(AbsoluteDate.J2000_EPOCH)
p  = pv.getPosition()
v  = pv.getVelocity()
self.x.append( np.array([p.getX(), p.getY(), p.getZ(),
v.getX(), v.getY(), v.getZ()]) )
self.t.append(t)

if is_last:
self.x = np.vstack(self.x)
self.t = np.array(self.t)
self.mapper.getStateJacobian(x0, self.jPhi.getDataRef())
print("Phi = {}".format(self.Phi))

@property
def Phi(self):
# convert to Numpy array
return orekit_matrix_to_ndarray(self.jPhi)
``````

The above is my handler. Below is where I create my propagator:

``````def create_propagator(t0, x0,
handler        = None,
min_step       = 0.001,
max_step       = 300.0,
rtol           = 1e-15,
atol           = 1e-9,
fixed_step     = 60.0,
dP             = None,
gravity_degree = 20,
gravity_order  = 20,
req            = None,
flattening     = None):
gravity_field  = GravityFieldFactory.getNormalizedProvider(gravity_degree, gravity_order)

j2000          = FramesFactory.getEME2000()
orbit          = CartesianOrbit(x0, j2000, t0, gravity_field.getMu())

if dP is not None: # Compute absolute and relative tolerances
tols = NumericalPropagator.tolerances(dP, orbit, OrbitType.CARTESIAN)
atol = orekit.JArray_double.cast_(tols[0])
rtol = orekit.JArray_double.cast_(tols[1])

integrator     = DormandPrince853Integrator(min_step, max_step, atol, rtol)
propagator     = NumericalPropagator(integrator)

pde = PartialDerivativesEquations("dYdY0", propagator)
initial_state = pde.setInitialJacobians(SpacecraftState(orbit))

propagator.setInitialState(initial_state)
propagator.setOrbitType(OrbitType.CARTESIAN)
if handler is not None:
handler.mapper = pde.getMapper()
propagator.setMasterMode(fixed_step, handler)

return propagator, pde

def propagate(t0, x0, tf, **kwargs):
handler                = MyHandler()
propagator, pde = create_propagator(t0, x0, handler, **kwargs)
final_state = propagator.propagate(tf)

return handler, final_state
``````

Finally, here’s my main():

``````if __name__ == '__main__':

moon = CelestialBodyFactory.getBody("MOON")
j2000 = FramesFactory.getEME2000()

# Initial state
t0 = orekit_time(708687952.5569172)
#tf = orekit_time(709099110.5780709)
tf = orekit_time(708689252.5569172)

x0_ = np.array([-6.45306258e+06, -1.19390257e+06, -8.56858164e+04,
1.83609046e+03, -9.56878337e+03, -4.95077925e+03])
dx1 = np.array([1.0, 0.0, 0.0, 0.0, 0.0, 0.0])
x01_ = x0_ + dx1

x0 = orekit_state(x0_)
x01 = orekit_state(x01_)

handler, xf = propagate(t0, x0, tf, dP = 0.001)
handler1, xf1 = propagate(t0, x01, tf, dP = 0.001)

dxf_pred = handler.Phi.dot(dx1)
dxf      = handler1.x[-1,:] - handler.x[-1,:]

print("dxf      = {}".format(dxf))
print("dxf pred = {}".format(dxf_pred))
``````

Now, this should produce a state transition matrix at the first step which resembles I + F dt. It does not.

It should also produce a dxf_pred which resembles dxf, but they’re really different.

Here is the output from the final step:

``````dxf      = [1.77740174e+00 1.11202185e+00 4.38334891e-01 5.37267846e-04
1.55869407e-03 6.73296128e-04]
dxf pred = [-6.12512929e-02 -7.83017986e-03  9.45141750e-04 -1.60447077e-05
-1.58988751e-05 -6.70553842e-06]
Phi      = [[-6.12512929e-02 -3.54320664e+08 -7.43803464e+07 -7.09772683e+05
-5.03871676e+06  3.87053544e+06]
[-7.83017986e-03 -4.14237366e+07 -1.11108263e+07  9.50666885e+06
5.15800605e+04 -2.68742312e+06]
[ 9.45141750e-04  7.41131148e+06  3.43213978e+05 -1.89151263e+07
-2.17965771e+05 -1.65755263e+06]
[-1.60447077e-05 -9.15390433e+04 -2.27524862e+04 -3.17616115e+02
-3.21761682e+03 -1.96602471e+01]
[-1.58988751e-05 -8.78025347e+04 -2.09868772e+04  4.76319673e+03
2.88483542e+03  1.70612129e+03]
[-6.70553842e-06 -3.67983823e+04 -8.72621438e+03 -8.46432817e+03
-1.21906677e+04  8.58227654e+02]]
``````

And here is from the first step:

``````Phi = [[-3.27572418e-02 -1.90349163e+08 -3.58331050e+07 -8.47356866e+04
3.72311932e+05  1.10261795e+06]
[-6.06052625e-03 -3.52171318e+07 -6.62960193e+06  6.20065411e+05
-2.88833030e+06 -5.74629223e+06]
[-4.34961070e-04 -2.52751670e+06 -4.75803359e+05 -2.25816835e+06
1.22054362e+07 -2.97306598e+06]
[-4.66021044e-06 -2.65279578e+04 -1.05537109e+04 -6.79131990e+02
-5.08109339e+03  5.46386936e+03]
[ 2.42866815e-05  1.43705299e+05  2.60237835e+04  9.23367695e+03
8.21816872e+02  1.01088864e+03]
[ 1.25656517e-05  7.42681435e+04  1.39071039e+04 -1.80985654e+04
-3.47281384e+03  7.25509944e+01]]
``````

Any ideas about what I might be doing wrong? Am I misunderstanding the PDE? (Does it produce the dynamics Jacobian instead of the state transition matrix?)

Hi @autumnsault,

This one was a little tricky. The issue is here:

You are initialiazing the PDE before setting orbit type to Cartesian.
Thus orbit type in the `JacobianMapper` is set to equinoctial. Indeed it takes its orbit type from the propagator, which, by default, is set to equinoctial.
Thus some conversions between your Cartesian initial state and the equinoctial type are done in the STM.

To fix your issue you need to put the line

``````propagator.setOrbitType(OrbitType.CARTESIAN)
``````

Before creating the PDE, so before

``````pde = PartialDerivativesEquations("dYdY0", propagator)
``````

With this I get the following result after 60s
(I cut matrix Phi at 60s because it is too long):

``````t = 0.0 / Phi = Array2DRowRealMatrix{{1.0,0.0,0.0,0.0,0.0,0.0},{0.0,1.0,0.0,0.0,0.0,0.0},{0.0,0.0,1.0,0.0,0.0,0.0},{0.0,0.0,0.0,1.0,0.0,0.0},{0.0,0.0,0.0,0.0,1.0,0.0},{0.0,0.0,0.0,0.0,0.0,1.0}}
t = 60.0 / Phi = Array2DRowRealMatrix{{1.0047255852,0.0015660757,0.0002083646,60.0933657016,0.033311656,0.0052610767},{0.0015664176,0.9978076847,0.000046714,0.0333150727,59.9571364923,0.0012489603}, ...
dpvff      = {P(1.0047255828976631, 0.0015664175152778625, 2.085363958030939E-4), V(1.5569024253636599E-4, 5.551503272727132E-5, 8.765430720814038E-6), A(2.4965587268610534E-6, 1.0863145574191435E-6, 2.3552176342533926E-7)}
dpvf pred  = {P(1.004725583190524, 0.0015664166804878776, 2.08536165873188E-4), V(1.5569020988084276E-4, 5.5515017480352316E-5, 8.765427476001655E-6), A(0.0, 0.0, 0.0)}
``````

I think it is closer to what you’re looking for.

Regards,
Maxime

3 Likes

This worked! Thank you so much for the tip!