.. user-manual:
******************
Cysgp4 user manual
******************
.. currentmodule:: cysgp4
Introduction
============
The `~cysgp4` package can be used to calculate the orbits of satellites around
Earth, based on so-called `Two-line element strings (TLE)`_. They not only
encapsulate not only the basic `orbital elements`_ that describe Keplarian
orbits (based on six free parameters), but also include more complicated
physical processes, such as drag caused by Earth's atmosphere and other
pertubations owing to the ellipsoidal gravitational field, influence of other
solar bodies etc. One of the more famous `simplified perturbations models`_
is SGP4, which is directly related to the TLEs published by NORAD and NASA,
and many tools exist to do the calculation.
Here we provide a Python wrapper around the `C++ SGP4 Satellite Library`_ by
Daniel Warner. It provides similar functionality as the well-known `sgp4
`_ Python package (by Brandon Rhodes), which
uses `Numba `_ internally to speed-up the
calculations. In contrast to `sgp4`_, `cysgp4` can work well with arrays of
TLEs and observing times and make use of multi-core platforms (via `OpenMP
`_) to boost processing times a lot.
Using cysgp4
============
The Python package `~cysgp4` provides a relatively thin wrapper around the
`C++ SGP4 Satellite Library`_ on the one hand, but on the other hand also
contains some higher-level functions, which make it possible to bulk-process
satellite orbits with few lines of code.
.. _wrapped-approach-label:
Option 1: C++ SGP4 library wrapper
----------------------------------
The package provides the following Classes.
PyDateTime
^^^^^^^^^^
`~cysgp4.PyDateTime` specifies a datetime. Various conversions are provided
for convenience, from and to `Modified Julian Date (MJD)
`_, Python's `~datetime`
class, so-called "ticks", which are the number of micro-seconds since 1.
January 0001, 00:00, but also the TLE-epoch format.
.. note::
TLE epochs have a very strange format: first two digits are the year,
next three digits are the day from beginning of year, then the fraction
of a day is given, e.g. 20180.25 would be 2020, day 180, 6 hours (
probably UT as no timezone is mentioned). See also `Two-line element
strings (TLE)`_ or `Celestrak `_.
There are several possibilities to create a `~cysgp4.PyDateTime` object,
e.g.::
>>> from datetime import datetime
>>> from cysgp4 import PyDateTime
>>> # from Python's datetime
>>> dt = datetime(2019, 1, 1, 12, 13, 14)
>>> pydt = PyDateTime(dt)
>>> pydt
>>> # from Modified Julian Date
>>> mjd = 56458.123
>>> pydt = PyDateTime.from_mjd(mjd)
>>> pydt
>>> # from "Ticks"
>>> ticks = 58628880000000000 # aka MJD zero point
>>> pydt = PyDateTime.from_ticks(ticks)
>>> pydt
Accessing and modifying the date and time is done via the following
properties and methods::
>>> print(pydt.datetime)
1858-11-17 00:00:00
>>> pydt.datetime = datetime(2010, 5, 1, 0, 10, 20)
>>> print(pydt.mjd)
55317.00717592592
>>> pydt.mjd = 55489.01358
>>> print(pydt.ticks)
63423130773311999
>>> pydt.ticks = 63681941594000000
>>> observer_longitude = 6.67 # degrees
>>> print('GMST: {:9.6f} hours'.format(pydt.gmst()))
GMST: 4.959715 hours
>>> print('LMST: {:9.6f} hours'.format(pydt.lmst(observer_longitude)))
LMST: 5.076129 hours
One can also set the date and time manually with the `~cysgp4.PyDateTime.set`
method.
PyTle
^^^^^
This helper class holds the `Two-line element strings (TLE)`_ and provides
access to (some of) its parameters. Although the name suggests that two lines
of strings are used, in practice often three lines are defined, the first
containing a satellite name (and/or) ID. It is important to note that for
many satellites, the corresponding TLEs get outdated quickly. Especially for
low earth orbit satellites, one should query up-to-date information at least
once per day.
`~cysgp4.PyTle` is only used to parse the TLE strings and store the orbital
parameters for use by other SGP4 routines. TLEs can be obtained from many
sources, such as `Celestrak`_::
>>> import requests
>>> import cysgp4
>>> url = 'https://celestrak.org/NORAD/elements/gp.php?GROUP=science&FORMAT=tle'
.. doctest-remote-data::
>>> ctrak_science = requests.get(url).text
>>> all_lines = ctrak_science.split('\r\n')
>>> all_lines[:3] # doctest: +SKIP
['AKEBONO (EXOS-D) ',
'1 19822U 89016A 19331.07700297 .00016037 93960-6 28103-3 0 9992',
'2 19822 75.0289 316.1359 1764871 261.4338 78.3806 12.01292738 15108']
>>> tle_list = list(zip(*tuple(
... all_lines[idx::3]
... for idx in range(3)
... )))
>>> len(tle_list) # doctest: +SKIP
65
>>> print(*tle_list[1], sep='\n') # doctest: +SKIP
HST
1 20580U 90037B 19321.38711875 .00000471 00000-0 17700-4 0 9991
2 20580 28.4699 288.8102 0002495 321.7771 171.5855 15.09299865423838
.. note::
As the above code example will download the latest (daily!) TLE from
`Celestrak`_, you'll be seeing a different set of numbers.
Because re-structuring the TLE text from the download or file into a list of
tuples is often needed, there is a utility routine to do just that
.. doctest-remote-data::
>>> tle_list = cysgp4.tle_tuples_from_text(ctrak_science)
And it will do the exact same thing. For testing, `~cysgp4` has one example
TLE file onboard::
>>> ctrak_science = cysgp4.get_example_tles()
>>> ctrak_science[:200]
'AKEBONO (EXOS-D) \r\n1 19822U 89016A 19321.49921565 .00016421 94291-6 28704-3 0 9992\r\n2 19822 75.0304 327.7460 1766579 276.4058 63.9734 12.00957719 13956\r\nHST \r\n1 2058'
>>> tle_list = cysgp4.tle_tuples_from_text(ctrak_science)
>>> tle_list # doctest: +SKIP
[('AKEBONO (EXOS-D) ',
'1 19822U 89016A 19321.49921565 .00016421 94291-6 28704-3 0 9992',
'2 19822 75.0304 327.7460 1766579 276.4058 63.9734 12.00957719 13956'),
('HST ',
'1 20580U 90037B 19321.38711875 .00000471 00000-0 17700-4 0 9991',
'2 20580 28.4699 288.8102 0002495 321.7771 171.5855 15.09299865423838'),
...
('ZHANGZHENG-1 (CSES) ',
'1 43194U 18015C 19322.55091545 .00001006 00000-0 48427-4 0 9993',
'2 43194 97.4135 85.6925 0016973 127.3937 1.4087 15.20935646 99453')]
Now, let's feed this into `~cysgp4.PyTle`::
>>> hst_tle = cysgp4.PyTle(*tle_list[1])
>>> hst_tle
>>> print(hst_tle)
Norad Number: 20580
Int. Designator: 90037B
Epoch: 2019-11-17 09:17:27.060000 UTC
Orbit Number: 42383
Mean Motion Dt2: 0.00000471
Mean Motion Ddt6: 0.00000000
Eccentricity: 0.00024950
BStar: 0.00001770
Inclination: 28.46990000
Right Ascending Node: 288.81020000
Argument Perigee: 321.77710000
Mean Anomaly: 171.58550000
Mean Motion: 15.09299865
PyCoordGeodetic
^^^^^^^^^^^^^^^
This is one of three coordinate frames, which are implemented. The `Geodetic
frame `_ is the same as the
well-known geographic longitude and latitude, plus an altitude above Earth's
surface. As such, it is fixed to the (rotating) Earth, in contrast to the
`ECI `_ system (see
below).
.. note::
The `Geodetic frame`_ goes usually along with a specific definition of
the Geoid (Earth's ellipsoid). For many application, the
`WGS84 `_ is used.
However, most TLEs, which can be downloaded, still use WGS72. The
underlying `C++ SGP4 Satellite Library`_ unfortunately has the WGS72
hardcoded.
Constructing and using a `~cysgp4.PyCoordGeodetic` object is straightforward::
>>> from cysgp4 import PyCoordGeodetic
>>> lon_deg, lat_deg = 6.88375, 50.525
>>> alt_km = 0.366
>>> geo = PyCoordGeodetic(lon_deg, lat_deg, alt_km)
>>> geo
>>> # Access is also possible via properties, e.g.:
>>> geo.lon
6.88375
>>> geo.alt = 0.4
PyEci
^^^^^
The `~cysgp4.PyEci` class holds an `ECI`_ location (latitude, longitude,
altitude) for a particular datetime. Note, internally, the coordinates (and
velocities) are stored in Cartesian form (read-only!). One can access these,
via the `~cysgp4.PyEci.loc` and `~cysgp4.PyEci.vel` properties. Setting
new parameters is only possible via a geographic location, i.e., via
`~cysgp4.PyCoordGeodetic`, which is then converted to Cartesian ECI based on
the provided `~cysgp4.PyDateTime`::
>>> from cysgp4 import PyEci, PyCoordGeodetic, PyDateTime
>>> pydt = PyDateTime.from_mjd(55555.)
>>> lon_deg, lat_deg = 6.88375, 50.525
>>> alt_km = 0.366
>>> geo = PyCoordGeodetic(lon_deg, lat_deg, alt_km)
>>> eci = PyEci(pydt, geo)
>>> eci
>>> # Access is also possible via properties, e.g.:
>>> eci.loc # read-only! # doctest: +FLOAT_CMP
(-725.3304166274728, 3997.924210010933, 4900.402205553537)
>>> eci.vel # read-only! # doctest: +FLOAT_CMP
(-0.2915332651982093, -0.05289193431368386, 0.0)
>>> eci.pydt = PyDateTime.from_mjd(55556.)
>>> # or the update method:
>>> eci.update(pydt, PyCoordGeodetic(0, 0, 0))
PyCoordTopocentric
^^^^^^^^^^^^^^^^^^
The third of the three coordinate frames is the `topocentric frame
`_ (or horizontal
frame), which is provided by `~cysgp4.PyCoordTopocentric`. It holds a
topocentric location (azimuth, elevation, range/distance and distance/range
rate).
.. note::
The topocentric position of a satellite is always relative to
an observer (in SGP4 the observer is defined in geographic coordinates).
The distance and distance change (aka range rate) is thus the distance
between the satellite and the observer. However,
`~cysgp4.PyCoordTopocentric` only holds the azimuth, elevation, distance
and distance rate parameters, but contains no information on the
observer. It is only useful in conjuction with a `~cysgp4.PyObserver` (
which holds a reference to a geographic location) and a datetime; see
also the discussion of `~cysgp4.Satellite` here: :ref:`satellite-label`.
Constructing and using a `~cysgp4.PyCoordTopocentric` object is
straightforward::
>>> from cysgp4 import PyCoordTopocentric
>>> az_deg, el_deg = 130.1, 10.53
>>> dist_km, dist_rate_kms = 1200., 0.03
>>> topo = PyCoordTopocentric(az_deg, el_deg, dist_km, dist_rate_kms)
>>> topo
>>> # Access is also possible via properties, e.g.:
>>> topo.az
130.1
>>> topo.dist = 1000.
PyObserver
^^^^^^^^^^
The `~cysgp4.PyObserver` class holds the location (as `~cysgp4.PyEci`) of
an observer.
.. note::
Usually, a `~cysgp4.PyDateTime` is attached to an ECI location in the
underlying `C++ SGP4 Satellite Library`_.
However, `~cysgp4.PyObserver` does not allow to specify the datetime
and internally it is always set to zero.
Using `~cysgp4.PyObserver` is similar to working with
`~cysgp4.PyCoordGeodetic`::
>>> from cysgp4 import PyObserver, PyCoordGeodetic
>>> lon_deg, lat_deg = 6.88375, 50.525
>>> alt_km = 0.366
>>> obs = PyObserver(lon_deg, lat_deg, alt_km)
>>> obs
>>> # Access is also possible via location property:
>>> obs.loc
>>> obs.loc = PyCoordGeodetic(1, 2, 3)
We will see in the next section, how the observer interplays with the
satellite orbit calculation.
.. _satellite-label:
Satellite interface
^^^^^^^^^^^^^^^^^^^
With the `~cysgp4.Satellite` class, one can calculate positions of a
satellite at given times. Here, the satellite is defined via an TLE (see
`~cysgp4.PyTle`). Furthermore, to be able to calculate apparent positions of
the satellite, an observer (see `~cysgp4.PyObserver`) needs to be defined.
The position calculations are lazy, i.e., only calculated if the newly
requested time differs by a certain amount from the last time at which a
calculation was performed. The granularity of this "cache" can be defined via
the `mjd_cache_resolution` parameter.
The following demonstrates how a typical use of the `~cysgp4.Satellite` class
would look like::
>>> from cysgp4 import *
>>> # Specify a datetime, for which the satellite position is desired
>>> pydt = PyDateTime.from_mjd(58805.57)
>>> # Define the observer...
>>> lon_deg, lat_deg = 6.88375, 50.525
>>> alt_km = 0.366
>>> obs = PyObserver(lon_deg, lat_deg, alt_km)
>>> # ... and TLE (i.e., satellite parameters)
>>> hst_tle = PyTle(
... 'HST',
... '1 20580U 90037B 19321.38711875 .00000471 00000-0 17700-4 0 9991',
... '2 20580 28.4699 288.8102 0002495 321.7771 171.5855 15.09299865423838',
... )
>>> # All of this goes into the Satellite constructor
>>> sat = Satellite(hst_tle, obs, pydt)
>>> # We can now query positions...
>>> sat.eci_pos().loc # ECI cartesian position # doctest: +FLOAT_CMP
(5879.5931344459295, 1545.7455647032068, 3287.4155452595)
>>> sat.eci_pos().vel # ECI cartesian velocity # doctest: +FLOAT_CMP
(-1.8205895517672226, 7.374044252723081, -0.20697960810978586)
>>> sat.geo_pos() # geographic position # doctest: +FLOAT_CMP
>>> sat.topo_pos() # topocentric position # doctest: +FLOAT_CMP
>>> # ... also for different times, also for different times
>>> sat.mjd += 1 / 720. # one minute later
>>> sat.topo_pos() # doctest: +FLOAT_CMP
>>> sat.topo_pos().az, sat.topo_pos().el # doctest: +FLOAT_CMP
(54.84463503781068, -38.274852915850126)
>>> # Change by less than cache resolution (1 ms) produces the same result
>>> sat.mjd += 0.0005 / 86400. # 0.5 ms
>>> sat.topo_pos().az, sat.topo_pos().el # doctest: +FLOAT_CMP
(54.84463503781068, -38.274852915850126)
>>> # Change by another 0.5 ms triggers re-calculation
>>> sat.mjd += 0.00051 / 86400.
>>> sat.topo_pos().az, sat.topo_pos().el # doctest: +FLOAT_CMP
(54.844568313870965, -38.274885794151324)
Option 2: Batch calculations
----------------------------
All of the above can of course be combined with for-loops to process many
different times, or several satellites. And this would be reasonably fast.
However, as the wrapper around `C++ SGP4 Satellite Library`_ is powered
with `Cython `_ it was possible to do implement the
for-loops directly in Cython to reach C-level speed. Furthermore, `OpenMP`_
can parallelize the computation such that it makes use of multi-core
platforms.
The batch processing is implemented with the `~cysgp4.propagate_many`
function. We start by defining datetimes (as MJDs), observers, and TLEs::
>>> import numpy as np
>>> from cysgp4 import PyTle, PyObserver, propagate_many
>>> from cysgp4 import get_example_tles, tles_from_text
>>> tle_text = get_example_tles()
>>> tles = np.array(tles_from_text(tle_text))
>>> observers = np.array([
... PyObserver(6.88375, 50.525, 0.366),
... PyObserver(16.88375, 50.525, 0.366),
... ])
>>> mjds = np.linspace(58805.5, 58806.5, 1000) # 1000 time steps
The `~cysgp4.propagate_many` will apply numpy's `broadcasting rules
`_. As in
this case all input parameters are arrays, we need to add two (length-1) axes
to each of them to make them "compatible" (think of this as an "outer
product". Here we choose to make `mjds` vary on the 1st axis, `tles` on the
3rd axis, and observers on the 2nd axis.::
>>> result = propagate_many(
... mjds[:, np.newaxis, np.newaxis],
... tles[np.newaxis, np.newaxis, :20], # use first 20 TLEs,
... observers[np.newaxis, :, np.newaxis],
... )
Here, `result` is a Python dictionary, with the following entries:
>>> print(sorted(result.keys()))
['eci_pos', 'eci_vel', 'geo', 'topo']
Each entry is an array containing the coordinates of the frame (on the last
axis). The broadcasted input arrays form the first axes of the output::
>>> # Shapes are as follows
>>> print(np.broadcast(
... mjds[:, np.newaxis, np.newaxis],
... tles[np.newaxis, np.newaxis, :20], # use first 20 TLEs,
... observers[np.newaxis, :, np.newaxis],
... ).shape)
(1000, 2, 20)
>>> print(result['eci_pos'].shape, result['topo'].shape)
(1000, 2, 20, 3) (1000, 2, 20, 4)
>>> # First mjd, first sat, first observer:
>>> print('{:.3f} {:.3f} {:.3f}'.format(*result['eci_pos'][0, 0, 0]))
6611.457 -4125.652 762.758
>>> print('{:.3f} {:.3f} {:.3f} {:.3f}'.format(*result['topo'][0, 0, 0]))
90.982 -34.119 9361.523 -2.074
>>> # One can also skip over coordinate frames.
>>> result = propagate_many(
... mjds[:, np.newaxis, np.newaxis],
... tles[np.newaxis, np.newaxis, :20], # use first 20 TLEs,
... observers[np.newaxis, :, np.newaxis],
... do_eci_pos=False, do_eci_vel=False, do_geo=False, do_topo=True
... )
>>> print(sorted(result.keys()))
['topo']
Benchmarks
==========
Because we make use of `Cython`_, `~cysgp4` turns out to be really fast. In
the following, a few benchmarks are run. First, it should be mentioned that
`~cysgp4.propagate_many` will utilize all CPU cores that are available on
your system. It is possible to change that by calling::
>>> num_cores = 1
>>> cysgp4.set_num_threads(num_cores) # set desired number of cores
With the Python `~timeit` package, it is possible to do some speed testing.
We can also try out, how well the parallelization is working::
>>> import numpy as np
>>> import cysgp4
>>> import timeit
>>> def prepare_arrays():
...
... tle_text = cysgp4.get_example_tles()
... tles = np.array(cysgp4.tles_from_text(tle_text))
... observers = np.array([
... cysgp4.PyObserver(6.88375, 50.525, 0.366),
... cysgp4.PyObserver(16.88375, 50.525, 0.366),
... ])
... mjds = np.linspace(58805.5, 58806.5, 1000) # 1000 time steps
...
... return (
... mjds[:, np.newaxis, np.newaxis],
... tles[np.newaxis, np.newaxis, :],
... observers[np.newaxis, :, np.newaxis],
... )
>>> mjds, tles, observers = prepare_arrays()
>>> for num_cores in [1, 2, 4, 8, 16]: # doctest: +SKIP
...
... cysgp4.set_num_threads(num_cores)
... run_times = timeit.Timer(
... 'cysgp4.propagate_many(mjds, tles, observers)',
... globals=globals()
... ).repeat(3, number=10) # do 3 repeats, 10 runs each
...
... print('{:2d} cores: {:.2f} seconds (best of three)'.format(
... num_cores, min(run_times)))
1 cores: 7.16 seconds (best of three)
2 cores: 3.63 seconds (best of three)
4 cores: 1.86 seconds (best of three)
8 cores: 0.99 seconds (best of three)
16 cores: 0.58 seconds (best of three)
How does this relate to using Python for-loops with the wrapped classes (see
:ref:`wrapped-approach-label`)? There is a function,
`~cysgp4.propagate_many_slow` that does exactly that and is merely included
in the package for benchmarking. You can see its `source code here <_modules/
cysgp4/helpers.html#propagate_many_slow>`_. We can compare this to the faster
routine like this::
>>> run_times = timeit.Timer( # doctest: +SKIP
... 'cysgp4.propagate_many_slow(mjds, tles, observers)',
... globals=globals()
... ).repeat(3, number=10) # do 3 repeats, 10 runs each
>>> print('slow version: {:.2f} seconds (best of three)'.format(
... min(run_times))) # doctest: +SKIP
slow version: 74.94 seconds (best of three)
Last but not least, we are interested in a comparison with the well-known
`sgp4`_ Python package (by Brandon Rhodes). Again, for the sole purpose of
benchmarking, a function was added, which uses `sgp4`_ but has the same
interface as `~cysgp4.propagate_many`. The source code of
`~cysgp4.propagate_many_sgp4` is available in the `manual <_modules/
cysgp4/helpers.html#propagate_many_sgp4>`_, as well. Here are the results::
>>> run_times = timeit.Timer( # doctest: +SKIP
... 'cysgp4.propagate_many_sgp4(mjds, tles, observers)',
... globals=globals()
... ).repeat(3, number=10) # do 3 repeats, 10 runs each
>>> print('sgp4 version: {:.2f} seconds (best of three)'.format(
... min(run_times))) # doctest: +SKIP
sgp4 version: 236.46 seconds (best of three)
As you can see, there is some significant speed-up to be gained.
See Also
========
- `C++ SGP4 Satellite Library `_ by Daniel Warner
- `sgp4 Python package `_ by Brandon Rhodes
- `Two-line element strings (TLE)
`_
- `Orbital elements `_
- `Simplified perturbations models
`_
Reference/API
=============
.. automodapi:: cysgp4
:inherited-members:
:no-inheritance-diagram:
:no-main-docstr:
.. :no-heading: