Space Science with Python
Space Science with Python — Around the Sun
Part 6 of the tutorial series will have a closer look at orbital mechanics. How can we calculate the distance to potentially hazardous asteroids?
Preface
This is the 6th part of my Python tutorial series “Space Science with Python”. All codes that are shown here are uploaded on GitHub. Enjoy!
Introduction
Space can be violent. You have probably heard of Supernovae, Black Holes with an everlasting appetite, or collisions of celestial bodies on a huge scale (even a small asteroid can cause a severe impact on a biosphere like our home planet).
Space can also be calm. Like our very cosmic neighbourhood: The Solar System. We have one star, 4 inner planets, 4 outer gas giants (who partly protect us from larger asteroids), comets, asteroids and dust. Sure there are collisions and other events, but in general, it is quite and we can assume that all bodies revolve around the Sun in a so-called 2-body problem: One large mass M is being revolved by a much smaller mass m. We can assume that in all cases the following inequation applies: m<<M. The resulting problem becomes even simpler.
And thanks to this assumption we can apply functions from the SPICE library, to compute simple orbital mechanics and positions with Python. Today, we will compute and work with the so-called orbital elements. Do not worry, we will work closely with Python and have a brief theory lesson before the programming part begins!
Orbits
Our home planet revolves the Sun in 365.25 days. During the year, the distance between the Earth and the Sun varies only slightly: Between around 147 and 152 Million kilometres; a variation of only a few percent. Our planet revolves the Sun almost on a perfect circle without any extreme variations, causing a friendly climate and biosphere over millions of years.
What is this “almost perfect circle”? It is an ellipse, described by different parameters that are explained in the following. We will not handle the theoretical derivation of the underlying equations since it would exceed the goal of this tutorial. Generally speaking, the mathematical fundament lies on Newton mechanics of motion. E.g., the semi-major axis (see below) can be derived from the Vis-Viva equation [Murray and Dermott, 2000, p. 30–31]. Chapter 3.1.1. of my thesis [Albin 2019] summarises briefly the derivation of all parameters. A little bit below you can see a sketch that shows the explained parameters:
- Semi-major axis (a) This parameter describes the “longest radius” of an ellipse. The smallest radius is called semi-minor axis (b) and is mostly not used in orbit dynamics.
- Eccentricity (e) The eccentricity describes the deviation of the orbit from a perfect circle. e has the following definition limits:
e = 0 (Circle)
0 < e < 1 (Ellipse)
e = 1 (Parabola)
e > 1 (Hyperbola)
Eccentricities that are larger or equal 1 describe open orbits. This applies, for example, for objects of interstellar origin, like 1I/ʻOumuamua [Meech et al. 2017, Micheli et al. 2018].
- Inclination (i) The inclination is the enclosed angle between the orbit’s plane and a reference frame, like e.g., the ecliptic plane in the Ecliptic Coordinate System (ECLIPJ2000). i describes different orbit cases:
0°≤i<90° (Prograde orbit. Example: The Earth revolves the Sun in a prograde orbit in ECLIPJ2000)
i=90° (Polar orbit. Example: Weather satellites, Earth observation satellites or reconnaissance satellites revolve the Earth on polar orbits. This allows the satellites to scan the Earth on changing longitudes due to Earth’s rotation).
90°<i≤180° (Retrograde orbit. Example: Some meteoroid streams of cometary origin are on a retrograde orbits around the Sun. These meteors enter the Earth’s atmosphere with higher velocities than prograde ones).
- Longitude of ascending node (Ω, ☊) This angle parameter describes where the ascending orbit crosses the reference plane with respect to a reference point (In ECLIPJ2000: the Vernal Equinox direction). At this node, the object moves from the southern to the northern hemisphere of the reference frame.
- Argument of periapsis (⍵) The argument of periapsis describes the angular distance between the periapsis (as seen from the focal point) and the longitude of ascending node. The periapsis is the closest point around the larger centre mass and can be computed with (1-e)a. The counterpart, the apoapsis, can be computed with (1+e)a.
- True Anomaly (v, ν) The true anomaly is the angle between the body’s position (as seen from the focal point) w.r.t. the argument of periapsis. This parameter is not needed to describe the orbit’s shape, but to determine the actual position of the object at a certain time called …
- Epoch (t)
- Mean Anomaly (M) The concept of the mean anomaly appears to be confusing is a parameter that is used e.g., by the SPICE library (instead of the true anomaly). The angular speed (degrees per second) of the true anomaly varies with time (higher at the periapsis, lower at the apoapsis). The mean anomaly M describes a virtual position of the object drawn on an imaginary circle around the focal point of the ellipse. The idea: In contrary to the true anomaly, the mean anomaly’s angular speed is constant.
Finally, before we dive into some use cases one needs to know that different systems have different suffixes. For example the periapsis is a general term for an ellipse. In our Solar System, the suffix -apsis is replaced with:
- -helion (objects that revolve the Sun)
- -gee (objects, like our Moon that revolve the Earth)
- -jove (objects that revolve around Jupiter, like Ganymede)
- -krone or -saturnium (objects that revolve around Saturn, like Titan)
- -astron (exoplanets around their stars)
The dwarf planet Ceres
Now let’s start with some programming fun. Today, we will use the data from SPICE’ spk kernels to derive orbital elements and vice versa. We begin with the dwarf planet / asteroid (1) Ceres. Ceres’ home is the asteroid main belt and with a diameter of over 900 km it is the largest representative in this crowded region. Between 2015 and 2018 NASA’s Dawn mission explored this world (see image above).
What are the orbital elements of Ceres? And can we derive values that are comparable with the literature? Let’s see:
First, as always, we need to import all necessary Python packages like datetime, numpy and spiceypy. From now on a few additional packages (pathlib, urllib, os) are needed to create an auxiliary download function.
The last tutorials required SPICE kernel files that have been uploaded on my GitHub repository. In future tutorials we need more and more kernels, depending on the tutorial’s task. Since the GitHub storage is limited, we create a helper function to download all future kernels in the _kernels directory. We define the function download_kernel in the following part. First, our function requires 2 input variables, a local storage path (dl_path) and the URL (dl_url) of the kernel (line 5). Line 6 to 18 are the docstring of the function. The URL contains the kernel’s file name. To obtain the name line 23 splits the URL at the “/”. Consequently, the very last entry of the resulting list has the file name. Now, the local download (sub-) directory is created, if it is not existing yet (line 27). If the kernel file is not present in this sub-directory (line 30) the kernel is downloaded and stored there (33).
For our Ceres project we need to search for the required spk kernels in the NAIF SPICE repository. Some asteroids data are stored under generic_kernels/spk/asteroids. There you find a readme file called aa_summaries.txt. It lists all spk asteroid data and NAIF IDs that are stored in the spk file codes_300ast_20100725.bsp.
But where do this data come from? Are they reliable? And since there are hundred thousands of asteroids: where are the others? To answer this question let’s take a look at the other readme file AAREADME_Asteroids_SPKs.txt. There you can read:
[...]We have left here a single SPK containing orbits for a set of 300 asteroids, the orbits for which were produced by Jim Baer in 2010 using his own CODES software. Some of these data may also be not very accurate. But the convenience of having data for 300 objects in a single SPK file could be useful for some purposes where high accuracy is not important. Read the *.cmt" file for details about this collection.[...]
So the kernel provides us spk information “for some purposes where high accuracy is not important”. Let’s assume that the accuracy is good enough to compare with the literature. Further, the file states that the “*.cmt file” contains more “details about this collection”. We take a look at codes_300ast_20100725.cmt. There you can find details about the underlying code that computed the spk file and other meta information. Doing science means doing literature research and reading and one detail in this cmt file is quite important:
Usage
--------------------------------------------------------
This file should be used together with the DE405 planetary
ephemeris since the integration was done using the DE405
ephemeris and with the frames kernel defining the
``ECLIPJ2000_DE405'' frame.Hold on! Another reference frame called ECLIPJ2000_DE405? How is this frame defined? The files states that:
This frames kernel defines the Ecliptic and Equinox of the J2000 frame using the DE405 obliquity angle value (84381.412 arcseconds).
[...]The corresponding frame kernel is defined as (a detailed explanation will follow in a future tutorial):
\begindata
FRAME_ECLIPJ2000_DE405 = 1900017
FRAME_1900017_NAME = 'ECLIPJ2000_DE405'
FRAME_1900017_CLASS = 4
FRAME_1900017_CLASS_ID = 1900017
FRAME_1900017_CENTER = 0
TKFRAME_1900017_SPEC = 'ANGLES'
TKFRAME_1900017_RELATIVE = 'J2000'
TKFRAME_1900017_ANGLES = ( 0.0, 0.0, -84381.412 )
TKFRAME_1900017_AXES = ( 1, 3, 1 )
TKFRAME_1900017_UNITS = 'ARCSECONDS'It appears to be a reference frame that is based on the Equatorial J2000 system and tilts one axis by 84381.412 arcseconds respectively 23,4°. This is the tilt of Earth’s rotation axis w.r.t. the ecliptic plane! Are ECLIPJ2000 and ECLIPJ2000_DE405 the same reference frame? We will find it out in a moment. First, we download the .bsp kernel and the kernel codes_300ast_20100725.tf. This kernel contains also the reference frame, but also the NAIF ID codes with the corresponding names. Since it contains miscellaneous information we create a new sub-directory called _kernels/_misc.
The SPICE meta file contains already the relative path names to the downloaded kernels, as well as the leap seconds kernel and planetary spk kernel. We load the meta file with furnsh. Our computations will be based on today’s date and the current time (line 5, line 8 converts the UTC date-time to a the Ephemeris time using the SPICE function utc2et).
Let’s take a look at the differences between ECLIPJ2000_DE405 and ECLIPJ2000. We compute in line 4 to 6 a 6x6 transformation matrix that transforms state vectors given in ECLIPJ2000_DE405 to ECLIPJ2000 at the current date-time (SPICE function sxform). Our theory: there is no difference between both systems, so we print (line 10 to 12) the matrix and expect the identity matrix (identity matrix times a vector results in the same vector):
The output is shown below. It is the identify matrix with 1s along the diagonal and 0s elsewhere. So both reference frame names represent the same frame.
Transformation matrix between ECLIPJ2000_DE405 and ECLIPJ2000
[1. 0. 0. 0. 0. 0.]
[0. 1. 0. 0. 0. 0.]
[0. 0. 1. 0. 0. 0.]
[0. 0. 0. 1. 0. 0.]
[0. 0. 0. 0. 1. 0.]
[0. 0. 0. 0. 0. 1.]Now let’s compute the state vector of Ceres in ECLIPJ2000. According to the .tf file the NAIF ID of Ceres is 2000001. We set the Sun (10) as the observer. Since we do not apply any light time corrections we can simply use the function spkgeo.
To compute the corresponding orbital elements we need to set the gravitational parameter (G) times the mass of the Sun (M). Line 2 extracts the G*M value from the needed kernel using the function bodvcd. The function requires the NAIF ID of the Sun (10), the requested parameter (GM) and the number of expected returns (1). The returned value is stored in an array and needs to be extracted in line 4.
SPICE provides two functions to determine the orbital elements based on a state vector:
oscltx requires the same input parameters as oscelt and returns the same values + 2 additional parameters. The input is:
- state The state vector of the object (the spatial components x, y, z in km and the corresponding velocity components vx, vy, vz given in km/s)
- et The Ephemeris time
- mu The gravitational parameter (here: G*M of the Sun)
The output is a single array that contains the following parameters:
- Periapsis in km
- Eccentricity
- Longitude of ascending node in radians
- Argument of periapsis in radians
- Mean anomaly at the epoch
- Epoch
- Gravitational parameter (same value as the input)
And the oscltx add-ons:
- Semi-major axis in km
- Orbital period in seconds
Line 2–4 compute the elements and the lines 7 to 24 assign the elements to individual constant names, where km are converted to AU (with the function convrt) and radians are converted to degrees (using the numpy function numpy.degrees). The very last line converts the orbital period given in seconds to years.
What is a reliable source or reference to compare the results with? It’s the (International Astronomical Union) Minor Planet Center (short: MPC). The MPC gathers all information regarding smaller objects of our Solar System (main belt asteroids, comets, Kuiper belt objects and so on). One section provides some general information for dwarf planets: MPC Dwarf Planets. There we find some data for the orbital elements and compare them with our results in the following coding part.
The results are:
Ceres' Orbital Elements
Semi-major axis in AU: 2.77 (MPC: 2.77)
Perihelion in AU: 2.55 (MPC: 2.56)
Eccentricity: 0.08 (MPC: 0.08)
Inclination in degrees: 10.6 (MPC: 10.6)
Long. of. asc. node in degrees: 80.3 (MPC: 80.3)
Argument of perih. in degrees: 73.7 (MPC: 73.6)
Orbit period in years: 4.61 (MPC: 4.61)All in all, the orbital elements from the spk kernel file appear to be feasible. Only minor differences can be seen like a 0.01 AU difference for the semi-major axis.
Can we convert the orbital elements now back to a state vector? Sure! For this purpose SPICE provides a function called conics. The function requires two inputs: an array with the orbital elements and the requested ET. The order of the orbital elements are …
- Periapsis in km
- Eccentricity
- Longitude of ascending node in radians
- Argument of periapsis in radians
- Mean anomaly at the epoch
- Epoch
- Gravitational parameter
… so basically the output order of the function oscelt. Line 2 to 9 convert the orbital elements to a state vector and line 11 to 14 print the state vectors from conics and spkgeo.
The output is shown below and you can see that the state vectors are similar!
State vector of Ceres from the kernel:
[ 3.07867363e+08 -3.13068701e+08 -6.66012287e+07 1.19236526e+01
1.14721598e+01 -1.83537545e+00]
State vector of Ceres based on the determined orbital elements:
[ 3.07867363e+08 -3.13068701e+08 -6.66012287e+07 1.19236526e+01
1.14721598e+01 -1.83537545e+00]Asteroid 1997BQ
Have you seen a meteor? A falling star? Probably. One can observe so-called sporadics or stream meteors (like the Perseids in August). These are small dust particles entering the Earth’s atmosphere with high velocities. The atmospheric friction and ionisation cause a nice glowing effect. The source of these particles: comets, asteroids, collisions of minor bodies, etc. But there are also bigger particles that cross our Earth’s path around the Sun (or at least they come quite close): Near-Earth Objects (NEOs).
On Spaceweather.com you can see a list (bottom of the webpage) where objects are listed that will have a close encounter with Earth. Today (beginning of May 2020) we can see several objects. So let’s pick a larger one, like:
Asteroid | Date(UT) | Miss Distance | Diameter (m)
---------|-------------|---------------|------------- 136795 | 2020-May-21 | 16.1 LD | 892That’s a huge asteroid! Almost 1 km in diameter. Such an impact would be devastating for our home planet. But luckily, the minimum distance between Earth and 136795, respectively 1997 BQ is around 16.1 Lunar Distances (LD; 1 LD = 384401 km).
Before we compute the orbital elements of this object let’s take a quick look whether our planet’s gravitational pull would alter the orbit of the object severely.
One could now develop an N-body numerical simulator that considers all perturbation effects. Or, we could do a first simple approach before we apply any mathematical “overkill”. For this approach we use a simple astrodynamical idea: The Sphere of Influence (SOI). Each object in our Solar System has a SOI: The Earth, the Moon, each asteroid and planet. Within each body’s SOI, the body itself is the major gravitational factor. This means that for simple computations e.g., within Earth’s SOI, we can neglect the gravitational pull of larger objects like the Sun. The SOI is theoretically an ellipsoid, but can be approximated by a sphere. The SOI radius r of a mass m (e.g., the Earth) depends on its semi-major axis a and the larger mass M (e.g., the Sun):
So first, we can compute the Earth’s SOI to determine whether our planet will alter the orbit or not. We set 1 AU as the semi-major axis. Line 10 converts 1 AU to km (SPICE function convrt). In line 13 and 14 we extract the G*M value of Earth using bodvcd. Finally, line 17 computes Earth’s SOI in km. We set 1 LD as 384401 km in line 20 and print the result in line 22.
Earth’s SOI is around 2.4 LD. Consequently, the orbit of 1997 BQ should not be altered to much by the encounter with our home planet.
SOI of the Earth in LD: 2.4054224328225597JPL’s Small Body Database provides the orbital elements at a certain epoch. However, we need to consider the provided errors (1 sigma deviation) of the values. In science, one has to round values, results or measurements based on the first 2 significant digits of the measurement error. E.g.,: x = 123.98765 with an error of xerr = 0.035 becomes x = 123.987 +/- 0.035.
The JPL database provides several digits for the error. So we build a small lambda function that rounds the values based on the position of the first 2 error digits. First, the log10 value of the error is taken. Example: 0.035 becomes -1.4559. Then the floor function is applied: -1.4559 becomes -2. We multiply the value with -1, and add +1 to the end. The result is 3. And that is what we need: we want to round 123.98765 at the 3rd position based on the error 0.035!
Let’s define all orbital elements with the data from the database. Please note: The database provides the epoch 2459000.5 (Julian Date). Thanks to the mighty utc2et function, this can easily be converted to ET using the string format “2459000.5 JD”.
Now we can compute the state vector of the asteroid after we create an array with orbital elements (line 2 to 9). We use the function conics to determine the state vector for the current date-time (line 12). Finally, the results are printed in line 14+15.
Current state vector of 1997BQ in km and km/s (2020-05-08T15:29:48):
[-1.01302335e+08 -9.80403762e+07 2.92324589e+06 2.10287249e+01
-2.97600309e+01 -6.83829911e+00]What is the current distance between 1997 BQ and the Earth? To find out, we compute Earth’s state vector, too (using spkgeo). The NAIF ID for our planet is 399 and the observer is the Sun (10). We use the same ET (DATETIME_ET) and set the reference frame ECLIPJ2000. The spatial components of both state vectors are then used to determine the length of the resulting vector between Earth and 1997 BQ (line 8+9, function vnorm).
For our timestamp the distance is around 38.69 LD.
Current distance between the Earth and 1997BQ in LD (2020-05-08T15:29:48):
38.69205689796575Conclusion
That was a lot of content to digest! Take your time and use the provided Jupyter Notebook and Python script from my GitHub repository to walk through this tutorial. Reference frames and orbital mechanics as shown in the last two tutorials are used in the future tutorials. Do not worry. The next tutorials will be smaller to give you some time to get familiar with everything. There are a lot of Python functions, concepts and technical terms that need to be understand. If you need support or have any questions do not hesitate to ask.
Thomas
References
- Albin, T. (2019). Machine learning and monte carlo based data analysis methods in cosmic dust research. University of Stuttgart.
- Meech, K., Weryk, R., Micheli, M. et al. A brief visit from a red and extremely elongated interstellar asteroid. Nature 552, 378–381 (2017). https://doi.org/10.1038/nature25020
- Micheli, M., Farnocchia, D., Meech, K.J. et al. Non-gravitational acceleration in the trajectory of 1I/2017 U1 (‘Oumuamua). Nature 559, 223–226 (2018). https://doi.org/10.1038/s41586-018-0254-4
- Carl D. Murray and Stanley F. Dermott. Solar System Dynamics. Cambridge University Press, 2000. ISBN 978–0–521–57597–3.
