From 4451c5a363800d3139eb14c4b65f9ce6e0bd9fd0 Mon Sep 17 00:00:00 2001
From: Travis Yeager <43413776+SuperdoerTrav@users.noreply.github.com>
Date: Mon, 28 Oct 2024 15:08:10 -0700
Subject: [PATCH 01/14] Update simple.py

added a condition for when t is not None, ensuring correct start time.
---
 ssapy/simple.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/ssapy/simple.py b/ssapy/simple.py
index f6a1e40..e7dd549 100644
--- a/ssapy/simple.py
+++ b/ssapy/simple.py
@@ -217,6 +217,8 @@ def ssapy_orbit(orbit=None, a=None, e=0, i=0, pa=0, raan=0, ta=0, r=None, v=None
     """
     if t is None:
         t = get_times(duration=duration, freq=freq, t=t0)
+    else:
+        t0 = t[0]
 
     if orbit is not None:
         pass

From 65b80e2145f3922c012cdff9f25b4abcf7af048f Mon Sep 17 00:00:00 2001
From: Travis Yeager <43413776+SuperdoerTrav@users.noreply.github.com>
Date: Fri, 6 Dec 2024 10:57:28 -0800
Subject: [PATCH 02/14] Update get_Times utils.py

added type annotation and changed t to t0 to better indicate it's an initial time, not a list of times.
---
 ssapy/utils.py | 25 ++++++++++++-------------
 1 file changed, 12 insertions(+), 13 deletions(-)

diff --git a/ssapy/utils.py b/ssapy/utils.py
index f7eaccd..35e8194 100644
--- a/ssapy/utils.py
+++ b/ssapy/utils.py
@@ -1978,27 +1978,26 @@ def dd_to_hms(degree_decimal):
     return f'{int(_h)}:{int(_m)}:{_s}'
 
 
-def get_times(duration, freq, t):
+def get_times(duration: Tuple[int, str], freq: Tuple[int, str], t0: Union[str, Time] = "2025-01-01") -> np.ndarray:
     """
     Calculate a list of times spaced equally apart over a specified duration.
 
     Parameters
     ----------
-    duration: int
-        The length of time to calculate times for.
-    freq: int, unit: str
-        frequency of time outputs in units provided
-    t: ssapy.utils.Time, optional
+    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".
-    example input:
-    duration=(30, 'day'), freq=(1, 'hr'), t=Time("2025-01-01", scale='utc')
+
     Returns
     -------
-    times: array-like
+    np.ndarray
         A list of times spaced equally apart over the specified duration.
     """
-    if isinstance(t, str):
-        t = Time(t, scale='utc')
+    if isinstance(t0, str):
+        t0 = Time(t0, scale='utc')
     unit_dict = {'second': 1, 'sec': 1, 's': 1, 'minute': 60, 'min': 60, 'hour': 3600, 'hr': 3600, 'h': 3600, 'day': 86400, 'd': 86400, 'week': 604800, 'month': 2630016, 'mo': 2630016, 'year': 31557600, 'yr': 31557600}
     dur_val, dur_unit = duration
     freq_val, freq_unit = freq
@@ -2014,7 +2013,7 @@ def get_times(duration, freq, t):
     freq_seconds = freq_val * unit_dict[freq_unit.lower()]
     timesteps = int(dur_seconds / freq_seconds) + 1
 
-    times = t + np.linspace(0, dur_seconds, timesteps) / unit_dict['day'] * u.day
+    times = t0 + np.linspace(0, dur_seconds, timesteps) / unit_dict['day'] * u.day
     return times
 
 
@@ -2149,4 +2148,4 @@ def points_on_circle(r, v, rad, num_points=4):
     # Translate the rotated points to the center point r
     points_rotated = rotated_points + r.reshape(1, 3)
 
-    return points_rotated
\ No newline at end of file
+    return points_rotated

From 4889a639dff1309d9fcc61bbe521ae162d0fa5e7 Mon Sep 17 00:00:00 2001
From: Travis Yeager <43413776+SuperdoerTrav@users.noreply.github.com>
Date: Fri, 6 Dec 2024 11:03:19 -0800
Subject: [PATCH 03/14] Update simple.py

Changed how t is None is handled, added further methods to make sure what is returned is correct dependent on inputs.
---
 ssapy/simple.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/ssapy/simple.py b/ssapy/simple.py
index e7dd549..0641b2f 100644
--- a/ssapy/simple.py
+++ b/ssapy/simple.py
@@ -215,10 +215,12 @@ def ssapy_orbit(orbit=None, a=None, e=0, i=0, pa=0, raan=0, ta=0, r=None, v=None
     - ValueError: If neither Keplerian elements nor position and velocity vectors are provided.
     - RuntimeError or ValueError: If an error occurs during computation.
     """
+    t0 = Time(t0, scale='utc')
     if t is None:
-        t = get_times(duration=duration, freq=freq, t=t0)
+        time_is_None = True
+        t = get_times(duration=duration, freq=freq, t0=t0)
     else:
-        t0 = t[0]
+        time_is_None = False
 
     if orbit is not None:
         pass
@@ -232,7 +234,10 @@ def ssapy_orbit(orbit=None, a=None, e=0, i=0, pa=0, raan=0, ta=0, r=None, v=None
 
     try:
         r, v = rv(orbit, t, prop)
-        return r, v, t
+        if time_is_None:
+            return r, v, t
+        else:
+            return r, v
     except (RuntimeError, ValueError) as err:
         print(err)
         return np.nan, np.nan, np.nan

From 5b963c90533ddf865ab68190cb2d8195a1f26905 Mon Sep 17 00:00:00 2001
From: Travis Yeager <43413776+SuperdoerTrav@users.noreply.github.com>
Date: Fri, 6 Dec 2024 11:06:37 -0800
Subject: [PATCH 04/14] including type annotation imports

---
 ssapy/simple.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/ssapy/simple.py b/ssapy/simple.py
index 0641b2f..177d132 100644
--- a/ssapy/simple.py
+++ b/ssapy/simple.py
@@ -14,6 +14,7 @@
 
 from astropy.time import Time
 import numpy as np
+from typing import Union, List, Tuple
 
 
 def keplerian_prop(integration_timestep: float = 10) -> RK78Propagator:

From d4d31fbb7a8f9b32e3c5eeff14e1aa6663f39dec Mon Sep 17 00:00:00 2001
From: Travis Yeager <43413776+SuperdoerTrav@users.noreply.github.com>
Date: Fri, 6 Dec 2024 11:09:01 -0800
Subject: [PATCH 05/14] Including typing imports for adding type annotations

---
 ssapy/utils.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/ssapy/utils.py b/ssapy/utils.py
index 35e8194..e2b7e78 100644
--- a/ssapy/utils.py
+++ b/ssapy/utils.py
@@ -3,6 +3,7 @@
 import numpy as np
 from astropy.time import Time
 import astropy.units as u
+from typing import Union, List, Tuple
 
 from . import datadir
 from .constants import RGEO, EARTH_RADIUS, MOON_RADIUS, WGS84_EARTH_OMEGA

From d1e05e9d1fe1b38930d16aeafd7382bff3c74fd8 Mon Sep 17 00:00:00 2001
From: Travis Yeager <43413776+SuperdoerTrav@users.noreply.github.com>
Date: Thu, 2 Jan 2025 15:14:16 -0800
Subject: [PATCH 06/14] Update orbit_plot

Removed the limits keyword and improved the automatic bounding box
Changed the Lagrange Point colors to white
Moved the second sub plots y axis labels to the right to prevent clipping issues.
---
 ssapy/plotUtils.py | 54 ++++++++++++++++++++++++----------------------
 1 file changed, 28 insertions(+), 26 deletions(-)

diff --git a/ssapy/plotUtils.py b/ssapy/plotUtils.py
index 337b67c..85dd32f 100644
--- a/ssapy/plotUtils.py
+++ b/ssapy/plotUtils.py
@@ -213,7 +213,7 @@ def check_numpy_array(variable):
         return "not numpy"
 
 
-def orbit_plot(r, t=[], limits=False, title='', figsize=(7, 7), save_path=False, frame="gcrf", show=False):
+def orbit_plot(r, t=[], title='', figsize=(7, 7), save_path=False, frame="gcrf", show=False):
     """
     Plot 2D and 3D projections of the orbit of one or more objects.
 
@@ -223,8 +223,6 @@ def orbit_plot(r, t=[], limits=False, title='', figsize=(7, 7), save_path=False,
         Position of orbiting object(s) in meters. If a list is provided, each element represents the orbit of a different object.
     t : array_like, optional
         Time array corresponding to the position vectors. Required for frames "itrf", "lunar", or "lunar fixed".
-    limits : float or tuple, optional
-        Limits for the x and y axes of the 2D plots. If not provided, it is calculated based on the data.
     title : str, optional
         Title of the plot. Default is an empty string.
     figsize : tuple of int, optional
@@ -259,9 +257,10 @@ def orbit_plot(r, t=[], limits=False, title='', figsize=(7, 7), save_path=False,
     r = np.array([[1e7, 2e7, 3e7], [4e7, 5e7, 6e7]])  # Replace with actual data
     t = np.linspace(0, 10, len(r))  # Replace with actual time data
 
-    fig, axes = orbit_plot(r, t, limits=1e8, title='Orbit Plot', frame='gcrf', show=True)
+    fig, axes = orbit_plot(r, t, title='Orbit Plot', frame='gcrf', show=True)
     """
-    def _make_scatter(fig, ax1, ax2, ax3, ax4, r, t, limits, title='', orbit_index='', num_orbits=1, frame=False):
+
+    def _make_scatter(fig, ax1, ax2, ax3, ax4, r, t, title='', orbit_index='', num_orbits=1, frame=False):
         if np.size(t) < 1:
             if frame in ["itrf", "lunar", "lunar_fixed"]:
                 raise("Need to provide t for itrf, lunar or lunar fixed frames")
@@ -312,7 +311,13 @@ def get_main_category(frame):
         xm = r_moon[:, 0] / RGEO
         ym = r_moon[:, 1] / RGEO
         zm = r_moon[:, 2] / RGEO
-            
+
+        if orbit_index == '' or orbit_index == 0:
+            lower_bound = np.array([np.inf, np.inf, np.inf])
+            upper_bound = np.array([-np.inf, -np.inf, -np.inf])
+        lower_bound_temp, upper_bound_temp = find_smallest_bounding_cube(r / RGEO * 1.2)
+        lower_bound = np.minimum(lower_bound, lower_bound_temp)
+        upper_bound = np.maximum(upper_bound, upper_bound_temp)
         if np.size(xm) > 1:
             gradient_colors = cm.Greys(np.linspace(0, .8, len(xm)))[::-1]
             blues = cm.Blues(np.linspace(.4, .9, len(xm)))[::-1]
@@ -329,11 +334,7 @@ def get_main_category(frame):
             stn = plot_settings[frame]
         except KeyError:
             raise ValueError("Unknown plot type provided. Accepted: 'gcrf', 'itrf', 'lunar', 'lunar fixed'")
-        if limits is False:
-            lower_bound, upper_bound = find_smallest_bounding_cube(r / RGEO)
-            lower_bound = lower_bound * 1.2
-            upper_bound = upper_bound * 1.2
-
+        
         if orbit_index == '':
             angle = 0
             scatter_dot_colors = cm.rainbow(np.linspace(0, 1, len(x)))
@@ -358,8 +359,8 @@ def get_main_category(frame):
                     pass
                 else:
                     pos[0] = pos[0] - LD / RGEO
-                ax1.scatter(pos[0], pos[1], color=color, label=point)
-                ax1.text(pos[0], pos[1], point, color=color)
+                ax1.scatter(pos[0], pos[1], color='w', label=point)
+                ax1.text(pos[0], pos[1], point, color='w')
 
         ax2.add_patch(plt.Circle((0, 0), stn[2], color='white', linestyle='dashed', fill=False))
         ax2.scatter(x, z, color=scatter_dot_colors, s=1)
@@ -369,6 +370,8 @@ def get_main_category(frame):
         ax2.set_aspect('equal')
         ax2.set_xlabel('x [GEO]')
         ax2.set_ylabel('z [GEO]')
+        ax2.yaxis.tick_right()  # Move y-axis ticks to the right
+        ax2.yaxis.set_label_position("right")  # Move y-axis label to the right
         ax2.set_xlim((lower_bound[0], upper_bound[0]))
         ax2.set_ylim((lower_bound[2], upper_bound[2]))
         ax2.set_title(f'{title}', color='white')
@@ -379,8 +382,8 @@ def get_main_category(frame):
                     pass
                 else:
                     pos[0] = pos[0] - LD / RGEO
-                ax2.scatter(pos[0], pos[2], color=color, label=point)
-                ax2.text(pos[0], pos[2], point, color=color)
+                ax2.scatter(pos[0], pos[2], color='w', label=point)
+                ax2.text(pos[0], pos[2], point, color='w')
 
         ax3.add_patch(plt.Circle((0, 0), stn[2], color='white', linestyle='dashed', fill=False))
         ax3.scatter(y, z, color=scatter_dot_colors, s=1)
@@ -399,9 +402,8 @@ def get_main_category(frame):
                     pass
                 else:
                     pos[0] = pos[0] - LD / RGEO
-                print(pos)
-                ax3.scatter(pos[1], pos[2], color=color, label=point)
-                ax3.text(pos[1], pos[2], point, color=color)
+                ax3.scatter(pos[1], pos[2], color='w', label=point)
+                ax3.text(pos[1], pos[2], point, color='w')
 
         ax4.scatter3D(x, y, z, color=scatter_dot_colors, s=1)
         ax4.scatter3D(0, 0, 0, color=stn[0], s=stn[1])
@@ -421,9 +423,9 @@ def get_main_category(frame):
                     pass
                 else:
                     pos[0] = pos[0] - LD / RGEO
-                ax4.scatter(pos[0], pos[1], pos[2], color=color, label=point)
-                ax4.text(pos[0], pos[1], pos[2], point, color=color)
-    
+                ax4.scatter(pos[0], pos[1], pos[2], color='w', label=point)
+                ax4.text(pos[0], pos[1], pos[2], point, color='w')
+
         return fig, ax1, ax2, ax3, ax4
     input_type = check_numpy_array(r)
 
@@ -432,13 +434,13 @@ def get_main_category(frame):
     ax2 = fig.add_subplot(2, 2, 2)
     ax3 = fig.add_subplot(2, 2, 3)
     ax4 = fig.add_subplot(2, 2, 4, projection='3d')
-    
+
     if input_type == "numpy array":
-        fig, ax1, ax2, ax3, ax4 = _make_scatter(fig, ax1, ax2, ax3, ax4, r=r, t=t, limits=limits, title=title, frame=frame)
+        fig, ax1, ax2, ax3, ax4 = _make_scatter(fig, ax1, ax2, ax3, ax4, r=r, t=t, title=title, frame=frame)
     if input_type == "list of numpy array":
         num_orbits = np.shape(r)[0]
         for i, row in enumerate(r):
-            fig, ax1, ax2, ax3, ax4 = _make_scatter(fig, ax1, ax2, ax3, ax4, r=row, t=t, limits=limits, title=title, orbit_index=i, num_orbits=num_orbits, frame=frame)
+            fig, ax1, ax2, ax3, ax4 = _make_scatter(fig, ax1, ax2, ax3, ax4, r=row, t=t, title=title, orbit_index=i, num_orbits=num_orbits, frame=frame)
 
     # Set axis color to white
     for i, ax in enumerate([ax1, ax2, ax3, ax4]):
@@ -458,8 +460,8 @@ def get_main_category(frame):
     for ax in [ax1, ax2, ax3, ax4]:
         for text in ax.get_xticklabels() + ax.get_yticklabels() + [ax.xaxis.label, ax.yaxis.label]:
             text.set_color('white')
-    
-    #Save the plot
+
+    # Save the plot
     fig.patch.set_facecolor('black')
     if show:
         plt.show()

From 0c1b51b9006f3835c1a5d966d99bf1e839f4f991 Mon Sep 17 00:00:00 2001
From: Travis Yeager <43413776+SuperdoerTrav@users.noreply.github.com>
Date: Mon, 6 Jan 2025 15:47:15 -0800
Subject: [PATCH 07/14] Added small fix to time initialization

In the case a time array was passed in, it was not ensured that the initial time was correctly assigned.
---
 ssapy/simple.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/ssapy/simple.py b/ssapy/simple.py
index 177d132..b3ffab5 100644
--- a/ssapy/simple.py
+++ b/ssapy/simple.py
@@ -221,6 +221,7 @@ def ssapy_orbit(orbit=None, a=None, e=0, i=0, pa=0, raan=0, ta=0, r=None, v=None
         time_is_None = True
         t = get_times(duration=duration, freq=freq, t0=t0)
     else:
+        t0 = t[0]
         time_is_None = False
 
     if orbit is not None:

From dcc594c48f0ef0aff3b732f5e75d6e703b2a686d Mon Sep 17 00:00:00 2001
From: Jason Bernstein <117031728+JasonBernstein1@users.noreply.github.com>
Date: Wed, 8 Jan 2025 19:45:13 -0700
Subject: [PATCH 08/14] Fix typos (#44)

Co-authored-by: Jason Bernstein <bernstein8@llnl.gov>
---
 ssapy/compute.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/ssapy/compute.py b/ssapy/compute.py
index f1f05cc..0cb6300 100644
--- a/ssapy/compute.py
+++ b/ssapy/compute.py
@@ -260,7 +260,7 @@ def __rv(orbit, time, propagator):
 
 def groundTrack(orbit, time, propagator=KeplerianPropagator(), format='geodetic'):
     """Calculate satellite ground track on the outer product of all supplied times and
-    statevectors or orbits.
+    state vectors or orbits.
 
     Parameters
     ----------
@@ -427,7 +427,7 @@ def dircos(
         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 accuration).
+        "linear" correction is almost always sufficiently accurate).
 
     Notes
     -----
@@ -503,7 +503,7 @@ def radec(
         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 accuration).
+        "linear" correction is almost always sufficiently accurate).
 
     Notes
     -----
@@ -590,7 +590,7 @@ def altaz(
         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 accuration).
+        "linear" correction is almost always sufficiently accurate).
 
     Notes
     -----
@@ -662,7 +662,7 @@ def quickAltAz(
         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 accuration).
+        "linear" correction is almost always sufficiently accurate).
 
     Notes
     -----
@@ -740,7 +740,7 @@ def radecRate(
         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 accuration).
+        "linear" correction is almost always sufficiently accurate).
 
     Notes
     -----
@@ -1436,7 +1436,7 @@ def earth_shine(r_sat, r_earth, r_sun, radius, albedo, albedo_earth, albedo_back
     """
     # https://amostech.com/TechnicalPapers/2013/POSTER/COGNION.pdf
     phase_angle = get_angle(r_sun, r_sat, r_earth)  # angle from Sun to object to Earth
-    earth_angle = np.pi - phase_angle  # Sun to Earth to oject.
+    earth_angle = np.pi - phase_angle  # Sun to Earth to object.
     r_earth_sat = np.linalg.norm(r_sat - r_earth, axis=-1)  # Earth is the observer.
     flux_earth_to_sat = 2 / 3 * albedo_earth * EARTH_RADIUS**2 / (np.pi * (r_earth_sat)**2) * (np.sin(earth_angle) + (np.pi - earth_angle) * np.cos(earth_angle))  # Fraction of sunlight reflected from the Earth to satellite
     # Fraction of light from back of solar panel

From dce5400b892b3c808f7767c403dcdc6bbf3704f2 Mon Sep 17 00:00:00 2001
From: Jason Bernstein <117031728+JasonBernstein1@users.noreply.github.com>
Date: Wed, 8 Jan 2025 19:47:04 -0700
Subject: [PATCH 09/14] Improve documentation in orbit_solver.py (#43)

* Add full citations for later references

* Add comma in comment

* Add docstring describing function and comment for context

* Clarify comment

* Replace comments with docstrings

* Replace comments with docstrings

* Replace comment with docstring and add details

* Replace comments with docstrings

* Replace comments with docstrings

* Tidy docstring

* Replace comment with docstring

* Replace comments with docstrings

---------

Co-authored-by: Jason Bernstein <bernstein8@llnl.gov>
---
 ssapy/orbit_solver.py | 71 ++++++++++++++++++++++++++++---------------
 1 file changed, 47 insertions(+), 24 deletions(-)

diff --git a/ssapy/orbit_solver.py b/ssapy/orbit_solver.py
index df50c14..92f76f0 100644
--- a/ssapy/orbit_solver.py
+++ b/ssapy/orbit_solver.py
@@ -12,6 +12,10 @@
 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.
 """
 
 
@@ -98,7 +102,7 @@ def __init__(self, r1, r2, t1, t2, mu=EARTH_MU,
         )
 
     def _finishOrbit(self, p):
-        # M&G (2.114)  good for both elliptical and hyperbolic orbits
+        # M&G (2.114), good for both elliptical and hyperbolic orbits
         ecosnu1 = p / self.r1mag - 1.0
         ecosnu2 = p / self.r2mag - 1.0  # M&G (2.114)
         # M&G (2.116)
@@ -129,6 +133,8 @@ def __init__(self, *args, **kwargs):
         TwoPosOrbitSolver.__init__(self, *args, **kwargs)
 
     def _getP(self):
+        """Compute p from Shefer (2)."""
+        # Solve for eta
         eta = 1
         d_eta = 1
         niter = 0
@@ -145,8 +151,7 @@ def _getP(self):
             d_eta = 1. + (self.ell + x) * X - eta
             eta += d_eta
             niter += 1
-
-        # Shefer (2)
+        # Plug eta into (2) to obtain p
         p = (0.5 * eta * self.kappa * self.sigma / self.tau)**2 / self.mu
         return p
 
@@ -157,11 +162,12 @@ def __init__(self, *args, **kwargs):
 
     @staticmethod
     def X(g):
-        return (2 * g - np.sin(2 * g)) / (np.sin(g)**3)  # Shefer (11)
+        """Compute X(g) from Shefer (11)."""
+        return (2 * g - np.sin(2 * g)) / (np.sin(g)**3)
 
     @staticmethod
     def dXdg(g):
-        # Shefer (12)
+        """Compute dX(g)/dg from Shefer (12)."""
         return (2 * (1 - np.cos(2 * g)) - 3 * (2 * g - np.sin(2 * g)) / np.tan(g)) / (np.sin(g)**3)
 
     def _getP(self):
@@ -213,18 +219,22 @@ def __init__(self, *args, **kwargs):
         self.nExam = kwargs.pop('nExam', 100)
         TwoPosOrbitSolver.__init__(self, *args, **kwargs)
 
-    def alpha(self, x):  # Shefer (18)
+    def alpha(self, x):
+        """Evaluate alpha(x) from Shefer (18)."""
         return (self.rbar + self.kappa * (x - 0.5)), self.kappa
 
-    def beta(self, xi):  # Shefer (A.4), (A.9)
+    def beta(self, xi):
+        """Evaluate beta(xi) and its derivative from Shefer (A.4) and (A.9)."""
         val = self.rbar - 0.5 * self.kappa + xi * (self.rbar + 0.5 * self.kappa)
         grad = self.rbar + 0.5 * self.kappa
         return val, grad
 
-    def Y(self, x):  # Shefer (17)
+    def Y(self, x):
+        """Evaluate Y(x) and dY(x)/dx from Shefer (17)."""
         XVal, XGrad = self.X(x)
         alphaVal, alphaGrad = self.alpha(x)
         val = self.kappa + alphaVal * XVal
+        # Evaluate dY(x)/dx using the chain rule
         grad = alphaVal * XGrad + alphaGrad * XVal
         return val, grad
 
@@ -233,7 +243,11 @@ def Yxi(self, xi):
         return self.kappa + betaVal * self.Z(xi)[0]
 
     @staticmethod
-    def X(x):  # Shefer (19), (20), (23)
+    def 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).
+        """
         import astropy.units as u
         if isinstance(x, u.Quantity):
             x = x.value
@@ -261,21 +275,26 @@ def X(x):  # Shefer (19), (20), (23)
             return np.inf  # X(x) actually asymptotes to Infinity as x -> 1
 
     @staticmethod
-    def Z(xi):  # Shefer (A.5)
+    def Z(xi):
+        """Compute Z(xi) function from Shefer (A.5)."""
         num = (1 + xi)**2 * np.arctanh(np.sqrt(-xi)) - (1 - xi) * np.sqrt(-xi)
         denom = 2 * xi * np.sqrt(-xi)
         val = num / denom
         grad = (4 - (3 - xi) * val) / (2 * xi * (1 + xi))
         return val, grad
 
-    def D(self, x):  # Shefer (43) and derivative
+    def D(self, x):
+        """Compute D(x) and its derivative, dD(x)/dx, from Shefer (43)."""
         aval, agrad = self.semiMajorAxis(x)
         val = self.tau * np.sqrt(self.mu) - 2 * np.pi * self.lam * aval**1.5
         grad = -3 * np.pi * self.lam * np.sqrt(aval) * agrad
         # print("----- D: ", x, aval, agrad, self.tau, self.mu, self.lam)
         return val, grad
 
-    def semiMajorAxis(self, x):  # Shefer (42) and derivative
+    def semiMajorAxis(self, x):
+        """Compute semi-major axis a(x) and its derivative, da(x)/dx, from
+        Shefer (42).
+        """
         assert x >= 0. and x <= 1., \
             "ERROR: invalid x in Shefer semi_major_axis: {}".format(x)
         alphaVal, alphaGrad = self.alpha(x)
@@ -293,7 +312,8 @@ def Flam0(self, x):  # Shefer (21), (22)
         grad = self.kappa * Ysq + 2 * alphaVal * Y * (self.kappa * XVal + alphaVal * XGrad)
         return val, grad
 
-    def F(self, x):  # Shefer (40), (41)
+    def F(self, x):
+        """Compute F(x) and dF(x)/dx from Shefer (40) and (41)."""
         if x >= 1.:
             # Should not get here
             val = 1.e16
@@ -307,7 +327,8 @@ def F(self, x):  # Shefer (40), (41)
             grad = alphaGrad * YVal**2 + alphaVal * 2 * YVal * YGrad - 2 * DVal * DGrad
         return val, grad
 
-    def G(self, xi):  # Shefer (A.7), (A.8)
+    def G(self, xi):
+        """Compute G(xi) and its derivative from Shefer (A.7) and (A.8)."""
         betaVal, betaGrad = self.beta(xi)
         Y = self.Yxi(xi)
         Ysq = Y**2
@@ -380,11 +401,7 @@ def _getInitialXiGuess(self):
         return xi
 
     def _getP(self):
-        """
-        Get the semi-latus rectum
-
-        This is the final result of the Shefer algorithm
-        """
+        """Get the semi-latus rectum. This is the final result of the Shefer algorithm."""
         x = self._getInitialXGuess()
         if x < -1. or x > 1.:
             xi = self._getInitialXiGuess()
@@ -413,7 +430,8 @@ def _getAllP(self):
         xs = find_all_zeros(lambda x: self.F(x)[0], xmin, xmax, n=self.nExam)
         return [self.sigma**2 / (4 * self.alpha(x)[0]) for x in xs]
 
-    def _getEta(self, p):  # Shefer (2)
+    def _getEta(self, p):
+        """Compute auxiliary value defined in Shefer (2)."""
         return 2 * np.sqrt(p * self.mu) * self.tau / (self.kappa * self.sigma)
 
     def solve(self):
@@ -571,25 +589,30 @@ def t3231(self):
     def t2131(self):
         return self.t21 / self.t31
 
-    def _getRho(self, n1, n3):  # M&G (2.129)
+    def _getRho(self, n1, n3):
+        """Compute three unknown station-satellite distances from Montenbruck and Gill (2.129).
+        """
         rho1 = -1 / (n1 * self.D) * (n1 * self.D11 - self.D12 + n3 * self.D13)
         rho2 = 1 / self.D * (n1 * self.D21 - self.D22 + n3 * self.D23)
         rho3 = -1 / (n3 * self.D) * (n1 * self.D31 - self.D32 + n3 * self.D33)
         return rho1, rho2, rho3
 
-    def _getR(self, rho1, rho2, rho3):  # M&G (2.122)
+    def _getR(self, rho1, rho2, rho3):
+        """Compute three Earth-satellite position vectors from Montenbruck and Gill (2.122).
+        """
         r1 = self.R1 + rho1 * self.e1
         r2 = self.R2 + rho2 * self.e2
         r3 = self.R3 + rho3 * self.e3
         return r1, r2, r3
 
-    def _getN(self, eta21, eta23):  # M&G (2.132)
+    def _getN(self, eta21, eta23):
+        """Compute n1 and n3 terms from Montenbruck and Gill (2.132)."""
         n1 = eta21 * self.t3231
         n3 = eta23 * self.t2131
         return n1, n3
 
-    # Use Shefer algorithm to improve eta estimates
     def _getEta(self, r1, r2, t1, t2):
+        """Use Shefer algorithm to improve eta estimates."""
         solver = SheferTwoPosOrbitSolver(r1, r2, t1, t2)
         p = solver._getP()
         eta = solver._getEta(p)

From 8cacc88a9138d4ab75f09ba7ca5c08c3b0bf5e22 Mon Sep 17 00:00:00 2001
From: Jason Bernstein <117031728+JasonBernstein1@users.noreply.github.com>
Date: Wed, 8 Jan 2025 19:51:53 -0700
Subject: [PATCH 10/14] Fix typos and uppercase acronym (#41)

Co-authored-by: Jason Bernstein <bernstein8@llnl.gov>
---
 tests/test_plots.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/test_plots.py b/tests/test_plots.py
index 8b21aa4..698933d 100644
--- a/tests/test_plots.py
+++ b/tests/test_plots.py
@@ -69,7 +69,7 @@ def initialize_DRO(t, delta_r=7.52064e7, delta_v=344):
 r, v, t = ssapy.simple.ssapy_orbit(orbit=dro_orbit, t=times)
 ssapy.plotUtils.orbit_plot(r=r, t=times, save_path=f"{save_folder}/DRO_orbit", frame='Lunar', show=False)
 r_lunar, v_lunar = ssapy.utils.gcrf_to_lunar_fixed(r, t=times, v=True)
-print("Succesfully converted gcrf to lunar frame.")
+print("Successfully converted GCRF to lunar frame.")
 ssapy.plotUtils.koe_plot(r, v, t=times, body='Earth', save_path=f"{save_folder}Keplerian_orbital_elements.png")
 
 ssapy.plotUtils.orbit_plot(r=r, t=times, save_path=f"{save_folder}/gcrf_plot.png", frame='gcrf', show=True)
@@ -165,7 +165,7 @@ def initialize_DRO(t, delta_r=7.52064e7, delta_v=344):
 
 print(f"Lagrange points were calculated correctly.")
 print(f"Rotate vector plot successfully created.")
-print(f"save_plot() executed succesfully.")
-print(f"save_animated_gif() executed succesfully.")
+print(f"save_plot() executed successfully.")
+print(f"save_animated_gif() executed successfully.")
 
 print(f"\nFinished plot testing!\n")

From 3a851b44709ad852995ade1f8182ab9192827f12 Mon Sep 17 00:00:00 2001
From: Jason Bernstein <117031728+JasonBernstein1@users.noreply.github.com>
Date: Wed, 8 Jan 2025 19:52:27 -0700
Subject: [PATCH 11/14] Add docstring explaining origin and purpose of test
 (#39)

Co-authored-by: Jason Bernstein <bernstein8@llnl.gov>
---
 tests/test_accel.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/tests/test_accel.py b/tests/test_accel.py
index 8dc1b49..ef372a8 100644
--- a/tests/test_accel.py
+++ b/tests/test_accel.py
@@ -564,6 +564,9 @@ def test_MG_3_4_scipy():
 
 @timer
 def test_MG_3_4_rk78():
+    """Exercise 3.4 from Montenbruck and Gill.  Tests SSAPy orbit propagation
+    with perturbations.
+    """
     t = Time("1999-03-01", scale='utc')
     # LEO sat
     kElements = [7178e3, 0.001, np.deg2rad(98.57), 0.0, 0.0, 0.0]

From 8c860e3bc12671ad7e650959e5422c6c781513cc Mon Sep 17 00:00:00 2001
From: Jason Bernstein <117031728+JasonBernstein1@users.noreply.github.com>
Date: Wed, 8 Jan 2025 19:53:17 -0700
Subject: [PATCH 12/14] Fix test of sun's gravitational parameter (#38)

Co-authored-by: Jason Bernstein <bernstein8@llnl.gov>
---
 tests/test_orbit.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_orbit.py b/tests/test_orbit.py
index 944a666..6a2644a 100755
--- a/tests/test_orbit.py
+++ b/tests/test_orbit.py
@@ -1886,7 +1886,7 @@ def test_musun():
     # failing example.
 
     a = u.AU.to(u.m)
-    mu = ssapy.constants.GM_SUN
+    mu = ssapy.constants.SUN_MU
     e = 0.001
     i = 0.001
     pa = 0.001

From 4713e356d8ca239be81d71137988f3f4e0d57b78 Mon Sep 17 00:00:00 2001
From: Travis Yeager <43413776+SuperdoerTrav@users.noreply.github.com>
Date: Fri, 10 Jan 2025 13:07:24 -0800
Subject: [PATCH 13/14] orbit_plot bounds updated

Bounds were not passed for the case of a single orbit.
Multiplication scaling on the bounds caused unequal axis scaling. Implemented a linear fix to the bounding boxes.
---
 ssapy/plotUtils.py | 172 +++++++++++++++++++++++++++++++--------------
 1 file changed, 118 insertions(+), 54 deletions(-)

diff --git a/ssapy/plotUtils.py b/ssapy/plotUtils.py
index 85dd32f..71e5343 100644
--- a/ssapy/plotUtils.py
+++ b/ssapy/plotUtils.py
@@ -260,7 +260,7 @@ def orbit_plot(r, t=[], title='', figsize=(7, 7), save_path=False, frame="gcrf",
     fig, axes = orbit_plot(r, t, title='Orbit Plot', frame='gcrf', show=True)
     """
 
-    def _make_scatter(fig, ax1, ax2, ax3, ax4, r, t, title='', orbit_index='', num_orbits=1, frame=False):
+    def _make_scatter(fig, ax1, ax2, ax3, ax4, r, t, bounds, title='', orbit_index='', num_orbits=1, frame=False):
         if np.size(t) < 1:
             if frame in ["itrf", "lunar", "lunar_fixed"]:
                 raise("Need to provide t for itrf, lunar or lunar fixed frames")
@@ -302,33 +302,74 @@ def get_main_category(frame):
             if transform_func:
                 r = transform_func(r, t)
                 r_moon = transform_func(r_moon, t)
+                r_earth = transform_func(np.zeros(np.shape(r_moon)), t)
         else:
             raise ValueError("Unknown plot type provided. Accepted: gcrf, itrf, lunar, lunar fixed")
 
-        x = r[:, 0] / RGEO
-        y = r[:, 1] / RGEO
-        z = r[:, 2] / RGEO
-        xm = r_moon[:, 0] / RGEO
-        ym = r_moon[:, 1] / RGEO
-        zm = r_moon[:, 2] / RGEO
-
-        if orbit_index == '' or orbit_index == 0:
-            lower_bound = np.array([np.inf, np.inf, np.inf])
-            upper_bound = np.array([-np.inf, -np.inf, -np.inf])
-        lower_bound_temp, upper_bound_temp = find_smallest_bounding_cube(r / RGEO * 1.2)
-        lower_bound = np.minimum(lower_bound, lower_bound_temp)
-        upper_bound = np.maximum(upper_bound, upper_bound_temp)
+        r = r / RGEO
+        r_moon = r_moon / RGEO
+        r_earth = r_earth / RGEO
+        x = r[:, 0]
+        y = r[:, 1]
+        z = r[:, 2]
+        xm = r_moon[:, 0]
+        ym = r_moon[:, 1]
+        zm = r_moon[:, 2]
+        xe = r_earth[:, 0]
+        ye = r_earth[:, 1]
+        ze = r_earth[:, 2]
+
+        lower_bound_temp, upper_bound_temp = find_smallest_bounding_cube(r)
+        bounds["lower"] = np.minimum(bounds["lower"], lower_bound_temp)
+        bounds["upper"] = np.maximum(bounds["upper"], upper_bound_temp)
+        print("Bounds: ", bounds)
+        print(f"Transformed positions (r): min={(r).min(axis=0)}, max={(r).max(axis=0)}")
+        print(f"Transformed positions (r_moon): min={(r_moon).min(axis=0)}, max={(r_moon).max(axis=0)}")
+        print(f"Transformed positions (r_earth): min={(r_earth).min(axis=0)}, max={(r_earth).max(axis=0)}")
+
         if np.size(xm) > 1:
-            gradient_colors = cm.Greys(np.linspace(0, .8, len(xm)))[::-1]
+            grey_colors = cm.Greys(np.linspace(0, .8, len(xm)))[::-1]
             blues = cm.Blues(np.linspace(.4, .9, len(xm)))[::-1]
         else:
-            gradient_colors = "grey"
+            grey_colors = "grey"
             blues = 'Blue'
         plot_settings = {
-            "gcrf": ("blue", 50, 1, xm, ym, zm, gradient_colors),
-            "itrf": ("blue", 50, 1, xm, ym, zm, gradient_colors),
-            "lunar": ("grey", 25, 1.3, xm, ym, zm, blues),
-            "lunar axis": ("blue", 50, 1, -xm, -ym, -zm, gradient_colors)
+            "gcrf": {
+                "primary_color": "blue",
+                "primary_size": 50,
+                "secondary_x": xm,
+                "secondary_y": ym,
+                "secondary_z": zm,
+                "secondary_color": grey_colors,
+                "secondary_size": 25
+            },
+            "itrf": {
+                "primary_color": "blue",
+                "primary_size": 50,
+                "secondary_x": xm,
+                "secondary_y": ym,
+                "secondary_z": zm,
+                "secondary_color": grey_colors,
+                "secondary_size": 25
+            },
+            "lunar": {
+                "primary_color": "grey",
+                "primary_size": 25,
+                "secondary_x": xe,
+                "secondary_y": ye,
+                "secondary_z": ze,
+                "secondary_color": blues,
+                "secondary_size": 50
+            },
+            "lunar axis": {
+                "primary_color": "blue",
+                "primary_size": 50,
+                "secondary_x": xm,
+                "secondary_y": ym,
+                "secondary_z": zm,
+                "secondary_color": grey_colors,
+                "secondary_size": 25
+            }
         }
         try:
             stn = plot_settings[frame]
@@ -341,16 +382,14 @@ def get_main_category(frame):
         else:
             angle = orbit_index * 10
             scatter_dot_colors = cm.rainbow(np.linspace(0, 1, num_orbits))[orbit_index]
-        ax1.add_patch(plt.Circle((0, 0), stn[2], color='white', linestyle='dashed', fill=False))
+        ax1.add_patch(plt.Circle((0, 0), 1, color='white', linestyle='dashed', fill=False))
         ax1.scatter(x, y, color=scatter_dot_colors, s=1)
-        ax1.scatter(0, 0, color=stn[0], s=stn[1])
+        ax1.scatter(0, 0, color=stn['primary_color'], s=stn['primary_size'])
         if xm is not False:
-            ax1.scatter(stn[3], stn[4], color=stn[6], s=5)
+            ax1.scatter(stn['secondary_x'], stn['secondary_y'], color=stn['secondary_color'], s=stn['secondary_size'])
         ax1.set_aspect('equal')
         ax1.set_xlabel('x [GEO]')
         ax1.set_ylabel('y [GEO]')
-        ax1.set_xlim((lower_bound[0], upper_bound[0]))
-        ax1.set_ylim((lower_bound[1], upper_bound[1]))
         ax1.set_title(f'Frame: {title2}', color='white')
         if 'lunar' in frame:
             colors = ['red', 'green', 'purple', 'orange', 'cyan']
@@ -359,21 +398,20 @@ def get_main_category(frame):
                     pass
                 else:
                     pos[0] = pos[0] - LD / RGEO
-                ax1.scatter(pos[0], pos[1], color='w', label=point)
-                ax1.text(pos[0], pos[1], point, color='w')
+                if bounds["lower"][0] <= pos[0] <= bounds["upper"][0] and bounds["lower"][1] <= pos[1] <= bounds["upper"][1]:
+                    ax1.scatter(pos[0], pos[1], color='w', label=point)
+                    ax1.text(pos[0], pos[1], point, color='w')
 
-        ax2.add_patch(plt.Circle((0, 0), stn[2], color='white', linestyle='dashed', fill=False))
+        ax2.add_patch(plt.Circle((0, 0), 1, color='white', linestyle='dashed', fill=False))
         ax2.scatter(x, z, color=scatter_dot_colors, s=1)
-        ax2.scatter(0, 0, color=stn[0], s=stn[1])
+        ax2.scatter(0, 0, color=stn['primary_color'], s=stn['primary_size'])
         if xm is not False:
-            ax2.scatter(stn[3], stn[5], color=stn[6], s=5)
+            ax2.scatter(stn['secondary_x'], stn['secondary_z'], color=stn['secondary_color'], s=stn['secondary_size'])
         ax2.set_aspect('equal')
         ax2.set_xlabel('x [GEO]')
         ax2.set_ylabel('z [GEO]')
         ax2.yaxis.tick_right()  # Move y-axis ticks to the right
         ax2.yaxis.set_label_position("right")  # Move y-axis label to the right
-        ax2.set_xlim((lower_bound[0], upper_bound[0]))
-        ax2.set_ylim((lower_bound[2], upper_bound[2]))
         ax2.set_title(f'{title}', color='white')
         if 'lunar' in frame:
             colors = ['red', 'green', 'purple', 'orange', 'cyan']
@@ -382,19 +420,18 @@ def get_main_category(frame):
                     pass
                 else:
                     pos[0] = pos[0] - LD / RGEO
-                ax2.scatter(pos[0], pos[2], color='w', label=point)
-                ax2.text(pos[0], pos[2], point, color='w')
+                if bounds["lower"][0] <= pos[0] <= bounds["upper"][0] and bounds["lower"][2] <= pos[2] <= bounds["upper"][2]:
+                    ax2.scatter(pos[0], pos[2], color='w', label=point)
+                    ax2.text(pos[0], pos[2], point, color='w')
 
-        ax3.add_patch(plt.Circle((0, 0), stn[2], color='white', linestyle='dashed', fill=False))
+        ax3.add_patch(plt.Circle((0, 0), 1, color='white', linestyle='dashed', fill=False))
         ax3.scatter(y, z, color=scatter_dot_colors, s=1)
-        ax3.scatter(0, 0, color=stn[0], s=stn[1])
+        ax3.scatter(0, 0, color=stn['primary_color'], s=stn['primary_size'])
         if xm is not False:
-            ax3.scatter(stn[4], stn[5], color=stn[6], s=5)
+            ax3.scatter(stn['secondary_y'], stn['secondary_z'], color=stn['secondary_color'], s=stn['secondary_size'])
         ax3.set_aspect('equal')
         ax3.set_xlabel('y [GEO]')
         ax3.set_ylabel('z [GEO]')
-        ax3.set_xlim((lower_bound[1], upper_bound[1]))
-        ax3.set_ylim((lower_bound[2], upper_bound[2]))
         if 'lunar' in frame:
             colors = ['red', 'green', 'purple', 'orange', 'cyan']
             for (point, pos), color in zip(lagrange_points_lunar_frame().items(), colors):
@@ -402,17 +439,14 @@ def get_main_category(frame):
                     pass
                 else:
                     pos[0] = pos[0] - LD / RGEO
-                ax3.scatter(pos[1], pos[2], color='w', label=point)
-                ax3.text(pos[1], pos[2], point, color='w')
+                if bounds["lower"][1] <= pos[1] <= bounds["upper"][1] and bounds["lower"][2] <= pos[2] <= bounds["upper"][2]:
+                    ax3.scatter(pos[1], pos[2], color='w', label=point)
+                    ax3.text(pos[1], pos[2], point, color='w')
 
         ax4.scatter3D(x, y, z, color=scatter_dot_colors, s=1)
-        ax4.scatter3D(0, 0, 0, color=stn[0], s=stn[1])
+        ax4.scatter3D(0, 0, 0, color=stn['primary_color'], s=stn['primary_size'])
         if xm is not False:
-            ax4.scatter3D(stn[3], stn[4], stn[5], color=stn[6], s=5)
-        ax4.set_xlim([lower_bound[0], upper_bound[0]])
-        ax4.set_ylim([lower_bound[1], upper_bound[1]])
-        ax4.set_zlim([lower_bound[2], upper_bound[2]])
-        ax4.set_aspect('equal')  # aspect ratio is 1:1:1 in data space
+            ax4.scatter3D(stn['secondary_x'], stn['secondary_y'], stn['secondary_z'], color=stn['secondary_color'], s=stn['secondary_size'])
         ax4.set_xlabel('x [GEO]')
         ax4.set_ylabel('y [GEO]')
         ax4.set_zlabel('z [GEO]')
@@ -423,10 +457,38 @@ def get_main_category(frame):
                     pass
                 else:
                     pos[0] = pos[0] - LD / RGEO
-                ax4.scatter(pos[0], pos[1], pos[2], color='w', label=point)
-                ax4.text(pos[0], pos[1], pos[2], point, color='w')
-
-        return fig, ax1, ax2, ax3, ax4
+                if bounds["lower"][0] <= pos[0] <= bounds["upper"][0] and bounds["lower"][1] <= pos[1] <= bounds["upper"][1] and bounds["lower"][2] <= pos[2] <= bounds["upper"][2]:
+                    ax4.scatter(pos[0], pos[1], pos[2], color='w', label=point)
+                    ax4.text(pos[0], pos[1], pos[2], point, color='w')
+
+        modifier = 1
+        # Update for ax1
+        ax1.set_xlim(bounds["lower"][0] + modifier * np.sign(bounds["lower"][0]),
+                    bounds["upper"][0] - modifier * np.sign(bounds["upper"][0]))
+        ax1.set_ylim(bounds["lower"][1] + modifier * np.sign(bounds["lower"][1]),
+                    bounds["upper"][1] - modifier * np.sign(bounds["upper"][1]))
+
+        # Update for ax2
+        ax2.set_xlim(bounds["lower"][0] + modifier * np.sign(bounds["lower"][0]),
+                    bounds["upper"][0] - modifier * np.sign(bounds["upper"][0]))
+        ax2.set_ylim(bounds["lower"][2] + modifier * np.sign(bounds["lower"][2]),
+                    bounds["upper"][2] - modifier * np.sign(bounds["upper"][2]))
+
+        # Update for ax3
+        ax3.set_xlim(bounds["lower"][1] + modifier * np.sign(bounds["lower"][1]),
+                    bounds["upper"][1] - modifier * np.sign(bounds["upper"][1]))
+        ax3.set_ylim(bounds["lower"][2] + modifier * np.sign(bounds["lower"][2]),
+                    bounds["upper"][2] - modifier * np.sign(bounds["upper"][2]))
+
+        # Update for ax4
+        ax4.set_xlim(bounds["lower"][0] + modifier * np.sign(bounds["lower"][0]),
+                    bounds["upper"][0] - modifier * np.sign(bounds["upper"][0]))
+        ax4.set_ylim(bounds["lower"][1] + modifier * np.sign(bounds["lower"][1]),
+                    bounds["upper"][1] - modifier * np.sign(bounds["upper"][1]))
+        ax4.set_zlim(bounds["lower"][2] + modifier * np.sign(bounds["lower"][2]),
+                    bounds["upper"][2] - modifier * np.sign(bounds["upper"][2]))
+
+        return fig, ax1, ax2, ax3, ax4, bounds
     input_type = check_numpy_array(r)
 
     fig = plt.figure(dpi=100, figsize=figsize, facecolor='black')
@@ -434,13 +496,15 @@ def get_main_category(frame):
     ax2 = fig.add_subplot(2, 2, 2)
     ax3 = fig.add_subplot(2, 2, 3)
     ax4 = fig.add_subplot(2, 2, 4, projection='3d')
-
+    
+    bounds = {"lower": np.array([np.inf, np.inf, np.inf]), "upper": np.array([-np.inf, -np.inf, -np.inf])}
+    
     if input_type == "numpy array":
-        fig, ax1, ax2, ax3, ax4 = _make_scatter(fig, ax1, ax2, ax3, ax4, r=r, t=t, title=title, frame=frame)
+        fig, ax1, ax2, ax3, ax4, bounds = _make_scatter(fig, ax1, ax2, ax3, ax4, r=r, t=t, title=title, frame=frame, bounds=bounds)
     if input_type == "list of numpy array":
         num_orbits = np.shape(r)[0]
         for i, row in enumerate(r):
-            fig, ax1, ax2, ax3, ax4 = _make_scatter(fig, ax1, ax2, ax3, ax4, r=row, t=t, title=title, orbit_index=i, num_orbits=num_orbits, frame=frame)
+            fig, ax1, ax2, ax3, ax4, bounds = _make_scatter(fig, ax1, ax2, ax3, ax4, r=row, t=t, title=title, orbit_index=i, num_orbits=num_orbits, frame=frame, bounds=bounds)
 
     # Set axis color to white
     for i, ax in enumerate([ax1, ax2, ax3, ax4]):

From 8b05b8985efba36aea0efcd99692f910bc963255 Mon Sep 17 00:00:00 2001
From: Peter McGill <18149966+astrophpeter@users.noreply.github.com>
Date: Mon, 13 Jan 2025 20:18:31 -0800
Subject: [PATCH 14/14] JOSS paper submission (#47)

* Add files needed for JOSS paper with the associated LLNL release number.

* Update the date on the paper.

---------

Co-authored-by: Alexx Perloff <perloff1@llnl.gov>
---
 JOSS/ground_track.png |   3 +
 JOSS/orbit_plot.png   |   3 +
 JOSS/paper.bib        | 404 ++++++++++++++++++++++++++++++++++++++++++
 JOSS/paper.md         | 132 ++++++++++++++
 4 files changed, 542 insertions(+)
 create mode 100644 JOSS/ground_track.png
 create mode 100644 JOSS/orbit_plot.png
 create mode 100644 JOSS/paper.bib
 create mode 100644 JOSS/paper.md

diff --git a/JOSS/ground_track.png b/JOSS/ground_track.png
new file mode 100644
index 0000000..ef09a69
--- /dev/null
+++ b/JOSS/ground_track.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7af1ebbca5db9702bb088200ef80afa4ada40f188cea7574fbfa37f897472a85
+size 1573431
diff --git a/JOSS/orbit_plot.png b/JOSS/orbit_plot.png
new file mode 100644
index 0000000..07981d9
--- /dev/null
+++ b/JOSS/orbit_plot.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:60b56a82fff649d7f9b705909219b7a1ba9176a1ae9ac968f53a1b0b1184b190
+size 126636
diff --git a/JOSS/paper.bib b/JOSS/paper.bib
new file mode 100644
index 0000000..64a1678
--- /dev/null
+++ b/JOSS/paper.bib
@@ -0,0 +1,404 @@
+@techreport{Bernstein2021,
+  author       = {Bernstein, Jason and Filippov, Andrey and Schneider, Michael and Miller, Caleb},
+  title        = {Quantifying Uncertainty in All-to-All Estimates of Space Object Conjunction Probabilities using U-Statistics},
+  institution  = {Lawrence Livermore National Lab. (LLNL), Livermore, CA (United States)},
+  annote       = {Predicting space object conjunctions is inherently probabilistic due to initial state and orbit model uncertainty. A commonly considered Monte Carlo estimator of the conjunction probability is the ’all-to-all’ estimator. Given independent random samples of the trajectories of both objects, the estimator is the percentage of all pairs of trajectories that result in a conjunction. Intuitively, the all-to-all estimator is the best possible estimator of the conjunction probability since it considers all pairs of Monte Carlo samples. However, its distribution is not available in closed-form, which limits its use in practice and makes this intuition difficult to make rigorous. In this paper, the all-to-all estimator is identified as a U-statistic, which implies that it has several favorable properties. Specifically, the estimator is the minimum variance unbiased estimator of the conjunction probability and is asymptotically Gaussian distributed. An approximate confidence interval for the conjunction probability is obtained from an estimate of the asymptotic Gaussian distribution. We show how to efficiently compute the confidence interval and demonstrate that the interval has the nominal coverage level. The confidence intervals are also seen to be narrower than those based on the commonly-used each-to-each estimator. Furthermore, the all-to-all estimator is shown to allow different Monte Carlo sample sizes, whereas the each-to-each estimator requires equal sample sizes.},
+  doi          = {10.2172/1825370},
+  url          = {https://www.osti.gov/biblio/1825370},
+  place        = {United States},
+  year         = {2021},
+  month        = {10}}
+
+@misc{MOU,
+  author       = "",
+  title        = "Memorandum of Understanding between the National Aeronautics and Space
+Administration and the United States Space Force",
+  howpublished = "",
+  month        = "",
+  year         = "2020",
+  note         = "",
+  annote       = ""
+}
+
+@INPROCEEDINGS{Duggan2019,
+
+  author={Duggan, Matthew and Simon, Xavier and Moseman, Travis},
+
+  booktitle={2019 IEEE Aerospace Conference}, 
+
+  title={Lander and Cislunar Gateway Architecture Concepts for Lunar Exploration}, 
+
+  year={2019},
+
+  volume={},
+
+  number={},
+
+  pages={1-9},
+
+  keywords={Logic gates;Moon;Space vehicles;Deep-space communications;Orbits;Propulsion;Robots},
+
+  doi={10.1109/AERO.2019.8741766}}
+
+}
+
+@inproceedings{Vallado2006,
+  title={Revisiting spacetrack report\# 3},
+  author={Vallado, David and Crawford, Paul and Hujsak, Ricahrd and Kelso, TS},
+  booktitle={AIAA/AAS Astrodynamics Specialist Conference and Exhibit},
+  pages={6753},
+  year={2006}
+}
+
+
+@misc{Rhodes2023,
+  author = {Brandon Rhodes},
+  title = {python-sgp4},
+  year = {2023},
+  publisher = {GitHub},
+  journal = {GitHub repository},
+  howpublished = {\url{https://github.com/brandon-rhodes/python-sgp4}},
+}
+
+@software{newville2024,
+  author       = {Matt Newville and
+                  Renee Otten and
+                  Andrew Nelson and
+                  Till Stensitzki and
+                  Antonino Ingargiola and
+                  Dan Allan and
+                  Austin Fox and
+                  Faustin Carter and
+                  Michał and
+                  Ray Osborn and
+                  Dima Pustakhod and
+                  Sebastian Weigand and
+                  lneuhaus and
+                  Andrey Aristov and
+                  Glenn and
+                  Mark and
+                  mgunyho and
+                  Christoph Deil and
+                  Allan L. R. Hansen and
+                  Gustavo Pasquevich and
+                  Leon Foks and
+                  Nicholas Zobrist and
+                  Oliver Frost and
+                  Stuermer and
+                  Jean-Christophe Jaskula and
+                  Shane Caldwell and
+                  Pieter Eendebak and
+                  Matteo Pompili and
+                  Jens Hedegaard Nielsen and
+                  Arun Persaud},
+  title        = {lmfit/lmfit-py: 1.3.2},
+  month        = jul,
+  year         = 2024,
+  publisher    = {Zenodo},
+  version      = {1.3.2},
+  doi          = {10.5281/zenodo.12785036},
+  url          = {https://doi.org/10.5281/zenodo.12785036}
+}
+
+@software{Kerkwijk2023,
+  author       = {Marten van Kerkwijk and
+                  Erik Tollerud and
+                  Antonio Valentino and
+                  Thomas Robitaille and
+                  Julien Woillez and
+                  E. M. Bray and
+                  Brigitta Sipőcz and
+                  Michael Droettboom and
+                  P. L. Lim and
+                  Christoph Deil and
+                  Michael Seifert and
+                  Simon Conseil and
+                  Tom Aldcroft and
+                  Adrian Price-Whelan and
+                  Hood Chatham and
+                  StuartLittlefair and
+                  Chris Beaumont and
+                  Chris Lamb and
+                  Daria Cara and
+                  Devin Crichton and
+                  Matt Davis and
+                  Miguel de Val-Borro and
+                  Sergio Pascual and
+                  Stefan Heimersheim and
+                  Stuart Mumford and
+                  Tomas Babej and
+                  Vatsala Swaroop and
+                  Vishnunarayan K I and
+                  Jani Šumak},
+  title        = {liberfa/pyerfa: v2.0.1.1},
+  month        = oct,
+  year         = 2023,
+  publisher    = {Zenodo},
+  version      = {v2.0.1.1},
+  doi          = {10.5281/zenodo.10023045},
+  url          = {https://doi.org/10.5281/zenodo.10023045}
+}
+
+@ARTICLE{Miller2022,
+
+  author={Miller, Caleb and Corcoran, Jem N. and Schneider, Michael D.},
+
+  journal={IEEE Signal Processing Letters}, 
+
+  title={Rare Events via Cross-Entropy Population Monte Carlo}, 
+
+  year={2022},
+
+  volume={29},
+
+  number={},
+
+  pages={439-443},
+
+  keywords={Proposals;Monte Carlo methods;Statistics;Sociology;Signal processing algorithms;Artificial intelligence;Optimization;Adaptive importance sampling;cross-entropy;rare events},
+
+  doi={10.1109/LSP.2021.3139572}}
+
+
+@ARTICLE{Rein2012,
+       author = {{Rein}, H. and {Liu}, S. -F.},
+        title = "{REBOUND: an open-source multi-purpose N-body code for collisional dynamics}",
+      journal = {\aap},
+     keywords = {methods: numerical, planets and satellites: rings, protoplanetary disks, Astrophysics - Earth and Planetary Astrophysics, Astrophysics - Instrumentation and Methods for Astrophysics, Mathematics - Dynamical Systems, Physics - Computational Physics},
+         year = 2012,
+        month = jan,
+       volume = {537},
+          eid = {A128},
+        pages = {A128},
+          doi = {10.1051/0004-6361/201118085},
+archivePrefix = {arXiv},
+       eprint = {1110.4876},
+ primaryClass = {astro-ph.EP},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2012A&A...537A.128R},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@inproceedings{Hughes2014,
+  title={Verification and validation of the general mission analysis tool (GMAT)},
+  author={Hughes, Steven P and Qureshi, Rizwan H and Cooley, Steven D and Parker, Joel J},
+  booktitle={AIAA/AAS astrodynamics specialist conference},
+  pages={4151},
+  year={2014}
+}
+
+
+@inproceedings{Yeager2023,
+  title = {{Long-term N-body Stability in Cislunar Space}},
+  author = {Travis Yeager and Kerianne Pruett and Michael Schneider},
+  booktitle = {Proceedings of the Advanced Maui Optical and Space Surveillance (AMOS) Technologies Conference},
+  year = {2023},
+  organization = {Lawrence Livermore National Laboratory},
+  note = {Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS)},
+  url = "https://amostech.com/TechnicalPapers/2023/Poster/Yeager.pdf"
+}
+
+@inproceedings{Higgins2024,
+  title = {{SOM-erizing Cislunar Orbits: Classification of Cislunar Orbits Using Self-Organizing Maps (SOMs)}},
+  author = {Denvir Higgins and Kerianne Pruett and Travis Yeager and Michael Schneider},
+  booktitle = {Proceedings of the Advanced Maui Optical and Space Surveillance (AMOS) Technologies Conference},
+  year = {2024},
+  organization = {Lawrence Livermore National Laboratory},
+  note = {Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS)},
+  url = "https://amostech.com/TechnicalPapers/2024/Poster/Higgins.pdf"
+}
+
+@inproceedings{Pruett2024,
+  title = {{Closely-Spaced Object Classification Using MuyGPyS}},
+  author = {Pruett, Kerianne and McNaughton, Nathan and Schneider, Michael},
+  booktitle = {Proceedings of the Advanced Maui Optical and Space Surveillance (AMOS) Technologies Conference},
+  year = {2024},
+  organization = {Lawrence Livermore National Laboratory},
+  note = {Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS)},
+  url = "https://amostech.com/TechnicalPapers/2023/Poster/Pruett.pdf"
+}
+
+@Article{Harris2020,
+ title         = {Array programming with {NumPy}},
+ author        = {Charles R. Harris and K. Jarrod Millman and St{\'{e}}fan J.
+                 van der Walt and Ralf Gommers and Pauli Virtanen and David
+                 Cournapeau and Eric Wieser and Julian Taylor and Sebastian
+                 Berg and Nathaniel J. Smith and Robert Kern and Matti Picus
+                 and Stephan Hoyer and Marten H. van Kerkwijk and Matthew
+                 Brett and Allan Haldane and Jaime Fern{\'{a}}ndez del
+                 R{\'{i}}o and Mark Wiebe and Pearu Peterson and Pierre
+                 G{\'{e}}rard-Marchant and Kevin Sheppard and Tyler Reddy and
+                 Warren Weckesser and Hameer Abbasi and Christoph Gohlke and
+                 Travis E. Oliphant},
+ year          = {2020},
+ month         = sep,
+ journal       = {Nature},
+ volume        = {585},
+ number        = {7825},
+ pages         = {357--362},
+ doi           = {10.1038/s41586-020-2649-2},
+ publisher     = {Springer Science and Business Media {LLC}},
+ url           = {https://doi.org/10.1038/s41586-020-2649-2}
+}
+
+@Article{Hunter2007,
+  Author    = {Hunter, J. D.},
+  Title     = {Matplotlib: A 2D graphics environment},
+  Journal   = {Computing in Science \& Engineering},
+  Volume    = {9},
+  Number    = {3},
+  Pages     = {90--95},
+  abstract  = {Matplotlib is a 2D graphics package used for Python for
+  application development, interactive scripting, and publication-quality
+  image generation across user interfaces and operating systems.},
+  publisher = {IEEE COMPUTER SOC},
+  doi       = {10.1109/MCSE.2007.55},
+  year      = 2007
+}
+
+@ARTICLE{Virtanen2020,
+       author = {{Virtanen}, Pauli and {Gommers}, Ralf and {Oliphant}, Travis E. and {Haberland}, Matt and {Reddy}, Tyler and {Cournapeau}, David and {Burovski}, Evgeni and {Peterson}, Pearu and {Weckesser}, Warren and {Bright}, Jonathan and {van der Walt}, St{\'e}fan J. and {Brett}, Matthew and {Wilson}, Joshua and {Millman}, K. Jarrod and {Mayorov}, Nikolay and {Nelson}, Andrew R.~J. and {Jones}, Eric and {Kern}, Robert and {Larson}, Eric and {Carey}, C.~J. and {Polat}, {\.I}lhan and {Feng}, Yu and {Moore}, Eric W. and {VanderPlas}, Jake and {Laxalde}, Denis and {Perktold}, Josef and {Cimrman}, Robert and {Henriksen}, Ian and {Quintero}, E.~A. and {Harris}, Charles R. and {Archibald}, Anne M. and {Ribeiro}, Ant{\^o}nio H. and {Pedregosa}, Fabian and {van Mulbregt}, Paul and {SciPy 1. 0 Contributors}},
+        title = "{SciPy 1.0: fundamental algorithms for scientific computing in Python}",
+      journal = {Nature Methods},
+     keywords = {Computer Science - Mathematical Software, Computer Science - Data Structures and Algorithms, Computer Science - Software Engineering, Physics - Computational Physics},
+         year = 2020,
+        month = feb,
+       volume = {17},
+        pages = {261-272},
+          doi = {10.1038/s41592-019-0686-2},
+archivePrefix = {arXiv},
+       eprint = {1907.10121},
+ primaryClass = {cs.MS},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2020NatMe..17..261V},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{ForemanMackey2013,
+       author = {{Foreman-Mackey}, Daniel and {Hogg}, David W. and {Lang}, Dustin and {Goodman}, Jonathan},
+        title = "{emcee: The MCMC Hammer}",
+      journal = {\pasp},
+     keywords = {Astrophysics - Instrumentation and Methods for Astrophysics, Physics - Computational Physics, Statistics - Computation},
+         year = 2013,
+        month = mar,
+       volume = {125},
+       number = {925},
+        pages = {306},
+          doi = {10.1086/670067},
+archivePrefix = {arXiv},
+       eprint = {1202.3665},
+ primaryClass = {astro-ph.IM},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2013PASP..125..306F},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@article{astropy2013,
+Adsnote = {Provided by the SAO/NASA Astrophysics Data System},
+Adsurl = {http://adsabs.harvard.edu/abs/2013A%26A...558A..33A},
+Archiveprefix = {arXiv},
+Author = {{Astropy Collaboration} and {Robitaille}, T.~P. and {Tollerud}, E.~J. and {Greenfield}, P. and {Droettboom}, M. and {Bray}, E. and {Aldcroft}, T. and {Davis}, M. and {Ginsburg}, A. and {Price-Whelan}, A.~M. and {Kerzendorf}, W.~E. and {Conley}, A. and {Crighton}, N. and {Barbary}, K. and {Muna}, D. and {Ferguson}, H. and {Grollier}, F. and {Parikh}, M.~M. and {Nair}, P.~H. and {Unther}, H.~M. and {Deil}, C. and {Woillez}, J. and {Conseil}, S. and {Kramer}, R. and {Turner}, J.~E.~H. and {Singer}, L. and {Fox}, R. and {Weaver}, B.~A. and {Zabalza}, V. and {Edwards}, Z.~I. and {Azalee Bostroem}, K. and {Burke}, D.~J. and {Casey}, A.~R. and {Crawford}, S.~M. and {Dencheva}, N. and {Ely}, J. and {Jenness}, T. and {Labrie}, K. and {Lim}, P.~L. and {Pierfederici}, F. and {Pontzen}, A. and {Ptak}, A. and {Refsdal}, B. and {Servillat}, M. and {Streicher}, O.},
+Doi = {10.1051/0004-6361/201322068},
+Eid = {A33},
+Eprint = {1307.6212},
+Journal = {\aap},
+Keywords = {methods: data analysis, methods: miscellaneous, virtual observatory tools},
+Month = oct,
+Pages = {A33},
+Primaryclass = {astro-ph.IM},
+Title = {{Astropy: A community Python package for astronomy}},
+Volume = 558,
+Year = 2013,
+Bdsk-Url-1 = {https://dx.doi.org/10.1051/0004-6361/201322068}}
+
+@ARTICLE{astropy2018,
+       author = {{Astropy Collaboration} and {Price-Whelan}, A.~M. and
+         {Sip{\H{o}}cz}, B.~M. and {G{\"u}nther}, H.~M. and {Lim}, P.~L. and
+         {Crawford}, S.~M. and {Conseil}, S. and {Shupe}, D.~L. and
+         {Craig}, M.~W. and {Dencheva}, N. and {Ginsburg}, A. and {Vand
+        erPlas}, J.~T. and {Bradley}, L.~D. and {P{\'e}rez-Su{\'a}rez}, D. and
+         {de Val-Borro}, M. and {Aldcroft}, T.~L. and {Cruz}, K.~L. and
+         {Robitaille}, T.~P. and {Tollerud}, E.~J. and {Ardelean}, C. and
+         {Babej}, T. and {Bach}, Y.~P. and {Bachetti}, M. and {Bakanov}, A.~V. and
+         {Bamford}, S.~P. and {Barentsen}, G. and {Barmby}, P. and
+         {Baumbach}, A. and {Berry}, K.~L. and {Biscani}, F. and {Boquien}, M. and
+         {Bostroem}, K.~A. and {Bouma}, L.~G. and {Brammer}, G.~B. and
+         {Bray}, E.~M. and {Breytenbach}, H. and {Buddelmeijer}, H. and
+         {Burke}, D.~J. and {Calderone}, G. and {Cano Rodr{\'\i}guez}, J.~L. and
+         {Cara}, M. and {Cardoso}, J.~V.~M. and {Cheedella}, S. and {Copin}, Y. and
+         {Corrales}, L. and {Crichton}, D. and {D'Avella}, D. and {Deil}, C. and
+         {Depagne}, {\'E}. and {Dietrich}, J.~P. and {Donath}, A. and
+         {Droettboom}, M. and {Earl}, N. and {Erben}, T. and {Fabbro}, S. and
+         {Ferreira}, L.~A. and {Finethy}, T. and {Fox}, R.~T. and
+         {Garrison}, L.~H. and {Gibbons}, S.~L.~J. and {Goldstein}, D.~A. and
+         {Gommers}, R. and {Greco}, J.~P. and {Greenfield}, P. and
+         {Groener}, A.~M. and {Grollier}, F. and {Hagen}, A. and {Hirst}, P. and
+         {Homeier}, D. and {Horton}, A.~J. and {Hosseinzadeh}, G. and {Hu}, L. and
+         {Hunkeler}, J.~S. and {Ivezi{\'c}}, {\v{Z}}. and {Jain}, A. and
+         {Jenness}, T. and {Kanarek}, G. and {Kendrew}, S. and {Kern}, N.~S. and
+         {Kerzendorf}, W.~E. and {Khvalko}, A. and {King}, J. and {Kirkby}, D. and
+         {Kulkarni}, A.~M. and {Kumar}, A. and {Lee}, A. and {Lenz}, D. and
+         {Littlefair}, S.~P. and {Ma}, Z. and {Macleod}, D.~M. and
+         {Mastropietro}, M. and {McCully}, C. and {Montagnac}, S. and
+         {Morris}, B.~M. and {Mueller}, M. and {Mumford}, S.~J. and {Muna}, D. and
+         {Murphy}, N.~A. and {Nelson}, S. and {Nguyen}, G.~H. and
+         {Ninan}, J.~P. and {N{\"o}the}, M. and {Ogaz}, S. and {Oh}, S. and
+         {Parejko}, J.~K. and {Parley}, N. and {Pascual}, S. and {Patil}, R. and
+         {Patil}, A.~A. and {Plunkett}, A.~L. and {Prochaska}, J.~X. and
+         {Rastogi}, T. and {Reddy Janga}, V. and {Sabater}, J. and
+         {Sakurikar}, P. and {Seifert}, M. and {Sherbert}, L.~E. and
+         {Sherwood-Taylor}, H. and {Shih}, A.~Y. and {Sick}, J. and
+         {Silbiger}, M.~T. and {Singanamalla}, S. and {Singer}, L.~P. and
+         {Sladen}, P.~H. and {Sooley}, K.~A. and {Sornarajah}, S. and
+         {Streicher}, O. and {Teuben}, P. and {Thomas}, S.~W. and
+         {Tremblay}, G.~R. and {Turner}, J.~E.~H. and {Terr{\'o}n}, V. and
+         {van Kerkwijk}, M.~H. and {de la Vega}, A. and {Watkins}, L.~L. and
+         {Weaver}, B.~A. and {Whitmore}, J.~B. and {Woillez}, J. and
+         {Zabalza}, V. and {Astropy Contributors}},
+        title = "{The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package}",
+      journal = {\aj},
+     keywords = {methods: data analysis, methods: miscellaneous, methods: statistical, reference systems, Astrophysics - Instrumentation and Methods for Astrophysics},
+         year = 2018,
+        month = sep,
+       volume = {156},
+       number = {3},
+          eid = {123},
+        pages = {123},
+          doi = {10.3847/1538-3881/aabc4f},
+archivePrefix = {arXiv},
+       eprint = {1801.02634},
+ primaryClass = {astro-ph.IM},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2018AJ....156..123A},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{astropy2022,
+       author = {{Astropy Collaboration} and {Price-Whelan}, Adrian M. and {Lim}, Pey Lian and {Earl}, Nicholas and {Starkman}, Nathaniel and {Bradley}, Larry and {Shupe}, David L. and {Patil}, Aarya A. and {Corrales}, Lia and {Brasseur}, C.~E. and {N{"o}the}, Maximilian and {Donath}, Axel and {Tollerud}, Erik and {Morris}, Brett M. and {Ginsburg}, Adam and {Vaher}, Eero and {Weaver}, Benjamin A. and {Tocknell}, James and {Jamieson}, William and {van Kerkwijk}, Marten H. and {Robitaille}, Thomas P. and {Merry}, Bruce and {Bachetti}, Matteo and {G{"u}nther}, H. Moritz and {Aldcroft}, Thomas L. and {Alvarado-Montes}, Jaime A. and {Archibald}, Anne M. and {B{'o}di}, Attila and {Bapat}, Shreyas and {Barentsen}, Geert and {Baz{'a}n}, Juanjo and {Biswas}, Manish and {Boquien}, M{'e}d{'e}ric and {Burke}, D.~J. and {Cara}, Daria and {Cara}, Mihai and {Conroy}, Kyle E. and {Conseil}, Simon and {Craig}, Matthew W. and {Cross}, Robert M. and {Cruz}, Kelle L. and {D'Eugenio}, Francesco and {Dencheva}, Nadia and {Devillepoix}, Hadrien A.~R. and {Dietrich}, J{"o}rg P. and {Eigenbrot}, Arthur Davis and {Erben}, Thomas and {Ferreira}, Leonardo and {Foreman-Mackey}, Daniel and {Fox}, Ryan and {Freij}, Nabil and {Garg}, Suyog and {Geda}, Robel and {Glattly}, Lauren and {Gondhalekar}, Yash and {Gordon}, Karl D. and {Grant}, David and {Greenfield}, Perry and {Groener}, Austen M. and {Guest}, Steve and {Gurovich}, Sebastian and {Handberg}, Rasmus and {Hart}, Akeem and {Hatfield-Dodds}, Zac and {Homeier}, Derek and {Hosseinzadeh}, Griffin and {Jenness}, Tim and {Jones}, Craig K. and {Joseph}, Prajwel and {Kalmbach}, J. Bryce and {Karamehmetoglu}, Emir and {Ka{l}uszy{'n}ski}, Miko{l}aj and {Kelley}, Michael S.~P. and {Kern}, Nicholas and {Kerzendorf}, Wolfgang E. and {Koch}, Eric W. and {Kulumani}, Shankar and {Lee}, Antony and {Ly}, Chun and {Ma}, Zhiyuan and {MacBride}, Conor and {Maljaars}, Jakob M. and {Muna}, Demitri and {Murphy}, N.~A. and {Norman}, Henrik and {O'Steen}, Richard and {Oman}, Kyle A. and {Pacifici}, Camilla and {Pascual}, Sergio and {Pascual-Granado}, J. and {Patil}, Rohit R. and {Perren}, Gabriel I. and {Pickering}, Timothy E. and {Rastogi}, Tanuj and {Roulston}, Benjamin R. and {Ryan}, Daniel F. and {Rykoff}, Eli S. and {Sabater}, Jose and {Sakurikar}, Parikshit and {Salgado}, Jes{'u}s and {Sanghi}, Aniket and {Saunders}, Nicholas and {Savchenko}, Volodymyr and {Schwardt}, Ludwig and {Seifert-Eckert}, Michael and {Shih}, Albert Y. and {Jain}, Anany Shrey and {Shukla}, Gyanendra and {Sick}, Jonathan and {Simpson}, Chris and {Singanamalla}, Sudheesh and {Singer}, Leo P. and {Singhal}, Jaladh and {Sinha}, Manodeep and {Sip{H{o}}cz}, Brigitta M. and {Spitler}, Lee R. and {Stansby}, David and {Streicher}, Ole and {{{S}}umak}, Jani and {Swinbank}, John D. and {Taranu}, Dan S. and {Tewary}, Nikita and {Tremblay}, Grant R. and {Val-Borro}, Miguel de and {Van Kooten}, Samuel J. and {Vasovi{'c}}, Zlatan and {Verma}, Shresth and {de Miranda Cardoso}, Jos{'e} Vin{'i}cius and {Williams}, Peter K.~G. and {Wilson}, Tom J. and {Winkel}, Benjamin and {Wood-Vasey}, W.~M. and {Xue}, Rui and {Yoachim}, Peter and {Zhang}, Chen and {Zonca}, Andrea and {Astropy Project Contributors}},
+        title = "{The Astropy Project: Sustaining and Growing a Community-oriented Open-source Project and the Latest Major Release (v5.0) of the Core Package}",
+      journal = {\apj},
+     keywords = {Astronomy software, Open source software, Astronomy data analysis, 1855, 1866, 1858, Astrophysics - Instrumentation and Methods for Astrophysics},
+         year = 2022,
+        month = aug,
+       volume = {935},
+       number = {2},
+          eid = {167},
+        pages = {167},
+          doi = {10.3847/1538-4357/ac7c74},
+archivePrefix = {arXiv},
+       eprint = {2206.14220},
+ primaryClass = {astro-ph.IM},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2022ApJ...935..167A},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@INPROCEEDINGS{2023amos.conf..208Y,
+       author = {{Yeager}, T. and {Pruett}, K. and {Schneider}, M.},
+        title = "{Long-term N-body Stability in Cislunar Space}",
+     keywords = {cislunar, n-body, dynamics, orbit, stability},
+    booktitle = {Proceedings of the Advanced Maui Optical and Space Surveillance (AMOS) Technologies Conference},
+         year = 2023,
+       editor = {{Ryan}, S.},
+        month = sep,
+          eid = {208},
+        pages = {208},
+          url = "https://amostech.com/TechnicalPapers/2023/Poster/Yeager.pdf",
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2023amos.conf..208Y},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
diff --git a/JOSS/paper.md b/JOSS/paper.md
new file mode 100644
index 0000000..a63492b
--- /dev/null
+++ b/JOSS/paper.md
@@ -0,0 +1,132 @@
+---
+title: 'SSAPy - Space Situational Awareness for Python'
+tags:
+  - Python
+  - space domain awareness
+  - orbits
+  - cislunar space
+authors:
+  - name: Joshua E. Meyers
+    affiliation: [1, 2]
+    orcid: 0000-0002-2308-4230
+  - name: Michael D. Schneider
+    affiliation: 3
+    orcid: 0000-0002-8505-7094
+  - name: Julia T. Ebert
+    affiliation: 4
+    orcid: 0000-0002-1975-772X
+  - name: Edward F. Schlafly
+    affiliation: 5
+    orcid: 0000-0002-3569-7421
+  - name: Travis Yeager
+    affiliation: 3
+    orcid: 0000-0002-2582-0190
+  - name: Alexx Perloff
+    affiliation: 3
+    orcid: 0000-0001-5230-0396
+  - name: Daniel Merl
+    affiliation: 3
+    orcid: 0000-0003-4196-5354
+  - name: Noah Lifset
+    affiliation: 6
+    orcid: 0000-0003-3397-7021
+  - name: Jason Bernstein
+    affiliation: 3
+    orcid: 0000-0002-3391-5931
+  - name: William A. Dawson
+    affiliation: 3
+    orcid: 0000-0003-0248-6123
+  - name: Nathan Golovich
+    affiliation: 3
+    orcid: 0000-0003-2632-572X
+  - name: Denvir Higgins
+    affiliation: 3
+    orcid: 0000-0002-7579-1092
+  - name: Peter McGill
+    affiliation: 3
+    orcid: 0000-0002-1052-6749
+    corresponding: true
+  - name: Caleb Miller
+    affiliation: 3
+    orcid: 0000-0001-6249-0031
+  - name: Kerianne Pruett
+    affiliation: 3
+    orcid: 0000-0002-2911-8657
+affiliations:
+ - name: SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA
+   index: 1
+ - name: Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA 94035, USA
+   index: 2
+ - name: Lawrence Livermore National Laboratory, 7000 East Ave., Livermore, CA 94550, USA
+   index: 3
+ - name: Fleet Robotics, 21 Properzi Way, Somerville, MA 02143, USA
+   index: 4
+ - name: Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
+   index: 5
+ - name: University of Texas at Austin, 2515 Speedway, Austin, TX 78712, USA
+   index: 6
+
+
+
+date: 13 January 2025
+bibliography: paper.bib
+
+aas-doi:
+aas-journal:
+---
+
+# Summary
+
+SSAPy is a fast and flexible orbit modeling and analysis tool for orbits spanning from
+low-Earth into the cislunar regime. Orbits can be flexibly specified from common
+input formats such as Keplerian elements or two-line
+element data files. SSAPy allows users to model satellites and specify parameters such
+as satellite area, mass, and drag coefficients. SSAPy includes a customizable force propagation
+with a range of Earth, Lunar, radiation, atmospheric, and maneuvering models. SSAPy makes
+use of various community integration methods and can calculate
+time-evolved orbital quantities, including satellite magnitudes and state vectors.
+Users can specify various space- and ground-based observation models with support for
+multiple coordinate and reference frames. SSAPy also supports orbit analysis and
+propagation methods such as multiple hypothesis tracking and has built-in uncertainty quantification.
+The majority of SSAPy's methods are vectorized and parallelizable, allowing effective use of
+high performance computer (HPC) systems. Finally, SSAPy has plotting functionality, allowing users to
+visualize orbits and trajectories, an example of which is shown in Figures 1 and 2.
+
+SSAPy has been used for the
+classification of cislunar [@Higgins2024], and closely-spaced [@Pruett2024], orbits as
+well as for studying the long-term stability of orbits in cislunar space [@Yeager2023]. SSAPy
+has also been used to build a case study for rare events analysis in the context of satellites
+passing close to each other in space [@Miller2022;@Bernstein2021].
+
+# Statement of need
+
+Cislunar space is a region between earth out to beyond the Moon's orbit that includes the
+Lagrange points. This region of space is of growing importance to scientific and other space exploration endeavors [e.g., @Duggan2019].
+Understanding, mapping, and modeling orbits through cislunar space is
+critical to all of these endeavors. The challenge for cislunar orbits is that n-body dynamics (e.g., gravitational forces
+from the Sun, Earth, Moon and other planets) are significant, leading to unpredictable and chaotic orbital motion.
+In this chaotic regime, orbits cannot be reduced to simple parametric descriptions making scalable orbit
+simulation and modeling a critical analysis tool [@Yeager2023]. Current orbit modeling software tools
+are predominantly used via graphical user interfaces (e.g., The General Mission Analysis Tool; @Hughes2014 or the Systems Tool Kit)
+and are not optimized for large scale simulation on HPC systems. Orbital modeling codes that
+can be run on HPC systems (e.g., REBOUND; @Rein2012) lack full observable generation and modeling capabilities
+with uncertainty quantification. SSAPy, with its full-featured modeling framework and scalable, parallelizable
+functionality, fills the gap in the orbital software landscape.
+
+
+![Example SSAPy visualization plot of an orbit ground track over the surface of the Earth. The 48 hour orbit has a semi-major axis of 0.75 GEO (27,000km), an eccentricity 0.2 and an inclination of 45 degrees.](ground_track.png)
+
+![Example SSAPy visualization plot. Simulated trajectory of the orbit shown in Figure 1. The color on this plot represents time.](orbit_plot.png){ width=50% }
+
+# Acknowledgements
+
+`SSAPy` depends on NumPy [@Harris2020], SciPy [@Virtanen2020], matplotlib [@Hunter2007], emcee [@ForemanMackey2013],
+astropy [@astropy2022], pyerfa [@Kerkwijk2023], lmfit [@newville2024], and sqp4 [@Vallado2006].
+We would like to thank Robert Armstrong and Iméne Goumiri for valuable contributions to this project.
+ This work was performed under the auspices of the U.S.
+Department of Energy by Lawrence Livermore National
+Laboratory (LLNL) under Contract DE-AC52-07NA27344.
+The document number is LLNL-JRNL-871602-DRAFT and the code number is LLNL-CODE-862420. SSAPy was developed with support
+from LLNL's Laboratory Directed Research and Development Program under projects 19-SI-004 and 22-ERD-054.
+
+# References