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()]) )
        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))

    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))
    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 =
    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


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.



This worked! Thank you so much for the tip!