Modified equinoctial elements in pykep

Modified equinoctial elements in pykep#

In pykep the modified equinoctial elements (MEE) are represented as:

\([p, f, g, h, k, L]\)

where \(p\) is the semi-latus rectum, \(f\) and \(g\) encode the eccentricity vector, \(h\) and \(k\) encode the orbital plane orientation and \(L\) is the true longitude. These elements are often preferable to Keplerian elements because they avoid some of the familiar singularities at low eccentricity and low inclination.

This notebook shows how to:

  • convert a Cartesian state to MEE with pykep.ic2mee(),

  • convert MEE back to Cartesian coordinates with pykep.mee2ic(),

  • use the retrograde convention,

  • access the symbolic heyoka expressions and Jacobians exposed by pykep.

import heyoka as hy
import numpy as np
import pykep as pk

From Cartesian coordinates to MEE#

Let us start from a Cartesian state around the Earth. As usual in pykep, units only need to be consistent. Here we work in SI units.

mu = pk.MU_EARTH
r = [7000e3, 0.0, 1200e3]
v = [0.0, 7.4e3, 1.0e3]

mee = pk.ic2mee([r, v], mu)
print(mee)

[7052411.762732841, -0.0067439609364834, -0.02255221909857627, 0.06678130484811101, -0.08471685529303227, 0.011283863497761658]

The returned vector is ordered as \([p, f, g, h, k, L]\). We can inspect it in a slightly more readable form.

labels = ["p", "f", "g", "h", "k", "L"]
for label, value in zip(labels, mee):
    print(f"{label} = {value}")

p = 7052411.762732841
f = -0.0067439609364834
g = -0.02255221909857627
h = 0.06678130484811101
k = -0.08471685529303227
L = 0.011283863497761658

Back to Cartesian coordinates#

The inverse transformation is available via pykep.mee2ic().

r_back, v_back = pk.mee2ic(mee, mu)
print("r_back =", r_back)
print("v_back =", v_back)

r_back = [7000000.000000001, 0.0, 1200000.0000000002]
v_back = [5.684341886080802e-14, 7399.999999999999, 999.9999999999997]

print("Position error (rel.):", np.linalg.norm(np.array(r_back) - np.array(r))/np.linalg.norm(np.array(r)))
print("Velocity error (rel.):", np.linalg.norm(np.array(v_back) - np.array(v))/np.linalg.norm(np.array(v)))

Position error (rel.): 1.35168987445e-16
Velocity error (rel.): 1.30302461055e-16

Retrograde convention#

pykep also supports the retrograde equinoctial convention via the retrogade flag. This is useful when one wants the non-singular formulation around inclinations close to \(i\). The same convention must be used consistently in both directions.

mee_retro = pk.ic2mee([r, v], mu, retrogade=True)
r_retro, v_retro = pk.mee2ic(mee_retro, mu, retrogade=True)

print("retrograde MEE =", mee_retro)
print("retrograde position error =", np.linalg.norm(np.array(r_retro) - np.array(r))/np.linalg.norm(np.array(r)))
print("retrograde velocity error =", np.linalg.norm(np.array(v_retro) - np.array(v))/np.linalg.norm(np.array(v)))

retrograde MEE = [7052411.762732841, 0.02350352264542442, -0.001291517959233943, 5.7388582924215195, -7.2801516623861575, 1.8177590265285917]
retrograde position error = 1.35526315396e-15
retrograde velocity error = 1.06355196642e-15

Symbolic expressions and Jacobians#

For use with heyoka, pykep exposes symbolic expressions for both transformations. The symbolic variables are:

  • ic2mee: \([x, y, z, v_x, v_y, v_z]\),

  • mee2ic: \([p, f, g, h, k, L]\).

The symbolic parameters are par[0] = mu and par[1] = I, where I = 1 for the standard convention and I = -1 for the retrograde one.

ic2mee_hy, J_ic2mee_hy = pk.ic2mee(jacobian=True)
mee2ic_hy, J_mee2ic_hy = pk.mee2ic(jacobian=True)

x, y, z, vx, vy, vz = hy.make_vars("x", "y", "z", "vx", "vy", "vz")
p, f, g, h, k, L = hy.make_vars("p", "f", "g", "h", "k", "L")

ic2mee_cf = hy.cfunc(ic2mee_hy, vars=[x, y, z, vx, vy, vz])
mee2ic_cf = hy.cfunc(mee2ic_hy, vars=[p, f, g, h, k, L])
J_ic2mee_cf = hy.cfunc(J_ic2mee_hy, vars=[x, y, z, vx, vy, vz])
J_mee2ic_cf = hy.cfunc(J_mee2ic_hy, vars=[p, f, g, h, k, L])

Let us verify two things:

  1. the compiled symbolic ic2mee gives the same result as the direct pykep.ic2mee,

  2. the Jacobians compose to the identity, as they should for inverse maps evaluated at matching points.

ic = [r[0], r[1], r[2], v[0], v[1], v[2]]
pars = [mu, 1.0]

mee_from_direct = np.array(pk.ic2mee([r, v], mu))
mee_from_symbolic = np.array(ic2mee_cf(ic, pars=pars))

print("max |direct - symbolic| =", np.max(np.abs(mee_from_direct - mee_from_symbolic)))

max |direct - symbolic| = 1.30104260698e-16

J1 = np.array(J_ic2mee_cf(ic, pars=pars)).reshape(6, 6)
J2 = np.array(J_mee2ic_cf(mee_from_symbolic, pars=pars)).reshape(6, 6)

print("J_ic2mee @ J_mee2ic =")
print(J1 @ J2)

J_ic2mee @ J_mee2ic =
[[  1.00000000e+00  -1.67410998e-09  -4.59340713e-11  -1.06989511e-09
    1.76674109e-09  -3.95778854e-11]
 [  6.63047005e-23   1.00000000e+00  -1.17738870e-17  -1.01275297e-16
    6.64003736e-17   1.73407323e-18]
 [ -4.35735939e-25  -5.75596865e-18   1.00000000e+00  -1.58440481e-17
   -2.57639672e-17  -2.29633630e-16]
 [  1.15776193e-24  -7.82308489e-19  -3.49012935e-19   1.00000000e+00
    1.24092691e-18   2.85987105e-18]
 [ -7.32111826e-25  -2.27350124e-18   9.89793830e-20  -8.48792944e-18
    1.00000000e+00   7.90294932e-18]
 [  3.75019563e-25  -2.42165385e-18  -9.87548881e-19   5.07060266e-17
   -3.56109756e-17   1.00000000e+00]]

This symbolic interface is useful when one needs differentiable coordinate transforms inside a larger heyoka workflow, for example in sensitivity analysis, indirect methods or custom trajectory optimization pipelines.