API

Acceleration

Classes for modeling accelerations.

class ssapy.accel.Accel(time_breakpoints=None)

Base class for accelerations.

class ssapy.accel.AccelConstNTW(accelntw, time_breakpoints=None)

Constant acceleration in NTW coordinates.

Intended to enable maneuvers. Semimajor axis changes are often done by accelerating in the in-track direction at perigee, while inclination change maneuvers are done by accelerating in the cross-track direction. So these maneuvers are both conveniently implemented in NTW space.

This class supports a constant acceleration burn in a fixed direction in NTW space at a specified time and duration.

Parameters:

• accelntw (array_like, shape(3,)) – Direction and magnitude of acceleration in NTW frame, m/s^2.

• time_breakpoints (array_like) – Times in GPS seconds when acceleration should be turned on and off. These alternate; first time is on, second time is off, third time is back on. must be sorted.

class ssapy.accel.AccelDrag(recalc_threshold=2592000, **defaultkw)

Acceleration due to atmospheric drag.

This class uses the Harris-Priester density model, which includes diurnal variation in the atmospheric bulge, but omits longer period seasonal variations.

The acceleration also depends on a drag coefficient, which is hard to determine a priori, but takes on typical values around ~2 to ~2.3 for most satellites.

See Section 3.5 of Montenbruck and Gill for more details.

Parameters:

• recalc_threshold (float, optional) – Number of seconds past which the code will recompute the precession/nutation matrix. Default: 86400*30 (30 days)

• defaultkw (dict) – default parameters for kwargs passed to __call__, (area, mass, CR)

class ssapy.accel.AccelEarthRad(**defaultkw)

Acceleration due to Earth radiation pressure.

This is a very simple model in which the direction of the acceleration is directly away from the Earth and the magnitude is modulated by a single solar radiation pressure coefficient CR. The coefficient is 1.0 for purely absorbed light, and 2.0 for purely reflected light. Rough typical values for a variety of different satellite components are

~ 1.2 for solar panels ~ 1.3 for a high gain antenna ~ 1.9 for a aluminum coated mylar solar sail.

The radiation pressure at the Earth’s surface is given as (230 + 459*k) W/m^2, where 230 is from the thermal radiation of the earth, and 459 is the reflected sunlight. k is the illuminated fraction of the earth as seen from the satellite, assuming the earth is point-like (i.e., neglecting that the satellite will see less than a full hemisphere for LEO objects). The radiation pressure goes down like 1/r^2 as an object moves away from the earth.

This is a simplification of the more complex model presented in MG 3.7.1, neglecting spatial variation in the emitted light and the different angles to different parts of the earth.

Parameters:

defaultkw (dict) – default parameters for kwargs passed to __call__, (area, mass, CR)

class ssapy.accel.AccelKepler(mu=398600441800000.0)

Keplerian acceleration. I.e., force is proportional to 1/|r|^2.

Parameters:

mu (float, optional) – Gravitational constant of central body in m^3/s^2. (Default: Earth’s gravitational constant in WGS84).

class ssapy.accel.AccelProd(accel, factor)

Acceleration defined as the product of an acceleration with a constant factor.

class ssapy.accel.AccelSolRad(**defaultkw)

Acceleration due to solar radiation pressure.

This is a relatively simple model in which the direction of the acceleration is directly away from the sun and the magnitude is modulated by a single solar radiation pressure coefficient CR. The coefficient is 1.0 for purely absorbed light, and 2.0 for purely reflected light. Rough typical values for a variety of different satellite components are

~ 1.2 for solar panels ~ 1.3 for a high gain antenna ~ 1.9 for a aluminum coated mylar solar sail.

Additionally, this class models the Earth’s shadow as a cylinder. Either the satellite is inside the cylinder, in which case the acceleration from this term is 0, or the satellite is outside the cylinder, in which case the full magnitude of the acceleration is computed.

More details can be found in Section 3.4 of Montenbruck and Gill.

Parameters:

defaultkw (dict) – default parameters for kwargs passed to __call__, (area, mass, CR)

class ssapy.accel.AccelSum(accels)

Acceleration defined as the sum of other accelerations.

Note, besides directly invoking this class, you can simply use the + operator on any two Accel subclasses too.

Parameters:

accels (list of Accel) – Accelerations to add

Body

Classes representing celestial bodies.

class ssapy.body.Body(mu, radius, position=<function Body.<lambda>>, orientation=<function Body.<lambda>>, harmonics=None)

A celestial body.

Parameters:

• mu (float) – Gravitational parameter of the body in m^3/s^2.

• radius (float) – Radius of the body in meters.

• position (callable, optional) – A callable that returns the position vector of the body in GCRF at a given time. [default: zero vector]

• orientation (callable, optional) – A callable that returns the orientation matrix of the body in GCRF at a given time. [default: identity matrix]

• harmonics (HarmonicCoefficients, optional) – Harmonic coefficients for the body. [default: None]

class ssapy.body.EarthOrientation(recalc_threshold=2592000)

Orientation of earth in GCRF. This is a callable class that returns the orientation matrix at a given time.

Parameters:

recalc_threshold (float) – Threshold for recomputing the orientation matrix. Default is 30 days.

class ssapy.body.MoonOrientation

Orientation of moon in GCRF. This is a callable class that returns the orientation matrix at a given time.

class ssapy.body.MoonPosition

Position of moon in GCRF. This is a callable class that returns the position vector at a given time.

class ssapy.body.PlanetPosition(planet_index)

Position of a planet in GCRF. This is a callable class that returns the position vector at a given time.

class ssapy.body.SunPosition

Position of sun in GCRF. This is a callable class that returns the position vector at a given time.

ssapy.body.get_body(name, model=None)

Get a Body object for a named body.

Parameters:

• name (str) – Name of the body. Must be one of “earth”, “moon”, “sun”, or other supported planets.

• model (str, optional only available for Earth) – Name of the Earth harmonic model to use. Default is EGM84. options: EGM96, EGM2008.

Returns:

body – Body object for the named body.

Return type:

Body

Compute

ssapy.compute.M_v_lambertian(r_sat, times, radius=1.0, albedo=0.2, sun_Mag=4.8, albedo_earth=0.3, albedo_moon=0.12, plot=False)

Calculate the apparent magnitude (M_v) of a satellite due to reflections from the Sun, Earth, and Moon.

This function computes the apparent magnitude of a satellite based on its reflected light from the Sun, Earth, and Moon, using the Lambertian reflection model. It optionally generates plots to visualize the results.

Parameters:

r_sat(n, 3) numpy.ndarray

Position of the satellite in meters.

timesTime or array_like

The times corresponding to the satellite and celestial body positions. Used to obtain the positions of the Sun and Moon.

radiusfloat, optional

Radius of the satellite in meters (default is 1.0 m).

albedofloat, optional

Albedo of the satellite’s surface, representing its reflectivity (default is 0.20).

sun_Magfloat, optional

Solar magnitude (apparent magnitude of the Sun) used in magnitude calculations (default is 4.80).

albedo_earthfloat, optional

Albedo of the Earth, representing its reflectivity (default is 0.30).

albedo_moonfloat, optional

Albedo of the Moon, representing its reflectivity (default is 0.12).

plotbool, optional

If True, generates plots to visualize solar phase angle and magnitudes (default is False).

Returns:

numpy.ndarray

The apparent magnitude (M_v) of the satellite as observed from Earth, considering reflections from the Sun, Earth, and Moon.

Notes:

• The function uses a Lambertian reflection model to compute the fraction of sunlight reflected by the satellite, Earth, and Moon.

• The apparent magnitude is calculated based on the distances between the satellite, Sun, Earth, and Moon, as well as their respective albedos.

• The function generates four subplots if plot is set to True: 1. Solar phase angle of the satellite. 2. Solar magnitude (M_v) of the satellite due to the Sun. 3. Magnitude (M_v) of the satellite due to reflections from the Earth. 4. Magnitude (M_v) of the satellite due to reflections from the Moon.

Example usage:

>>> r_sat = np.array([[1e7, 2e7, 3e7]])
>>> times = Time("2024-01-01")
>>> M_v = M_v_lambertian(r_sat, times, plot=True)
>>> M_v
array([15.63])

ssapy.compute.altaz(orbit, time, observer, propagator=KeplerianPropagator(), obsAngleCorrection=None)

Calculate observed altitude and azimuth of orbiting objects as viewed at specified times and locations.

Parameters:

• orbit (Orbit or list of Orbit (n,)) – Orbits for which to calculate direction cosines.

• time (array_like or astropy.time.Time (m,)) – If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• observer (EarthObserver) – Where on Earth to compute alt and az.

• propagator (Propagator, optional) – The propagator instance to use.

• obsAngleCorrection ({None, "linear", "exact"}, optional) – Correct actual angle to observed angle, meaning account for light-time delay, and aberration due to the observer’s velocity. None means don’t do any correction. “linear” means do an aberration correction and first order light-time correction. “exact” means do an aberration correction and iteratively solve for the exact light-time correction. (The “linear” correction is almost always sufficiently accurate).

Notes

If orbit is scalar-valued (an Orbit instead of a list of Orbit), then that dimension will be squeezed out in the return value. Likewise, if time is scalar, that dimension will be squeezed out.

For Keplerian orbit propagation it is more efficient to use a “vector Orbit” instead of a list of single scalar Orbits.

When doing light time corrections, the time argument is the arrival time of the photons at the observer, as opposed to the emission time at the satellite.

Returns:

alt, az – Altitude and azimuth in radians.

Return type:

array_like (n, m)

ssapy.compute.calc_M_v(r_sat, r_earth, r_sun, r_moon=False, radius=0.4, albedo=0.2, sun_Mag=4.8, albedo_earth=0.3, albedo_moon=0.12, albedo_back=0.5, albedo_front=0.05, area_panels=100, return_components=False)

Calculate the apparent magnitude (M_v) of a satellite due to reflections from the Sun, Earth, and optionally the Moon.

This function computes the apparent magnitude of a satellite based on its reflected light from the Sun, Earth, and optionally the Moon. It uses separate functions to calculate the flux contributions from each of these sources and combines them to determine the overall apparent magnitude.

Parameters:

r_sat(n, 3) numpy.ndarray

Position of the satellite in meters.

r_earth(n, 3) numpy.ndarray

Position of the Earth in meters.

r_sun(n, 3) numpy.ndarray

Position of the Sun in meters.

r_moon(n, 3) numpy.ndarray or False, optional

Position of the Moon in meters. If False, the Moon’s contribution is ignored (default is False).

radiusfloat, optional

Radius of the satellite in meters (default is 0.4 m).

albedofloat, optional

Albedo of the satellite’s surface, representing its reflectivity (default is 0.20).

sun_Magfloat, optional

Solar magnitude (apparent magnitude of the Sun) used in magnitude calculations (default is 4.80).

albedo_earthfloat, optional

Albedo of the Earth, representing its reflectivity (default is 0.30).

albedo_moonfloat, optional

Albedo of the Moon, representing its reflectivity (default is 0.12).

albedo_backfloat, optional

Albedo of the back surface of the satellite (default is 0.50).

albedo_frontfloat, optional

Albedo of the front surface of the satellite (default is 0.05).

area_panelsfloat, optional

Area of the satellite’s panels in square meters (default is 100 m^2).

return_componentsbool, optional

If True, returns the magnitude as well as the flux components from the Sun, Earth, and Moon (default is False).

Returns:

float

The apparent magnitude (M_v) of the satellite as observed from Earth.

dict, optional

If return_components is True, a dictionary containing the flux components from the Sun, Earth, and Moon.

Notes:

• The function uses separate calculations for flux contributions from the Sun, Earth, and Moon:

  • sun_shine calculates the flux from the Sun.

  • earth_shine calculates the flux from the Earth.

  • moon_shine calculates the flux from the Moon (if applicable).

• The apparent magnitude is calculated based on the distances between the satellite, Sun, Earth, and optionally the Moon, as well as their respective albedos and other parameters.

Example usage:

>>> r_sat = np.array([[1e7, 2e7, 3e7]])
>>> r_earth = np.array([1.496e11, 0, 0])
>>> r_sun = np.array([0, 0, 0])
>>> Mag_v = calc_M_v(r_sat, r_earth, r_sun, return_components=True)
>>> Mag_v
(15.63, {'sun_bus': 0.1, 'sun_panels': 0.2, 'earth_bus': 0.05, 'earth_panels': 0.1, 'moon_bus': 0.03, 'moon_panels': 0.07})

ssapy.compute.calc_gamma(r, t)

Calculate the gamma angle between position and velocity vectors in the ITRF frame.

Parameters: r (numpy.ndarray): The position vectors in the GCRF frame, shaped (n, 3), where n is the number of vectors. t (numpy.ndarray or astropy.time.Time): The times corresponding to the position vectors. Can be an array of GPS seconds or an Astropy Time object.

Returns: numpy.ndarray: An array of gamma angles in degrees between the position and velocity vectors for each time point.

Notes: - This function first converts the given position vectors from the GCRF frame to the ITRF frame. - It then calculates the angle between the position and velocity vectors at each time point in the ITRF frame. - If the input time array is an Astropy Time object, it converts it to GPS time before processing.

ssapy.compute.calculate_orbital_elements(r_, v_, mu_barycenter=398600441800000.0)

Calculate the Keplerian orbital elements from position and velocity vectors.

This function computes the Keplerian orbital elements (semi-major axis, eccentricity, inclination, true longitude, argument of periapsis, longitude of ascending node, true anomaly, and specific angular momentum) for one or more celestial objects given their position and velocity vectors.

Parameters:

r_(n, 3) numpy.ndarray

Array of position vectors (in meters) of the celestial objects. Each row represents a position vector.

v_(n, 3) numpy.ndarray

Array of velocity vectors (in meters per second) of the celestial objects. Each row represents a velocity vector.

mu_barycenterfloat, optional

Gravitational parameter (standard gravitational constant times mass) of the central body (default is EARTH_MU). This parameter defines the gravitational influence of the central body on the orbiting object.

Returns:

dict

A dictionary containing the orbital elements for each celestial object: - ‘a’: Semi-major axis (in meters). - ‘e’: Eccentricity (dimensionless). - ‘i’: Inclination (in radians). - ‘tl’: True longitude (in radians). - ‘ap’: Argument of periapsis (in radians). - ‘raan’: Longitude of ascending node (in radians). - ‘ta’: True anomaly (in radians). - ‘L’: Specific angular momentum (in meters squared per second).

Notes:

• Position and velocity vectors should be provided in the same units.

• The function assumes that the input vectors are provided in an array where each row corresponds to a different celestial object.

• Orbital elements are computed using standard orbital mechanics formulas.

• The inclination is measured from the reference plane, and the argument of periapsis and true anomaly are measured in the orbital plane.

Example:

>>> r = np.array([[1e7, 1e7, 1e7], [1e8, 1e8, 1e8]])
>>> v = np.array([[1e3, 2e3, 3e3], [4e3, 5e3, 6e3]])
>>> calculate_orbital_elements(r, v, mu_barycenter=3.986e14)
{'a': [1.5707e7, 2.234e8], 'e': [0.123, 0.456], 'i': [0.785, 0.654], 'tl': [2.345, 3.456], 'ap': [0.123, 0.456], 'raan': [1.234, 2.345], 'ta': [0.567, 1.678], 'L': [1.234e10, 2.345e11]}

ssapy.compute.dircos(orbit, time, obsPos=None, obsVel=None, observer=None, propagator=KeplerianPropagator(), obsAngleCorrection=None)

Calculate observed direction-cosines of orbiting objects as viewed at specified times and positions.

The direction cosines are the cosines of the angles between the vector pointing towards the orbiting object and the x, y, z axes. An equivalent description is that they are the components of the unit vector pointing towards the orbiting object.

Parameters:

• orbit (Orbit or list of Orbit (n,)) – Orbit(s) for which to calculate direction cosines.

• time (array_like or astropy.time.Time (m,)) – If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• obsPos (array_like (m, 3), optional) – Position of observer at given time.

• obsVel (array_like (m, 3), optional) – Velocity of observer at given time. Only required if correcting for diurnal aberration.

• observer (Observer or list of Observers (m,), optional) – Observer(s) used to calculate obsPos and obsVel.

• propagator (Propagator, optional) – The propagator instance to use.

• obsAngleCorrection ({None, "linear", "exact"}, optional) – Correct actual angle to observed angle, meaning account for light-time delay, and aberration due to the observer’s velocity. None means don’t do any correction. “linear” means do an aberration correction and first order light-time correction. “exact” means do an aberration correction and iteratively solve for the exact light-time correction. (The “linear” correction is almost always sufficiently accurate).

Notes

Exactly 1 of obsPos and observer must be supplied. observer and obsPos follow similar broadcasting rules as detailed below explicitly only for obsPos.

The length of time and obsPos must match or be broadcastable to match. If orbit is scalar-valued (an Orbit instead of a list of Orbit), then that dimension will be squeezed out in the return value. Likewise, if both time and obsPos are scalar, that dimension will be squeezed out.

For Keplerian orbit propagation it is more efficient to use a “vector Orbit” instead of a list of single scalar Orbits.

When doing light time corrections, the time argument is the arrival time of the photons at the observer, as opposed to the emission time at the satellite.

Returns:

dircos – Direction cosines on the outer product of orbit(s) and time/obsPos.

Return type:

array_like (n, m, 3)

ssapy.compute.earthShadowCoords(r, time)

Determine components of position r parallel and perpendicular to sun unit vector.

The sun unit vector points from the center of the sun through the center of the Earth. Decomposing a satellite position into these coordinates yields a simple model of whether or not the satellite is in Earth’s shadow:

r_par, r_perp = earthShadowCoords(r, time) inShadow = r_par > 0 and r_perp < EARTH_RADIUS

Parameters:

• r (ndarray (3,)) – Position (GCRF) in meters

• time (float or astropy.time.Time) – If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

Returns:

• r_par (float) – Position of satellite projected onto the sun unit vector in meters.

• r_perp (float) – Distance of satellite from Earth-Sun line in meters.

ssapy.compute.earth_shine(r_sat, r_earth, r_sun, radius, albedo, albedo_earth, albedo_back, area_panels)

Calculate the fractional flux of sunlight reflected from the Earth to the satellite.

This function computes the flux of sunlight reflected from the Earth to the satellite, including contributions from the back surface of the satellite’s solar panels.

Parameters:

r_sat(n, 3) numpy.ndarray

Array of coordinates representing the position of the satellite.

r_earth(n, 3) numpy.ndarray

Array of coordinates representing the position of the Earth.

r_sun(n, 3) numpy.ndarray

Array of coordinates representing the position of the Sun.

radiusfloat

Radius of the satellite in meters.

albedofloat

Albedo of the satellite’s surface.

albedo_earthfloat

Albedo of the Earth.

albedo_backfloat

Albedo of the back surface of the satellite’s solar panels.

area_panelsfloat

Area of the satellite’s solar panels in square meters.

Returns:

dict

Dictionary containing the flux contributions from the Earth to the satellite: - ‘earth_bus’: Fraction of light reflected off the satellite’s bus from the Earth. - ‘earth_panels’: Fraction of light reflected off the satellite’s panels from the Earth.

Notes:

• The function assumes that the solar panels are always facing the Sun and calculates flux based on the phase angle.

• Flux contributions from the back surface of the solar panels are computed.

Example:

>>> r_sat = np.array([[1e7, 1e7, 1e7]])
>>> r_earth = np.array([[1.496e11, 0, 0]])
>>> r_sun = np.array([[0, 0, 0]])
>>> earth_shine(r_sat, r_earth, r_sun, radius=0.4, albedo=0.20, albedo_earth=0.30, albedo_back=0.50, area_panels=100)
{'earth_bus': array([...]), 'earth_panels': array([...])}

ssapy.compute.find_passes(orbit, observers, tStart, tSpan, dt, propagator=KeplerianPropagator(), horizon=np.float64(0.3490658503988659))

Find satellite overhead passes for a collection of observers.

Uses a brute force test of a grid of time points from tStart to tStart+tSpan separated by dt.

Returns passes even if they occur during the daytime or if the satellite is not illuminated by the sun. The only criterion for a successful “pass” is for the topocentric altitude of the satellite to be above the input horizon. More details about a pass can subsequently be obtained by running the refine_passes function.

Note this function is only suitable for `EarthObserver`s and not `OrbitalObserver`s.

Parameters:

• orbit (Orbit) – Satellite orbit in question.

• observers (List of EarthObserver.) – Earth observers for which to check satellite visibility.

• tStart (float or astropy.time.Time) – Beginning of search window. If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• tSpan (float or Quantity) – Time span in which to search for passes. If float, then seconds.

• dt (float or Quantity) – Time increment to use during search. Satellite visibility will be computed every dt increment. Smaller values will decrease the probability that a short duration pass will be missed, but will make the search take longer to complete. If float, then seconds.

• propagator (Propagator, optional) – The propagator instance to use.

• horizon (float or Quantity, optional) – Minimum altitude for which to consider a satellite “visible”. If float, then should be in radians.

Returns:

passDict – keys are EarthObserver`s. values are lists (possibly empty) of `astropy.Time corresponding to visible passes of the satellite. Only one time is returned per pass, so multiple times in the return list indicate multiple distinct passes.

Return type:

dict

ssapy.compute.groundTrack(orbit, time, propagator=KeplerianPropagator(), format='geodetic')

Calculate satellite ground track on the outer product of all supplied times and state vectors or orbits.

Parameters:

• r (array_like (n,3) Position (m)) –

• time (array_like or astropy.time.Time (m,)) – If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• or –

• orbit (r : array_like (n,3) Position (m) or) – Orbit or list of Orbit (n,) Desired orbit(s)

• propagator (Propagator, optional) – The propagator instance to use.

• format ('geodetic' or 'cartesian') – If ‘geodetic’, then returns longitude, latitude, height. If ‘cartesian’, then returns xyz in ITRF frame.

Notes

If orbit or time is scalar valued (as opposed to a list of Orbit, e.g.) then the corresponding dimension of the output will be squeezed out.

For Keplerian orbit propagation it is more efficient to use a “vector Orbit” instead of a list of single scalar Orbits.

Returns:

• lon, lat, height (array_like (n, m, 3)) – Radians and meters.

• or

• x, y, z (array_like(n, m, 3)) – Meters.

ssapy.compute.lagrange_points_lunar_frame()

Calculate the positions of the lunar Lagrange points in the lunar frame, This frame is defined by the coordinate transformation in utils.py gcrf_to_lunar().

Returns:

dict

A dictionary containing the coordinates of the five Lagrange points: - “L1”: Position of the first Lagrange point between the Earth and the Moon. - “L2”: Position of the second Lagrange point beyond the Moon. - “L3”: Position of the third Lagrange point directly opposite the Moon, relative to the Earth. - “L4”: Position of the fourth Lagrange point, forming an equilateral triangle with the Earth and Moon. - “L5”: Position of the fifth Lagrange point, forming an equilateral triangle with the Earth and Moon, but on the opposite side.

Notes:

• The function assumes that the Earth and Moon are the two primary bodies, with the Earth-Moon distance denoted as LD.

• The gravitational parameters of the Earth and Moon are denoted as EARTH_MU and MOON_MU, respectively.

• The positions of L4 and L5 are calculated using the fact that these points form an equilateral triangle with the Earth and Moon.

Example usage:

>>> lagrange_points = lagrange_points_lunar_frame()
>>> lagrange_points["L1"]
array([1.01e6, 0.0, 0.0])
>>> lagrange_points["L4"]
array([1.5e6, 1.5e6, 0.0])

ssapy.compute.lunar_lagrange_points(t)

Calculate the positions of the lunar Lagrange points in the GCRF frame at a given time.

This function computes the positions of the five Lagrange points (L1, L2, L3, L4, and L5) in the Earth-Moon system at a specific time t. It considers the positions of the Earth and Moon and uses the orbital period of the Moon to find the Lagrange points.

Parameters:

tTime

The time at which to calculate the Lagrange points. The position of the Moon at this time is used to compute the Lagrange points.

Returns:

dict

A dictionary containing the coordinates of the five Lagrange points: - “L1”: Position of the first Lagrange point between the Earth and the Moon. - “L2”: Position of the second Lagrange point beyond the Moon. - “L3”: Position of the third Lagrange point directly opposite the Moon, relative to the Earth. - “L4”: Position of the fourth Lagrange point, calculated as the Moon’s position offset by one-sixth of the lunar period. - “L5”: Position of the fifth Lagrange point, calculated as the Moon’s position offset by negative one-sixth of the lunar period.

Notes:

• The function assumes a circular orbit for the Moon.

• The lunar period used is approximately 2.36 million seconds.

• The gravitational parameters of the Earth and Moon are denoted as EARTH_MU and MOON_MU, respectively.

Example usage:

>>> t = Time("2024-01-01")
>>> lagrange_points = lunar_lagrange_points(t)
>>> lagrange_points["L1"]
array([1.02e6, 0.0, 0.0])
>>> lagrange_points["L4"]
array([1.5e6, 1.5e6, 0.0])

ssapy.compute.lunar_lagrange_points_circular(t)

Calculate the positions of the lunar Lagrange points in the GCRF frame for a given time.

This function calculates the positions of the five Lagrange points (L1, L2, L3, L4, and L5) in the Earth-Moon system at a specific time t. It accounts for the rotation of the Moon’s orbit around the Earth, providing the positions in a circular approximation of the Earth-Moon system.

Parameters:

tTime

The time at which to calculate the Lagrange points. The position of the Moon at this time is used to compute the Lagrange points.

Returns:

dict

A dictionary containing the coordinates of the five Lagrange points: - “L1”: Position of the first Lagrange point between the Earth and the Moon. - “L2”: Position of the second Lagrange point beyond the Moon. - “L3”: Position of the third Lagrange point directly opposite the Moon, relative to the Earth. - “L4”: Position of the fourth Lagrange point, forming an equilateral triangle with the Earth and Moon. - “L5”: Position of the fifth Lagrange point, forming an equilateral triangle with the Earth and Moon, but on the opposite side.

Notes:

• The function assumes a circular orbit for the Moon and uses a rotation matrix to align the z-axis with the Moon’s normal vector.

• The positions of L4 and L5 are calculated using a rotation matrix to align with the Moon’s orientation.

• The gravitational parameters of the Earth and Moon are denoted as EARTH_MU and MOON_MU, respectively.

Example usage:

>>> t = Time("2024-01-01")
>>> lagrange_points = lunar_lagrange_points_circular(t)
>>> lagrange_points["L1"]
array([1.02e6, 0.0, 0.0])
>>> lagrange_points["L4"]
array([1.5e6, 1.5e6, 0.0])

ssapy.compute.moon_normal_vector(t)

Calculate the normal vector to the Moon’s orbital plane at a given time.

This function computes the normal vector to the Moon’s orbital plane by taking the cross product of the Moon’s position vectors at two different times: the given time t and one week later. The resulting vector is normalized to provide the direction of the orbital plane’s normal.

Parameters:

tTime

The time at which to calculate the normal vector to the Moon’s orbital plane.

Returns:

numpy.ndarray

A 3-element array representing the normal vector to the Moon’s orbital plane at the given time. The vector is normalized to have a unit length.

Notes:

• The function assumes a circular orbit for the Moon and uses a time step of one week (604800 seconds) to calculate the normal vector.

• The normal vector is perpendicular to the plane defined by the Moon’s position vectors at the given time and one week later.

Example usage:

>>> t = Time("2024-01-01")
>>> normal_vector = moon_normal_vector(t)
>>> normal_vector
array([-0.093,  0.014,  0.995])

ssapy.compute.moon_shine(r_moon, r_sat, r_earth, r_sun, radius, albedo, albedo_moon, albedo_back, albedo_front, area_panels)

Calculate the fractional flux of sunlight reflected from the Moon to the satellite.

This function computes the flux of sunlight reflected from the Moon to the satellite, including contributions from both the front and back surfaces of the satellite’s solar panels.

Parameters:

r_moon(n, 3) numpy.ndarray

Array of coordinates representing the position of the Moon.

r_sat(n, 3) numpy.ndarray

Array of coordinates representing the position of the satellite.

r_earth(n, 3) numpy.ndarray

Array of coordinates representing the position of the Earth.

r_sun(n, 3) numpy.ndarray

Array of coordinates representing the position of the Sun.

radiusfloat

Radius of the satellite in meters.

albedofloat

Albedo of the satellite’s surface.

albedo_moonfloat

Albedo of the Moon.

albedo_backfloat

Albedo of the back surface of the satellite’s solar panels.

albedo_frontfloat

Albedo of the front surface of the satellite’s solar panels.

area_panelsfloat

Area of the satellite’s solar panels in square meters.

Returns:

dict

Dictionary containing the flux contributions from the Moon to the satellite: - ‘moon_bus’: Fraction of light reflected off the satellite’s bus from the Moon. - ‘moon_panels’: Fraction of light reflected off the satellite’s panels from the Moon.

Notes:

• The function assumes that the solar panels are always facing the Sun and calculates flux based on the phase angles.

• Flux contributions from both the front and back surfaces of the solar panels are computed.

Example:

>>> r_moon = np.array([[1e8, 1e8, 1e8]])
>>> r_sat = np.array([[1e7, 1e7, 1e7]])
>>> r_earth = np.array([[0, 0, 0]])
>>> r_sun = np.array([[1e11, 0, 0]])
>>> moon_shine(r_moon, r_sat, r_earth, r_sun, radius=0.4, albedo=0.20, albedo_moon=0.12, albedo_back=0.50, albedo_front=0.05, area_panels=100)
{'moon_bus': array([...]), 'moon_panels': array([...])}

ssapy.compute.quickAltAz(orbit, time, observer, propagator=KeplerianPropagator(), obsAngleCorrection=None)

Quickly estimate observed altitude and azimuth of orbiting objects as viewed at specified times and locations.

This algorithm approximates “up” as pointing directly away from the center of the Earth, instead of normal to the reference ellipsoid. Use altAz if you want values wrt the reference ellipsoid.

Parameters:

• orbit (Orbit or list of Orbit (n,)) – Orbits for which to calculate direction cosines.

• time (array_like or astropy.time.Time (m,)) – If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• observer (EarthObserver) – Where on Earth to compute alt and az.

• propagator (Propagator, optional) – The propagator instance to use.

• obsAngleCorrection ({None, "linear", "exact"}, optional) – Correct actual angle to observed angle, meaning account for light-time delay, and aberration due to the observer’s velocity. None means don’t do any correction. “linear” means do an aberration correction and first order light-time correction. “exact” means do an aberration correction and iteratively solve for the exact light-time correction. (The “linear” correction is almost always sufficiently accurate).

Notes

If orbit is scalar-valued (an Orbit instead of a list of Orbit), then that dimension will be squeezed out in the return value. Likewise, if time is scalar, that dimension will be squeezed out.

For Keplerian orbit propagation it is more efficient to use a “vector Orbit” instead of a list of single scalar Orbits.

When doing light time corrections, the time argument is the arrival time of the photons at the observer, as opposed to the emission time at the satellite.

Returns:

alt, az – Altitude and azimuth in radians.

Return type:

array_like (n, m)

ssapy.compute.radec(orbit, time, obsPos=None, obsVel=None, observer=None, propagator=KeplerianPropagator(), obsAngleCorrection=None, rate=False)

Calculate observed right ascension, declination, and slant range of orbiting objects as viewed at specified times and positions.

Parameters:

• orbit (Orbit or list of Orbit (n,)) – Orbit(s) for which to calculate ra and dec.

• time (array_like or astropy.time.Time (m,)) – If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• obsPos (array_like (m, 3), optional) – Positions of observers at given times.

• obsVel (array_like (m, 3), optional) – Velocity of observers at given times.

• observer (Observer or list of Observers (m,), optional) – Observer(s) used to calculate obsPos.

• propagator (Propagator, optional) – The propagator instance to use.

• rate (bool) – If True, return additionally the time derivatives of the quantity, times cos(dec) in the case of right ascension.

• obsAngleCorrection ({None, "linear", "exact"}, optional) – Correct actual angle to observed angle, meaning account for light-time delay, and aberration due to the observer’s velocity. None means don’t do any correction. “linear” means do an aberration correction and first order light-time correction. “exact” means do an aberration correction and iteratively solve for the exact light-time correction. (The “linear” correction is almost always sufficiently accurate).

Notes

Exactly 1 of obsPos and observer must be supplied. observer and obsPos follow similar broadcasting rules as detailed below explicitly only for obsPos.

The length of time and obsPos must match or be broadcastable to match. If orbit is scalar-valued (an Orbit instead of a list of Orbit), then that dimension will be squeezed out in the return value. Likewise, if both time and obsPos are scalar, that dimension will be squeezed out.

For Keplerian orbit propagation it is more efficient to use a “vector Orbit” instead of a list of single scalar Orbits.

When doing light time corrections, the time argument is the arrival time of the photons at the observer, as opposed to the emission time at the satellite.

Returns:

• ra, dec (array_like (n, m)) – Right ascension and declination in radians.

• range (array_like (n, m)) – (Slant) range in meters.

• If rate, also

• raRate (array_like (n, m)) – Time derivatives of right ascension times cos(dec), rad / sec.

• decRate (array_like (n, m)) – Time derivative of declination, rad / sec.

• rangeRate (array_like (n, m)) – Time derivative of range, meter / s.

ssapy.compute.radecRate(orbit, time, obsPos=None, obsVel=None, observer=None, propagator=KeplerianPropagator(), obsAngleCorrection=None)

Calculate ra/dec rate and slant range rate of orbit at specified times and observer positions and velocities.

DEPRECATED. Use radec(…, rate=True) in new code.

Parameters:

• orbit (Orbit or list of Orbit (n,)) – Orbit(s) for which to calculate slant range rate.

• time (array_like or astropy.time.Time (m,)) – If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• obsPos (array_like (m, 3), optional) – Positions of observers at given times.

• obsVel (array_like (m, 3), optional) – Velocity of observers at given times.

• observer (Observer or list of Observers (m,), optional) – Observer(s) used to calculate obsPos and obsVel.

• propagator (Propagator, optional) – The propagator instance to use.

• obsAngleCorrection ({None, "linear", "exact"}, optional) – Correct actual angle to observed angle, meaning account for light-time delay, and aberration due to the observer’s velocity. None means don’t do any correction. “linear” means do an aberration correction and first order light-time correction. “exact” means do an aberration correction and iteratively solve for the exact light-time correction. (The “linear” correction is almost always sufficiently accurate).

Notes

Exactly 1 of obsPos and observer must be supplied. observer and obsPos follow similar broadcasting rules as detailed below explicitly only for obsPos. If obsPos is specified, obsVel must also be specified and congruent to obsPos.

The length of time and obsPos must match or be broadcastable to match. If orbit is scalar-valued (an Orbit instead of a list of Orbit), then that dimension will be squeezed out in the return value. Likewise, if both time and obsPos are scalar, that dimension will be squeezed out.

For Keplerian orbit propagation it is more efficient to use a “vector Orbit” instead of a list of single scalar Orbits.

When doing light time corrections, the time argument is the arrival time of the photons at the observer, as opposed to the emission time at the satellite.

Returns:

• ra (array_like (n, m)) – right ascension in radians

• raRate (array_like (n, m)) – Rate of change of right ascension*cos(dec) in radians per second.

• dec (array_link (n, m)) – declination in radians

• decRate (array_like (n, m)) – Rate of change of declination in radians per second.

• slantRange (array_like (n, m)) – Range in meters

• slantRangeRate (array_like (n, m)) – Slant range rate in meters per second.

ssapy.compute.radecRateObsToRV(ra, dec, slantRange, raRate=None, decRate=None, slantRangeRate=None, obsPos=None, obsVel=None)

Convert object angles and observer position to 3D observer position

This only does the geometric part; it ignores light travel time. This is the inverse of rvObsToRaDecRate.

If obsVel is None, then the returned velocity will also be None.

Parameters:

• ra (array_like (...)) – right ascension in radians

• dec (array_like (...)) – declination in radians

• slantRange (array_like (...)) – Range in meters

• raRate (array_like (...)) – Rate of change of right ascension*cos(dec) in radians per second.

• decRate (array_like (...)) – Rate of change of declination in radians per second.

• slantRangeRate (array_like (...)) – Slant range rate in meters per second.

• obsPos (array_like (..., 3)) – Observer position in meters

• obsVel (array_like (..., 3)) – Observer velocity in meters

Returns:

• r (array_like (…, 3)) – object position in meters

• v (array_like (…, 3)) – object velocity in meters per second observer velocity in meters per second

• v is None if obsVel is None.

ssapy.compute.refine_pass(orbit, observer, time, propagator=KeplerianPropagator(), horizon=np.float64(0.3490658503988659), maxSpan=86400.0)

Refine a satellite overhead pass.

Parameters:

• orbit (Orbit) – Orbit in question.

• observer (EarthObserver.) – Observer for which to refine pass.

• time (float or astropy.time.Time) – A time when satellite is visible to observer. (Found using find_passes, for instance) If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• propagator (Propagator, optional) – The propagator instance to use.

• horizon (float or Quantity, optional) – Minimum altitude for which to consider a satellite “visible”. If float, then should be in radians.

• maxSpan (float or Quantity, optional) – Maximum amount of time before or after time to search for rise/set times, or time of max altitude. If float, then seconds.

Returns:

Key/values are: tStart : Time

Time at which Orbit rises above horizon

tEndTime

Time at which Orbit rises above horizon

tMaxAltTime

Time Orbit passes through maximum altitude

maxAltQuantity

Maximum altitude

durationQuantity

Duration of pass.

illumAtStartbool

Is satellite illuminated at tStart?

illumAtEndbool

Is satellite illuminated at tEnd?

tTerminatorTime or None

If illumAtStart != illumAtEnd, then time satellite passes through cylindrical terminator shadow. Otherwise, None.

sunAltStartQuantity

Altitude of sun at tStart

sunAltEndQuantity

Altitude of sun at tEnd

Return type:

dict

ssapy.compute.rv(orbit, time, propagator=KeplerianPropagator())

Calculate positions and velocities on the outer product of all supplied orbits and times.

Parameters:

• orbit (Orbit or list of Orbit (n,)) – Desired orbit(s)

• time (array_like or astropy.time.Time (m,)) – If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• propagator (Propagator, optional) – The propagator instance to use.

Notes

If orbit or time is scalar valued (as opposed to a list of Orbit, e.g.) then the corresponding dimension of the output will be squeezed out.

For Keplerian orbit propagation it is more efficient to use a “vector Orbit” instead of a list of single scalar Orbits.

Returns:

• r (array_like (n, m, 3)) – Position in meters.

• v (array_like (n, m, 3)) – Velocity in meters per second.

ssapy.compute.rvObsToRaDecRate(r, v, obsPos=None, obsVel=None)

Convert object and observer position and velocity to angles.

This only does the geometric part; it ignores light travel time, which may be applied to the object location before input. Assumes that r, v, obsPos, and obsVel all have common shape.

Parameters:

• r (array_like (..., 3)) – object position in meters

• v (array_like (..., 3)) – object velocity in meters per second

• obsPos (array_like (..., 3)) – observer position in meters

• obsVel (array_like (..., 3), optional) – observer velocity in meters per second

Returns:

• ra (array_like (…)) – right ascension in radians

• dec (array_link (…)) – declination in radians

• slantRange (array_like (…)) – Range in meters

• raRate (array_like (…)) – Rate of change of right ascension*cos(dec) in radians per second.

• decRate (array_like (…)) – Rate of change of declination in radians per second.

• slantRangeRate (array_like (…)) – Slant range rate in meters per second.

ssapy.compute.sun_shine(r_sat, r_earth, r_sun, radius, albedo, albedo_front, area_panels)

Calculate the fractional flux of sunlight reflected from the Sun to the satellite.

This function computes the flux of sunlight reflected from the Sun to the satellite, including contributions from the front surface of the satellite’s solar panels.

Parameters:

r_sat(n, 3) numpy.ndarray

Array of coordinates representing the position of the satellite.

r_earth(n, 3) numpy.ndarray

Array of coordinates representing the position of the Earth.

r_sun(n, 3) numpy.ndarray

Array of coordinates representing the position of the Sun.

radiusfloat

Radius of the satellite in meters.

albedofloat

Albedo of the satellite’s surface.

albedo_frontfloat

Albedo of the front surface of the satellite’s solar panels.

area_panelsfloat

Area of the satellite’s solar panels in square meters.

Returns:

dict

Dictionary containing the flux contributions from the Sun to the satellite: - ‘sun_bus’: Fraction of light reflected off the satellite’s bus from the Sun. - ‘sun_panels’: Fraction of light reflected off the satellite’s panels from the Sun.

Notes:

• The function assumes that the solar panels are always facing the Sun and calculates flux based on the phase angle.

• Flux contributions from the front surface of the solar panels are computed.

Example:

>>> r_sat = np.array([[1e7, 1e7, 1e7]])
>>> r_earth = np.array([[0, 0, 0]])
>>> r_sun = np.array([[1e11, 0, 0]])
>>> sun_shine(r_sat, r_earth, r_sun, radius=0.4, albedo=0.20, albedo_front=0.05, area_panels=100)
{'sun_bus': array([...]), 'sun_panels': array([...])}

Constants

A collection of physical constants.

ssapy.constants.WGS84_EARTH_MU

Earth gravitational constant from WGS84 model [m³/s²]

ssapy.constants.WGS72_EARTH_MU

Earth gravitational constant from WGS72 model [m³/s²]

ssapy.constants.WGS84_EARTH_OMEGA

Earth angular velocity from WGS84 model [rad/s]

ssapy.constants.WGS72_EARTH_OMEGA

Earth angular velocity from WGS72 model [rad/s]

ssapy.constants.WGS84_EARTH_RADIUS

Earth radius at equator from WGS84 model [m]

ssapy.constants.WGS72_EARTH_RADIUS

Earth radius at equator from WGS72 model [m]

ssapy.constants.WGS84_EARTH_FLATTENING

Earth flattening f = (a - b) / a where a and b are the major and minor axes respectively, from WGS84 model [unitless]

ssapy.constants.WGS72_EARTH_FLATTENING

Earth flattening f = (a - b) / a where a and b are the major and minor axes respectively, from WGS72 model [unitless]

ssapy.constants.WGS84_EARTH_POLAR_RADIUS

Earth polar radius, calculated by multiplying the equatorial radius by (1 - flattening), from WGS84 model [m]

ssapy.constants.WGS72_EARTH_POLAR_RADIUS

Earth polar radius, calculated by multiplying the equatorial radius by (1 - flattening), from WGS72 model [m]

ssapy.constants.RGEO

GEO-synchronous radius [m]

ssapy.constants.VGEO

GEO-synchronous velocity [m/s]

ssapy.constants.RGEOALT

GEO-synchronous altitude [m]

ssapy.constants.VLEO

Approximate orbital velocity for low earth orbit (altitude of 500km) [m/s]

ssapy.constants.LD

Lunar semi-major axis [m]

Gravitational constants

ssapy.constants.SUN_MU

Sun gravitational constant, from IAU 1976 model [m³/s²]

ssapy.constants.MOON_MU

Moon gravitational constant, from the DE200 ephemeris [m³/s²]

ssapy.constants.MERCURY_MU

Mercury gravitational constant [m³/s²]

ssapy.constants.VENUS_MU

Venus gravitational constant [m³/s²]

ssapy.constants.EARTH_MU

Earth gravitational constant, from WGS84 model [m³/s²]

ssapy.constants.MARS_MU

Mars gravitational constant [m³/s²]

ssapy.constants.JUPITER_MU

Jupiter gravitational constant [m³/s²]

ssapy.constants.SATURN_MU

Saturn gravitational constant [m³/s²]

ssapy.constants.URANUS_MU

Uranus gravitational constant [m³/s²]

ssapy.constants.NEPTUNE_MU

Neptune gravitational constant [m³/s²]

Mass

ssapy.constants.SUN_MASS

Sun mass, from the DE405 ephemeris [kg]

ssapy.constants.MOON_MASS

Moon mass, from the DE405 ephemeris [kg]

ssapy.constants.MERCURY_MASS

Mercury mass, from the DE405 ephemeris [kg]

ssapy.constants.VENUS_MASS

Venus mass, from the DE405 ephemeris [kg]

ssapy.constants.EARTH_MASS

Earth mass, from the DE405 ephemeris [kg]

ssapy.constants.MARS_MASS

Mars mass, from the DE405 ephemeris [kg]

ssapy.constants.JUPITER_MASS

Jupiter mass, from the DE405 ephemeris [kg]

ssapy.constants.SATURN_MASS

Saturn mass, from the DE405 ephemeris [kg]

ssapy.constants.URANUS_MASS

Uranus mass, from the DE405 ephemeris [kg]

ssapy.constants.NEPTUNE_MASS

Neptune mass, from the DE405 ephemeris [kg]

Radius

ssapy.constants.MOON_RADIUS

Moon radius (source: 10.2138/rmg.2006.60.3) [m]

ssapy.constants.MERCURY_RADIUS

Mercury radius [m]

ssapy.constants.VENUS_RADIUS

Venus radius [m]

ssapy.constants.EARTH_RADIUS

Earth radius [m]

ssapy.constants.MARS_RADIUS

Mars radius [m]

ssapy.constants.JUPITER_RADIUS

Jupiter radius [m]

ssapy.constants.SATURN_RADIUS

Saturn radius [m]

ssapy.constants.URANUS_RADIUS

Uranus radius [m]

ssapy.constants.NEPTUNE_RADIUS

Neptune radius [m]

Correlate Tracks

class ssapy.correlate_tracks.CircVelocityPrior(sigma=np.float64(0.19804205158855578))

Gaussian prior that log(GM/r/v^2) = 0.

Parameters:

sigma (float) – Prior standard deviation on log(GM/r/v^2), dimensionless.

sigma

__call__(orbit)

Returns log prior probability of given orbit.

class ssapy.correlate_tracks.GaussPrior(mu, cinvcholfac, translator)

Gaussian prior on parameters.

Parameters:

• mu (array_like (6)) – Mean of prior

• cinvcholfac (array_like (6, 6)) – Inverse Cholesky factor for covariance matrix, chi = cinvcholfac*(param-mu)

• translator (function(orbit) -> param) – function that translates from orbits to parameters

mu, cinvcholfac, tranlator

__call__(orbit)

Returns log prior probability of given orbit.

class ssapy.correlate_tracks.Hypothesis(tracks, nsat=1000)

Assignment of detections to multiple tracks.

This class represents an assignment of detections to multiple tracks. The goal of Multiple Hypothesis Tracking to to find the Hypothesis (track assignment) that has the highest likelihood.

Parameters:

tracks (list of Tracks) – The tracks in the hypothesis. Each track must correspond to different observations, and all observations under consideration must be accounted for

tracks

the tracks in this hypothesis

Type:

list

lnprob

the log probability of this hypothesis

Type:

float

nsat

the total number of satellites in the sky that could be observed

Type:

int

static addto(hypothesis, track, oldtrack=None, **kw)

Add a new track to a Hypothesis.

Parameters:

• hypothesis (Hypothesis) – hypothesis to which to add

• track (Track) – track to add to hypothesis

• oldtrack (Track) – track to remove from hypothesis, if the new track updates the old track.

Returns:

• new Hypothesis, with track added, replacing oldtrack if oldtrack is

• not None.

difference(hypothesis)

Print the difference between two hypotheses.

Often two hypotheses are very similar. This makes it easier to inspect only the differences between them.

Parameters:

hypothesis (Hypothesis) – hypothesis to which to compare this hypothesis

ntracklet()

Number of observations in tracks in the hypothesis.

summarize(verbose=False)

Summarize the Hypothesis as a string.

Parameters:

verbose (bool) – Include information about each track in the hypothesis.

Return type:

A string summarizing the hypothesis.

class ssapy.correlate_tracks.MHT(data, nsat=1000, truth=None, hypotheses=None, propagator=None, mode='rv', fitonly=None, approximate=False, orbitattr=None, priors=None)

Multiple Hypothesis Tracking of many hypotheses explaining data.

This class manages the assessment of which hypothesis are most likely to explain data, building a set of tracks observation by observation.

Initialize the MHT with the data, do mht.run(), and inspect the most likely hypotheses using [h.lnprob for h in mht.hypotheses].

Parameters:

• data (array_like[N]) – data to model. Must contain angles of observations, observer locations and velocities, times, detection IDs, etc.

• nsat (int) – total number of satellites in the sky. Affects the overall prior and the preference for new tracks vs. additional assignments to existing tracks.

• truth (dict) – when true assignments are known, this argument can provide debug information about when the tracker has lost the true assignment

• hypotheses (list of Hypothesis) – initialize the MHT with this existing list of Hypotheses. Default None.

• propagator (instance of ssa Propagator) – propagator to use. Default to None (Keplerian)

• fitonly (ndarray[N] bool) – entries in data to fit

add_tracklet(satid, **kw)

Add an observation from data to the MHT analysis.

Parameters:

satid (int) – the ID of the observation to add.

static flag_inconsistency(track2hyp, hyp)

Debug method checking to see if there are any inconsistencies between the tracks tracked by the MHT and the hypotheses tracked by the MHT. Should only be needed for debugging.

Every hypothesis in track2hyp[track] should include the Track track. Every hypothesis in hyp should be in the list of track2hyp[track] for each of its tracks. Every track should be in a hypothesis.

Parameters:

• track2hyp (dict[Track -> list(Hypothesis)]) –

• hyp (list(Hypothesis)) –

Return type:

True if everything is consistent.

prune(satid, nkeepmax=10000, pkeep=1e-09, keeponlytrue=False, nconfirm=6)

Prune unlikely hypotheses from the MHT.

Parameters:

• satid (int) – the detection ID most recently added to the MHT

• nkeepmax (int) – keep no more than nkeepmax hypotheses

• pkeep (float) – keep no hypotheses more than pkeep times less likely than the most likely hypothesis

• keeponlytrue (bool) – keep only the single hypothesis known to be true for speed comparison purposes.

• nconfirm (int) – only keep one variant of tracks that agree in the last nconfirm detections

prune_stale_hypotheses(newdeadtracks)

Prune hypotheses that are different from better hypotheses only in dead tracks.

Note: implentation relies on identical tracks always having the same id. The MHT code tries to guarantee this—otherwise it must do extra work matching identical tracks to new detections, etc. But using id(Track) feels horrible to me.

Parameters:

newdeadtracks (list of Track) – list of tracks that died since the last pruning.

Return type:

boolean mask indicating hypotheses to keep

prune_tracks(satid, nconfirm=6)

Prune tracks that are now nearly identical from the MHT analysis.

As the MHT adds observations to tracks, eventually some tracks have a number of confirming observations. For these long tracks, we really only need to keep different Tracks and Hypotheses for cases where there are recent disagreements with how those tracks continue; if hypotheses agree that the track continues in the same way, we don’t need to continue tracking the past differences. So we prune out cases of former disagreement.

This tries to identify such overlapping tracks and keep only one of them.

Parameters:

• satid (int) – the satID of the most-recently added detection to the MHT

• nconfirm (int) – if two tracks are identical in their last nconfirm observations delete the tracks not in the most likely hypothesis.

Return type:

boolean array indicating hypotheses to keep

run(first=None, last=None, verbose=False, order='forward', **kw)

Run the MHT analysis.

Parameters:

• first (int) – First tracklet to consider

• last (int) – Last tracklet to consider

• verbose (bool) – print extra progress information while running

• **kw (dict) – extra keyword arguments to pass to MHT.prune.

class ssapy.correlate_tracks.Track(satIDs, data, guess=None, initial_lnprob=0, priors=None, mode='rv', propagator=None, orbitattr=None)

Set of observations to be fit as an object moving through space.

Subclasses TrackBase, which could have other implementations (for example, a subclass that doesn’t fit the whole track simultaneously, but instead approximates past information with a single Gaussian prior on the parameters).

Parameters:

• satIDs (see TrackBase) –

• data (see TrackBase) –

• guess (array_like, float) – guess of initial parameters fitting orbit

• initial_lnprob (float) – adjust probability of this track by initial_lnprob (default: 0)

• priors (rvsampler.Prior instance) – list of priors to use when fitting this track

• mode (str) – fitting mode to use. One of ‘rv’, ‘angle’, or ‘equinoctial’ parameters have different meanings in these modes

• propagator (Propagator instance) – ssa Propagator instance to use for propagation. Default None (Keplerian)

addto(satid)

Add a detection to this track.

Parameters:

satID (int or list(int)) – the detection ID of the observation to add to this track

Return type:

A new Track, including the additional observation.

gaussian_approximation(propagator=None)

Get Gaussian approximation of this track.

Return type:

TrackGauss, giving Gaussian approximation of this track.

class ssapy.correlate_tracks.TrackBase(satIDs, data, volume=None, mode='rv', propagator=None, orbitattr=None)

Set of observations fit together as an orbital track.

This class represents a set of observations that are being modeled as the motion of a single object through space.

Parameters:

• satIDs (array_like) – list of detection IDs considered part of this track

• data (array_like) – observations that could be part of this track. Observations that are actually part of this track are specified in satIDs. data is must contain a number of fields used in fitting, like ra, dec, obsPos, obsVel, and time.

• volume (float) – the volume in 6-D orbital element space in which objects could be found, needed to normalize false-positive vs. new satellite odds default to a very rough number only applicable for RV parameterization

• propagator (ssa Propagator instance) – the propagator to use when propagating the orbit. Default to None (Keplerian)

satIDs

the detection IDs that are part of this Track

Type:

array_like, int

data

observations that could be part of this track; see constructor

Type:

array_like

volume

the prior volume in 6-D parameter where satellites may be observed

Type:

float

times

the GPS times at which this object was observed

Type:

array_like, float

mode

the type of parameterization to use to fit the track one of ‘rv’, ‘angle’, ‘equinoctial’

Type:

str

lnprob

the log probability of this association of observations to an orbit contains contributions both from priors and likelihood.

Type:

float

gate(arc, return_nwrap=False)

Compute chi^2 that this track could be associated with arc.

Parameters:

• arc (array_like, structure array) – must contain fields: time: future time (GPS) (float) rStation_GCRF: GCRF position of observer at future time (3, float) vStation_GCRF: GCRF velocity of observer at future time (3, float) ra: observed ra (float) pmra: observed proper motion in ra (float) dec: observed dec (float) pmdec: observed proper motion in dec (float) dra: uncertainty in ra (float) ddec: uncertainty in dec (float) dpmra: uncertainty in proper motion in ra (float) dpmdec uncertainty in proper motion in dec (float)

• return_nwrap (bool) – if True, also return the standard deviation of the number of orbits the object could have made at the time of observation.

Returns:

• chi2, nwrapsig

• chi2 (float) – The chi-squared implied by associating this observation with this track.

• nwrapsig (float) – if return_nwrap, the standard deviation of the number of orbits the object could have made at the time of observation

predict(arc0, return_sigma=False, return_nwrap=True)

Predict the location and uncertainty of this object at future times.

Parameters:

• arc0 (array_like, structure array) – must contain fields: time: future times (GPS) (float) rStation_GCRF: GCRF position of observer at future time (3, float) vStation_GCRF: GCRF velocity of observer at future time (3, float)

• return_sigma (bool) – if True, also return the sigma points used to compute the future covariance

• return_nwrap (bool) – if True, also return the number of orbital wrappings of the sigma points.

Returns:

• mean, covar, fsigma

• mean (array_like (4, n)) – ra, dec, dra/dt*cos(dec), ddec/dt at future times (rad, rad/s)

• covar (array_like (4, 4, n)) – covariance of mean parameters at future times

• fsigma (array_like (4, 13, n)) – 13 sigma points used to compute covar fsigma[:, 0, :] is identical to mean.

propagaterdz(param, arc0=None, return_nwrap=False)

Propagate param to times and locations specified in arc0.

Parameters:

• param (array_like) – parameters corresponding to self.mode for orbit to propagate

• arc0 (arary_like, must contain observation fields) – array containing fields specifying times and observer locations to which track fit should be propagated.

• return_nwrap (bool) – if True, additionally return nwrap

Returns:

• array([rr, pmrr, dd, pmdd, nwrap])

• rr (ra at future times, radians)

• dd (dec at future times, radians)

• pmrr (proper motion in ra at future times, rad / s)

• pmdd (proper motion in dec at future times, rad / s)

• nwrap (number of orbits completed. Only provided if return_nwrap) – is True.

update(time, rStation=None, vStation=None)

Shift this track to a new time.

Only implemented for numerical propagators where time-shifting can provide a big speed up.

class ssapy.correlate_tracks.TrackGauss(satIDs, data, param, covar, chi2, mode='rv', propagator=None, sigma_points=None, priors=None, orbitattr=None)

Set of observations to be fit as an object moving through space.

Subclasses TrackBase. This implementation does not fit the entire track simultaneously, but instead approximates past tracklets with a single Gaussian prior on the parameters. This information is then updated with each new observation in a Kalman-filter approach.

Parameters:

• satIDs (see TrackBase) –

• data (see TrackBase) –

• param (array_like, float) – mean of Gaussian

• covar (array_like, float) – covariance of Gaussian

• chi2 (float) – chi2 of fit at center of Gaussian

• mode (str) – fitting mode to use. One of ‘rv’, ‘angle’, or ‘equinoctial’ parameters have different meanings in these modes

• propagator (Propagator instance) – ssa Propagator instance to use for propagation. Default None (Keplerian)

• priors (list of Priors instances) – priors used in initially fitting this orbit. Their influence should now be completely recorded in param + covar, but they are recorded here for reference.

addto(satid)

Add a detection to this track.

This also shifts the track’s epoch to the epoch of the added observation.

Parameters:

satID (int) – the detection ID of the observation to add to this track

Return type:

A new TrackGauss, including the additional observation.

at(t, rStation=None, vStation=None)

Get Gaussian approximation of this track at different epoch.

This just copies the track and then updates it to be at the new time.

Parameters:

• t (astropy.time.Time) – epoch for new TrackGauss

• rStation (array_like (3)) – position of station at new time (m, GCRF); required if mode == ‘angle’

• vStation (array_like (3)) – velocity of station at new time (m, GCRF); required if mode == ‘angle’

Return type:

Gaussian-approximated track (TrackGauss) at new time.

gaussian_approximation(propagator)

Get Gaussian approximation of this track.

Return type:

self. No op; this track already is a Gaussian approximation.

update(t, rStation=None, vStation=None)

Move this track to a different epoch.

Does not bother updating if t is different from the current track time by less than one microsecond.

Parameters:

• t (astropy.time.Time) – epoch for new TrackGauss

• rStation (array_like (3)) – position of station at new time (m, GCRF); required if mode == ‘angle’

• vStation (array_like (3)) – velocity of station at new time (m, GCRF); required if mode == ‘angle’

class ssapy.correlate_tracks.VolumeDistancePrior(scale=np.float64(42164172.36450565))

Prior on range like r^2*exp(-r/scale), preferring larger distances where there is more volume.

__call__(orbit)

Returns log prior probability of given orbit.

class ssapy.correlate_tracks.ZeroRadialVelocityPrior(sigma=0.25)

Gaussian prior that v_R = 0.

Parameters:

sigma (float) – Prior standard deviation on np.dot(v, r)/|v|/|r|, dimensionless.

sigma

__call__(orbit)

Returns log prior probability of given orbit.

ssapy.correlate_tracks.combinatoric_lnprior(nsat, ntrack, ndet)

Prior probability of assigning detections to tracks of satellites.

This function intends to give the prior probability that ndet detections of nsat total satellites correspond to ntrack satellites.

Parameters:

• nsat (int) – number of satellites that could possibly be observed

• ntrack (int) – number of tracks in current solution

• ndet (int) – number of detections of these satellites

ssapy.correlate_tracks.data_for_satellite(data, satIDList)

Extract data for specific satID.

This helper routine takes a set of observations containing a column ‘satID’, and returns the subset of those observations for which ‘satID’ is in the list satIDList.

Parameters:

• data (array_like (n), including column 'satID') – data from which to select matching rows

• satIDList (list) – list of satIDs to select from data

ssapy.correlate_tracks.fit_arc(arc, guess, verbose=False, propagator=KeplerianPropagator(), mode='rv', priors=None, damp=-1, optimizerkw={}, lsq=True, orbitattr=None, **kw)

Fit an orbit to an arc.

See documentation for fit_arc_blind. fit_arc implements the same interface, but takes a additional guess & epoch arguments. These specify the initial guess in appropriate units given the chosen mode. epoch specifies the time at which this initial guess applies.

Parameters:

• arc (array_like (n)) – numpy ndarray containing several fields necessary for describing observations. rStation_GCRF (m), vStation_GCRF (m), time (astropy.time.Time or gps seconds), ra, dec, pmra, pmdec are all required fields.

• guess (array_like (n_par)) – parameters for initial guess guess[nparam] should be the epoch in GPS guess[nparam:] should be extra parameters not optimized over

• verbose (bool) – True for more verbose output

• propagator (ssapy.propagator.Propagator instance) – propagator to use in fitting

• mode (str) – Mode for fitting, must be one of ‘rv’, ‘equinoctial’, or ‘angle’. In the first case, the parameters defining the orbit are taken to be the position and velocity of the object. In the second case, the parameters are taken to be the equinoctial elements. In the third case, the parameters are taken to be the angular position and proper motion of the object, together with its line-of-sight range and velocity.

• priors (list of objects implementing ssa Prior interface) – any priors to be applied in fitting

• damp (float) – Damping used in pseudo-Huber loss function. Default of -1 means no damping.

• orbitattr (list(str)) – Names of orbit attributes used for propagation (mass, area, …)

• optimizerkw (dict) – any extra keywords to pass to the optimizer

• **kw (dict) – any extra keywords to pass to optimizer.optimize()

ssapy.correlate_tracks.fit_arc_blind(arc, verbose=False, mode='rv', priors=None, propagator=None, damp=-1, orbitattr=None, optimizerkw={}, lsq=True, factor=2, **kw)

Fit an orbit to an arc.

The arc can be arbitrarily long. Orbits are initially fit to the initial two observations and extended using a geomtrically increasing time baseline. If no new observations are inside of a given baseline, the next obs is added. The fit is least-squares, Levenberg-Marquardt, given specified priors and propagator. If the factor is set to 1, the orbit is built out iteratively by satID.

Parameters:

• arc (array_like (n)) – numpy ndarray containing several fields necessary for describing observations. rStation_GCRF (m), vStation_GCRF (m), time (astropy.time.Time or gps seconds), ra, dec, pmra, pmdec are all required fields.

• verbose (bool) – True for more verbose output

• mode (str) – Mode for fitting, must be one of ‘rv’, ‘equinoctial’, or ‘angle’. In the first case, the parameters defining the orbit are taken to be the position and velocity of the object. In the second case, the parameters are taken to be the equinoctial elements. In the third case, the parameters are taken to be the angular position and proper motion of the object, together with its line-of-sight range and velocity.

• priors (list of objects implementing ssa Prior interface) – any priors to be applied in fitting

• propagator (ssapy.propagator.Propagator instance) – propagator to use in fitting

• damp (float) – damp chi residuals according to pseudo-Huber loss function. Default of -1 means do not damp

• orbitattr (list(str)) – names of additional orbit propagation keywords to guess

• optimizerkw (dict) – any extra keywords to pass to the optimizer

• factor (float) – factor for geometrically increasing time baseline (default 2)

• **kw (dict) – any extra keywords to pass to optimizer.optimize()

ssapy.correlate_tracks.fit_arc_blind_via_track(arc, propagator=None, verbose=False, reset_if_too_uncertain=False, mode='rv', priors=None, order='forward', approximate=False, orbitattr=None, factor=1)

Fit a series of detections using the Track interface.

Parameters:

• arc (ndarry with many fields) – The observations to fit. Must contain many fields (ra, dec, …)

• propagator (Propagator instance) – Propagator to use for orbit propagation. Default None (Keplerian).

• verbose (bool) – print additional information during fitting if True

• reset_if_too_uncertain (bool) – if the uncertainty in the number of wraps grows to be greater than pi, start track fitting over, ignoring old tracklets.

• factor (float) – factor for geometrically increasing time baseline (default 1)

• order (str) – order to process observations; default ‘forward’. ‘backward’ processes observations in reverse chronological order.

Returns:

• tracklist

• A list of Tracks, giving the fits to the track as each observation

• in added sequentially in time to the fit.

ssapy.correlate_tracks.fit_arc_with_gaussian_prior(arc, mu, cinvcholfac, verbose=False, propagator=KeplerianPropagator(), mode='rv', lsq=True, optimizerkw={}, orbitattr=None, **kw)

Fit an orbit to an arc.

See documentation for fit_arc_blind. fit_arc_with_gaussian_prior implements the same interface, but takes an additional mu, cinvcholfac, and epoch arguments. These specify a gaussian prior on the parameters. The mean and covariance of the Gaussian are given by mu and cinvcholfac, the inverse Cholesky matrix corresponding to the covariance matrix. The initial guess is taken to be the mean of the Gaussian prior. Presently, any other priors are assumed to have been folded into the Gaussian. epoch specifies the time at which Gaussian prior parameters are applicable. For greatest speed, the arc should contain a single observation and the epoch should be equal to the time of that observation; then no propagations are needed in this function.

Parameters:

• arc (array_like (n)) – numpy ndarray containing several fields necessary for describing observations. rStation_GCRF (m), vStation_GCRF (m), time (astropy.time.Time or gps seconds), ra, dec, pmra, pmdec are all required fields.

• mu (array_like (n_par)) – parameters for initial guess

• cinvcholfac (array_like (n_par, n_par)) – inverse cholesky matrix for covariance matrix

• verbose (bool) – True for more verbose output

• propagator (ssapy.propagator.Propagator instance) – propagator to use in fitting

• mode (str) – Mode for fitting, must be one of ‘rv’, ‘equinoctial’, or ‘angle’. In the first case, the parameters defining the orbit are taken to be the position and velocity of the object. In the second case, the parameters are taken to be the equinoctial elements. In the third case, the parameters are taken to be the angular position and proper motion of the object, together with its line-of-sight range and velocity.

• orbitattr (list(str)) – list of strings of names of additional propagation attributes to fit (mass, area, cr, cd, …)

• optimizerkw (dict) – any extra keywords to pass to the optimizer

• **kw (dict) – any extra keywords to pass to optimizer.optimize()

ssapy.correlate_tracks.orbit_to_param(orbit, mode='rv', orbitattr=None, rStation=None, vStation=None, fitonly=False)

Convert Orbits to parameters.

Note: Orbits are converted so that they are all bound orbits to avoid challenges with hyperbolic orbit propagation.

Parameters:

• param (array_like (n, n_param)) – parameters for orbits

• mode (str) – parameter type, one of ‘rv’, ‘angle’, ‘orbit’

ssapy.correlate_tracks.param_to_orbit(param, mode='rv', orbitattr=None)

Convert parameters to Orbits.

Note: Orbits are converted so that they are all bound orbits to avoid challenges with hyperbolic orbit propagation.

Parameters:

• param (array_like (n, n_param)) – parameters for orbits

• mode (str) – parameter type, one of ‘rv’, ‘angle’, ‘orbit’

ssapy.correlate_tracks.summarize_tracklet(arc)

Compute mean ra/dec, uncertainty in mean, proper motion, uncertainty for a short arc.

Parameters:

arc (the short arc of detections) –

Returns:

• [(ra, dec), (dra, ddec), (pmra, pmdec), (dpmra, dpmdec)]

• ra, dec (the mean position in right ascension and declination)

• dra, ddec (the uncertainty in the mean RA and declination)

• pmra, pmdec (the proper motion in right ascension and declination)

• dpmra, dpmdec (the uncertainty in the proper motion in RA and declination)

ssapy.correlate_tracks.summarize_tracklets(data, posuncfloor=None, pmuncfloor=None)

Add fields to default data set to accommodate MHT fitting.

Parameters:

• data (array_like) – data to add fields to

• posuncfloor (float or None) – add a positional uncertainty floor to the ra/dec uncertainties (deg)

• pmuncfloor (float or None) – add a proper motion uncertainty floor to the ra/dec proper motion uncertainties (deg/s)

Returns:

• array parallel to data with additional fields dra, ddec,

• pmra, pmdec, dpmra, dpmdec

ssapy.correlate_tracks.time_ordered_satIDs(data, with_time=False, order='forward')

Give satIDs in data in order of observation. Optionally includes first times for each ID.

Parameters:

• data (array_like) – the data to consider

• with_time (bool) – Option to also output list of corresponding gps times

• order ('forward' or 'backward') –

Returns:

• list of satIDs in data, ordered by observation time

• list of corresponding gps times of first obs per ID

Ellipsoid

Class to handle transformations between ECEF x,y,z coords and geodetic longitude, latitude, and height.

Technically, only handles a one-axis ellipsoid, defined via a flattening parameter f, but that’s good enough for simple Earth models.

Gravity

Classes for gravity-related accelerations.

class ssapy.gravity.AccelHarmonic(body, n_max=None, m_max=None)

Acceleration due to a harmonic potential.

Parameters:

• body (ssapy.Body) – The body.

• n_max (int) – The maximum degree of the potential.

• m_max (int) – The maximum order of the potential.

class ssapy.gravity.AccelThirdBody(body)

Acceleration due to a third body.

Parameters:

body (ssapy.Body) – The third body.

class ssapy.gravity.HarmonicCoefficients

Class to hold coefficients for a spherical harmonic expansion of a gravitational potential.

classmethod fromEGM(filename)

Construct a HarmonicCoefficients object from a .egm/.cof file pair as available from https://geographiclib.sourceforge.io/html/gravity.html

SSAPy comes with four models for Earth gravity of varying degrees and orders:

name | degree | order

WGS84 | 20 | 0 EGM84 | 180 | 180 EGM96 | 360 | 360 EGM2008 | 2190 | 2159

Note that many coefficients with large degree and/or order underflow double precision floats when denormalized for use in computing accelerations and are therefore effectively ignored.

param filename:

Either the name of the .egm file or one of the names listed above.

type filename:

str

classmethod fromTAB(filename, n_max=40, m_max=40)

Construct a HarmonicCoefficients object from a .tab file as available from https://pgda.gsfc.nasa.gov/products/50

Io

Collection of functions to read from and write to various file formats.

ssapy.io.append_csv(file_names, save_path='combined_data.csv', sep=None, dtypes=False, progress=None)

Appends multiple CSV files into a single CSV file.

Parameters:

• file_names (list) – A list of CSV file names.

• save_path (str) – The path to the output CSV file. If not specified, the output will be saved to the current working directory.

• sep (str) – The delimiter used in the CSV files.

• dtypes (dict) – A dictionary specifying data types for columns.

Returns:

None

ssapy.io.append_csv_on_disk(csv_files, output_file)

Append multiple CSV files into a single CSV file.

This function merges multiple CSV files into one output CSV file. The output file will contain the header row from the first CSV file and data rows from all input CSV files.

Parameters:

csv_fileslist of str

List of file paths to the CSV files to be merged. All CSV files should have the same delimiter and structure.

output_filestr

Path to the output CSV file where the merged data will be written.

Notes:

• The function assumes all input CSV files have the same delimiter. It determines the delimiter from the first CSV file using the guess_csv_delimiter function.

• Only the header row from the first CSV file is included in the output file. Headers from subsequent files are ignored.

• This function overwrites the output file if it already exists.

Example:

>>> csv_files = ['file1.csv', 'file2.csv', 'file3.csv']
>>> output_file = 'merged_output.csv'
>>> append_csv_on_disk(csv_files, output_file)
Completed appending of: merged_output.csv.

Dependencies:

• guess_csv_delimiter function: A utility function used to guess the delimiter of the CSV files.

• csv module: Standard library module used for reading and writing CSV files.

ssapy.io.append_dict_to_csv(file_name, data_dict, delimiter='\t')

Append data from a dictionary to a CSV file.

This function appends rows of data to a CSV file, where each key-value pair in the dictionary represents a column. If the CSV file does not already exist, it creates the file and writes the header row using the dictionary keys.

Parameters:

file_namestr

Path to the CSV file where data will be appended.

data_dictdict

Dictionary where keys are column headers and values are lists of data to be written to the CSV file. All lists should be of the same length.

delimiterstr, optional

The delimiter used in the CSV file (default is tab ` `).

Notes:

• The function assumes that all lists in the dictionary data_dict have the same length.

• If the CSV file already exists, only the data rows are appended. If it doesn’t exist, a new file is created with the header row based on the dictionary keys.

• The delimiter parameter allows specifying the delimiter used in the CSV file. Common values are , for commas and ` ` for tabs.

Example:

>>> data_dict = {
>>>     'Name': ['Alice', 'Bob', 'Charlie'],
>>>     'Age': [25, 30, 35],
>>>     'City': ['New York', 'Los Angeles', 'Chicago']
>>> }
>>> append_dict_to_csv('people.csv', data_dict, delimiter=',')
This will append data to 'people.csv', creating it if it does not exist, with columns 'Name', 'Age', 'City'.

Dependencies:

• os.path.exists: Used to check if the file already exists.

• csv: Standard library module used for reading and writing CSV files.

ssapy.io.append_h5(filename, pathname, append_data)

Append data to key in HDF5 file.

Parameters:

• filename (str) – The filename of the HDF5 file.

• pathname (str) – The path to the key in the HDF5 file.

• append_data (any) – The data to be appended.

Returns:

None

ssapy.io.b3obs2pos(b3line)

Return an SGP4 Satellite imported from B3OBS data.

Intended to mimic the sgp4 twoline2rv function

Data format reference: http://help.agi.com/odtk/index.html?page=source%2Fod%2FODObjectsMeasurementFormatsandMeasurementTypes.htm

ssapy.io.combine_h5(filename, files, verbose=False, overwrite=False)

Combine multiple HDF5 files into a single HDF5 file.

This function reads datasets from a list of HDF5 files and writes them to a specified output HDF5 file. If overwrite is True, it will remove any existing file at the specified filename before combining the files. The verbose parameter, if set to True, will display progress bars during the process.

Parameters:

filenamestr

The path to the output HDF5 file where the combined datasets will be stored.

fileslist of str

A list of paths to the HDF5 files to be combined.

verbosebool, optional

If True, progress bars will be displayed for the file and key processing. Default is False.

overwritebool, optional

If True, any existing file at filename will be removed before writing the new combined file. Default is False.

Returns:

None

The function performs file operations and does not return any value.

Examples:

>>> combine_h5('combined.h5', ['file1.h5', 'file2.h5'], verbose=True, overwrite=True)

ssapy.io.exists(pathname)

Check if a file or directory exists.

Parameters:

pathnamestr

The path to the file or directory.

Returns:

bool

True if the path exists as either a file or a directory, False otherwise.

ssapy.io.exists_in_csv(csv_file, column, number, sep='\t')

Checks if a number exists in a specific column of a CSV file.

This function reads a specified column from a CSV file and checks if a given number is present in that column.

Parameters:

csv_filestr

Path to the CSV file.

columnstr or int

The column to search in.

numberint or float

The number to check for existence in the column.

sepstr, default=’ ‘

Delimiter used in the CSV file.

Returns:

bool

True if the number exists in the column, False otherwise.

ssapy.io.file_exists_extension_agnostic(filename)

Check if a file with the given name and any extension exists.

Parameters:

filenamestr

The name of the file to check, without extension.

Returns:

bool

True if a file with the given name and any extension exists, False otherwise.

ssapy.io.get_tel_pos_itrf_to_gcrs(time, tel_label='511')

Convert telescope locations in ITRF (i.e., fixed to the earth) to GCRS (i.e., geocentric celestial frame)

Parameters:

time (astropy.time.Time) – Time at which to evaluate the position

ssapy.io.guess_csv_delimiter(file_name)

Guess the delimiter used in a CSV file.

Parameters:

file_name (str) – The path to the CSV file.

Returns:

Guessed delimiter (one of ‘,’, ‘ ‘, ‘;’)

Return type:

str

ssapy.io.h5_key_exists(filename, key)

Checks if a key exists in an HDF5 file.

Parameters:

• filename (str) – The filename of the HDF5 file.

• key (str) – The key to check.

Returns:

True if the key exists, False otherwise.

ssapy.io.h5_keys(file_path)

List all groups in HDF5 file.

Parameters:

file_path (str) – The file_path of the HDF5 file.

Returns:

A list of group keys in the HDF5 file.

ssapy.io.h5_root_keys(file_path)

Retrieve the keys in the root group of an HDF5 file.

This function opens an HDF5 file and returns a list of keys (dataset or group names) located in the root group of the file.

Parameters:

file_pathstr

The path to the HDF5 file from which the root group keys are to be retrieved.

Returns:

list of str

A list of keys in the root group of the HDF5 file. These keys represent the names of datasets or groups present at the root level of the file.

ssapy.io.listdir(dir_path='*', files_only=False, exclude=None, sorted=False, index=0)

Lists files and directories in a specified path with optional filtering and sorting.

Parameters:

dir_pathstr, default=’*’

The directory path or pattern to match files and directories.

files_onlybool, default=False

If True, only returns files, excluding directories.

excludestr or None, optional

If specified, excludes files and directories whose base name contains this string.

sortedbool, default=False

If True, sorts the resulting list by numeric values in filenames.

indexint, default=0

sorted required to be true. Index of the digit used for sorting.

Returns:

list

A list of file or directory paths based on the specified filters and sorting.

ssapy.io.load_b3obs_file(file_name)

Convenience function to load all entries in a B3OBS file

ssapy.io.load_np(filename_)

Load a NumPy array from a binary file.

This function loads a NumPy array from a file in .npy format. If the file cannot be read, it handles common exceptions and prints an error message. If loading fails, it returns an empty list.

Parameters:

filename_str

The path to the file from which the NumPy array will be loaded.

Returns:

numpy.ndarray or list

The loaded NumPy array. If an error occurs during loading, returns an empty list.

Examples:

>>> arr = load_np('array.npy')
>>> print(arr)
[1 2 3 4 5]

ssapy.io.make_tle(a, e, i, pa, raan, trueAnomaly, t)

Create a TLE from Kozai mean orbital elements

Parameters:

• a (float) – Kozai mean semi-major axis in meters

• e (float) – Kozai mean eccentricity

• i (float) – Kozai mean inclination in radians

• pa (float) – Kozai mean periapsis argument in radians

• raan (float) – Kozai mean right ascension of the ascending node in radians

• trueAnomaly (float) – Kozai mean true anomaly in radians

• t (float or astropy.time.Time) – If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

Notes

Dynamic TLE terms, including the drag coefficient and ballistic coefficient, are ignored in this function.

ssapy.io.makedf(df)

Convert an input into a pandas DataFrame.

This function takes an input which can be a list or a dictionary and converts it into a pandas DataFrame. If the input is already a DataFrame, it returns it unchanged.

Parameters:

dflist, dict, or pd.DataFrame

The input data to be converted into a DataFrame. This can be a list or dictionary to be transformed into a DataFrame, or an existing DataFrame which will be returned as is.

Returns:

pd.DataFrame

A DataFrame created from the input data if the input is a list or dictionary. If the input is already a DataFrame, the original DataFrame is returned unchanged.

ssapy.io.mkdir(pathname)

Creates a directory if it does not exist.

Parameters:

pathnamestr

The path to the directory to be created.

ssapy.io.overwrite_h5(filename, pathname, new_data)

Overwrite key in HDF5 file.

Parameters:

• filename (str) – The filename of the HDF5 file.

• pathname (str) – The path to the key in the HDF5 file.

• new_data (any) – The data to be overwritten.

Returns:

None

ssapy.io.parseB3(filename)

Load data from a B3 observation file

Parameters:

filename (string) – Name of the B3 obs file to load

Returns:

A catalog of observations

Return type:

astropy.table.Table

Note that angles and positions are output in TEME frame.

ssapy.io.parseB3Line(line)

Read one line of a B3 file and parse into distinct catalog components

ssapy.io.parse_tle(tle)

Parse a TLE returning Kozai mean orbital elements.

Parameters:

tle (2-tuple of str) – Line1 and Line2 of TLE as strings

Returns:

• a (float) – Kozai mean semi-major axis in meters

• e (float) – Kozai mean eccentricity

• i (float) – Kozai mean inclination in radians

• pa (float) – Kozai mean periapsis argument in radians

• raan (float) – Kozai mean right ascension of the ascending node in radians

• trueAnomaly (float) – Kozai mean true anomaly in radians

• t (float) – GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

Notes

Dynamic TLE terms, including the drag coefficient and ballistic coefficient, are ignored in this function.

ssapy.io.read_csv(file_name, sep=None, dtypes=None, col=False, to_np=False, drop_nan=False, skiprows=[])

Read a CSV file with options.

Parameters:

• file_name (str) – The path to the CSV file.

• sep (str, optional) – The delimiter used in the CSV file. If None, delimiter will be guessed.

• dtypes (dict, optional) – Dictionary specifying data types for columns.

• col (bool or list of str, optional) – Specify columns to read. If False, read all columns.

• to_np (bool, optional) – Convert the loaded data to a NumPy array.

• drop_nan (bool, optional) – Drop rows with missing values (NaNs) from the loaded DataFrame.

• skiprows (list of int, optional) – Rows to skip while reading the CSV file.

Returns:

The loaded data in either a DataFrame or a NumPy array format.

Return type:

DataFrame or NumPy array

ssapy.io.read_csv_header(file_name, sep=None)

Get the header of a CSV file.

Parameters:

• file_name (str) – The filename of the CSV file.

• sep (str) – The delimiter used in the CSV file.

Returns:

A list of the header fields.

ssapy.io.read_h5(filename, pathname)

Load data from HDF5 file.

Parameters:

• filename (str) – The filename of the HDF5 file.

• pathname (str) – The path to the data in the HDF5 file.

Returns:

The data loaded from the HDF5 file.

ssapy.io.read_h5_all(file_path)

Read all datasets from an HDF5 file into a dictionary.

This function recursively traverses an HDF5 file and extracts all datasets into a dictionary. The keys of the dictionary are the paths to the datasets, and the values are the dataset contents.

Parameters:

file_pathstr

The path to the HDF5 file from which datasets will be read.

Returns:

dict

A dictionary where keys are the paths to datasets within the HDF5 file, and values are the contents of these datasets.

Examples:

>>> data = read_h5_all('example.h5')
>>> print(data.keys())
dict_keys(['/group1/dataset1', '/group2/dataset2'])
>>> print(data['/group1/dataset1'])
[1, 2, 3, 4, 5]

ssapy.io.read_tle(sat_name, tle_filename)

Get the TLE data from the file for the satellite with the given name

Parameters:

• sat_name (str) – NORAD name of the satellite

• tle_filename (str) – Path and name of file where TLE is

Returns:

line1, line2 – Both lines of the satellite TLE

Return type:

str

ssapy.io.read_tle_catalog(fname, n_lines=2)

Read in a TLE catalog file

Parameters:

• fname (string) – The filename

• n_lines (int) – number of lines per tle, usually 2 but sometimes 3

Returns:

catalog – lists containing TLEs

Return type:

list

ssapy.io.rmdir(source_)

Deletes a directory and its contents if it exists.

Parameters:

source_str

The path to the directory to be deleted.

ssapy.io.rmfile(pathname)

Deletes a file if it exists.

Parameters:

pathnamestr

The path to the file to be deleted.

ssapy.io.save_csv(file_name, df, sep='\t', dtypes=None)

Save a Pandas DataFrame to a CSV file.

Parameters:

• file_name (str) – The path to the CSV file.

• df (DataFrame) – The Pandas DataFrame to save.

• sep (str) – The delimiter used in the CSV file.

• dtypes (dict) – A dictionary specifying data types for columns.

Returns:

None

ssapy.io.save_csv_array_to_line(filename, array, delimiter='\t')

Appends a single row of data to a CSV file with a specified delimiter.

Parameters: filename (str): The name of the file to which the row will be appended. array (list): A list of values representing a single row of data to be appended to the CSV file. delimiter (str, optional): The delimiter to use between columns in the CSV file.

Default is tab (’ ‘).

Example: save_csv_array_to_line(‘output.csv’, [‘Alice’, 30, ‘New York’], delimiter=’,’)

ssapy.io.save_csv_header(filename, header, delimiter='\t')

Saves a header row to a CSV file with a specified delimiter.

Parameters: filename (str): The name of the file where the header will be saved. header (list): A list of strings representing the column names. delimiter (str, optional): The delimiter to use between columns in the CSV file.

Default is tab (’ ‘).

Example: save_csv_header(‘output.csv’, [‘Name’, ‘Age’, ‘City’], delimiter=’,’)

ssapy.io.save_csv_line(file_name, df, sep='\t', dtypes=None)

Save a Pandas DataFrame to a CSV file, appending the DataFrame to the file if it exists.

Parameters:

• file_name (str) – The path to the CSV file.

• df (DataFrame) – The Pandas DataFrame to save.

• sep (str) – The delimiter used in the CSV file.

Returns:

None

ssapy.io.save_h5(filename, pathname, data)

Save data to HDF5 file with recursive attempt in case of write errors.

Parameters:

• filename (str) – The filename of the HDF5 file.

• pathname (str) – The path to the data in the HDF5 file.

• data (any) – The data to be saved.

• max_retries (int) – Maximum number of recursive retries in case of write errors.

• retry_delay (tuple) – A tuple representing the range of delay (in seconds) between retries.

Returns:

None

ssapy.io.save_np(filename_, data_)

Save a NumPy array to a binary file.

This function saves a NumPy array to a file in .npy format. If the file cannot be created or written to, it handles common exceptions and prints an error message.

Parameters:

filename_str

The path to the file where the NumPy array will be saved.

data_numpy.ndarray

The NumPy array to be saved.

Returns:

None

The function does not return any value. It handles exceptions internally and prints error messages if any issues occur.

Examples:

>>> arr = np.array([1, 2, 3, 4, 5])
>>> save_np('array.npy', arr)

Linker

class ssapy.linker.BinarySelectorParams(num_tracks, num_orbits, init_value=1, dtype=<class 'numpy.float64'>)

Special methods for the binary model selection variables

TODO: Use a masked array here?

get_linked_track_indices(model_index)

Select those tracks that are linked to the given model_index

get_unlinked_track_indices()

Select those models that are not linked to any tracks

class ssapy.linker.Linker(iods, num_orbits=None, a_orbit=0.01)

Monte Carlo linking of uncorrelated orbit tracks

Parameters:

• iods (list) – List of Particles class instances summarizing initial orbit determinations (iods)

• num_orbits (int, optional) – Number of orbit models to fit. If not specified, set num_orbits = num_tracks

• a_orbit (float, optional) – Clustering hyperparameter in the Dirichlet prior on P_orbit (Default: 0.01)

iods

List of Particles class instances summarizing initial orbit determinations (iods)

Type:

list

num_tracks

Number of observed tracks. Equal to the length of the iods list.

Type:

int

num_orbits

Number of orbit models to fit. If not specified, set num_orbits = num_tracks

Type:

int, optional

orbit_selectors

The orbit selector parameters

Type:

BinarySelectorParams

p_orbit

The prior parameters on the orbit selectors

Type:

ModelSelectorParams

sample_Porbit_conditional_dist(track_ndx)

Sample from the conditional distribution for p_orbit at track_index

sample_orbit_selectors_from_data_conditional(track_ndx)

Sample from the conditional distribution for orbit_selector at track_index

update_orbit_parameters()

Find the groups of linked tracks and update all particles across linked groups

update_params_using_carlin_chib()

Perform MCMC step using algorithm of Carlin & Chib (1995).

save_step(outfile_head)

Save MCMC step to file

sample(nStep=100, outfile_head='linker')

Run the top-level MCMC to link tracks and refine orbit parameters

Parameters:

• nStep (int) – Number of MCMC steps

• outfile_head (str, optional) – Head (i.e., first part) of the output file name for chain steps

sample_Porbit_conditional_dist(track_ndx)

Sample from the conditional distribution for p_orbit at track_index

Parameters:

track_index (int) – The observation list index

Returns:

p_orbit – Draw from the conditional distribution for the p.orbit parameters

Return type:

array_like

sample_orbit_selectors_from_data_conditional(track_ndx, verbose=True)

Sample from the conditional distribution for orbit_selector at track_index ignoring the dependence orbit selector dependence of the conditional priors for the orbital elements.

Draw from :math:P(y | k, theta_k) P(k | P_orbit):

Parameters:

track_index (int) – The observation list index

Returns:

orbit_selectors – Draw from the conditional distribution for the orbit_selectors parameters

Return type:

array_like

save_step(outfile_head)

Save MCMC step to file

Parameters:

outfile_head (str) – Head (i.e., first part) of the output file name for chain steps

update_orbit_parameters()

Find the groups of linked tracks and update all particles across linked groups

update_params_using_carlin_chib(verbose=True)

Perform MCMC step using algorithm of Carlin & Chib (1995).

Sample the “active” orbit parameters from their conditional distribution (in a Gibbs update), the “inactive” orbit parameters from the pseudo-prior, and Gibbs updates for the orbit selectors and their prior parameters. This method is distinct from the Reversible Jump MCMC (RJMCMC) updates that have no need of the pseudo-prior at the expense of finding a suitable proposal distribution for the model updates.

The Gibbs updates for each class of parameters in turn defines a blocking procedure that can, in principle, give more efficient sampling than the joint updates in RJMCMC if a good proposal distribution for the orbit selectors is not known (see Godsill 2001).

Here the pseudo-prior for the “inactive” orbit parameters is not conjugate to the orbit selector prior, so there is no direct method to perform Gibbs samples of the orbit selectors. Instead, this function performs an MH update of orbit selectors with a proposal from the distribution determined by the product of the orbit likelihoods and the orbit selector prior (which would be the orbit selector conditional distribution if the orbit priors were independent of the orbit selectors).

Returns:

lnP – Current log-posterior value

Return type:

float

References

Carlin, B., Chib, S., 1995. Bayesian model choice via Markov chain Monte Carlo methods. Journal of the Royal Statistical Society. Series B (Method- ological) 57 (3), 473?484.

Godsill, S., 2001. On the relationship between Markov chain Monte Carlo methods for model uncertainty. Journal of Computational and Graphical Statistics 10 (2), 230?248.

class ssapy.linker.ModelSelectorParams(num_tracks, num_orbits, init_val=1.0, dtype=<class 'numpy.float64'>)

Container to hold model selection parameters or hyperparameters

Orbit Solver

Classes to solve for Keplerian orbital parameters from different initial inputs

# TODO: Provide overview of different solvers

TwoPosOrbitSolver:

GaussTwoPosOrbitSolver:

DanchickTwoPosOrbitSolver:

SheferTwoPosOrbitSolver:

ThreeAngleOrbitSolver

References: - Shefer, V. A. (2010). New method of orbit determination from two position vectors based on solving Gauss’s equations. Solar System Research, 44, 252-266. - Montenbruck, O., Gill, E., & Lutze, F. H. (2002). Satellite orbits: models, methods, and applications. Appl. Mech. Rev., 55(2), B27-B28.

class ssapy.orbit_solver.DanchickTwoPosOrbitSolver(*args, **kwargs)

static X(g)

Compute X(g) from Shefer (11).

static dXdg(g)

Compute dX(g)/dg from Shefer (12).

class ssapy.orbit_solver.GaussTwoPosOrbitSolver(*args, **kwargs)

class ssapy.orbit_solver.SheferTwoPosOrbitSolver(*args, **kwargs)

D(x)

Compute D(x) and its derivative, dD(x)/dx, from Shefer (43).

F(x)

Compute F(x) and dF(x)/dx from Shefer (40) and (41).

G(xi)

Compute G(xi) and its derivative from Shefer (A.7) and (A.8).

static X(x)

Evaluate X(x) function from Shefer (19) for elliptical orbits and (20) for hyperbolic orbits. The derivative of X(x) is given in (23).

Y(x)

Evaluate Y(x) and dY(x)/dx from Shefer (17).

static Z(xi)

Compute Z(xi) function from Shefer (A.5).

alpha(x)

Evaluate alpha(x) from Shefer (18).

beta(xi)

Evaluate beta(xi) and its derivative from Shefer (A.4) and (A.9).

semiMajorAxis(x)

Compute semi-major axis a(x) and its derivative, da(x)/dx, from Shefer (42).

class ssapy.orbit_solver.ThreeAngleOrbitSolver(e1, e2, e3, R1, R2, R3, t1, t2, t3, mu=398600441800000.0, tol=1e-08, maxiter=100)

Determine orbit of satellite for set of three angle-only observations.

Might only work well for smallish changes in sector position of observations.

Parameters:

• e1 ((3,) array_like) – Unit vectors indicating observed directions (dimensionless).

• e2 ((3,) array_like) – Unit vectors indicating observed directions (dimensionless).

• e3 ((3,) array_like) – Unit vectors indicating observed directions (dimensionless).

• R1 ((3,) array_like) – Vectors indicating the positions of the observers (in meters).

• R2 ((3,) array_like) – Vectors indicating the positions of the observers (in meters).

• R3 ((3,) array_like) – Vectors indicating the positions of the observers (in meters).

• t1 (float or astropy.time.Time) – Times of observations. If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• t2 (float or astropy.time.Time) – Times of observations. If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• t3 (float or astropy.time.Time) – Times of observations. If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• mu (float, optional) – Gravitational constant of central body in m^3/s^2. (Default: Earth’s gravitational constant in WGS84).

• tol (float, optional) – Tolerance to use to stop iteratively improving solution. (default 1e-8).

• maxiter (int, optional) – Maximum number of iterations to use. (default: 100)

class ssapy.orbit_solver.TwoPosOrbitSolver(r1, r2, t1, t2, mu=398600441800000.0, kappaSign=1, lam=0, eps=3e-16, maxiter=100)

Parameters:

• r1 ((3,) array_like) – Positions at t1 and t2 in meters.

• r2 ((3,) array_like) – Positions at t1 and t2 in meters.

• t1 (float or astropy.time.Time) – Times of observations. If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• t2 (float or astropy.time.Time) – Times of observations. If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• mu (float, optional) – Gravitational constant of central body in m^3/s^2. (Default: Earth’s gravitational constant in WGS84).

• kappaSign (int, optional) – +1 for prograde, -1 for retrograde

• lam (int, optional) – Number of complete orbits between observations. (Default: 0)

• eps (float, optional) – Iteration tolerance.

• maxiter (int, optional) – Maximum number of iterations.

Orbit

Module to handle Earth satellite state vector manipulation and propagation.

Notes

Although much of the Orbit class here is agnostic with respect to the coordinate reference frame, we strongly recommend that positions and velocities be specified in the GCRF frame. Some classes in this module do insist on this frame, (e.g., EarthObserver). We will try to indicate where the GCRF frame is required, but make no guarantees that we don’t miss any such instances.

class ssapy.orbit.EarthObserver(lon, lat, elevation=0, fast=False)

An earth-bound observer.

Parameters:

• lon (float) – Geodetic longitude in degrees (increasing to the East).

• lat (float) – Geodetic latitude in degrees.

• elevation (float, optional) – Elevation in meters.

• fast (bool, optional) – Use fast lookup tables for Earth Orientation parameters. ~ meter and ~10 micron/sec accurate for dates between approximately 1973 and the present. Less accurate outside this range.

lon

lat

elevation

itrs

ITRS coordinates in meters.

Type:

array_like (3,)

getRV(time)

Get position and velocity at specified time(s) in the GCRF frame.

getRV(time)

Get position and velocity at specified time(s).

Parameters:

time (array_like or astropy.time.Time (n,)) – If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

Returns:

• r (array_like (n, 3)) – Position in meters.

• v (array_like (n, 3)) – Velocity in meters per second.

Notes

The returned position and velocity are in the GCRF frame.

sunAlt(time)

Get Sun altitude for observer at time.

Parameters:

time (array_like or astropy.time.Time (n,)) – If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

Returns:

alt – Altitude of sun in radians.

Return type:

array_like(n,)

class ssapy.orbit.Orbit(r, v, t, mu=398600441800000.0, propkw=None)

Orbital state of one or more objects.

This class represents one or more instantaneous (osculating) states of orbiting bodies. Representations in Cartesian, Keplerian, and equinoctial elements are available. The default initializer requires Cartesian elements, though class methods are also available to initialize with Keplerian or equinoctial elements.

Note that generally this class can represent either a single scalar orbit, in which case attributes will generally be scalar floats, or a vector of orbits, in which case attributes will be ndarrays of floats. For attributes that are intrinsically vectors even for a single orbit (position, velocity, …), a 2d array will be used for a “vector of orbits”, with the first dimension representing the different orbits.

For simplicity, we will only indicate the “single scalar Orbit” dimensions in docstrings. For a vector-Orbit, most scalar input arguments can be supplied as broadcastable arrays, and scalar attributes become array attributes. When a multi-dimensional array is required, the first dimension is the one over different orbits.

Parameters:

• r ((3,) array_like) – Position of orbiting object in meters.

• v ((3,) array_like) – Velocity of orbiting objects in meters per second.

• t (float or astropy.time.Time) – If float or array of float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• mu (float, optional) – Gravitational constant of central body in m^3/s^2. (Default: Earth’s gravitational constant in WGS84).

r

Position of orbiting object in meters.

Type:

(3,) array_like

v

Velocity of orbiting object in meters per second.

Type:

(3,) array_like

t

GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

Type:

float

mu

Gravitational constant of central body in m^3/s^2

Type:

float

a

Semimajor axis in meters.

Type:

float

hx, hy

Components of the equinoctial inclination vector.

Type:

float

ex, ey

Components of the equinoctial eccentricity vector.

Type:

float

lv

True longitude in radians.

Type:

float

lE

Eccentric longitude in radians.

Type:

float

lM

Mean longitude in radians.

Type:

float

e

Keplerian eccentricity.

Type:

float

i

Keplerian inclination in radians.

Type:

float

pa

Keplerian periapsis argument in radians.

Type:

float

raan

Keplerian right ascension of the ascending node in radians.

Type:

float

trueAnomaly

Keplerian true anomaly in radians.

Type:

float

eccentricAnomaly

Keplerian eccentric anomaly in radians.

Type:

float

meanAnomaly

Keplerian mean anomaly in radians.

Type:

float

period

Orbital period in seconds.

Type:

float

meanMotion

Keplerian mean motion in radians per second.

Type:

float

p

Semi-latus rectum in meters.

Type:

float

angularMomentum

(Specific) angular momentum in m^2/s.

Type:

(3,) array_like

energy

(Specific) orbital energy in m^2/s^2.

Type:

float

LRL

Laplace-Runge-Lenz vector in m^3/s^2.

Type:

(3,) array_like

periapsis

Periapsis coordinate of orbit in meters.

Type:

(3,) array_like

apoapsis

Apoapsis coordinate of orbit in meters.

Type:

(3,) array_like

equinoctialElements

keplerianElements

fromKeplerianElements(a, e, i, pa, raan, trueAnomaly, t, mu)

Construct an Orbit from Keplerian elements.

fromEquinoctialElements(a, hx, hy, ex, ey, lv, t, mu)

Construct an Orbit from Equinoctial elements.

at(t, propagator)

Propagate orbit to time t and return new Orbit.

Notes

Although much of the Orbit class here is agnostic with respect to the coordinate reference frame, we strongly recommend that positions and velocities be specified in the GCRF frame. Some other classes in this module do insist on this frame, (e.g., EarthObserver). We will try to indicate where the GCRF frame is required, but make no guarantees that we don’t miss any such instances.

at(t, propagator=KeplerianPropagator())

Propagate this orbit to time t.

Parameters:

• t (float or astropy.time.Time) – If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• propagator (Propagator, optional) – The propagator instance to use.

Returns:

Propagated Orbit.

Return type:

Orbit

property equinoctialElements

Equinoctial elements (a, hx, hy, ex, ey, lv).

classmethod fromEquinoctialElements(a, hx, hy, ex, ey, lv, t, mu=398600441800000.0, propkw=None)

Construct an Orbit from equinoctial elements.

Parameters:

• a (float) – Semimajor axis in meters.

• hx (float) – Components of the equinoctial inclination vector.

• hy (float) – Components of the equinoctial inclination vector.

• ex (float) – Components of the equinoctial eccentricity vector.

• ey (float) – Components of the equinoctial eccentricity vector.

• lv (float) – True longitude in radians.

• t (float or astropy.time.Time) – If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• mu (float, optional) – Gravitational constant of central body in m^3/s^2. (Default: Earth’s gravitational constant in WGS84).

Returns:

The Orbit with given parameters.

Return type:

Orbit

classmethod fromKeplerianElements(a, e, i, pa, raan, trueAnomaly, t, mu=398600441800000.0, propkw=None)

Construct an Orbit from Keplerian elements.

Parameters:

• a (float) – Semimajor axis in meters.

• e (float) – Keplerian eccentricity.

• i (float) – Keplerian inclination in radians.

• pa (float) – Keplerian periapsis argument in radians.

• raan (float) – Keplerian right ascension of the ascending node in radians.

• trueAnomaly (float) – Keplerian true anomaly in radians.

• t (float or astropy.time.Time) – If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• mu (float, optional) – Gravitational constant of central body in m^3/s^2. (Default: Earth’s gravitational constant in WGS84).

Returns:

The Orbit with given parameters.

Return type:

Orbit

classmethod fromKozaiMeanKeplerianElements(a, e, i, pa, raan, trueAnomaly, t, mu=398600441800000.0, _useTEME=False, propkw=None)

Construct an Orbit from Kozai mean Keplerian elements. By default, this method also converts from the TEME reference frame to the GCRF frame. This is mainly to support TLEs and SGP4, so the default gravitational parameter is WGS84.

Parameters:

• a (float) – Semimajor axis in meters.

• e (float) – Keplerian eccentricity.

• i (float) – Keplerian inclination in radians.

• pa (float) – Keplerian periapsis argument in radians.

• raan (float) – Keplerian right ascension of the ascending node in radians.

• trueAnomaly (float) – Keplerian true anomaly in radians.

• t (float or astropy.time.Time) – If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• mu (float, optional) – Gravitational constant of central body in m^3/s^2. (Default: Earth’s gravitational constant in WGS84).

Returns:

The Orbit with given parameters.

Return type:

Orbit

classmethod fromTLE(sat_name, tle_filename, propkw=None)

Construct an Orbit from a TLE file

Parameters:

• sat_name (str) – NORAD name of the satellite

• tle_filename (str) – Path and name of file where TLE is

Return type:

Orbit

classmethod fromTLETuple(tle, propkw=None)

Construct an Orbit from a TLE tuple

Parameters:

tle (2-tuple of str) – Line1 and Line2 of TLE as strings

Return type:

Orbit

property keplerianElements

Keplerian elements (a, e, i, pa, raan, trueAnomaly).

class ssapy.orbit.OrbitalObserver(orbit, propagator=KeplerianPropagator())

An observer in orbit.

Parameters:

• orbit (Orbit) – The orbit of the observer.

• propagator (Propagator, optional) – The propagator instance to use.

orbit

propagator

getRV(time)

Get position and velocity at specified time(s).

getRV(time)

Get position and velocity at specified time(s).

Parameters:

time (array_like or astropy.time.Time (n,)) – If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

Returns:

• r (array_like (n, 3)) – Position in meters.

• v (array_like (n, 3)) – Velocity in meters per second.

Particles

Classes for sampling orbit model parameters.

class ssapy.particles.Particles(particles, rvprobability, lnpriors=None, ln_weights=None)

A class for importance (re-)sampling orbit model parameter samples

Parameters:

• particles ((num_particles, 6) array_like) – Positions and velocities of orbiting object at epoch in meters and meters / second. The ‘chain’ part of the output from rvsampler.XXSampler.sample

• rvprobability (class instance) – An instance of an RVProbability for this epoch. The Particles class wraps the lnlike method of the supplied RVProbability object.

• ln_weights ((num_particles,) array_like, optional) – Weights for each particle in the input particles

particles

Positions and velocities of orbiting object at epoch in meters and meters / second. The ‘chain’ part of the output from RVSampler.sample

Type:

(num_particles, ) array_like

rvprobability

An instance of an RVProbability for this epoch. The Particles class wraps the lnlike method of the supplied RVProbability object.

Type:

class instance

ln_wts

Log weights for each particle.

Type:

(num_particles,) array_like

num_particles

Number of particles

Type:

int

initial_particles

Copy of input particles to enable reset_to_pseudo_prior()

Type:

(num_particles, 6) array_like

initial_ln_wts

Copy of input log weights to enable reset_to_pseudo_prior()

Type:

(num_particles,) array_like

epoch

The time at which the model parameters (position and velocity) are specified. If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC.

Type:

float or astropy.time.Time

orbits

List of Orbit class instances derived from each particle

Type:

list

lnpriors

Log prior probabilities for each particle. The priors are as set in the input RVProbability object.

Type:

(num_particles,) array_like

reset_to_pseudo_prior()

Reset the particles and weights to their values at instantiation

move(epoch)

Move particles to the specified epoch

lnlike(epoch_particles)

Evaluate the log-likelihood of the input particles given internally stored measurements

reweight(epoch_particles)

Reweight particles using cross-epoch likelihoods

resample(num_particles)

Resample particles to achieve near-uniform weighting per particle

fuse(epoch_particles, verbose)

Fuse particles from a different epoch and fuse internal particles and weights

draw_orbit()

Draw a single orbit model from the internal list of particles

mean()

Evaluate the weighted mean of the particle values

draw_orbit()

Draw a single orbit model from the internal particle collection with probability proportional to the internal weights.

Returns:

particle – Position and velocity for one particle

Return type:

array (1, 6)

property epoch

Get the epoch time at which the particles (i.e., position and velocity) are defined

fuse(epoch_particles, verbose=False)

Add particles from a new epoch and fuse all epoch particles and weights accordingly

Parameters:

• epoch_particles (Particles) – The particles object from another observation

• verbose (bool, optional) – Enable verbose output to stdout

Return type:

None. Side effect is to fuse the internal particles state.

lnlike(orbits)

Evaluate this epoch likelihood given one or more particles

Parameters:

orbits (List of Orbit class instances) – A list of the Orbit models for each particle, e.g., as generated by Particles.orbits

Returns:

lnlike – Array of log-likelihood values for each particle

Return type:

array (num_particles, )

property lnpriors

Get an array of log prior probabilities for each particle

mean()

Evaluate the weighted mean of all particles

Returns:

mean – Weighted mean of the particle values (3*m, 3*m/s)

Return type:

(6,) array_like

move(epoch)

Move particles to the given epoch

Parameters:

epoch (astropy.time.Time) –

Returns:

particles – Position and velocity for each particle at the new epoch

Return type:

array (num_particles, 6)

property orbits

Get a list of orbits corresponding to each stored particle

resample(num_particles)

Resample particles to achieve near-uniform weighting per particle

Parameters:

num_particles (int) – The number of particles to keep

Return type:

None. Side effect is to update the internal particles state.

reweight(epoch_particles)

Reweight particles using cross-epoch likelihoods

The weights for the kth particle are,

w_k = log(Prod_i L_(d_i | theta_k) / L(d_j | theta_k))

= Sum_{i /= j} log(L(d_i | theta_k))

where j is the current epoch, theta are the particle parameters, and d_i is the data from epoch i.

Parameters:

epoch_particles (Particles) – The particles object from another observation

Returns:

status – Status (True means success). Side effect is to update the internal particles state.

Return type:

bool

Plot Utilities

Utility functions for plotting.

ssapy.plotUtils.check_numpy_array(variable: ndarray | list) → str

Checks if the input variable is a NumPy array, a list of NumPy arrays, or neither.

Parameters:

variable (Union[np.ndarray, list]) – The variable to check. It can either be a NumPy array or a list of NumPy arrays.

Returns:

Returns a string indicating the type of the variable: - “numpy array” if the variable is a single NumPy array, - “list of numpy array” if it is a list of NumPy arrays, - “not numpy” if it is neither.

Return type:

str

ssapy.plotUtils.draw_dashed_circle(ax, normal_vector, radius, dashes, dash_length=0.1, label='Dashed Circle')

Draw a dashed circle on a 3D axis with a given normal vector.

Parameters: - ax (matplotlib.axes._subplots.Axes3DSubplot): The 3D axis on which to draw the circle. - normal_vector (array-like): A 3-element array representing the normal vector to the plane of the circle. - radius (float): The radius of the circle. - dashes (int): The number of dashes to be used in drawing the circle. - dash_length (float, optional): The relative length of each dash, as a fraction of the circle’s circumference. Default is 0.1. - label (str, optional): The label for the circle. Default is ‘Dashed Circle’.

Returns: None

This function draws a dashed circle on a 3D axis. The circle is defined in the xy-plane, then rotated to align with the given normal vector. The circle is divided into dashes to create the dashed effect.

Example usage: ``` import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from your_module import draw_dashed_circle

fig = plt.figure() ax = fig.add_subplot(111, projection=’3d’) normal_vector = [0, 0, 1] radius = 5 dashes = 20

draw_dashed_circle(ax, normal_vector, radius, dashes)

plt.show() ```

ssapy.plotUtils.draw_earth(time, ngrid=100, R=6378137.0, rfactor=1)

Parameters:

• time (array_like or astropy.time.Time (n,)) – If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• ngrid (int) – Number of grid points in Earth model.

• R (float) – Earth radius in meters. Default is WGS84 value.

• rfactor (float) – Factor by which to enlarge Earth (for visualization purposes)

ssapy.plotUtils.draw_moon(time, ngrid=100, R=1738100.0, rfactor=1)

Parameters:

• time (array_like or astropy.time.Time (n,)) – If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• ngrid (int) – Number of grid points in Earth model.

• R (float) – Earth radius in meters. Default is WGS84 value.

• rfactor (float) – Factor by which to enlarge Earth (for visualization purposes)

ssapy.plotUtils.format_date_axis(time_array, ax)

Format the x-axis of a plot with time-based labels depending on the span of the time array.

Parameters: - time_array (array-like): An array of time objects (e.g., astropy.time.Time) to be used for the x-axis labels. - ax (matplotlib.axes.Axes): The matplotlib axes object on which to set the x-axis labels.

Returns: None

This function formats the x-axis labels of a plot based on the span of the provided time array. The labels are set to show either hours and day-month or month-year formats, depending on the time span.

The function performs the following steps: 1. If the time span is less than one month:

• If the time span is less than a day, the labels show ‘HH:MM dd-Mon’.

• Otherwise, the labels show ‘dd-Mon-YYYY’.

2. If the time span is more than one month, the labels show ‘Mon-YYYY’.

The function selects six nearly evenly spaced points in the time array to set the x-axis labels.

Example usage: ``` import matplotlib.pyplot as plt from astropy.time import Time import numpy as np

# Example time array time_array = Time([‘2024-07-01T00:00:00’, ‘2024-07-01T06:00:00’, ‘2024-07-01T12:00:00’,

‘2024-07-01T18:00:00’, ‘2024-07-02T00:00:00’])

fig, ax = plt.subplots() ax.plot(time_array.decimalyear, np.random.rand(len(time_array))) format_date_axis(time_array, ax) plt.show() ```

ssapy.plotUtils.globe_plot(r, t, limits=False, title='', figsize=(7, 8), save_path=False, el=30, az=0, scale=1)

Plot a 3D scatter plot of position vectors on a globe representation.

Parameters: - r (array-like): Position vectors with shape (n, 3), where n is the number of points. - t (array-like): Time array corresponding to the position vectors. This parameter is not used in the current function implementation but is included for consistency. - limits (float, optional): The limit for the plot axes. If not provided, it is calculated based on the data. Default is False. - title (str, optional): Title of the plot. Default is an empty string. - figsize (tuple of int, optional): Figure size (width, height) in inches. Default is (7, 8). - save_path (str, optional): Path to save the generated plot. If not provided, the plot will not be saved. Default is False. - el (int, optional): Elevation angle (in degrees) for the view of the plot. Default is 30. - az (int, optional): Azimuth angle (in degrees) for the view of the plot. Default is 0. - scale (int, optional): Scale factor for resizing the Earth image. Default is 1.

Returns: - fig (matplotlib.figure.Figure): The figure object containing the plot. - ax (matplotlib.axes._subplots.Axes3DSubplot): The 3D axis object used in the plot.

The function creates a 3D scatter plot of the position vectors on a globe. The globe is represented using a textured Earth image, and the scatter points are colored using a rainbow colormap. The plot’s background is set to black, and the plot is displayed with customizable elevation and azimuth angles.

Example usage: ``` import numpy as np from your_module import globe_plot

# Example data r = np.array([[1, 2, 3], [4, 5, 6]]) # Replace with actual data t = np.arange(len(r)) # Replace with actual time data

globe_plot(r, t, save_path=’globe_plot.png’) ```

ssapy.plotUtils.groundTrackVideo(r, time)

Parameters:

• r ((3,) array_like) – Position of orbiting object in meters.

• t (float or astropy.time.Time) – If float or array of float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

ssapy.plotUtils.ground_track_plot(r, t, ground_stations=None, save_path=False)

Parameters:

• r ((3,) array_like - Orbit positions in meters.) –

• t ((n,) array_like - array of Astropy Time objects or time in gps seconds.) –

• ground_stations (optional -) –

ssapy.plotUtils.koe_hist_2d(stable_data, title='Initial orbital elements of\n1 year stable cislunar orbits', limits=[1, 50], bins=200, logscale=False, cmap='coolwarm', save_path=False)

Create a 2D histogram plot for various Keplerian orbital elements of stable cislunar orbits.

Parameters: - stable_data (object): An object with attributes a, e, i, and ta, which are arrays of semi-major axis, eccentricity, inclination, and true anomaly, respectively. - title (str, optional): Title of the figure. Default is “Initial orbital elements of

1 year stable cislunar orbits”.

• limits (list, optional): Color scale limits for the histogram. Default is [1, 50].

• bins (int, optional): Number of bins for the 2D histograms. Default is 200.

• logscale (bool or str, optional): Whether to use logarithmic scaling for the color bar. Default is False. Can also be ‘log’ to apply logarithmic scaling.

• cmap (str, optional): Colormap to use for the histograms. Default is ‘coolwarm’.

• save_path (str, optional): Path to save the generated plot. If not provided, the plot will not be saved. Default is False.

Returns: - fig (matplotlib.figure.Figure): The figure object containing the 2D histograms.

This function creates a 3x3 grid of 2D histograms showing the relationships between various orbital elements, including semi-major axis, eccentricity, inclination, and true anomaly. The color scale of the histograms can be adjusted with a logarithmic or linear normalization. The plot is customized with labels and a color bar.

Example usage: ``` import numpy as np from your_module import koe_hist_2d

# Example data class StableData:

def __init__(self):

self.a = np.random.uniform(1, 20, 1000) self.e = np.random.uniform(0, 1, 1000) self.i = np.radians(np.random.uniform(0, 90, 1000)) self.ta = np.radians(np.random.uniform(0, 360, 1000))

stable_data = StableData() koe_hist_2d(stable_data, save_path=’orbit_histograms.pdf’) ```

ssapy.plotUtils.koe_plot(r, v, t=<Time object: scale='utc' format='iso' value=['2025-01-01 00:00:00.000' '2025-01-01 00:59:57.946'  '2025-01-01 01:59:55.893' ... '2025-12-31 22:00:04.107'  '2025-12-31 23:00:02.054' '2026-01-01 00:00:00.000']>, elements=['a', 'e', 'i'], save_path=False, body='Earth')

Plot Keplerian orbital elements over time for a given trajectory.

Parameters: - r (array-like): Position vectors for the orbit. - v (array-like): Velocity vectors for the orbit. - t (array-like, optional): Time array for the plot, given as a sequence of astropy.time.Time objects or a Time object with np.linspace. Default is one year of hourly intervals starting from “2025-01-01”. - elements (list of str, optional): List of orbital elements to plot. Options include ‘a’ (semi-major axis), ‘e’ (eccentricity), and ‘i’ (inclination). Default is [‘a’, ‘e’, ‘i’]. - save_path (str, optional): Path to save the generated plot. If not provided, the plot will not be saved. Default is False. - body (str, optional): The celestial body for which to calculate the orbital elements. Options are ‘Earth’ or ‘Moon’. Default is ‘Earth’.

Returns: - fig (matplotlib.figure.Figure): The figure object containing the plot. - ax1 (matplotlib.axes.Axes): The primary axis object used in the plot.

The function calculates orbital elements for the given position and velocity vectors, and plots these elements over time. It creates a plot with two y-axes: one for the eccentricity and inclination, and the other for the semi-major axis. The x-axis represents time in decimal years.

Example usage: ``` import numpy as np from astropy.time import Time from your_module import koe_plot

# Example data r = np.array([[[1, 0, 0], [0, 1, 0]]]) # Replace with actual data v = np.array([[[0, 1, 0], [-1, 0, 0]]]) # Replace with actual data t = Time(“2025-01-01”, scale=’utc’) + np.linspace(0, int(1 * 365.25), int(365.25 * 24))

koe_plot(r, v, t, save_path=’orbital_elements_plot.png’) ```

ssapy.plotUtils.orbit_divergence_plot(rs, r_moon=[], t=False, limits=False, title='', save_path=False)

Plot multiple cislunar orbits in the GCRF frame with respect to the Earth and Moon.

Parameters: - rs (numpy.ndarray): A 3D array of shape (n, 3, m) where n is the number of time steps,

3 represents the x, y, z coordinates, and m is the number of orbits.

• r_moon (numpy.ndarray, optional): A 2D array of shape (3, n) representing the Moon’s position at each time step.

  If not provided, it is calculated based on the time t.

• t (astropy.time.Time, optional): The time at which to calculate the Moon’s position if r_moon is not provided. Default is False.

• limits (float, optional): The plot limits in units of Earth’s radius (GEO). If not provided, it is calculated as 1.2 times the maximum norm of rs. Default is False.

• title (str, optional): The title of the plot. Default is an empty string.

• save_path (str, optional): The file path to save the plot. If not provided, the plot is not saved. Default is False.

Returns: None

This function creates a 3-panel plot of multiple cislunar orbits in the GCRF frame. Each panel represents a different plane (xy, xz, yz) with Earth at the center. The orbits are plotted with color gradients to indicate progression. The Moon’s position is also plotted if provided or calculated.

Example usage: ``` import numpy as np from astropy.time import Time from your_module import orbit_divergence_plot

# Example data rs = np.random.randn(100, 3, 5) # 5 orbits with 100 time steps each t = Time(“2025-01-01”)

orbit_divergence_plot(rs, t=t, title=’Cislunar Orbits’) ```

ssapy.plotUtils.save_animated_gif(gif_name, frames, fps=30)

Create a GIF from a sequence of image frames.

Parameters: - gif_name (str): The name of the output GIF file, including the .gif extension. - frames (list of str): A list of file paths to the image frames to be included in the GIF. - fps (int, optional): Frames per second for the GIF. Default is 30.

Returns: None

This function uses the imageio library to write a GIF file. It prints messages indicating the start and completion of the GIF writing process. Each frame is read from the provided file paths and appended to the GIF.

Example usage: frames = [‘frame1.png’, ‘frame2.png’, ‘frame3.png’] write_gif(‘output.gif’, frames, fps=24)

ssapy.plotUtils.save_plot(figure, save_path, dpi=200)

Save a Python figure as a PNG/JPG/PDF/ect. image. If no extension is given in the save_path, a .png is defaulted.

Parameters:

• figure (matplotlib.figure.Figure) – The figure object to be saved.

• save_path (str) – The file path where the image will be saved.

Returns:

None

ssapy.plotUtils.save_plot_to_pdf(figure, pdf_path)

Save a Matplotlib figure to a PDF file, with support for merging with existing PDFs.

Parameters: - figure (matplotlib.figure.Figure): The Matplotlib figure to be saved. - pdf_path (str): The path to the PDF file. If the file exists, the figure will be appended to it.

Returns: None

This function saves a Matplotlib figure as a PNG in-memory and then converts it to a PDF. If the specified PDF file already exists, the new figure is appended to it. Otherwise, a new PDF file is created. The function also keeps track of how many times it has been called using a global variable save_plot_to_pdf_call_count.

The function performs the following steps: 1. Expands the user directory if the path starts with ~. 2. Generates a temporary PDF path by appending “_temp.pdf” to the original path. 3. Saves the figure as a PNG in-memory using a BytesIO buffer. 4. Opens the in-memory PNG using PIL and creates a new figure to display the image. 5. Saves the new figure with the image into a temporary PDF. 6. If the specified PDF file exists, merges the temporary PDF with the existing one.

Otherwise, renames the temporary PDF to the specified path.

7. Closes the original and temporary figures and prints a message indicating the save location.

Example usage: ``` import matplotlib.pyplot as plt

fig, ax = plt.subplots() ax.plot([1, 2, 3], [4, 5, 6]) save_plot_to_pdf(fig, ‘~/Desktop/my_plot.pdf’) ```

ssapy.plotUtils.scatter_2d(x, y, cs, xlabel='x', ylabel='y', title='', cbar_label='', dotsize=1, colorsMap='jet', colorscale='linear', colormin=False, colormax=False, save_path=False)

Create a 2D scatter plot with optional color mapping.

Parameters: - x (numpy.ndarray): Array of x-coordinates. - y (numpy.ndarray): Array of y-coordinates. - cs (numpy.ndarray): Array of values for color mapping. - xlabel (str, optional): Label for the x-axis. Default is ‘x’. - ylabel (str, optional): Label for the y-axis. Default is ‘y’. - title (str, optional): Title of the plot. Default is an empty string. - cbar_label (str, optional): Label for the color bar. Default is an empty string. - dotsize (int, optional): Size of the dots in the scatter plot. Default is 1. - colorsMap (str, optional): Colormap to use for the color mapping. Default is ‘jet’. - colorscale (str, optional): Scale for the color mapping, either ‘linear’ or ‘log’. Default is ‘linear’. - colormin (float, optional): Minimum value for color scaling. If False, it is set to the minimum value of cs. Default is False. - colormax (float, optional): Maximum value for color scaling. If False, it is set to the maximum value of cs. Default is False. - save_path (str, optional): File path to save the plot. If not provided, the plot is not saved. Default is False.

Returns: - fig (matplotlib.figure.Figure): The figure object. - ax (matplotlib.axes._subplots.AxesSubplot): The 2D axis object.

This function creates a 2D scatter plot with optional color mapping based on the values provided in cs. The color mapping can be adjusted using either a linear or logarithmic scale. The plot can be customized with axis labels, title, and colormap. The plot can also be saved to a specified file path.

Example usage: ``` import numpy as np from your_module import scatter_2d

# Example data x = np.random.rand(100) y = np.random.rand(100) cs = np.random.rand(100)

scatter_2d(x, y, cs, xlabel=’X-axis’, ylabel=’Y-axis’, cbar_label=’Color Scale’, title=’2D Scatter Plot’) ```

ssapy.plotUtils.scatter_3d(x, y=None, z=None, cs=None, xlabel='x', ylabel='y', zlabel='z', cbar_label='', dotsize=1, colorsMap='jet', title='', save_path=False)

Create a 3D scatter plot with optional color mapping.

Parameters: - x (numpy.ndarray): Array of x-coordinates or a 2D array with shape (n, 3) representing the x, y, z coordinates. - y (numpy.ndarray, optional): Array of y-coordinates. Required if x is not a 2D array with shape (n, 3). Default is None. - z (numpy.ndarray, optional): Array of z-coordinates. Required if x is not a 2D array with shape (n, 3). Default is None. - cs (numpy.ndarray, optional): Array of values for color mapping. Default is None. - xlabel (str, optional): Label for the x-axis. Default is ‘x’. - ylabel (str, optional): Label for the y-axis. Default is ‘y’. - zlabel (str, optional): Label for the z-axis. Default is ‘z’. - cbar_label (str, optional): Label for the color bar. Default is an empty string. - dotsize (int, optional): Size of the dots in the scatter plot. Default is 1. - colorsMap (str, optional): Colormap to use for the color mapping. Default is ‘jet’. - title (str, optional): Title of the plot. Default is an empty string. - save_path (str, optional): File path to save the plot. If not provided, the plot is not saved. Default is False.

Returns: - fig (matplotlib.figure.Figure): The figure object. - ax (matplotlib.axes._subplots.Axes3DSubplot): The 3D axis object.

This function creates a 3D scatter plot with optional color mapping based on the values provided in cs. The plot can be customized with axis labels, title, and colormap. The plot can also be saved to a specified file path.

Example usage: ``` import numpy as np from your_module import scatter_3d

# Example data x = np.random.rand(100) y = np.random.rand(100) z = np.random.rand(100) cs = np.random.rand(100)

scatter_3d(x, y, z, cs, xlabel=’X-axis’, ylabel=’Y-axis’, zlabel=’Z-axis’, cbar_label=’Color Scale’, title=’3D Scatter Plot’) ```

ssapy.plotUtils.set_color_theme(fig, *axes, theme)

Set the color theme of the figure and axes to white or black and the text color to white or black.

Parameters: - fig (matplotlib.figure.Figure): The figure to modify. - axes (list of matplotlib.axes._subplots.AxesSubplot): One or more axes to modify. - theme (str) either black/dark or white.

Returns: - fig (matplotlib.figure.Figure): The modified figure. - axes (tuple of matplotlib.axes._subplots.AxesSubplot): The modified axes.

This function changes the background color of the given figure and its axes to black or white. It also sets the color of all text items (title, labels, tick labels) to white or black.

Example usage: ``` import matplotlib.pyplot as plt

fig, ax = plt.subplots() ax.plot([1, 2, 3], [4, 5, 6]) set_color_theme(fig, ax, theme=’black’) plt.show() ```

Propagator

Classes for propagating orbits.

class ssapy.propagator.KeplerianPropagator

A basic Keplerian propagator for finding the position and velocity of an orbiting object at some future or past time.

class ssapy.propagator.Propagator

Abstract base class for orbit propagators.

Note, the interface for propagators is through functions in ssapy.compute.

class ssapy.propagator.RK4Propagator(accel, h)

Runge-Kutta 4th order numerical integrator.

Parameters:

• accel (ssapy.Accel) – Accel object containing the acceleration model by which to propagate.

• h (float) – Step size in seconds. Reasonable values are ~50s for LEO propagations over a ~day for ~meter accuracy, or 1000s for GEO propagations over a few days with ~meter accuracy. For best results, check for convergence.

class ssapy.propagator.RK78Propagator(accel, h, tol=(1e-06, 1e-06, 1e-06, 1e-09, 1e-09, 1e-09))

Runge-Kutta 8th order numerical integrator with adaptive step size computed from embedded 7th order integrator error estimate.

Parameters:

• accel (ssapy.Accel) – Accel object containing the acceleration model by which to propagate.

• h (float) – Initial step size in seconds. A few 10s of seconds is usually a good starting point here; it’ll automatically be adjusted by the algorithm.

• tol (float or array of float.) – Tolerance for a single integrator step. Used to adaptively change the integrator step size. Broadcasts to 6-dimensions. A good target is usually ~[1e-6, 1e-6, 1e-6, 1e-9, 1e-9, 1e-9] for cm accuracy at GEO over a few days, or around LEO over a few hours.

class ssapy.propagator.RK8Propagator(accel, h)

Runge-Kutta 8th order numerical integrator.

Parameters:

• accel (ssapy.Accel) – Accel object containing the acceleration model by which to propagate.

• h (float) – Step size in seconds. ~70s yields accuracy of ~1e-6 meters at GEO over a couple of days. ~20s yields accuracy of ~1e-5 meters at LEO over a few hours.

class ssapy.propagator.RKPropagator

class ssapy.propagator.SGP4Propagator(t=None, truncate=False)

Propagate using simplified perturbation model SGP4.

Parameters:

• t (float or astropy.time.Time, optional) –

  Reference time at which to compute frame transformation between GCRF and TEME. SGP4 calculations occur in the TEME frame, but useful input and output is in the GCRF frame. In principle, one could do the transformation at every instant in time for which the orbit is queried. However, the rate of change in the transformation is small, ~0.15 arcsec per day, so here we just use a single transformation.

  If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

  If None, then use the time of the orbit being propagated.

• truncate (bool, optional) – Truncate elements to precision of TLE ASCII format? This may be required in order to reproduce the results of running sgp4 directly from a TLE.

class ssapy.propagator.SciPyPropagator(accel, ode_kwargs=None)

Propagate using the scipy.integrate.solve_ivp ODE solver.

Parameters:

• accel (ssapy.Accel) – Accel object containing the acceleration model by which to propagate.

• ode_kwargs (dict) – Keyword arguments to pass to scipy.integrate.solve_ivp. Of particular interest may be the kwarg rtol, which usually yields reasonable results when set ~1e-7. For best results, check for convergence.

class ssapy.propagator.SeriesPropagator(order=2)

Propagate using Taylor series expansion of Keplerian motion. This is quite a bit faster than the full Keplerian propagator, but only valid for short time intervals. Using order=2 for a few seconds of time though is definitely reasonable.

Parameters:

order (int, optional) – Order of series expansion. 1 = constant velocity, 2 = constant acceleration

ssapy.propagator.default_numerical(*args, cls=None, accel=None, extra_accel=None)

Construct a numerical propagator with sensible default acceleration.

Parameters:

• *args (list) – Arguments to Propagator

• cls (Propagator) – class to use. Default of None means SciPyPropagator.

• accel (Accel) – acceleration model to use. Default of None means Earth(4, 4), sun, moon.

• extra_accel (Accel or list of Accel, optional) – Additional accelerations to add to accel.

Return type:

Instance of Propagator with desired Accel model.

Position and Velocity Sampler

class ssapy.rvsampler.APrior(amean, asigma)

Gaussian prior on orbit semimajor axis a.

Parameters:

• amean (float) – Prior mean on a in meters.

• asigma (float) – Prior standard deviation on a in meters.

aMean

aSigma

__call__(orbit)

Returns log prior probability of given orbit.

class ssapy.rvsampler.AreaPrior(mean=0, sigma=2)

Gaussian prior on log10(area) of object.

Parameters:

• mean (float) – mean of log10(area)

• sigma (float) – standard deviation of log10(area)

mean, sigma

__call__(orbit)

Returns log prior probability of given orbit.

class ssapy.rvsampler.DirectInitializer(samples, replace=False)

Directly specify initial position and velocity samples. This can be useful if restarting sampling after some previous sampling step has completed, for instance, if adding new observations to a previously sampled arc.

Parameters:

• samples (array_like (nSample, 6)) – Samples in position (in meters) and velocity (in meters per second). Columns are [x, y, z, vx, vy, vz].

• replace (bool, optional) – Whether to sample with or without replacement.

samples

replace

__call__(nSample)

Returns nSample initial samples.

class ssapy.rvsampler.DistanceProjectionInitializer(arc, rho, indices=[0, 1], rSigma=np.float64(4216417.236450565), vSigma=np.float64(307.46599995993046))

Initialize position and velocity from two closely spaced (in time) observations of ra/dec. The position is initialized by projecting along the direction of the first observation a given slant range. The velocity is then initialized assuming that the slant range derivative is zero, and that the motion in between observations is inertial.

Parameters:

• arc (QTable with columns for 'ra', 'dec', 'rStation_GCRF'.) – The set of observations used to initialize samples. Note, only the first row of the given QTable will be used.

• rho (Distance in meters to project, measured from the position of the) – observer (not the center of the Earth).

• indices (list of indices indicating which two observations in arc to use) – for initialization. Default: [0, 1]; i.e., use the first two observations.

• rSigma (float, optional) – Standard deviation of sample positions dispersion in meters.

• vSigma (float, optional) – Standard deviation of sample velocity dispersion in meters per second.

observations

Indicated two rows of the input arc.

Type:

QTable

rho

rSigma

vSigma

__call__(nSample)

Returns nSample initial samples.

class ssapy.rvsampler.EPrior(emean, esigma)

Gaussian prior on orbit eccentricity e.

Parameters:

• emean (float) – Prior mean on e.

• esigma (float) – Prior standard deviation on e.

eMean

eSigma

__call__(orbit)

Returns log prior probability of given orbit.

class ssapy.rvsampler.EmceeSampler(probfn, initializer, nWalker=50)

A sampler built on the emcee package.

The emcee packages implements the Goodman and Weare (2010) affine-invariant sampler. This is often an efficient sampler to use when selecting a proposal distribution is not simple.

Parameters:

• probfn (Callable) – Callable that accepts sampling parameters p and returns posterior probability.

• initializer (Callable) – Callable that accepts number of desired initial samples and returns samples. Note, the initial samples should be (at least mostly) unique.

• nWalker (int) – Number of ‘walkers’ to use in the Goodman & Weare algorithm. This should generally be at least 12 for the 6-dimensional problem.

probfn

initializer

nWalker

sample(nBurn, nStep)

Generate samples, first discarding nBurn steps, and then keeping nStep steps.

sample(nBurn=1000, nStep=500)

Generate samples.

Parameters:

• nBurn (int) – Number of initial steps to take but discard.

• nStep (int) – Number of subsequent steps to keep and return.

Returns:

• chain (array (nStep, nWalker, 6)) – Generated samples. Columns are [x, y, z, vx, vy, vz].

• lnprob (array (nStep, nWalker)) – The log posterior probabilities of the samples.

• lnprior (array(nStep, nWalker)) – The log prior values for each step.

class ssapy.rvsampler.EquinoctialExEyPrior(sigma)

Gaussian prior on equinoctial ex and ey, favoring ex = ey = 0.

Parameters:

sigma (float) – standard deviation on ex and ey

sigma

__call__(orbit)

Returns log prior probability of given orbit.

class ssapy.rvsampler.GEOProjectionInitializer(arc, rSigma=np.float64(4216417.236450565), vSigma=np.float64(307.46599995993046))

Initialize position and velocity samples by projecting an observed ra/dec into the equatorial plane and assuming a GEO statationary orbital velocity. Samples are dispersed from this given position and velocity using an isotropic Gaussian distribution.

Parameters:

• arc (QTable with columns for 'ra', 'dec', 'rStation_GCRF'.) – The set of observations used to initialize samples. Note, only the first row of the given QTable will be used.

• rSigma (float, optional) – Standard deviation of sample positions dispersion in meters.

• vSigma (float, optional) – Standard deviation of sample velocity dispersion in meters per second.

rSigma

vSigma

observation

First row of input arc.

Type:

QTable

__call__(nSample)

Returns nSample initial samples.

class ssapy.rvsampler.GaussianRVInitializer(rMean, vMean, rSigma=np.float64(4216417.236450565), vSigma=np.float64(307.46599995993046))

Generate position and velocity samples as an isotropic Gaussian distribution around an asserted position and velocity.

Parameters:

• rMean (array_like (3,)) – Mean of sample positions in meters.

• vMean (array_like (3,)) – Mean of sample velocities in meters per second.

• rSigma (float, optional) – Standard deviation of sample positions dispersion in meters.

• vSigma (float, optional) – Standard deviation of sample velocity dispersion in meters per second.

rMean

vMean

rSigma

vSigma

__call__(nSample)

Returns nSample initial samples.

class ssapy.rvsampler.LeastSquaresOptimizer(probfn, initparam, translatorcls, absstep=None, fracstep=None, **kw)

Base class for LeastSquaresOptimizers

optimize(**fit_kws)

Run the optimizer and return the resulting fit parameters.

Returns:

fit – Least-squares fit as [x, y, z, vx, vy, vz] in meters, and meters/second.

Return type:

(6,) array_like

class ssapy.rvsampler.Log10AreaPrior(mean=-0.5, sigma=1.5)

Gaussian prior on log10(area) of object.

Parameters:

• mean (float) – mean of log10(area)

• sigma (float) – standard deviation of log10(area)

mean, sigma

__call__(orbit)

Returns log prior probability of given orbit.

class ssapy.rvsampler.MHSampler(probfn, initializer, proposer, nChain)

Generate MCMC samples using a Metropolis-Hastings sampler.

Parameters:

• probfn (Callable) – Callable that accepts sampling parameters p and returns posterior probability.

• initializer (Callable) – Callable that accepts number of desired initial samples and returns samples. Note, the initial samples should be (at least mostly) unique.

• proposer (Callable) – Callable that accepts a “current” sample and returns a “proposed” sample to either accept or reject at each step.

• nChain (int) – Number of independent chains to use.

probfn

initializer

proposer

nChain

chain

MCMC sample chain.

Type:

array_like (nStep, nChain, 6)

lnprob

Log posterior probability of sample chain.

Type:

array_like (nStep, nChain)

nAccept

Total number of accepted proposals.

Type:

int

nStep

Total number of proposal steps.

Type:

int

reset()

Reset chains.

sample(nBurn, nStep)

Generate samples, first discarding nBurn steps, and then keeping nStep steps.

property acceptanceRatio

Ratio of accepted to proposed steps.

reset()

Reset chain, including chain, lnprob, nAccept, and nStep attributes.

sample(nBurn=1000, nStep=500)

Generate samples.

Parameters:

• nBurn (int) – Number of initial steps to take but discard.

• nStep (int) – Number of subsequent steps to keep and return.

Returns:

• chain (array (nStep, nChain, 6)) – Generated samples. Columns are [x, y, z, vx, vy, vz].

• lnprob (array (nStep, nChain)) – The log posterior probabilities of the samples.

• lnprior (array (nStep, nChain)) – The log prior probabilities of the samples.

class ssapy.rvsampler.MVNormalProposal(cov)

A multivariate normal proposal distribution.

Params

covarray_like(6, 6)

Covariance of multivariate normal distribution from which to sample. The order of variables is [x, y, z, vx, vy, vz] in meters and meters per second.

cov

propose(p)

Generate a new proposal.

propose(p)

Generate a proposal.

Parameters:

p (array (6,)) – Mean around which to generate a proposal.

Returns:

newp – Generated proposal.

Return type:

array (6,)

class ssapy.rvsampler.ParamOrbitAngle(initparam, epoch, initObsPos, initObsVel, **kwargs)

class ssapy.rvsampler.ParamOrbitEquinoctial(*args, **kwargs)

class ssapy.rvsampler.ParamOrbitRV(*args, **kwargs)

class ssapy.rvsampler.ParamOrbitTranslator(initparam, epoch, fixed=None, orbitattr=None)

Class for making parameters into Orbits and vice-versa.

class ssapy.rvsampler.RPrior(rmean, rsigma)

Gaussian prior on distance from origin |r|.

Parameters:

• rmean (float) – Prior mean on |r| in meters.

• rsigma (float) – Prior standard deviation on |r| in meters.

rMean

rSigma

__call__(orbit)

Returns log prior probability of given orbit.

class ssapy.rvsampler.RVProbability(arc, epoch, priors=[<ssapy.rvsampler.RPrior object>, <ssapy.rvsampler.APrior object>], propagator=KeplerianPropagator(), meanpm=False, damp=-1)

A class to manage MCMC sampling of orbits (parameterized by position and velocity at a given epoch) given some angular observations.

Parameters:

• arc (QTable with one row per observation and columns:) –

  ‘ra’, ‘dec’Angle

  Observed angular position in ICRS (topocentric).

  ’rStation_GCRF’Quantity

  Position of observer in GCRF.

  ’sigma’Quantity

  Isotropic uncertainty in observed angular position.

  ’time’float or astropy.time.Time

  Time of observation. If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

  The arc is a linked set of observations assumed to be of the same object.

• epoch (float or astropy.time.Time) – The time at which the model parameters (position and velocity) are specified. It probably makes sense to make this the same as one of the observation times, but it’s not required to be. If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• priors (list of priors, optional) – A list of class instances similar to RPrior() or APrior() to apply jointly to the posterior probability.

• propagator (class, optional) – The propagator class to use.

• meanpm (bool) – fit model to mean positions and proper motions rather than to individual epochs along each streak.

• damp (float) – damp chi residuals with pseudo-Huber loss function. Default -1: do not damp.

arc

epoch

priors

propagator

lnlike(orbit)

Calculate log likelihood of observations given orbit.

lnprob(p)

Calculate log posterior probability of orbit parameters p given observations (up to a constant). Parameters p are [x, y, z, vx, vy, vz] at epoch.

__call__(p)

Alias for lnprob(p)

chi(orbit)

Calculate chi residuals for use with lmfit.

This is essentially the input to a chi^2 statistic, just without the squaring and without the summing over observations. I.e., it’s

[(data_i - model_i)/err_i for i in 0..nobs]

Parameters:

orbit (Orbit) – Orbit in question. The model.

Returns:

• chilike (array_like (nobs,)) – chi residual array for likelihoods

• chiprior (array_like (nprior,)) – chi residual array for priors

lnlike(orbit)

Calculate the log likelihood of the observations given an orbit.

Parameters:

orbit (Orbit) – Orbit in question.

Returns:

lnlike – Log likelihood.

Return type:

float

lnprior(orbit)

Calculate the log prior probability of the observations given an orbit.

Parameters:

orbit (Orbit) – Orbit in question.

Returns:

lnprior – Log prior probability.

Return type:

float

lnprob(p)

Calculate the log posterior probability of parameters p given observations.

Parameters:

p (array (6,)) – Sampling parameters. [x, y, z, vx, vy, vz] in meters and meters per second. Understood to be the parameters at epoch.

Returns:

• lnprob (float) – Log posterior probability.

• lnprior (float) – Log prior probability.

class ssapy.rvsampler.RVSigmaProposal(rSigma, vSigma)

Isotropic multivariate normal proposal distribution.

Params

rSigmafloat

Isotropic standard deviation in position to use.

vSigmafloat

Isotropic standard deviation in velocity to use.

cov

propose(p)

Generate a new proposal.

class ssapy.rvsampler.VPrior(vmean, vsigma)

Gaussian prior on velocity magnitude |v|.

Parameters:

• vmean (float) – Prior mean on |v| in meters per second.

• vsigma (float) – Prior standard deviation on |v| in meters per second.

vMean

vSigma

__call__(orbit)

Returns log prior probability of given orbit.

ssapy.rvsampler.circular_guess(arc, indices=[0, 1])

Guess position and velocity from two closely spaced (in time) observations of ra/dec.

The state is inferred by finding the circular orbit passing through the two points, assuming the motion in between the observations is inertial.

Parameters:

• arc (QTable with columns for 'ra', 'dec', 'rStation_GCRF'.) – The set of observations used to make the guess. Note, only indices rows of the given QTable will be used.

• indices (list of indices indicating which two observations in arc to use) – for initialization. Default: [0, 1]; i.e., use the first two observations.

Returns:

state – state corresponding to orbit passing through points

Return type:

array_like, (6,)

ssapy.rvsampler.damper(chi, damp)

Pseudo-Huber loss function.

ssapy.rvsampler.damper_deriv(chi, damp, derivnum=1)

Derivative of the pseudo-Huber loss function.

ssapy.rvsampler.sample_ball(p0, std, nSample=1)

Produce a ball of walkers around an initial parameter value.

Parameters:

• p0 (array_like, (d,)) – Central parameter values.

• std (array_like, (d,)) – Axis-aligned standard deviations.

• nSample (int, optional) – The number of samples to produce.

Returns:

samples – The samples.

Return type:

array_like (nSample, d)

Simple

Functions to simplify certain aspects of SSAPy. For when you want a quick answer and do not need want to use high fidelity defaults.

ssapy.simple.best_prop(integration_timestep=10, kwargs={'CD': 2.3, 'CR': 1.3, 'area': 0.022, 'mass': 250})

Create and return an RK78Propagator with a comprehensive set of accelerations.

Parameters:

integration_timestepfloat, optional

The time step for the RK78Propagator integration (default is 10 seconds).

kwargsdict, optional

Dictionary of parameters for non-conservative forces (mass, area, CD, CR). If not provided, defaults are used.

Returns:

RK78Propagator

An instance of RK78Propagator configured with cached accelerations.

ssapy.simple.fourbody_prop(integration_timestep: float = 10) → RK78Propagator

Create and return an RK78Propagator with a set of accelerations for a four-body problem.

The four bodies considered are Earth (Keplerian effect), Moon, Sun, and the Earth itself.

Parameters:

integration_timestepfloat, optional

The time step for the RK78Propagator integration (default is 10 seconds).

Returns:

RK78Propagator

An instance of RK78Propagator configured with the four-body accelerations.

ssapy.simple.get_similar_orbits(r0, v0, rad=100000.0, num_orbits=4, duration=(90, 'days'), freq=(1, 'hour'), start_date='2025-1-1', mass=250)

Generate similar orbits by varying the initial position.

Parameters:

r0array_like

Initial position vector of shape (3,).

v0array_like

Initial velocity vector of shape (3,).

radfloat

Radius of the circle around the initial position to generate similar orbits.

num_orbitsint

Number of similar orbits to generate.

durationtuple

Duration of the orbit simulation.

freqtuple

Frequency of output data.

start_datestr

Start date for the simulation.

massfloat

Mass of the satellite.

Returns:

trajectoriesndarray

Stacked array of shape (3, n_times, num_orbits) containing the trajectories.

ssapy.simple.keplerian_prop(integration_timestep: float = 10) → RK78Propagator

Create and return an RK78Propagator for a Keplerian orbit.

This propagator uses only the Keplerian acceleration for a two-body problem.

Parameters:

integration_timestepfloat, optional

The time step for the RK78Propagator integration (default is 10 seconds).

Returns:

RK78Propagator

An instance of RK78Propagator configured with Keplerian acceleration.

ssapy.simple.ssapy_kwargs(mass=250, area=0.022, CD=2.3, CR=1.3)

Generate a dictionary of default parameters for a space object used in simulations.

Parameters:

massfloat, optional

Mass of the object in kilograms (default is 250 kg).

areafloat, optional

Cross-sectional area of the object in square meters (default is 0.022 m^2).

CDfloat, optional

Drag coefficient of the object (default is 2.3).

CRfloat, optional

Radiation pressure coefficient of the object (default is 1.3).

Returns:

dict

A dictionary containing the parameters for the space object.

ssapy.simple.ssapy_orbit(orbit=None, a=None, e=0, i=0, pa=0, raan=0, ta=0, r=None, v=None, duration=(30, 'day'), freq=(1, 'hr'), t0=<Time object: scale='utc' format='iso' value=2025-01-01 00:00:00.000>, t=None, prop=RK78Propagator(<ssapy.accel.AccelSum object>, 60, (1e-06, 1e-06, 1e-06, 1e-09, 1e-09, 1e-09)))

Compute the orbit of a spacecraft given either Keplerian elements or position and velocity vectors.

Parameters: - orbit (Orbit, optional): An Orbit object if you already have an orbit defined. - a (float, optional): Semi-major axis of the orbit in meters. - e (float, optional): Eccentricity of the orbit (default is 0, i.e., circular orbit). - i (float, optional): Inclination of the orbit in degrees. - pa (float, optional): Argument of perigee in degrees. - raan (float, optional): Right ascension of the ascending node in degrees. - ta (float, optional): True anomaly in degrees. - r (array-like, optional): Position vector in meters. - v (array-like, optional): Velocity vector in meters per second. - duration (tuple, optional): Duration of the simulation as a tuple (value, unit), where unit is ‘day’, ‘hour’, etc. Default is 30 days. - freq (tuple, optional): Frequency of the output as a tuple (value, unit), where unit is ‘day’, ‘hour’, etc. Default is 1 hour. - t0 (str, optional): Start date of the simulation in ‘YYYY-MM-DD’ format. Default is “2025-01-01”. - t (array-like, optional): Specific times at which to compute the orbit. If None, times will be generated based on duration and frequency. - prop (function, optional): A function to compute the perturbation effects. Default is ssapy_prop().

Returns: - r (array-like): Position vectors of the spacecraft at the specified times. - v (array-like): Velocity vectors of the spacecraft at the specified times. - t (array-like): Times at which the orbit was computed. Returned only if t was None.

Raises: - ValueError: If neither Keplerian elements nor position and velocity vectors are provided. - RuntimeError or ValueError: If an error occurs during computation.

ssapy.simple.ssapy_prop(integration_timestep=60, propkw={'CD': 2.3, 'CR': 1.3, 'area': 0.022, 'mass': 250})

Setup and return an RK78 propagator with specified accelerations and radiation pressure effects.

Parameters:

integration_timestepint

Time step for the numerical integration (in seconds).

propkwdict, optional

Keyword arguments for radiation pressure accelerations. If None, default arguments are used.

Returns:

RK78Propagator

An RK78 propagator configured with the specified accelerations and time step.

ssapy.simple.threebody_prop(integration_timestep: float = 10) → RK78Propagator

Create and return an RK78Propagator with a set of accelerations for a three-body problem.

The three bodies considered are Earth (Keplerian effect), Moon, and the Earth itself.

Parameters:

integration_timestepfloat, optional

The time step for the RK78Propagator integration (default is 10 seconds).

Returns:

RK78Propagator

An instance of RK78Propagator configured with the three-body accelerations.

Utilities

class ssapy.utils.LRU_Cache(user_function, maxsize=1024)

Simplified Least Recently Used Cache.

Mostly stolen from http://code.activestate.com/recipes/577970-simplified-lru-cache/, but added a method for dynamic resizing. The least recently used cached item is overwritten on a cache miss.

Parameters:

• user_function – A python function to cache.

• maxsize – Maximum number of inputs to cache. [Default: 1024]

Example:

>>> def slow_function(*args) # A slow-to-evaluate python function
>>>    ...
>>>
>>> v1 = slow_function(*k1)  # Calling function is slow
>>> v1 = slow_function(*k1)  # Calling again with same args is still slow
>>> cache = LRU_Cache(slow_function)
>>> v1 = cache(*k1)  # Returns slow_function(*k1), slowly the first time
>>> v1 = cache(*k1)  # Returns slow_function(*k1) again, but fast this time.

resize(maxsize)

Resize the cache.

Increasing the size of the cache is non-destructive, i.e., previously cached inputs remain in the cache. Decreasing the size of the cache will necessarily remove items from the cache if the cache is already filled. Items are removed in least recently used order.

Parameters:

maxsize – The new maximum number of inputs to cache.

ssapy.utils.catalog_to_apparent(ra, dec, t, observer=None, pmra=0, pmdec=0, parallax=0, skipAberration=False)

Convert ra/dec of stars from catalog positions (J2000/ICRS) to apparent positions, correcting for proper motion, parallax, annual and diurnal aberration.

Parameters:

• ra (array_like) – J2000 right ascension in radians.

• dec (array_like) – J2000 declination in radians.

• t (float or astropy.time.Time) – If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• observer (Observer, optional) – Observer to use for diurnal aberration correction. If None, then no diurnal aberration correction is performed.

• pmra (array_like, optional) – proper motion in right ascension / declination, in milliarcsec per year

• pmdec (array_like, optional) – proper motion in right ascension / declination, in milliarcsec per year

• parallax (array_like, optional) – annual parallax in arcseconds

• skipAberration (bool, optional) – Don’t apply aberration correction. Mostly useful during testing…

Returns:

ara, adec – Apparent RA and dec in radians

Return type:

array_like

ssapy.utils.cluster_emcee_walkers(chain, lnprob, lnprior, thresh_multiplier=1, verbose=False)

Down-select emcee walkers to those with the largest mean posteriors

Follows the algorithm of Hou, Goodman, Hogg et al. (2012)

ssapy.utils.continueClass(cls)

Re-open the decorated class, adding any new definitions into the original. For example: .. code-block:: python

class Foo:

pass

@continueClass class Foo:

def run(self):

return None

is equivalent to: .. code-block:: python

class Foo:

def run(self):

return None

Warning

Python’s built-in super function does not behave properly in classes decorated with continueClass. Base class methods must be invoked directly using their explicit types instead.

ssapy.utils.find_all_zeros(f, low, high, n=100)

Attempt to find all zeros of a function between given bounds. Uses n initial samples to characterize the zero-landscape of the interval.

ssapy.utils.find_extrema_brackets(fs)

Find triplets bracketing extrema from an ordered set of function evaluations.

ssapy.utils.find_file(filename, ext=None)

Find a file in the current directory or the ssapy datadir. If ext is not None, also try appending ext to the filename.

ssapy.utils.find_smallest_bounding_cube(r: ndarray, pad: float = 0.0) → Tuple[ndarray, ndarray]

Find the smallest bounding cube for a set of 3D coordinates, with optional padding.

Parameters: r (np.ndarray): An array of shape (n, 3) containing the 3D coordinates. pad (float): Amount to increase the bounding cube in all dimensions.

Returns: tuple: A tuple containing the lower and upper bounds of the bounding cube.

ssapy.utils.gcrf_to_teme(t)

Return the rotation matrix that converts GCRS cartesian coordinates to TEME cartesian coordinates.

Parameters:

t (float or astropy.time.Time) – If float or array of float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

Returns:

rot – Rotation matrix to apply to GCRS to yield TEME.

Return type:

array (3,3)

ssapy.utils.get_angle(a, b, c)

Calculate the angle between two vectors where b is the vertex of the angle.

This function computes the angle between vectors ba and bc, where b is the vertex and a and c are the endpoints of the angle.

Parameters:

a(n, 3) numpy.ndarray

Array of coordinates representing the first vector.

b(n, 3) numpy.ndarray

Array of coordinates representing the vertex of the angle.

c(n, 3) numpy.ndarray

Array of coordinates representing the second vector.

Returns:

numpy.ndarray

Array of angles (in radians) between the vectors ba and bc.

Notes:

• The function handles multiple vectors by using broadcasting.

• The angle is calculated using the dot product formula and the arccosine function.

Example:

>>> a = np.array([[1, 0, 0]])
>>> b = np.array([[0, 0, 0]])
>>> c = np.array([[0, 1, 0]])
>>> get_angle(a, b, c)
array([1.57079633])

ssapy.utils.get_kernel_cov(kernel_mat, weights)

Get the covariance matrix of kernel_mat. This a wrapper for numpy’s cov

Parameters:

• kernel_mat (numpy.ndarray) – Data matrix

• weights (numpy.ndarray) – 1D array of weights of the data

Returns:

numpy.ndarray

Return type:

Covariance matrix

ssapy.utils.get_times(duration: Tuple[int, str], freq: Tuple[int, str], t0: str | Time = '2025-01-01') → ndarray

Calculate a list of times spaced equally apart over a specified duration.

Parameters:

• duration (tuple) – A tuple containing the length of time and the unit (e.g., (30, ‘day’)).

• freq (tuple) – A tuple containing the frequency value and its unit (e.g., (1, ‘hr’)).

• t0 (str or Time, optional) – The starting time. Default is “2025-01-01”.

Returns:

A list of times spaced equally apart over the specified duration.

Return type:

np.ndarray

ssapy.utils.iers_interp(t)

Interpolate IERS values

Parameters:

t (float or astropy.time.Time) – If float or array of float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

Returns:

• dut1 (array_like (n)) – Time difference UT1 - TT in days.

• pmx, pmy (array_like (n)) – Polar motion values in radians.

ssapy.utils.interpolate_points_between(r, m)

Interpolates points between the given points.

Parameters:

• r – An (n, 3) numpy array of the original points.

• m – The number of points to interpolate between each pair of points.

Returns:

An (n * m) numpy array of the interpolated points.

ssapy.utils.isAttributeSafeToTransfer(name, value)

Return True if an attribute is safe to monkeypatch-transfer to another class. This rejects special methods that are defined automatically for all classes, leaving only those explicitly defined in a class decorated by continueClass or registered with an instance of TemplateMeta.

class ssapy.utils.lazy_property(fget)

meant to be used for lazy evaluation of an object attribute. property should represent non-mutable data, as it replaces itself.

ssapy.utils.lb_to_tan(lb, b, mul=None, mub=None, lcen=None, bcen=None)

Convert lb-like coordinates & proper motions to orthographic tangent plane.

Everything is in radians. If mul is None (default), transformed proper motions will not be returned. The tangent plane is always chosen so that +Y is towards b = 90, and +X is towards +lb.

Parameters:

• lb (array_like (n)) – right ascension of point

• b (array_like (n)) – declination of point

• mul (array_like (n)) – proper motion in ra of point (arbitrary units) rate of change in lb is mul / np.cos(b); i.e., length of proper motion vector on sphere is np.hypot(mul, mub)

• mub (array_like (n)) – proper motion in dec of point (arbitrary units)

• lcen (array_like (n)) – right ascension to use for center of tangent plane if None, use spherical mean of (lb, b)

• bcen (array_like (n)) – declination to use for center of tangent plane

Returns:

• if mul is None, (x, y) otherwise (x, y, vx, vy)

• x (array_like (n)) – x coordinate of tangent plane projection of lb, b

• y (array_like (n)) – y coordinate of tangent plane projection of lb, b

• vx (array_like (n)) – x coordinate of tangent plane projection of (mul, mub)

• vy (array_like (n)) – y coordinate of tangent plane projection of (mul, mub)

ssapy.utils.lb_to_tp(lb, b)

Convert lb-like coordinates to theta-phi like coordinates.

Here ‘theta-phi’ coordinates refers to the system where theta is the angle between zenith and the point in question, and phi is the corresponding azimuthal angle.

This just sets theta = pi - b and renames lb -> phi. Everything is in radians.

ssapy.utils.lb_to_unit(r, d)

Convert lb-like coordinates to unit vectors.

Everything is in radians.

ssapy.utils.moonPos(t)

Compute GCRF position of the moon.

Parameters:

t (float or astropy.time.Time) – If float or array of float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

Returns:

r – position in meters

Return type:

array_like (n)

ssapy.utils.newton_raphson(guess, f, fprime=None, eps=3e-16, maxiter=100, **kwargs)

Find a root of a univariate function with known gradient using Newton-Raphson iterations.

Parameters:

• guess (float) – Initial guess.

• f (callable function of 1 argument.) – Function for which to find a zero. If fprime is None, then the return value of f should be a 2-tuple with the value of the function and the first derivative. If fprime is given separately, then f should just return the value of the function.

• fprime (callable function of 1 argument, optional) – Derivative of f. Default None.

• eps (float, optional) – Absolute tolerance for finding a zero. Default 3e-16.

• maxiter (int, optional) – Maximum number of iterations to try. Default 100.

Returns:

solution

Return type:

float

ssapy.utils.ntw_to_r(r, v, ntw, relative=False)

Convert NTW coordinates to cartesian coordinates, using r, v to define NTW system.

T gives the projection of (rcoord - r) along V (tangent to track) W gives the projection of (rcoord - r) along (V cross r) (normal to plane) N gives the projection of (rcoord - r) along (V cross (V cross r))

(in plane, perpendicular to T)

Parameters:

• r (array_like (n, 3)) – central positions defining coordinate system

• v (array_like (n, 3)) – velocity defining coordinate system

• ntw (array_like (n, 3)) – ntw coordinates to transform to cartesian coordinates

• relative (bool) – if True, just rotate the NTW coordinates to Cartesian; do not offset the origin so that NTW = 0 -> Cartesian r.

Returns:

r – cartesian x, y, z coordinates

Return type:

array_like (n, 3)

ssapy.utils.num_wraps(ang)

Return number of times angle ang has wrapped around.

Returns:

Number of wraps

Return type:

int

ssapy.utils.perpendicular_vectors(v)

Returns two vectors that are perpendicular to v and each other.

ssapy.utils.regularize_default(particles, weights, num_particles_out=None, dimension=6)

Perform particle regularization. This generates a perturbation from the particles’ original values to prevent singularity issues. From “Bayesian Multiple Target Tracking” (2nd ed.) p. 101-102

Parameters:

• particles (numpy.ndarray) – Particles to reqularize. Each row is an nD particle

• weights (numpy.ndarray) – 1D array of particle weight. Should be same length as number of particles.

• num_particles_out (int | None, optional) – Number of particles to return, defaults to None. If not specified, will return the same number of particles as in the input particles.

• dimension (int, optional) – Dimension of the parameter space for resampling. Assumes the first dimension columns of particles are the parameters to use. Any remaining columns are carried through without modification (e.g., time columns). Default: 6

Returns:

Deltas from original particles and their weights

Return type:

(numpy.ndarray, numpy.ndarry)

ssapy.utils.resample(particles, ln_weights, obs_times=None, pod=False)

Resample particles to achieve more uniform weights. This produces the same number of particles as are input (i.n. dimension matches “particles”) From “Bayesian Multiple Target Tracking” (2nd ed.) p. 100

Parameters:

• particles (numpy.ndarray) – 2D array of particles. Each row is an nD particles

• ln_weights (numpy.ndarray) – 1D array of particle log-weights. Length should match number of particles

• obs_times ([type], optional) – Time at which observation occured. This is only used if pod=True, defaults to None

• pod (bool, optional) – Specify whether this is the special case of regularizing orbit parameters, defaults to False

Returns:

Resampled particles and weights (same number as input)

Return type:

(numpy.ndarray, numpy.ndarray)

ssapy.utils.rotate_points_3d(points, axis=array([0, 0, 1]), theta=-1.5707963267948966)

Rotate a set of 3D points about a 3D axis by an angle theta in radians.

Parameters:

• points (np.ndarray) – The set of 3D points to rotate, as an Nx3 array.

• axis (np.ndarray) – The 3D axis to rotate about, as a length-3 array. Default is the z-axis.

• theta (float) – The angle to rotate by, in radians. Default is pi/2.

Returns:

The rotated set of 3D points, as an Nx3 array.

Return type:

np.ndarray

ssapy.utils.rotation_matrix_from_vectors(vec1, vec2)

Find the rotation matrix that aligns vec1 to vec2 :param vec1: A 3d “source” vector :param vec2: A 3d “destination” vector :return mat: A transform matrix (3x3) which when applied to vec1, aligns it with vec2.

ssapy.utils.rv_to_ntw(r, v, rcoord)

Convert coordinates to NTW coordinates, using r, v to define NTW system.

T gives the projection of (rcoord - r) along V (tangent to track) W gives the projection of (rcoord - r) along (V cross r) (normal to plane) N gives the projection of (rcoord - r) along (V cross (V cross r))

(in plane, perpendicular to T)

Parameters:

• r (array_like (n, 3)) – central positions defining coordinate system

• v (array_like (n, 3)) – velocity defining coordinate system

• rcoord (array_like (n, 3)) – positions to transform to NTW coordinates

Returns:

ntw – n, t, w positions

Return type:

array_like (n, 3)

ssapy.utils.sample_points(x, C, npts, sqrt=False)

Sample points around x according to covariance matrix.

Parameters:

• x (array_like (n)) – point to sample around

• C (array_like (n, n)) – covariance matrix corresponding to x

• npts (int) – number of points to sample

• sqrt (bool) – use sqrt(C) rather than an SVD. The SVD is often more stable.

Return type:

Gaussian samples around x corresponding to covariance matrix C

ssapy.utils.sigma_points(f, x, C, scale=1, fixed_dimensions=None)

Compute f(sigma points) for sigma points of C around x.

There are many possible definitions of sigma points. This one takes the singular value decomposition of C and uses the eigenvectors times sqrt(dimension)*scale*(+/-1) as the sigma points. It then evaluates the given function f at x plus those sigma points.

Parameters:

• f (function) – the function to evaluate at the sigma points

• x (array_like (n)) – the central value to evaluate the function around

• C (array_like (n, n)) – the covariance matrix corresponding to x

• scale (float) – return scale-sigma points rather than n-sigma points. e.g., for 5 sigma, set scale = 5.

• fixed_dimensions (array_like, (n), bool) – boolean array specifying dimensions of x that are fixed and not specified in C

ssapy.utils.subsample_high_lnprob(chain, lnprob, lnprior, nSample, thresh=-10)

Select MCMC samples with probabilities above some relative threshold

Parameters:

• chain (array_like, (nWalker, nStep, 6)) – Input MCMC chain from EmceeSampler.sample()

• lnprob (array_like, (nWalker, nStep)) – Input MCMC lnprob from EmceeSampler.sample()

• lnprior (array_like, (nWalker, nStep)) – Input MCMC lnprior from EmceeSampler.sample()

• nSample (int) – Number of samples to return. If fewer than nSample samples exceed threshold, then return the nSample most probable samples.

• thresh (float, optional) – Threshold w.r.t. max(lnprob) below which to exclude samples.

Returns:

samples – Output samples

Return type:

array_like, (nSample, 6)

ssapy.utils.sunPos(t, fast=True)

Compute GCRF position of the sun.

Parameters:

• t (float or astropy.time.Time) – If float or array of float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

• fast (bool) – Use fast approximation?

Returns:

r – position in meters

Return type:

array_like (n)

ssapy.utils.tan_to_lb(xx, yy, lcen, bcen)

Convert orthographic tangent plane coordinates to lb coordinates.

Parameters:

• xx (array_like (n)) – tangent plane x coordinate (radians)

• yy (array_like (n)) – targent plane y coordinate (radians

• lcen (float, array_like (n)) – right ascension of center of tangent plane

• bcen (float, array_linke (n)) – declination of center of tangent plane

Returns:

• lb (array_like (n)) – right ascension corresponding to xx, yy

• b (array_like (n)) – declination corresponding to xx, yy

ssapy.utils.teme_to_gcrf(t)

Return the rotation matrix that converts TEME cartesian coordinates to GCRS cartesian coordinates.

Parameters:

t (float or astropy.time.Time) – If float or array of float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC

Returns:

rot – Rotation matrix to apply to TEME to yield GCRS.

Return type:

array (3,3)

ssapy.utils.tp_to_lb(t, p)

Convert theta-phi-like coordinates to lb-like coordinates.

Here ‘theta-phi’ coordinates refers to the system where theta is the angle between zenith and the point in question, and phi is the corresponding azimuthal angle.

This just sets b = pi - theta and renames phi -> lb. Everything is in radians.

ssapy.utils.tp_to_unit(t, p)

Convert theta-phi-like coordinates to unit vectors.

Here ‘theta-phi’ coordinates refers to the system where theta is the angle between zenith and the point in question, and phi is the corresponding azimuthal angle.

Everything is in radians.

ssapy.utils.unitAngle3(r1, r2)

Robustly compute angle between unit vectors r1 and r2. Vectorized for multiple triplets.

ssapy.utils.unit_to_lb(unit)

Convert unit vectors to lb-like coordinates.

Everything is in radians.

ssapy.utils.unit_to_tp(unit)

Convert unit vectors to theta-phi-like coordinates.

Here ‘theta-phi’ coordinates refers to the system where theta is the angle between zenith and the point in question, and phi is the corresponding azimuthal angle.

Everything is in radians.

ssapy.utils.unit_vector(vector)

Returns the unit vector of the vector.

ssapy.utils.unscented_transform_mean_covar(f, x, C, scale=1)

Compute mean and covar using unscented transform given a transformation f, a point x, and a covariance C.

This uses the sigma point convention from sigma_points. It assumes that f(sigma_points)[i] is f evaluated at the ith sigma point. If f does not obey this convention, this function will produce undefined results.

Parameters:

• f (function) – the function to evaluate at the sigma points.

• x (array_like (n)) – the central value to evaluate the function around

• C (array_like (n, n)) – the covariance matrix corresponding to x

• scale (float) – return scale-sigma points rather than n-sigma points. e.g., for 5 sigma, set scale = 5.

ssapy.utils.xyz_to_lb(x, y, z)

Convert x, y, z vectors to lb-like coordinates.

Everything is in radians.

ssapy.utils.xyz_to_tp(x, y, z)

Convert x, y, z vectors to theta-phi-like coordinates.

Here ‘theta-phi’ coordinates refers to the system where theta is the angle between zenith and the point in question, and phi is the corresponding azimuthal angle.

Everything is in radians.