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)
propagator.addForceModel(NewtonianAttraction(gravity_field.getMu()))
#propagator.addForceModel(ThirdBodyAttraction(CelestialBodyFactory.getMoon()))
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?)