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. We assume non dimensional units for the state.
class cr3bp_udpla:
def __init__(self, ref_mjd2000, state_nd, mu, T_units, name = "station", tol = 1e-12):
import heyoka as hy
self.ref_mjd2000 = ref_mjd2000
self.ref_state = state_nd
self.mu = mu
self.T_units = T_units
self.name = name
self.ta = pk.ta.get_cr3bp(tol)
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):
dt = (mjd2000 - self.ref_mjd2000) * pk.DAY2SEC / self.T_units
self.ta.state[:] = self.ref_state
self.ta.time = self.ref_mjd2000 * pk.DAY2SEC / self.T_units
self.ta.propagate_for(dt)
return [self.ta.state[:3], self.ta.state[3:]] # these will be casted to list by the pykep planet iface.
def acc(self, mjd2000):
epoch = mjd2000 * pk.DAY2SEC / self.T_units
self.ta.propagate_until(epoch) # We use propagate until so that if eph was called just before we avoid propagating.
state = self.ta.state
return self.acc_cfunc(state, pars =self.ta.pars)
def get_name(self):
return "A body under nd CR3BP dynamics"
def get_extra_info(self):
retval = "Reference state (nd): " + str(self.ref_state)
retval += "\nReference epoch (mjd2000): " + str(self.ref_mjd2000)
return retval
L_Units = 384400000
T_Units = 4.35 * pk.DAY2SEC
state = [ 0.99263187, -0.15446297, 0. , -0.74481059, -0.53599621, 0.]
station = pk.planet(cr3bp_udpla(0., state, pk.CR3BP_MU_EARTH_MOON, T_Units, name = "station"))
r,v = station.eph(epoch)
print(station)
print(station.extra_info)
print("Ephemerides at epoch: ", epoch)
print("r:", r)
print("v:", v)
A body under nd CR3BP dynamics
Reference state (nd): [0.99263187, -0.15446297, 0.0, -0.74481059, -0.53599621, 0.0]
Reference epoch (mjd2000): 0.0
Ephemerides at epoch: 2026-10-22T00:00:00.000000
r: [0.07233978305614301, -0.8336255008104227, 0.0]
v: [-0.47833877065155095, 0.8387520292044038, 0.0]
And thats it! Hope you enjoyed.