ENH: Update usage of sgp4

[jklenzing]

Description

Addresses #61, and #55 (partial); Replaces #8

Adds Keplerian inputs to the sgp4 object to streamline the user interface and bypass TLEs.

• Includes several method functions to convert between altitude of perigee / apogee and eccentricity / mean_motion
• Includes array calculation of orbits to speed things up
• Adds version limits on sgp4 (>=2.7) [Current version is 2.20]
• Improves generation of metadata for missions_sgp4
• Removed unused kwargs from missions_sgp4

NOTE: the array generation routine produces slight differences in the orbital propagation, on the order of 20 cm, versus the previous version.

Type of change

• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality
  to not work as expected)
• This change requires a documentation update

How Has This Been Tested?

Assuming module is registered:

import numpy as np
import pysat
from pysatMissions.instruments.methods import orbits

from sgp4.earth_gravity import wgs72
from sgp4.io import twoline2rv

# TLEs (Two Line Elements for ISS)
# format of TLEs is fixed and available from wikipedia...
line1 = '1 25544U 98067A   18001.00000000  .00002728  00000-0  48567-4 0  9998'
line2 = '2 25544  51.6402 181.0633 0004018  88.8954  22.2246 15.54059185113452'

# Use sgp4 to extract Keplerian elements from TLEs
satellite = twoline2rv(line1, line2, wgs72)

# Define object old school, using TLEs
old = pysat.Instrument('missions', 'sgp4', TLE1=line1, TLE2=line2,
                       num_samples=86400)

# Grab altitude of perigee, apogee from new conversion routines
alt_periapsis, alt_apoapsis = orbits.convert_from_keplerian(satellite.ecco,
                                                            satellite.no_kozai)

# Calculate using new routines.  Note that input should be in degrees
new = pysat.Instrument('missions', 'sgp4',
                       alt_periapsis=alt_periapsis,
                       alt_apoapsis=alt_apoapsis,
                       inclination=np.degrees(satellite.inclo),
                       bstar=satellite.bstar,
                       arg_periapsis=np.degrees(satellite.argpo),
                       raan=np.degrees(satellite.nodeo),
                       mean_anomaly=np.degrees(satellite.mo),
                       num_samples=86400)

# Load a day of each, inspect the difference
old.load(2018, 1)
new.load(2018, 1)
diff = new.data - old.data

Difference in position for test cases better than microns in each dimension.

Test Configuration:

• Operating system: Mac OS X 10.15.7
• Version number: python 3.8.2
• Any details about your local setup that are relevant: sgp4 2.20

Checklist:

• Make sure you are merging into the develop (not main) branch
• My code follows the style guidelines of this project
• I have performed a self-review of my own code
• I have commented my code, particularly in hard-to-understand areas
• I have made corresponding changes to the documentation
• My changes generate no new warnings
• I have added tests that prove my fix is effective or that my feature works
• New and existing unit tests pass locally with my changes
• Any dependent changes have been merged and published in downstream modules
• Add a note to CHANGELOG.md, summarizing the changes

[jklenzing]

Pinging @rstoneback and @JonathonMSmith. Not ready for "review", but testing and discussion is welcome.

[jklenzing]

Additional test: compare the mean motion to the expected orbital period

import numpy as np
from pysatMissions.instruments.methods import orbits
ecc1, mm1 = orbits.convert_to_keplerian(400)
period1 = 2 * np.pi / mm1
ecc2, mm2 = orbits.convert_to_keplerian(400, 850)
period2 = 2 * np.pi / mm2

yields orbital periods of 92.4 min and 97 min respectively. The eccentricity for the second should match C/NOFS (0.032).

[jklenzing]

Note: orbital conversion routines here assume a spherical Earth. This is resulting in a ~10km error when converting to altitude above an ellipsoid. Needs to be fixed here or in the next pull (geospacepy).

The version of this with geospacepy (with ECEF and lat/lon/alt data) can be found in the enh/geospacepy branch. No pull yet as this doesn't work with github actions.

[rstoneback]

Thanks @jklenzing. Being able to specify an orbit from a more intuitive physical base rather than sort out the format of TLEs is a nice user improvement.

[rstoneback]

This could use a bit more description.

[rstoneback]

This could use a bit more description.

[rstoneback]

Note: orbital conversion routines here assume a spherical Earth. This is resulting in a ~10km error when converting to altitude above an ellipsoid. Needs to be fixed here or in the next pull (geospacepy).

The version of this with geospacepy (with ECEF and lat/lon/alt data) can be found in the enh/geospacepy branch. No pull yet as this doesn't work with github actions.

So the altitudes that come out with respect to a mean Earth aren't correct when converted into a geodetic altitude?

[rstoneback]

Missed this the first time around,

the array generation routine produces slight differences in the orbital propagation, on the order of 20 cm, versus the previous version.

What do you think this is from?

[jklenzing]

Note: orbital conversion routines here assume a spherical Earth. This is resulting in a ~10km error when converting to altitude above an ellipsoid. Needs to be fixed here or in the next pull (geospacepy).
The version of this with geospacepy (with ECEF and lat/lon/alt data) can be found in the enh/geospacepy branch. No pull yet as this doesn't work with github actions.

So the altitudes that come out with respect to a mean Earth aren't correct when converted into a geodetic altitude?

In the other branch, I tried to make a C/NOFS like satellite, 400km x 850km. The resulting orbit was 390 by 840. The Earth radius at the equator is about 7 km greater than the mean radius. This gets us close, but there still could be an error in my math. I need to check my assumptions with Liam once I pull together some summary plots from the other branch.

[jklenzing]

Missed this the first time around,

the array generation routine produces slight differences in the orbital propagation, on the order of 20 cm, versus the previous version.

What do you think this is from?

The 20 cm comes from sgp4. The sgp4 1.x routines we were using previously are now stored in sgp4.io for backwards compatibility. The new 2.x routines (including the array calculations) are in sgp4.api. Using only the new api, the difference between initializing with Keplerian elements and with TLEs is on the order of microns, whereas comparing with the old 1.x methods is about 20 cm. Sorry, I updated code above and didn't quite adjust all the comments for clarity.

[jklenzing]

To compare with latitudes, etc, I've added a draft of the experimental code at #66.