# satkit: Satellite Mission Analysis and Design for Python
#
# Copyright (C) 2023 Egemen Imre
#
# Licensed under GNU GPL v3.0. See LICENSE.rst for more info.
"""
TLE helper functions and factory.
"""
from enum import Enum
import numpy as np
from org.orekit.propagation.analytical.tle import TLE
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import Constants, IERSConventions
from pint import Quantity
from satkit import u
from satkit.propagation.orbits import OrbitUtils
[docs]class TleDefaultUnits(Enum):
"""TLE Default Units."""
EPOCH = None
DATE = None
SAT_NR = None
LAUNCH_NR = None
LAUNCH_PIECE = None
LAUNCH_YR = None
CLASSIFICATION = None
EPHEMERIS_TYPE = None
INCL = "rad"
INCLINATION = "rad"
E = None
ECCENTRICITY = None
MEAN_ANOMALY = "rad"
ARG_OF_PERIGEE = "rad"
RAAN = "rad"
MEAN_MOTION = "rad/s"
N = "rad/s"
N_DOT = "rad/s**2"
N_DOT_DOT = "rad/s**3"
BSTAR = None
ELEMENT_NR = None
REV_NR = None
[docs]class TleFactory:
"""
This class generates TLEs for common orbit types.
A two-line element set (TLE) is a data format encoding a list of TEME
(True Equator, Mean Equinox) mean orbital elements
of an Earth-orbiting object for a given point in time, called the Epoch Time.
These orbital elements are solely for use with the SGP4 propagator due to the
analytical orbit theory used in its derivation.
See the `TLE page in Wikipedia <https://en.wikipedia.org/wiki/Two-line_element_set>`_
or `Space-Track definition <https://www.space-track.org/documentation#tle>`_
for more information.
The TLEs are usually generated by external sources and are used to propagate
the orbits with the initial condition encapsulated by the TLE.
"""
[docs] @staticmethod
@u.wraps(None, (None, "rad", None, None, None, None, None, None, None), False)
def init_geo(
epoch: type[AbsoluteDate],
longitude: float | Quantity,
launch_year: int = 2000,
launch_nr: int = 1,
launch_piece: str = "A",
sat_num: int = 99999,
classification: str = "U",
rev_nr: int = 0,
el_nr: int = 1,
) -> TLE:
"""
Initialises a geostationary satellite TLE.
Due to the nature of the Earth's geopotential, the orbit may drift off
by some degrees within weeks.
Parameters
----------
epoch : type[AbsoluteDate]
Epoch Time corresponding to the orbital elements (nominally very near
the time of true ascending node passage)
longitude : float or Quantity
initial longitude of the satellite [rad]
launch_year : int
launch year (e.g., '2014')
launch_nr : int
launch number within the year (e.g., '17')
launch_piece : str
launch piece (3 chars max)
sat_num : int
satellite catalog number (e.g., '39552'; see TLE class documentation)
classification : str
Classification (`U` for Unclassified, `C` for Classified, `S` for Secret)
rev_nr : int
Revolution number of the object at Epoch Time (should be in range 0 <= `rev_nr` < 10000) [revolutions]
el_nr : int
Element set number (should be in range 0 <= `el_nr` < 10000). Incremented when a new TLE is generated for this object.
Returns
-------
TLE
`TLE` object initialised with the required GEO parameters.
"""
# init GEO specific values - period is one sidereal day
mean_motion = 2 * np.pi / (1.0 * u.sidereal_day).m_as("sec")
# sidereal time hour angles in radians
ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, True)
gmst = IERSConventions.IERS_2010.getGMSTFunction(ut1).value(epoch)
sidereal_time = gmst % (2 * np.pi)
# RAAN in rad
raan = longitude + sidereal_time
# print(f"sidereal time: {np.degrees(sidereal_time)} deg")
bstar = 0.0 # no drag
n_dot = 0.0 # mean motion assumed constant
n_dotdot = 0.0 # mean motion assumed constant
# init standard circular orbit (angles in rad)
inclination = 0.0
arg_perigee = 0.0
mean_anomaly = 0.0
eccentricity = 0.0
eph_type = 0
# all Quantity values converted to values by now
tle = TLE(
sat_num,
classification,
launch_year,
launch_nr,
launch_piece,
eph_type,
el_nr,
epoch,
mean_motion,
n_dot,
n_dotdot,
eccentricity,
inclination,
arg_perigee,
raan,
mean_anomaly,
rev_nr,
bstar,
)
return tle
[docs] @staticmethod
@u.wraps(
None,
(
None,
"m",
None,
None,
"rad",
"rad",
None,
None,
None,
None,
None,
None,
None,
None,
),
False,
)
def init_sso(
epoch: type[AbsoluteDate],
altitude: float | Quantity,
ltan: float,
eccentricity: float = 0.0,
arg_perigee: float | Quantity = 0.0 * u.rad,
mean_anomaly: float | Quantity = 0.0 * u.rad,
bstar: float = 0.0,
launch_year: int = 2000,
launch_nr: int = 1,
launch_piece: str = "A",
sat_num: int = 99999,
classification: str = "U",
rev_nr: int = 0,
el_nr: int = 1,
) -> TLE:
"""
Initialises a sun-synchronous Earth orbit TLE.
Sets the inclination in accordance with the target node rotation rate that
satisfies sun-synchronous conditions up to J4. Sets the RAAN in accordance
with the requested Local Time of the Ascending Node (LTAN).
Target node rotation rate is 360 deg / 365.2421897 days (See
Fundamentals of Astrodynamics Vallado 4th ed pg. 863).
The node rotation rate is calculated via a J4 secular drift rate (See
Fundamentals of Astrodynamics Vallado 4th ed pg. 650).
Equations use WGS84 parameters for Earth equator radius and mu.
J2 and J4 are from EGM96.
Parameters
----------
epoch : type[AbsoluteDate]
Epoch Time corresponding to the orbital elements (nominally very near
the time of true ascending node passage)
altitude : float or Quantity
altitude above `earth_radius` around Equator [m]
ltan : float
Local Time of the Ascending Node defined as [0, 24)
eccentricity : float
mean eccentricity of the orbit [dimensionless]
arg_perigee : float or Quantity
TEME mean argument of perigee [rad]
mean_anomaly : float or Quantity
mean anomaly of the orbit [rad]
bstar : float or Quantity
sgp4 type drag coefficient [1 / earth radius] (see TLE class documentation)
launch_year : int
launch year (e.g., '2014')
launch_nr : int
launch number within the year (e.g., '17')
launch_piece : str
launch piece (3 chars max)
sat_num : int
satellite catalog number (e.g., '39552'; see TLE class documentation)
classification : str
Classification (`U` for Unclassified, `C` for Classified, `S` for Secret)
rev_nr : int
Revolution number of the object at Epoch Time (should be in range 0 <= `rev_nr` < 10000) [revolutions]
el_nr : int
Element set number (should be in range 0 <= `el_nr` < 10000). Incremented when a new TLE is generated for this object.
Returns
-------
TLE
`TLE` object initialised with the required SSO parameters.
"""
# MU in m3/s2
mu = Constants.WGS84_EARTH_MU
# R_E in m
earth_radius = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
# Set J2 and J4 (EGM96)
j2 = 0.0010826266835531513
# j4 = -1.619621591367e-06
# j6 = 5.406812391070849e-07
# check Quantity objects and add units in m
# sma = (
# altitude + earth_radius
# if isinstance(altitude, u.Quantity)
# else altitude * u.m + earth_radius
# )
# sma = altitude [m] + earth radius [m]
sma = altitude + earth_radius
# mean motion [rad/s]
mean_motion = np.sqrt(mu / sma**3)
# mean motion assumed constant
n_dot = 0.0 # rad/s2
n_dotdot = 0.0 # rad/s3
# Set ephemeris type
eph_type = 0
# target omega_dot value for sun sync
om_dot_sun_sync = (360.0 / 365.2421897 * u["deg/day"]).m_as("rad/s")
# Start inclination procedure
# ---------------------------
# Inclination initial guess from J2 only
inclination = np.arccos(
(-2 * sma**3.5 * om_dot_sun_sync * (1 - eccentricity**2) ** 2)
/ (3 * earth_radius**2 * j2 * np.sqrt(mu))
)
# Set up iteration to compute the J2^2 and J4 correction
om_dot_next = 0.0
i = inclination
e = eccentricity
om_dot = om_dot_sun_sync
om_dot_threshold = (1e-8 * u["deg/day"]).m_as("rad/s")
while (om_dot_sun_sync - om_dot_next) > om_dot_threshold:
# in rad/s
om_dot_next = OrbitUtils.compute_raan_drift_rate(sma, e, i).m_as("rad/s")
d_om = om_dot_next - om_dot
di = -d_om / (om_dot * np.tan(i))
i = i + di
om_dot = om_dot_next
# Set final Inclination value
inclination = i
# Start RAAN procedure
# ---------------------------
# sidereal time hour angles in radians
ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, True)
gmst = IERSConventions.IERS_2010.getGMSTFunction(ut1).value(epoch)
sidereal_time = gmst % (2 * np.pi)
frac_day = (
epoch.getComponents(TimeScalesFactory.getUTC())
.getTime()
.getSecondsInUTCDay()
* u.sec
).m_as("day")
# compute the RAAN value and normalise (in rad)
raan = sidereal_time + ((ltan / 24.0) - frac_day) * 2 * np.pi
raan = np.mod(raan, 2 * np.pi)
# all Quantity values converted to values by now
# numpy uses float64 but Orekit requires float
tle = TLE(
sat_num,
classification,
launch_year,
launch_nr,
launch_piece,
eph_type,
el_nr,
epoch, # user defined
float(mean_motion), # in 1/s (or rad/s)
n_dot, # defined 0.
n_dotdot, # defined 0.
eccentricity, # user defined float
float(inclination), # in rad
arg_perigee, # in rad
float(raan), # in rad
mean_anomaly, # in rad
rev_nr,
bstar, # user defined
)
return tle
[docs]class TLEUtils:
"""
Utilities for TLE based orbits.
Basic conversions and functions of orbit related operations for convenience.
"""
[docs] @staticmethod
def compute_sma(tle: TLE) -> Quantity:
"""
Computes the (mean) semimajor axis from the TLE.
The orbit plane rotation rate is calculated via a J4 secular drift rate (See
Fundamentals of Astrodynamics Vallado 4th ed pg. 650).
Equations use WGS84 parameters for mu.
Parameters
----------
tle
TLE
Returns
-------
sma : Quantity
Semimajor axis with units [m]
"""
return OrbitUtils.compute_sma(tle.getMeanMotion())
[docs] @staticmethod
def compute_raan_drift_rate(
tle: TLE,
include_j4: bool = True,
) -> Quantity:
"""
Computes the RAAN (or orbital plane) drift rate from the TLE.
The orbit plane rotation rate is calculated via a J4 secular drift rate (See
Fundamentals of Astrodynamics Vallado 4th ed pg. 650).
Equations use WGS84 parameters for Earth equator radius and mu.
J2 and J4 are from EGM96.
Parameters
----------
tle
TLE
include_j4
True if J2^2 and J4 effects are to be included, False for J2 only.
Returns
-------
Quantity
RAAN drift rate in angles/time (e.g. deg/day)
"""
sma = TLEUtils.compute_sma(tle)
return OrbitUtils.compute_raan_drift_rate(
sma, tle.getE(), tle.getI(), include_j4
)