Demonstration of coordinate transformations between orbital elements and Cartesian coordinates
The goal here is to demonstrate how these coordinate transformations work, and what they do. The goal is to from a set of orbital elements in the usual Keplerian basis \((a, e, i, \Omega, \omega, \mathcal{M})\) or the cometary basis \((q, e, i, \Omega, \omega, T_p)\) into a set of state vectors \((x, y, z, v_x, v_y, v_z)\).
Under the hood, the transformations are using the universal variable formulation of the two body problem, so that eccentric (\(0 \leq e <1\)), parabolic (\(e=1\)) and hyperbolic (\(e>1\)) orbits can be handled correclty and simultaneously.
For the reader interested into getting in the details, Danby’s Fundamental of Celestial Mechanics is a good reference. The implementation inside sorcha follows Everhart & Pitkin (1983).
A simple example
Let’s start with the simplest possible example. Let’s define the reduced mass of the system to \(\mu = 1\), and the reference time for the elements (i.e. the epoch where they are defined) to \(t = 0\), and unitless distances and times.
Let’s start with a simple series of orbits:
\(q = 10\), \(e = 0\), \(i = 0^\circ\), \(\Omega = 10^\circ\), \(\omega = 0^\circ\), \(T_p = 0\)
\(q = 10\), \(e = 0.1\), \(i = 0^\circ\), \(\Omega = 10^\circ\), \(\omega = 0^\circ\), \(T_p = 0\)
\(q = 10\), \(e = 0.9999\), \(i = 0^\circ\), \(\Omega = 10^\circ\), \(\omega = 0^\circ\), \(T_p = 0\)
\(q = 10\), \(e = 1\), \(i = 0^\circ\), \(\Omega = 10^\circ\), \(\omega = 0^\circ\), \(T_p = 0\)
\(q = 10\), \(e = 1.0001\), \(i = 0^\circ\), \(\Omega = 10^\circ\), \(\omega = 0^\circ\), \(T_p = 0\)
\(q = 10\), \(e = 6.\), \(i = 0^\circ\), \(\Omega = 10^\circ\), \(\omega = 0^\circ\), \(T_p = 0\)
(In all of these cases, the object is at perihelion at the reference time)
[1]:
import numpy as np
from sorcha.ephemeris.orbit_conversion_utilities import universal_cartesian, universal_cometary
[2]:
# define orbits (no e)
q = 10
i = 0 * np.pi/180
Omega = 10 * np.pi/180
omega = 0
Tp = 0
#define constants
epochMJD_TDB = 0
mu = 1
[3]:
for e in [0, 0.1, 0.9999, 1.0, 1.0001, 6.]:
x,y,z,vx,vy,vz = universal_cartesian(mu, q, e, i, Omega, omega, Tp, epochMJD_TDB)
print(f'{e}: {x}, {y}, {z}, {vx}, {vy}, {vz}')
0: 9.84807753012208, 1.7364817766693033, 0.0, -0.05491237529650835, 0.31142355569111246, 0.0
0.1: 9.84807753012208, 1.7364817766693033, 0.0, -0.05759258508501801, 0.32662378073744874, 0.0
0.9999: 9.84807753012208, 1.7364817766693033, 0.0, -0.07765588441652757, 0.4404084054777871, 0.0
1.0: 9.84807753012208, 1.7364817766693033, 0.0, -0.0776578258864434, 0.4404194161008241, 0.0
1.0001: 9.84807753012208, 1.7364817766693033, 0.0, -0.0776597673078231, 0.440430426448599, 0.0
6.0: 9.84807753012208, 1.7364817766693033, 0.0, -0.1452844889344078, 0.8239492807661573, 0.0
Unsurprisingly (in other words, intentionally), all orbits are at the same \((x,y,z)\) position, just different velocities (note \(v_z = 0\) in all cases because \(i = 0^\circ\)).
We can also invert these transformations:
[4]:
for e in [0, 0.1, 0.9999, 1.0, 1.0001, 6.]:
x,y,z,vx,vy,vz = universal_cartesian(mu, q, e, i, Omega, omega, Tp, epochMJD_TDB)
q_out,e_out,i_out,Omega_out,omega_out,Tp_out = universal_cometary(mu, x, y, z, vx, vy, vz, epochMJD_TDB)
print(f'{e}: {q_out}, {e_out}, {i_out}, {Omega_out}, {omega_out}, {Tp_out}')
0: 9.999999999999998, 1.5969451812864107e-18, 0.0, 0.0, 0.17453292519943295, 0.0
0.1: 10.0, 0.10000000000000009, 0.0, 0.0, 0.17453292519943273, -6.305290618946865e-15
0.9999: 9.999999999999998, 0.9998999999999998, 0.0, 0.0, 0.17453292519943286, -1.6393726832476508e-15
1.0: 9.999999999999998, 1.0, 0.0, 0.0, 0.1745329251994329, -9.946086886476325e-16
1.0001: 9.999999999999998, 1.0000999999999993, 0.0, 0.0, 0.17453292519943292, -7.98005785532596e-16
6.0: 9.999999999999998, 5.999999999999998, 0.0, 0.0, 0.17453292519943292, -2.182462553244727e-16
Note that any differences are of order \(10^{-14}\) or less - this only slightly above the machine precision level for doubles - ensuring higher precision would sacrifice runtimes and convergence for the solutions. For the all applications of sorcha, this is more than enough.
Real world example #1
Let’s also go through these transformations for a real asteroid. Let’s take asteroid 3666 Holman. Using JPL Horizons, we have that:
2457545.5000000000 = A.D. 2016-Jun-06 00:00:00.0000 TDB
EC= .1273098034941495 QR= 2.719440725689577 TP= 2457934.5526586706
OM= 120.3869311657135 W= 55.06308036878693 IN= 2.363582123711951
At the same epoch,
X =-7.569545429706993E-02 Y = 3.024083648650882E+00 Z =-6.044399403284755E-02
VX=-9.914117209213893E-03 VY=-1.485136186100886E-03 VZ= 3.840061650310168E-04
in these units (distances are in au and times are in days), \(\mu = 2.9591220828559115 \cdot 10^{-4} \, \mathrm{au}^3/\mathrm{d}^2\).
Let’s copy these numbers directly:
[5]:
EC= .1273098034941495
QR= 2.719440725689577
TP= 2457934.5526586706
OM= 120.3869311657135 * np.pi/180
W= 55.06308036878693 * np.pi/180
IN= 2.363582123711951 * np.pi/180
epochMJD_TDB = 2457545.5
[6]:
mu = 2.9591220828559115e-04
x, y, z, vx, vy, vz = universal_cartesian(mu, QR, EC, IN, OM, W, TP, epochMJD_TDB)
print(x,y,z,vx,vy,vz)
-0.0756954542915965 3.0240836486516987 -0.06044399403305996 -0.009914117209238994 -0.0014851361860867482 0.0003840061650316154
[7]:
X =-7.569545429706993E-02
Y = 3.024083648650882E+00
Z =-6.044399403284755E-02
VX=-9.914117209213893E-03
VY=-1.485136186100886E-03
VZ= 3.840061650310168E-04
print(x- X, y- Y, z - Z, vx - VX, vy - VY, vz - VZ)
5.473427266977637e-12 8.166800569142652e-13 -2.1240648129250417e-13 -2.51014486973844e-14 1.413777948877648e-14 5.98588019429247e-16
The differences are within reasonable tolerance limits.
So finally, if we want to rotate these from the ecliptic coordinate frame to the equatorial frame, we have that:
[8]:
from sorcha.ephemeris.simulation_constants import ECL_TO_EQ_ROTATION_MATRIX
[9]:
pos_ecl = np.array([x,y,z])
vel_ecl = np.array([vx,vy,vz])
pos_eq = ECL_TO_EQ_ROTATION_MATRIX @ pos_ecl
vel_eq = ECL_TO_EQ_ROTATION_MATRIX @ vel_ecl
print(pos_eq, vel_eq)
[-0.07569545 2.75049926 -1.25836767] [-0.00991412 -0.00120984 0.00094307]
Real world example #2
As a second test, let’s compare against 2I/Borisov, doing the inverse transformations:
From JPL, again:
2460188.500000000 = A.D. 2023-Sep-01 00:00:00.0000 TDB
EC= 3.345550605771202E+00 QR= 1.998821321255714E+00 IN= 4.412305585291305E+01
OM= 3.080238378883400E+02 W = 2.091328388499988E+02 Tp= 2458825.978725019377
X =-6.658710421827730E-01 Y =-2.313385391700901E+01 Z =-1.432926077530441E+01
VX= 1.089064793811738E-03 VY=-1.681872446865637E-02 VZ=-9.215726774314793E-03
[10]:
X =-6.658710421827730E-01
Y =-2.313385391700901E+01
Z =-1.432926077530441E+01
VX= 1.089064793811738E-03
VY=-1.681872446865637E-02
VZ=-9.215726774314793E-03
EC= 3.345550605771202E+00
QR= 1.998821321255714E+00
IN= 4.412305585291305E+01 * np.pi/180
OM= 3.080238378883400E+02 *np.pi/180
W = 2.091328388499988E+02 * np.pi/180
Tp= 2458825.978725019377
epochMJD_TDB = 2460188.500000000
q, e, incl, longnode, argperi, tp = universal_cometary(mu, X, Y, Z, VX, VY, VZ, epochMJD_TDB)
print(q, e, incl, longnode, argperi, Tp)
print(q-QR, e-EC, incl-IN, longnode-OM, argperi-W, tp - Tp)
1.9988213212529575 3.345550605755576 0.7700926006746867 -0.9071551613987933 -2.633128695205852 2458825.9787250194
-2.7564617255393387e-12 -1.5625722937784303e-11 -8.881784197001252e-16 -6.283185307179584 -6.283185307181043 0.0
Note that the differences in \(\Omega\) and \(\omega\) are factors of \(2 \pi\) coming from choices about the ranges of these angles:
[11]:
print(longnode-OM + 2 * np.pi, argperi-W + 2 * np.pi)
2.6645352591003757e-15 -1.4566126083082054e-12
But what about the mean anomaly?
So far, we have used \((q, e, i, \Omega, \omega, T_p)\) instead of \((a, e, i, \Omega, \omega, \mathcal{M})\). For a bound orbit, the second representation is also possible (and oftentimes simpler/easier to interpret). In these cases, the usual conversions apply:
Going back to our first example:
[12]:
a = 10
e = 0.1
i = 0 * np.pi/180
Omega = 10 * np.pi/180
omega = 0
M = 0 * np.pi/180 #note that we need radians here as well!
#define constants
epochMJD_TDB = 0
mu = 1
print(universal_cartesian(mu, a * (1-e), e, i, Omega, omega, epochMJD_TDB - M * np.sqrt(a**3/mu), epochMJD_TDB))
(8.863269777109872, 1.562833599002373, 0.0, -0.06070791506856635, 0.34429169503525636, 0.0)
The orbit is still \(q = a(1-e) = 9\) units from the center at perihelion (\(\mathcal{M} = 0^\circ\)), as desired.