Planets and User Defined Planets (UDPLAs)#

The ephemerides of a body are the position end velocity Cartesian vectors of a moving object. In pykep the class planet offers a common interface to access various ephemerides, regardless on how they are computed. Whether the underlying computations are simply Keplerian, or based on more advanced trajectory propagations, interpolations or predictions, a unified interface is offered by this type-erasing class.

The user can code his own python class following the mandatory interface of a planet and his own coded objects will be treated uniformly as any other planet in pykep. These User Defined Planets (or UPLAs) can implement heterogeneous techniques and interface to any third-party API, but once they are used to construct a planet, they will appear as any other planet in pykep.

For convenience, a number of already coded UDPLAs are offered so that, to start with, users can compute the positions of planet, satellites comets and spacecraft without the need to code their own UDPLA.

import numpy as np
import pykep as pk

UDPLAs in pykep#

Let us start to use the UDPLA keplerian which describes the motion of an object in a Keplerian orbit. All UDPLAs provided by pykep are in the module udpla

SPICE kernels can be used to create the UDPLA spice, without needing to download kernels, the DE440 ephemerides are conveniently provided in de440s. For TLEs the UDPLA tle can be used to build a planet.

# We instantiate the udpla
udpla = pk.udpla.keplerian(pk.epoch(0.), [1., 0, 0, 0, 0, 0], mu_central_body = 1., name = "A Circular Orbit")
# And use it to construct a pykep planet, hence erasing its type.
pla = pk.planet(udpla)

Note

Non-dimensional units are often useful in astrodynamics and most pykep class allow their use. Of course only if everything is consistent the computations will make sense. pykep does not check this consistency and thus the user must be careful when using anything which is non SI, and to not mix units.

Let us print on screen this planet.

print("Short repr: \n", pla)
print("Long repr: \n", pla.extra_info)

Short repr: 
 A Circular Orbit
Long repr: 
 Keplerian planet elements: 
Semi major axis: 1
Semi major axis (AU): 6.684587122268445e-12
Eccentricity: 0
Inclination (deg.): 0
Big Omega (deg.): 0
Small omega (deg.): 0
True anomly (deg.): 0
Mean anomly (deg.): 0
Elements reference epoch (MJD2000): 0
Elements reference epoch (UTC): 2000-01-01T00:00:00.000000
r at ref. = [1, 0, 0]
v at ref. = [0, 1, 0]

The textual representation of a planet can be extracted in two parts.

The first part is common to all objects of type planet and reports the body name only.

The second part is instead original with the udpla type, and is essentially whatever is returned by the planet interface optional method extra_info. In this case, the udpla is a c++ based class exposed to python as :class:pykep.udpla.keplerian and its extra info report the orbital parameters as well as the Cartesian components of the body state at the reference epoch.

Let see where the body is:

epoch = pk.epoch("2026-10-22 00:00:00:00")
pla.eph(epoch)

[[-0.9810645314764517, -0.19368114280665175, -0.0],
 [0.19368114280665175, -0.9810645314764517, 0.0]]

Custom UDPLAs: a tutorial#

Lets see how users can define their own custom udplas and thus leverage the common planet interface. A minimal version includes the method eph only. Here we will code a very simple and strange object moving in a rectilinear constant motion along the x-axis. Its a ficticious object and it only serves the purpose of this tutorial.

class my_udpla:
    # This is the only mandatory method in the planet interface
    def eph(self, mjd2000):
        return [mjd2000*np.array([1., 0., 0.]), np.array([1., 0., 0.])]

# Lets instantiate a planet
pla = pk.planet(my_udpla())

# And test
print("call to *eph* works, calls user implementation:", pla.eph(2.))
print("call to *eph_v* also works, calls default implementation:\n", pla.eph_v([1., 2.]))
print("call to *mu_central_body* returns a default -1 value signalling its not implemented:", pla.mu_central_body)
print("call to *name* returns a default value with little info:", pla.name)
# etc....

call to *eph* works, calls user implementation: [[2.0, 0.0, 0.0], [1.0, 0.0, 0.0]]
call to *eph_v* also works, calls default implementation:
 [[ 1.  0.  0.  1.  0.  0.]
 [ 2.  0.  0.  1.  0.  0.]]
call to *mu_central_body* returns a default -1 value signalling its not implemented: -1.0
call to *name* returns a default value with little info: <class '__main__.my_udpla'>

We can then add a bunch of optional method unlocking various functionalities. A rather complete final udpla could look like this:

class my_udpla:
    # This is the only mandatory method in the planet interface
    def eph(self, mjd2000):
        return [mjd2000*np.array([1., 0., 0.]), np.array([1., 0., 0.])]

    # Often used to compute gradients in some optimal control, the acc method (optional)
    # returns a (acc), while eph are r,v (pos, vel). Its default implementation
    # will return -mu/r^3 if mu is present, or throw otherwise.
    def acc(self, mjd2000):
        return np.array([4., 5., 4.])*mjd2000
    
    # Vectorized versions for eph and acc have default naive implementations, which can be (optionally) overidden.
    def acc_v(self, mjd2000s):
        # We implement a naive vectorization
        retval = np.zeros((len(mjd2000s), 3))
        for i, mjd2000 in enumerate(mjd2000s):
            acc = self.acc(mjd2000)
            retval[i,:3] = acc
        return retval
    def eph_v(self, mjd2000s):
        # We implement a naive vectorization
        retval = np.zeros((len(mjd2000s), 6))
        for i, mjd2000 in enumerate(mjd2000s):
            r, v = self.eph(mjd2000)
            retval[i,:3] = r
            retval[i,3:] = v
        return retval
    
    # A bunch of methods allow to add information to
    # the udpla (and exposed in the planet interface) which can be used in some contexts and 
    # also affect the availability of default implementations of, for example, acc and elements.
    def get_mu_central_body(self):
        return 3.
    
    def get_mu_self(self):
        return 1.1
    
    def get_radius(self):
        return 0.345
    
    def get_safe_radius(self):
        return 0.124
    
    def get_name(self):
        return self.whoami
    
    def get_extra_info(self):
        return "Extra info"
    
    # We can also have a constructor to have some internal state.
    def __init__(self, whoami):
        self.whoami = whoami
    
pla2 = pk.planet(my_udpla("I have all methods! I feel complete."))

print("call to *eph* works, calls user implementation:", pla2.eph(2.))
print("call to *eph_v* calls user implementation:\n", pla2.eph_v([1., 2.]))
print("call to *mu_central_body* calls user implementation:", pla2.mu_central_body)
print("call to *name* calls user implementation:", pla2.name)
print("call to *period* makes use of the user provided mu:", pla2.period(2.))
print("call to *elements* makes use of the user provided mu:", pla2.elements(2., pk.el_type.MEE))

call to *eph* works, calls user implementation: [[2.0, 0.0, 0.0], [1.0, 0.0, 0.0]]
call to *eph_v* calls user implementation:
 [[ 1.  0.  0.  1.  0.  0.]
 [ 2.  0.  0.  1.  0.  0.]]
call to *mu_central_body* calls user implementation: 3.0
call to *name* calls user implementation: I have all methods! I feel complete.
call to *period* makes use of the user provided mu: 6.664324407237548
call to *elements* makes use of the user provided mu: [0.0, nan, nan, nan, nan, nan]

Note the nans deriving from the impossibility to define osculating equinoctial elemnts for this weird body moving along a line passing throught eh origin … all good. Equinoctial elements are non-singular …. up to a certain extent!

Custom UDPLAs: an example in the CR3BP#

Lets say we want to define a planet in the CR3BP. For example we may have an object station moving along a periodic orbit around L1. We need to define an udpla able to compute \(r,v\) and possibly \(a\) at some given epoch.

class cr3bp_udpla:
    def __init__(self, when, state_nd, mu_cr3bp, TIME, L, name=None, tol=1e-16):
        if name is None:
            name = "An unkown body in the CR3BP with mu = " + str(mu_cr3bp)
        import heyoka as hy

        # When can be a mjd2000 or a pykep epoch, consistently with other pykep planet interfaces.
        if isinstance(when, pk.epoch):
            self.ref_mjd2000 = when.mjd2000
        else:
            self.ref_mjd2000 = when
        self.ref_state = state_nd
        self.mu_cr3bp = mu_cr3bp
        # These are the nd units declared by the user.
        self.TIME = TIME
        self.L = L
        self.V = self.L / self.TIME
        self.ACC = self.V / self.TIME
        self.name = name
        # We use the nd numerical integrator of the cr3bp given by pykep.
        self.ta = pk.ta.get_cr3bp(tol)
        # Set the parameters and the initial state of the integrator.
        self.ta.pars[0] = self.mu_cr3bp
        self.ta.state[:] = self.ref_state
        self.ta.time = self.ref_mjd2000 * pk.DAY2SEC / self.TIME
        # The cfunc for the acceleration to implement the acc method of the pykep planet interface.
        dyn = pk.ta.cr3bp_dyn()
        x, y, z, vx, vy, vz = hy.make_vars("x", "y", "z", "vx", "vy", "vz")
        self.acc_cfunc = hy.cfunc([e[1] for e in dyn[3:]], vars=[x, y, z, vx, vy, vz])

    def eph(self, mjd2000):
        epoch_nd = mjd2000 * pk.DAY2SEC / self.TIME
        self.ta.propagate_until(epoch_nd)
        return [self.ta.state[:3] * self.L, self.ta.state[3:] * self.V]

    def acc(self, mjd2000):
        epoch_nd = mjd2000 * pk.DAY2SEC / self.TIME
        self.ta.propagate_until(epoch_nd)  # avoid double propagation when possible
        return self.acc_cfunc(self.ta.state, pars=self.ta.pars) * self.ACC

    def get_name(self):
        return self.name

    def get_extra_info(self):
        retval = "Reference state (nd): " + str(self.ref_state)
        retval += "\nReference epoch (mjd2000): " + str(self.ref_mjd2000)
        return retval

# We instantiate it as a UDPla
udpla = cr3bp_udpla(
    0.0,
    [1.0809931218390707, 0.0, -0.20235953267405354, 0.0, -0.19895001215078018, 0.0],
    pk.CR3BP_MU_EARTH_MOON,
    L=384400000,
    TIME=375000,
    tol=1e-16,
)
# And cancel its type to a pk.planet, which is the only thing that pykep functions will see.
station = pk.planet(udpla)

# Let's call the ephemerides at the reference epoch, which should be the same as the reference state, but in dimensional units.
r,v = station.eph(epoch)
print(station)
print(station.extra_info)
print("Ephemerides at epoch: ", epoch)
print("r:", r)
print("v:", v)

An unkown body in the CR3BP with mu = 0.01215058439470971
Reference state (nd): [1.0809931218390707, 0.0, -0.20235953267405354, 0.0, -0.19895001215078018, 0.0]
Reference epoch (mjd2000): 0.0
Ephemerides at epoch:  2026-10-22T00:00:00.000000
r: [391975414.64036447, -40215343.553091094, -26478840.474551987]
v: [-119.30550871978843, 49.025077996573835, 337.59994073443704]

# The above udpla is actually provided for your convenience already in pykep as pk.udpla.cr3bp
station = pk.udpla.cr3bp(
    0.0,
    [1.0809931218390707, 0.0, -0.20235953267405354, 0.0, -0.19895001215078018, 0.0],
    pk.CR3BP_MU_EARTH_MOON,
    L=384400000,
    TIME=375000,
    tol=1e-16,
)

And thats it! Hope you enjoyed.

Planets and User Defined Planets (UDPLAs)

Contents

Planets and User Defined Planets (UDPLAs)#

UDPLAs in pykep#

Custom UDPLAs: a tutorial#

Custom UDPLAs: an example in the CR3BP#