From 6bf7134408daab9fca9b771675b02a1e222d15fd Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Fri, 12 Sep 2025 12:27:20 -0400
Subject: [PATCH 01/43] ENH: added Gussenhoven model

Added the Gussenhoven model for the EAB to the list of models.
---
 ocbpy/boundaries/models.py | 235 +++++++++++++++++++++++++++++++++++++
 1 file changed, 235 insertions(+)

diff --git a/ocbpy/boundaries/models.py b/ocbpy/boundaries/models.py
index 4c36d7fa..751312e9 100644
--- a/ocbpy/boundaries/models.py
+++ b/ocbpy/boundaries/models.py
@@ -9,11 +9,17 @@
 ----------
 .. [8] Starkov, G. V. (1994) Mathematical model of the auroral boundaries,
    Geomagnetism and Aeronomy, English Translation, 34(3), 331-336.
+.. [9] Gussenhoven, M. S. et al. (1983) Systematics of the Equatorward Diffuse
+   Auroral Boundary, J. Geophys. Res, 88(A7), 5692-5708.
+.. [10] Umbach and Jones (2003) A Few Methods for Fitting Circles to Data,
+   IEEE Transactions on Instruments and Measurements, 52(6), pp 1881-1885
 
 """
 
 import numpy as np
 
+from ocbpy import ocb_time
+
 
 def starkov_auroral_boundary(mlt, al=-1, bnd='ocb'):
     """Calculate the location of the auroral boundaries.
@@ -121,3 +127,232 @@ def starkov_coefficient_values(al, coeff_name, bnd):
             log_al**2) + coeff_terms[coeff_name][bnd][3] * (log_al**3)
 
     return coeff
+
+
+def gussenhoven_equatorward_auroral_boundary(mlt, kp=0, model='circle'):
+    """Calculate the location of the diffuse EAB.
+
+    Parameters
+    ----------
+    mlt : float or array-like
+        Magnetic local time in hours
+    kp : float or int
+        Decimal Kp geomagnetic index, ranges from 0 to 9 with "+" translating
+        to plus one third and "-" translating to minus one third (default=0)
+    model : str
+        Model flavour, with 'binned' using the value from the available hourly
+        solutions or NaN if that hour is absent, 'closest' to use the closest
+        hourly bin, or 'circle' to use the value of the circle fit to the
+        available hourly solutions (default='circle')
+
+    Returns
+    -------
+    bnd_lat : float or array-like
+        Location of the boundary in degrees away from the pole in corrected
+        geomagnetic coordinates for the specified magnetic local times.
+
+    Raises
+    ------
+    ValueError
+        If an unknown model method is requested.
+
+    References
+    ----------
+    [9]_
+
+    """
+    closest = True if model.lower() == 'closest' else False
+    
+    # If desired, calculate the integer hour
+    if model.lower() in ['binned', 'closest']:
+        imlt = np.floor(mlt).astype(int)
+
+        if imlt.shape == ():
+            imlt = np.array([imlt])
+    else:
+        imlt = None
+
+    # Get the desired latitude boundaries
+    colats, mlts = gussenhoven_colatitudes(kp, mlt_inds=imlt, closest=closest)
+
+    # If desired, fit the co-latitude boundaries and return the locations
+    # at the exact MLT values
+    if model.lower() in ['binned', 'closest']:
+        bnd_lat = colats
+    elif model.lower() == "circle":
+        # Fit a circle to the boundaries at this Kp
+        phi_cent, r_cent, radius, _ = circle_fit(mlts, colats)
+        
+        # Calculate the circle location at the desired MLTs. Use the positive
+        # solution of the quadratic equation, as radius >= r_cent
+        theta = ocb_time.hr2rad(mlt) - phi_cent
+        bnd_lat = r_cent * np.cos(theta) + np.sqrt(
+            radius**2 - r_cent**2 * (np.sin(theta))**2)
+    else:
+        raise ValueError('unknown model method: {:}'.format(model))
+
+    return bnd_lat
+
+
+def gussenhoven_colatitudes(kp, mlt_inds=None, closest=False):
+    """Calculate the Gussenhoven EAB model co-latitude values.
+
+    Parameters
+    ----------
+    kp : float or int
+        Decimal Kp geomagnetic index, ranges from 0 to 9 with "+" translating
+        to plus one third and "-" translating to minus one third
+    mlt_inds : array-like or NoneType
+        MLT hourly bins at which boundary locations will be returned or None
+        to return all available values (default=None)
+    closest : bool
+        If False, will return NaN for the request for an MLT value that does not
+        exist. If True, will find the closest binned value and return that
+        instead (default=False)
+
+    Returns
+    -------
+    colats : array-like
+        Co-latitude values in degrees away from the pole for the requested
+        MLT bins
+    mlts : array-like
+        MLT bin values for each co-latitude
+
+    References
+    ----------
+    [9]_
+
+    """
+    # Set the function values for each MLT bin
+    offset = {0: 66.1, 1: 65.1, 4: 67.7, 5: 67.8, 6: 68.2, 7: 68.9, 8: 69.3,
+              9: 69.5, 10: 69.6, 11: 70.1, 12: 69.4, 15: 70.9, 16: 71.6,
+              17: 71.1, 18: 71.2, 19: 70.4, 20: 69.4, 21: 68.6, 22: 67.9,
+              23: 67.8}
+    slope = {0: -1.99, 1: -1.55, 4: -1.48, 5: -1.87, 6: -1.90, 7: -1.91,
+             8: -1.87, 9: -1.67, 10: -1.41, 11: -1.25, 12: -0.84, 15: -0.81,
+             16: -1.28, 17: -1.31, 18: -1.74, 19: -1.83, 20: -1.89, 21: -1.86,
+             22: -1.78, 23: -2.07}
+
+    # Ensure the MLT indices are set
+    mlt_keys = np.array(list(offset.keys()))
+    if mlt_inds is None:
+        mlt_inds = mlt_keys
+
+    # Initialize the output
+    colats = np.full(shape=mlt_inds.shape, fill_value=np.nan)
+    mlts = np.asarray(mlt_inds)
+
+    # Cycle through each MLT index
+    for i, imlt in enumerate(mlt_inds):
+        if imlt in mlt_keys:
+            # Set the calculation key
+            colat_key = imlt
+        elif closest:
+            # Find the closest time with a boundary solution
+            iclose = abs(imlt - mlt_keys).argmin()
+            colat_key = mlt_keys[iclose]
+            mlts[i] = colat_key
+        else:
+            colat_key = None
+
+        if colat_key is not None:
+            # Find the closest co-latitude solution
+            colat = 90 - (offset[colat_key] + slope[colat_key] * kp)
+
+            # Ensure there are no values outside of the allowed range
+            if colat < 0.0:
+                colat = 0.0
+            elif colat > 90.0:
+                colat = 90.0
+
+            # Save the output
+            colats[i] = colat
+
+    return colats, mlts
+
+
+def circle_fit(mlt, colat):
+    """Fit a circle to boundary estimates.
+
+    Parameters
+    ----------
+    mlt : array-like
+        Array of MLT values in hours for the boundary locations
+    colat : array-like
+        Array of magnetic co-latitude values in degrees denoting the boundary
+        locations, paired to the `mlt` inputs
+
+    Returns
+    -------
+    phi_cent : float
+        MLT location of the circle centers in radians
+    r_cent : float
+        Co-latitude of the circle center in degrees
+    radius : float
+        Radius of the circle in degrees latitude
+    r_err : float
+        RMS error of the fit
+
+    Raises
+    ------
+    ValueError
+        If the inputs are the wrong shape
+
+    References
+    ----------
+    Section D in [10]_
+
+    """
+    # Ensure the arrays are of equal length
+    phi = ocb_time.hr2rad(np.asarray(mlt))
+    colat = np.asarray(colat)
+
+    if phi.shape != colat.shape:
+        raise ValueError("input MLT and MLat arrays must be the same shape")
+
+    if len(phi.shape) != 1:
+        raise ValueError('routine only supports 1D arrays')
+
+    # Get the cartesian values of the boundary points.  Define the x-axis to
+    # lie along 0 MLT and the y-axis to lie along 6 MLT, with the origin at the
+    # magnetic pole.
+    x_bound = colat * np.cos(phi)
+    y_bound = colat * np.sin(phi)
+
+    # Apply the circle fitting algorithm, a modified least-squares method,
+    # from Umbach and Jones section D
+    sum_x = np.sum(x_bound)
+    sum_y = np.sum(y_bound)
+    sum_x2 = np.sum(x_bound**2)
+    sum_y2 = np.sum(y_bound**2)
+    sum_xy = np.sum((x_bound * y_bound))
+    sum_xy2 = np.sum(x_bound * (y_bound**2))
+    sum_yx2 = np.sum(y_bound * (x_bound**2))
+    sum_x3 = np.sum(x_bound * x_bound * x_bound)
+    sum_y3 = np.sum(y_bound * y_bound * y_bound)
+
+    fit_a = (colat.shape[0] * sum_x2) - sum_x**2
+    fit_b = (colat.shape[0] * sum_xy) - (sum_x * sum_y)
+    fit_c = (colat.shape[0] * sum_y2) - sum_y**2
+    fit_d = 0.5 * ((colat.shape[0] * sum_xy2) - (sum_x * sum_y2)
+                   + (colat.shape[0] * sum_x3) - (sum_x * sum_x2))
+    fit_e = 0.5 * ((colat.shape[0] * sum_yx2) - (sum_y * sum_x2)
+                   + (colat.shape[0] * sum_y3) - (sum_y * sum_y2))
+
+    # Determine the centre of the fitted circle in Cartesian coordinates
+    circle_denom = (fit_a * fit_c) - fit_b**2
+    major_a = ((fit_d * fit_c) - (fit_b * fit_e)) / circle_denom
+    minor_b = ((fit_a * fit_e) - (fit_b * fit_d)) / circle_denom
+
+    # Convert the circle centre to co-latitude in degrees
+    r_cent = np.sqrt(major_a**2 + minor_b**2)
+    phi_cent = np.arctan2(minor_b, major_a)
+
+    # Determine the radius of the fitted circle
+    rad_bound = np.sqrt((x_bound - major_a)**2 + (y_bound - minor_b)**2)
+    radius = rad_bound.mean()
+
+    # Calculate the error of the radius
+    r_err = np.sqrt(np.mean((rad_bound - radius)**2))
+
+    return phi_cent, r_cent, radius, r_err

From 96a47b8cf7b98cc04ed4bc110cd2911dee8282a1 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Fri, 12 Sep 2025 18:08:23 -0400
Subject: [PATCH 02/43] DOC: updated changelog

Added a summary of changes to the changelog.
---
 Changelog.rst | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/Changelog.rst b/Changelog.rst
index 09b9bbb5..66184879 100644
--- a/Changelog.rst
+++ b/Changelog.rst
@@ -3,6 +3,10 @@ Changelog
 
 Summary of all changes made since the first stable release
 
+0.7.0 (XX-XX-2025)
+------------------
+* ENH: Added the Gussenhoven (1983) model for the EAB
+
 0.6.0 (07-07-2025)
 ------------------
 * ENH: Updated the AMPERE boundaries to include 2022-2024, inclusive

From 9879d47c83b44a399ab373857901a63d07d87901 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Fri, 12 Sep 2025 18:08:50 -0400
Subject: [PATCH 03/43] BUG: fixed for list input

Cast data as an array in case a list is input.
---
 ocbpy/boundaries/models.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/ocbpy/boundaries/models.py b/ocbpy/boundaries/models.py
index 751312e9..d873231d 100644
--- a/ocbpy/boundaries/models.py
+++ b/ocbpy/boundaries/models.py
@@ -237,6 +237,8 @@ def gussenhoven_colatitudes(kp, mlt_inds=None, closest=False):
     mlt_keys = np.array(list(offset.keys()))
     if mlt_inds is None:
         mlt_inds = mlt_keys
+    else:
+        mlt_inds = np.asarray(mlt_inds)
 
     # Initialize the output
     colats = np.full(shape=mlt_inds.shape, fill_value=np.nan)

From 73f979392aecab65f50b97dd7b70229cf3c2d79d Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Fri, 12 Sep 2025 18:09:10 -0400
Subject: [PATCH 04/43] TST: added Gussenhoven unit tests

Added unit tests for the new model functions.
---
 ocbpy/tests/test_models.py | 199 ++++++++++++++++++++++++++++++++++++-
 1 file changed, 198 insertions(+), 1 deletion(-)

diff --git a/ocbpy/tests/test_models.py b/ocbpy/tests/test_models.py
index a437d792..d1ea6323 100644
--- a/ocbpy/tests/test_models.py
+++ b/ocbpy/tests/test_models.py
@@ -15,7 +15,7 @@ class TestStarkovModel(unittest.TestCase):
     """"Unit tests for the Starkov 1994 routines."""
 
     def setUp(self):
-        """Initialize the test case by copying over necessary files."""
+        """Initialize the test case by setting some values to test against."""
         self.mlt = np.arange(0, 24, 1)
         self.coeff_out = {
             'A0': {'ocb': -.07, 'eab': 1.16, 'diffuse': 3.44},
@@ -97,3 +97,200 @@ def test_bound_loc_float(self):
                     self.assertLessEqual(lat, self.max_lat[bnd][ia])
 
         return
+
+
+class TestGussenhovelModel(unittest.TestCase):
+    """"Unit tests for the Gussenhoven 1983 routines."""
+
+    def setUp(self):
+        """Initialize the test case."""
+        self.mlt = np.arange(0, 24, 1)
+        self.bad_mlt = [2, 3, 13, 14]
+        self.colat = {0: [23.9, 24.9, 22.3, 22.2, 21.8, 21.1, 20.7, 20.5, 20.4,
+                          19.9, 20.6, 19.1, 18.4, 18.9, 18.8, 19.6, 20.6, 21.4,
+                          22.1, 22.2],
+                      9: [41.81, 38.85, 35.62, 39.03, 38.9 , 38.29, 37.53,
+                          35.53, 33.09, 31.15, 28.16, 26.39, 29.92, 30.69,
+                          34.46, 36.07, 37.61, 38.14, 38.12, 40.83]}
+        return
+
+    def tearDown(self):
+        """Clean up the test environment."""
+        del self.mlt, self.bad_mlt, self.colat
+        return
+
+    def test_gussenhoven_colatitudes_good(self):
+        """Test the colat calculation for MLTs with solutions"""
+
+        for kp in self.colat.keys():
+            with self.subTest(kp=kp):
+                # Calculate the colatitude values with default kwargs
+                out_lat, _ = models.gussenhoven_colatitudes(kp)
+
+                # Compare the output
+                self.assertEqual(len(out_lat), len(self.colat[kp]))
+                diff_lat = out_lat - np.asarray(self.colat[kp])
+                self.assertLess(sum(diff_lat), 1.0e-3)
+        return
+
+    def test_gussenhoven_colatitudes_bad(self):
+        """Test the colat calculation for MLTs without solutions."""
+        # Cycle through each Kp
+        for kp in self.colat.keys():
+            with self.subTest(kp=kp):
+                # Calculate the colatitude values with default kwargs
+                out_lat, out_mlt = models.gussenhoven_colatitudes(
+                    kp, mlt_inds=self.bad_mlt)
+
+                # Compare the output
+                self.assertTrue(np.isnan(out_lat).all())
+                self.assertListEqual(list(out_mlt), self.bad_mlt)
+        return
+
+    def test_gussenhoven_colatitudes_bad_close(self):
+        """Test the colat calculation for MLTs at nearest solutions."""
+        # Reshape the MLT to exclude the bad MLT values
+        self.mlt = list(self.mlt)
+        for imlt in self.bad_mlt:
+            self.mlt.pop(self.mlt.index(imlt))
+
+        # Cycle through each Kp
+        for kp in self.colat.keys():
+            with self.subTest(kp=kp):
+                # Calculate the colatitude values with default kwargs
+                out_lat, out_mlt = models.gussenhoven_colatitudes(
+                    kp, mlt_inds=self.bad_mlt, closest=True)
+
+                # Compare the output
+                self.assertFalse(np.isnan(out_lat).any())
+                self.assertEqual(len(out_mlt), len(self.bad_mlt))
+
+                for i, imlt in enumerate(out_mlt):
+                    self.assertFalse(imlt in self.bad_mlt)
+                    self.assertAlmostEqual(
+                        out_lat[i], self.colat[kp][self.mlt.index(imlt)])
+        return
+
+    def test_bad_model(self):
+        """Test a ValueError is raised for an unknown model type."""
+        model = "not a model"
+        with self.assertRaisesRegex(ValueError, model):
+            models.gussenhoven_equatorward_auroral_boundary(
+                self.mlt, model=model)
+        return
+
+    def test_gussenhoven_eab_binned(self):
+        """Test the EAB for MLTs with binned solutions."""
+        # Cycle through each Kp
+        for kp in self.colat.keys():
+            with self.subTest(kp=kp):
+                # Calculate the colatitude values with default kwargs
+                out_lat = models.gussenhoven_equatorward_auroral_boundary(
+                    self.mlt, kp=kp, model='binned')
+
+                # Compare the output
+                self.assertEqual(len(out_lat[~np.isnan(out_lat)]),
+                                 len(self.colat[kp]))
+                self.assertEqual(len(out_lat[np.isnan(out_lat)]),
+                                 len(self.bad_mlt))
+                diff_lat = out_lat[~np.isnan(out_lat)] - np.asarray(
+                    self.colat[kp])
+                self.assertLess(sum(diff_lat), 1.0e-3)
+        return
+
+    def test_gussenhoven_eab_closest(self):
+        """Test the EAB for MLTs with closest value solutions."""
+        # Cycle through each Kp
+        for kp in self.colat.keys():
+            with self.subTest(kp=kp):
+                # Calculate the colatitude values with default kwargs
+                out_lat = models.gussenhoven_equatorward_auroral_boundary(
+                    self.mlt, kp=kp, model='closest')
+
+                # Compare the output
+                self.assertEqual(len(out_lat), len(self.mlt))
+
+                j = 0
+                for i, imlt in enumerate(self.mlt):
+                    if imlt in self.bad_mlt:
+                        try:
+                            j -= 1
+                            self.assertAlmostEqual(
+                                out_lat[i], self.colat[kp][j])
+                            j += 1
+                        except AssertionError:
+                            j += 1
+                            self.assertAlmostEqual(
+                                out_lat[i], self.colat[kp][j],
+                                msg="failed at {:} MLT with j={:}".format(
+                                    imlt, j))
+                    else:
+                        self.assertAlmostEqual(
+                            out_lat[i], self.colat[kp][j],
+                            msg="failed at {:} MLT with j={:}".format(imlt, j))
+                        j += 1
+        return
+
+    def test_gussenhoven_eab_circle(self):
+        """Test the EAB for MLTs with circle fit solutions."""
+        # Cycle through each Kp
+        for kp in self.colat.keys():
+            with self.subTest(kp=kp):
+                # Calculate the colatitude values with default kwargs
+                out_lat = models.gussenhoven_equatorward_auroral_boundary(
+                    self.mlt, kp=kp, model='circle')
+
+                # Compare the output
+                self.assertEqual(len(out_lat), len(self.mlt))
+
+                j = 0
+                for i, imlt in enumerate(self.mlt):
+                    if imlt in self.bad_mlt:
+                        try:
+                            j -= 1
+                            self.assertLessEqual(
+                                abs(out_lat[i] - self.colat[kp][j]), 5.0)
+                            j += 1
+                        except AssertionError:
+                            j += 1
+                            self.assertLessEqual(
+                                abs(out_lat[i] - self.colat[kp][j]), 5.0,
+                                msg="failed at {:} MLT with j={:}".format(
+                                    imlt, j))
+                    else:
+                        self.assertLessEqual(
+                            abs(out_lat[i] - self.colat[kp][j]), 5.0,
+                            msg="failed at {:} MLT with j={:}".format(imlt, j))
+                        j += 1
+        return
+
+
+
+class TestFits(unittest.TestCase):
+    """"Unit tests for the fitting routines."""
+
+    def setUp(self):
+        """Initialize the test case by setting some values to test against."""
+        self.mlt = np.arange(0, 24, 1)
+        self.rvals = np.ones(shape=self.mlt.shape)
+        return
+
+    def tearDown(self):
+        """Clean up the test environment."""
+        del self.mlt, self.rvals
+        return
+
+    def test_circle_fit(self):
+        """Test the circle fitting to a unit circle."""
+        # Run the fitting routine
+        phi_cent, r_cent, radius, r_err = models.circle_fit(
+            self.mlt, self.rvals)
+
+        # Test the output.  The anglular offset can be any value with a
+        # radial offset of zero.  It should be constrained within 0-2pi
+        self.assertGreaterEqual(phi_cent, 0.0)
+        self.assertLessEqual(phi_cent, 2.0 * np.pi)
+        self.assertAlmostEqual(r_cent, 0.0)
+        self.assertAlmostEqual(radius, 1.0)
+        self.assertAlmostEqual(r_err, 0.0)
+        return

From d6dfd62a22dad6539826a909367d6e6e5f8778ed Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 17 Sep 2025 15:14:57 -0400
Subject: [PATCH 05/43] STY: removed whitespace

Removed extra whitespace.
---
 ocbpy/boundaries/models.py | 4 ++--
 ocbpy/tests/test_models.py | 3 +--
 2 files changed, 3 insertions(+), 4 deletions(-)

diff --git a/ocbpy/boundaries/models.py b/ocbpy/boundaries/models.py
index d873231d..fafc9457 100644
--- a/ocbpy/boundaries/models.py
+++ b/ocbpy/boundaries/models.py
@@ -162,7 +162,7 @@ def gussenhoven_equatorward_auroral_boundary(mlt, kp=0, model='circle'):
 
     """
     closest = True if model.lower() == 'closest' else False
-    
+
     # If desired, calculate the integer hour
     if model.lower() in ['binned', 'closest']:
         imlt = np.floor(mlt).astype(int)
@@ -182,7 +182,7 @@ def gussenhoven_equatorward_auroral_boundary(mlt, kp=0, model='circle'):
     elif model.lower() == "circle":
         # Fit a circle to the boundaries at this Kp
         phi_cent, r_cent, radius, _ = circle_fit(mlts, colats)
-        
+
         # Calculate the circle location at the desired MLTs. Use the positive
         # solution of the quadratic equation, as radius >= r_cent
         theta = ocb_time.hr2rad(mlt) - phi_cent
diff --git a/ocbpy/tests/test_models.py b/ocbpy/tests/test_models.py
index d1ea6323..12189d0d 100644
--- a/ocbpy/tests/test_models.py
+++ b/ocbpy/tests/test_models.py
@@ -109,7 +109,7 @@ def setUp(self):
         self.colat = {0: [23.9, 24.9, 22.3, 22.2, 21.8, 21.1, 20.7, 20.5, 20.4,
                           19.9, 20.6, 19.1, 18.4, 18.9, 18.8, 19.6, 20.6, 21.4,
                           22.1, 22.2],
-                      9: [41.81, 38.85, 35.62, 39.03, 38.9 , 38.29, 37.53,
+                      9: [41.81, 38.85, 35.62, 39.03, 38.90, 38.29, 37.53,
                           35.53, 33.09, 31.15, 28.16, 26.39, 29.92, 30.69,
                           34.46, 36.07, 37.61, 38.14, 38.12, 40.83]}
         return
@@ -265,7 +265,6 @@ def test_gussenhoven_eab_circle(self):
         return
 
 
-
 class TestFits(unittest.TestCase):
     """"Unit tests for the fitting routines."""
 

From b5adeb0f0f9b27da991b1b1086e39cd3f49f87fe Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 17 Sep 2025 15:21:50 -0400
Subject: [PATCH 06/43] DOC: added Gussenhoven model to docs

Added the citation for the Gussenhoven EAB to the Citation section and the model description to the Models section.
---
 docs/citing.rst     | 14 ++++++++++++++
 docs/ocb_models.rst | 13 +++++++++++++
 2 files changed, 27 insertions(+)

diff --git a/docs/citing.rst b/docs/citing.rst
index fee63c5d..1e4fd674 100644
--- a/docs/citing.rst
+++ b/docs/citing.rst
@@ -112,3 +112,17 @@ when publishing your work.
 
 * Starkov, G. V. (1994) Mathematical model of the auroral boundaries,
   Geomagnetism and Aeronomy, English Translation, 34(3), 331-336.
+
+
+.. _cite-guss:
+
+Gussenhoven Model
+-----------------
+
+The Gussenhoven mathematical diffuse auroral boundary model is described in the
+following article. If you use this model please cite both it and your Kp data
+source when publishing your work.
+
+* Gussenhoven, M. S. et al. (1983) Systematics of the Equatorward Diffuse
+  Auroral Boundary, J. Geophys. Res, 88(A7), 5692-5708.
+
diff --git a/docs/ocb_models.rst b/docs/ocb_models.rst
index fc4b3652..e8e5be27 100644
--- a/docs/ocb_models.rst
+++ b/docs/ocb_models.rst
@@ -27,6 +27,19 @@ in the code here using the boundary keyworkds: 'ocb', 'eab', and 'diffuse',
 respectively.
 
 
+.. _bound-model-guss:
+
+Gussenhoven
+-----------
+
+The Gussenhoven 1983 model (see :ref:`cite-guss`) uses a mathematical
+formulation based on DMSP data and the Kp index. They specify a single boundary:
+the equatorward edge of the diffuse aurora. As this model defines the boundary
+using separate linear fits for different MLT bins, the results may be returned
+at the binned times ('binned'), at the nearest binned time ('closest'), or at
+the requested times using a circle fit to all of binned times ('circle').
+
+
 .. _bound-model-module:
 
 Boundary Models Module

From 026cd27fbca468281e8133e2438164ef6b0506e7 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 17 Sep 2025 15:55:17 -0400
Subject: [PATCH 07/43] BUG: fixed doc reference

Fixed a reference in the doc examples.
---
 docs/examples/ex_model_boundaries.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/examples/ex_model_boundaries.rst b/docs/examples/ex_model_boundaries.rst
index a0731001..f4925b5c 100644
--- a/docs/examples/ex_model_boundaries.rst
+++ b/docs/examples/ex_model_boundaries.rst
@@ -32,7 +32,7 @@ be calculated for a set of MLTs using the
 
 This provides the AACGMV2 magnetic latitude for the entire range of MLT at the
 specified AL values. We can plot these boundaries, using the formatting function
-previously defined in Example :ref:`format-polar-axes`.
+previously defined in :ref:`set_up_polar_plot <format-polar-axes>`.
 
 ::
 

From aa81268876c0db6c361049e1d94b62a7773b73dd Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 17 Sep 2025 15:55:31 -0400
Subject: [PATCH 08/43] BUG: fixed angle range

Fixed the test angle range.
---
 ocbpy/tests/test_models.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/ocbpy/tests/test_models.py b/ocbpy/tests/test_models.py
index 12189d0d..a34b9e6f 100644
--- a/ocbpy/tests/test_models.py
+++ b/ocbpy/tests/test_models.py
@@ -286,9 +286,9 @@ def test_circle_fit(self):
             self.mlt, self.rvals)
 
         # Test the output.  The anglular offset can be any value with a
-        # radial offset of zero.  It should be constrained within 0-2pi
-        self.assertGreaterEqual(phi_cent, 0.0)
-        self.assertLessEqual(phi_cent, 2.0 * np.pi)
+        # radial offset of zero.  It should be constrained within +/-pi
+        self.assertGreaterEqual(phi_cent, -np.pi)
+        self.assertLessEqual(phi_cent, np.pi)
         self.assertAlmostEqual(r_cent, 0.0)
         self.assertAlmostEqual(radius, 1.0)
         self.assertAlmostEqual(r_err, 0.0)

From a89d51014dc296bab846f477a352c20cd3b339c1 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Thu, 18 Sep 2025 09:53:52 -0400
Subject: [PATCH 09/43] BUG: coveralls implementation

The coveralls implementation broke on the mac OS.  Updated to a simpler method that will hopefully work now.
---
 .github/workflows/main.yml | 35 ++++++++---------------------------
 1 file changed, 8 insertions(+), 27 deletions(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 36f417fc..1d0bb489 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -50,33 +50,14 @@ jobs:
         coverage report
         coverage xml --rcfile=pyproject.toml
 
-    - name: Install and run Coveralls Reporter(Linux)
-      if: startsWith(matrix.os, 'ubuntu')
-      env:
-        COVERALLS_REPO_TOKEN: ${{ secrets.COVERALLS_REPO_TOKEN }}
-        COVERALLS_PARALLEL: true
-      run: |
-        curl -sL https://coveralls.io/coveralls-linux.tar.gz | tar -xz
-        ./coveralls report -f coverage.xml --parallel --repo-token=${{ secrets.COVERALLS_REPO_TOKEN }} --build-number ${{ github.run_number }}
-
-    - name: Install and run Coveralls Reporter (Windows)
-      if: startsWith(matrix.os, 'windows')
-      env:
-        COVERALLS_REPO_TOKEN: ${{ secrets.COVERALLS_REPO_TOKEN }}
-        COVERALLS_PARALLEL: true
-      run: |
-        curl -L https://github.com/coverallsapp/coverage-reporter/releases/latest/download/coveralls-windows.exe -o coveralls.exe
-        ./coveralls.exe report -f coverage.xml --parallel --repo-token=${{ secrets.COVERALLS_REPO_TOKEN }} --build-number ${{ github.run_number }}
-
-    - name: Report  and run Coveralls (macOS)
-      if: startsWith(matrix.os, 'macos')
-      env:
-        COVERALLS_REPO_TOKEN: ${{ secrets.COVERALLS_REPO_TOKEN }}
-        COVERALLS_PARALLEL: true
-      run: |
-        brew tap coverallsapp/coveralls --quiet
-        brew install coveralls --quiet
-        coveralls report -f coverage.xml --parallel --repo-token=${{ secrets.COVERALLS_REPO_TOKEN }} --build-number ${{ github.run_number }}
+    - name:  Coveralls Parallel
+      uses: coverallsapp/github-action@v2
+      with:
+        github-token: ${{ secrets.GITHUB_TOKEN }}
+        flag-name: run=${{ join(matrix.*, '-') }}
+        parallel: true
+        format: cobertura
+        debug: true
 
   finish:
     name: Finish Coverage Analysis

From cdefb2428394e056831dc94abb0945bc85b2cc98 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 1 Apr 2026 14:37:13 -0400
Subject: [PATCH 10/43] MAINT: removed support for old zenodo_get

Removed the support for the old zenodo_get downloading formulation.
---
 ocbpy/boundaries/dmsp_ssj_files.py | 9 ++-------
 1 file changed, 2 insertions(+), 7 deletions(-)

diff --git a/ocbpy/boundaries/dmsp_ssj_files.py b/ocbpy/boundaries/dmsp_ssj_files.py
index 74279af7..b171f4e0 100644
--- a/ocbpy/boundaries/dmsp_ssj_files.py
+++ b/ocbpy/boundaries/dmsp_ssj_files.py
@@ -113,13 +113,8 @@ def fetch_ssj_boundary_files(stime=None, etime=None, out_dir=None,
     sys.stdout = zenodo_io
     sys.stderr = zenodo_io
 
-    # TODO(#151): remove the old (second) way of calling zenodo_get
-    if hasattr(zenodo_get, "download"):
-        zenodo_get.download(doi=doi, output_dir=out_dir)
-        zenodo_checksum = None
-    else:
-        zenodo_get.zenodo_get([doi, '-o', out_dir])
-        zenodo_checksum = os.path.join(out_dir, 'md5sums.txt')
+    zenodo_get.download(doi=doi, output_dir=out_dir)
+    zenodo_checksum = None
 
     # Parse the output and retrieve files from the zip archive
     sys.stdout = sys.__stdout__

From 57cb632ba621bbdf2a6bf5e0b95a87584a789219 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 1 Apr 2026 14:37:45 -0400
Subject: [PATCH 11/43] MAINT: added a version limit

Added a lower limit for zenodo_get.
---
 pyproject.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index e843883a..57023e4c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -55,14 +55,14 @@ classifiers = [
 
 [project.optional-dependencies]
 pysat_instruments = [ "pysat>=3.2.1" ]
-dmsp_ssj = [ "zenodo-get" ]
+dmsp_ssj = [ "zenodo-get>=2.0.0" ]
 doc = [
     "numpydoc",
     "pyproject-parser",
     "pysat>=3.2.1",
     "sphinx>=1.3",
     "sphinx-rtd-theme",
-    "zenodo-get",
+    "zenodo-get>=2.0.0",
 ]
 test = [
      "coverage[toml]",

From 4e91faaf6bdbfa805b91526a533b9190af5ae900 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 1 Apr 2026 14:43:16 -0400
Subject: [PATCH 12/43] MAINT: added Python 3.14

Added support for Python 3.14
---
 pyproject.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/pyproject.toml b/pyproject.toml
index 57023e4c..d9ed3035 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -46,6 +46,7 @@ classifiers = [
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
     "Programming Language :: Python :: 3.13",
+    "Programming Language :: Python :: 3.14",
     'Operating System :: Unix',
     'Operating System :: POSIX',
     'Operating System :: POSIX :: Linux',

From 2c3d75c7fce8af1b552db97fd15e70bf437f36fd Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 1 Apr 2026 14:43:45 -0400
Subject: [PATCH 13/43] MAINT: added support for Python 3.14

Added testing for Python 3.14 and updated GitHub action versions.
---
 .github/workflows/main.yml | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 1d0bb489..1a07e83d 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -13,15 +13,15 @@ jobs:
       fail-fast: false
       matrix:
         os: ["ubuntu-latest", "macos-latest", "windows-latest"]
-        python-version: ["3.10", "3.11", "3.12", "3.13"]
+        python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"]
         install-extras: ["base", "pysat_instruments", "dmsp_ssj"]
 
     name: Python ${{ matrix.python-version }} on ${{ matrix.os }} with ${{ matrix.install-extras }}
     runs-on: ${{ matrix.os }}
     steps:
-    - uses: actions/checkout@v4
+    - uses: actions/checkout@v6
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v5
+      uses: actions/setup-python@v6
       with:
         python-version: ${{ matrix.python-version }}
 

From 3b0e1d103cd6e56a0b7e8cbb19acf5de4f4a9a26 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 1 Apr 2026 14:44:21 -0400
Subject: [PATCH 14/43] MAINT: updated GitHub Action versions

Updated the general and Python GitHub action versions to the most recent ones.
---
 .github/workflows/docs.yml           | 6 +++---
 .github/workflows/pip_rc_install.yml | 6 +++---
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 78e36ed9..c7a1634c 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -15,13 +15,13 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.11"]
+        python-version: ["3.14"]
 
     name: Documentation tests
     steps:
-    - uses: actions/checkout@v4
+    - uses: actions/checkout@v6
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v5
+      uses: actions/setup-python@v6
       with:
         python-version: ${{ matrix.python-version }}
 
diff --git a/.github/workflows/pip_rc_install.yml b/.github/workflows/pip_rc_install.yml
index 1124b9ba..d54da2ba 100644
--- a/.github/workflows/pip_rc_install.yml
+++ b/.github/workflows/pip_rc_install.yml
@@ -13,14 +13,14 @@ jobs:
       fail-fast: false
       matrix:
         os: ["ubuntu-latest", "macos-latest", "windows-latest"]
-        python-version: ["3.11"]  # Keep this version at the highest supported Python version
+        python-version: ["3.14"]  # Keep this version at the highest supported Python version
 
     name: Python ${{ matrix.python-version }} on ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
     steps:
-    - uses: actions/checkout@v4
+    - uses: actions/checkout@v6
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v5
+      uses: actions/setup-python@v6
       with:
         python-version: ${{ matrix.python-version }}
 

From 0ab8802581edefa3d6e988a24d7aa82f500e3e78 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 1 Apr 2026 14:45:52 -0400
Subject: [PATCH 15/43] DOC: update changelog

Add a summary of changes to the changelog.
---
 Changelog.rst | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/Changelog.rst b/Changelog.rst
index 66184879..a066ae82 100644
--- a/Changelog.rst
+++ b/Changelog.rst
@@ -6,6 +6,9 @@ Summary of all changes made since the first stable release
 0.7.0 (XX-XX-2025)
 ------------------
 * ENH: Added the Gussenhoven (1983) model for the EAB
+* DEP: Removed support for older versions of `zenodo_get`
+* MAINT: Added support for Python 3.14
+* MAINT: Updated GitHub CI actions
 
 0.6.0 (07-07-2025)
 ------------------

From cc6c36870f54a48374d3be9e55d4aa6fa84046b4 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 1 Apr 2026 14:48:00 -0400
Subject: [PATCH 16/43] DOC: update supported versions

Update the supported versions for Python and dependencies in the docs.
---
 README.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index 1abf5c72..4e767f25 100644
--- a/README.md
+++ b/README.md
@@ -45,7 +45,7 @@ These routines may be used as a guide to write routines for other datasets.
 
 # Python versions
 
-This module currently supports Python version 3.10 - 3.13.
+This module currently supports Python version 3.10 - 3.14.
 
 # Dependencies
 
@@ -53,7 +53,7 @@ The listed dependecies were tested with the following versions:
   * numpy
   * aacgmv2
   * pysat (3.2.1+)
-  * zenodo_get
+  * zenodo_get (2.0.0+)
 
 Testing is performed using the python module, unittest.  To limit dependency
 issues, the pysat and zenodo_get dependencies are optional.

From 351e7c2c01940c5192047306376026ce801ccc16 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 1 Apr 2026 16:51:37 -0400
Subject: [PATCH 17/43] MAINT: numpy deprecation warnings

Addressed numpy deprecation warnings for:
- not providing output location in greater/less evaluations, and
- extracting an item from a zero-sized array.
---
 ocbpy/ocb_scaling.py              |  6 +--
 ocbpy/ocb_time.py                 | 12 ++---
 ocbpy/tests/test_boundary_dual.py | 12 ++---
 ocbpy/tests/test_boundary_ocb.py  | 16 +++----
 ocbpy/tests/test_ocb_scaling.py   | 18 ++++----
 ocbpy/vectors.py                  | 74 +++++++++++++++++--------------
 6 files changed, 74 insertions(+), 64 deletions(-)

diff --git a/ocbpy/ocb_scaling.py b/ocbpy/ocb_scaling.py
index 0378157a..965a7f32 100644
--- a/ocbpy/ocb_scaling.py
+++ b/ocbpy/ocb_scaling.py
@@ -473,7 +473,7 @@ def vect_mag(self, vect_mag):
             self._vect_mag = vect_sqrt
         else:
             if np.any(np.greater(abs(vect_mag - vect_sqrt), 1.0e-3,
-                                 where=~np.isnan(vect_mag))):
+                                 where=~np.isnan(vect_mag), out=None)):
                 ocbpy.logger.warning("".join([
                     "inconsistent vector components with a maximum difference ",
                     "of {:} > 1.0e-3".format(abs(vect_mag - vect_sqrt).max())]))
@@ -1221,7 +1221,7 @@ def archav(hav):
     alpha = np.full(shape=hav.shape, fill_value=np.nan)
 
     # If the number is positive, calculate the angle
-    norm_mask = (np.greater_equal(hav, 1.0e-16, where=~np.isnan(hav))
+    norm_mask = (np.greater_equal(hav, 1.0e-16, where=~np.isnan(hav), out=None)
                  & ~np.isnan(hav))
     if np.any(norm_mask):
         if hav.shape == ():
@@ -1231,7 +1231,7 @@ def archav(hav):
 
     #  The number is small enough that machine precision may have changed
     # the sign, but it's a single-precission zero
-    small_mask = (np.less(abs(hav), 1.0e-16, where=~np.isnan(hav))
+    small_mask = (np.less(abs(hav), 1.0e-16, where=~np.isnan(hav), out=None)
                   & ~np.isnan(hav))
     if np.any(small_mask):
         if hav.shape == ():
diff --git a/ocbpy/ocb_time.py b/ocbpy/ocb_time.py
index db915b27..9d0b46b7 100644
--- a/ocbpy/ocb_time.py
+++ b/ocbpy/ocb_time.py
@@ -434,19 +434,19 @@ def fix_range(values, min_val, max_val, val_range=None):
 
     # Fix the values, allowing for deviations that are multiples of the
     # value range.  Also propagate NaNs
-    ibad = (np.greater_equal(fixed_vals, max_val, where=~np.isnan(fixed_vals))
-            & ~np.isnan(fixed_vals))
+    ibad = np.greater_equal(fixed_vals, max_val, where=~np.isnan(fixed_vals),
+                            out=None) & ~np.isnan(fixed_vals)
     while np.any(ibad):
         fixed_vals[ibad] -= val_range
         ibad = (np.greater_equal(fixed_vals, max_val,
-                                 where=~np.isnan(fixed_vals))
+                                 where=~np.isnan(fixed_vals), out=None)
                 & ~np.isnan(fixed_vals))
 
-    ibad = (np.less(fixed_vals, min_val, where=~np.isnan(fixed_vals))
+    ibad = (np.less(fixed_vals, min_val, where=~np.isnan(fixed_vals), out=None)
             & ~np.isnan(fixed_vals))
     while np.any(ibad):
         fixed_vals[ibad] += val_range
-        ibad = (np.less(fixed_vals, min_val, where=~np.isnan(fixed_vals))
-                & ~np.isnan(fixed_vals))
+        ibad = np.less(fixed_vals, min_val, where=~np.isnan(fixed_vals),
+                       out=None) & ~np.isnan(fixed_vals)
 
     return fixed_vals
diff --git a/ocbpy/tests/test_boundary_dual.py b/ocbpy/tests/test_boundary_dual.py
index 5982e5ae..b64c8da3 100644
--- a/ocbpy/tests/test_boundary_dual.py
+++ b/ocbpy/tests/test_boundary_dual.py
@@ -586,7 +586,7 @@ def eval_coords(self, hemisphere=1, tol=1.0e-7, ind=None, revert=False,
                              len(numpy.isnan(rlat)))
             if not numpy.isnan(rlat).all():
                 self.assertTrue(
-                    numpy.less(abs(self.out[0] - rlat), tol,
+                    numpy.less(abs(self.out[0] - rlat), tol, out=None,
                                where=~numpy.isnan(rlat)).all(),
                     msg="unequal {:s}: {:} != {:}".format(lat_str, self.out[0],
                                                           rlat))
@@ -594,7 +594,7 @@ def eval_coords(self, hemisphere=1, tol=1.0e-7, ind=None, revert=False,
                              len(numpy.isnan(rmlt)))
             if not numpy.isnan(rmlt).all():
                 self.assertTrue(
-                    numpy.less(abs(self.out[1] - rmlt), tol,
+                    numpy.less(abs(self.out[1] - rmlt), tol, out=None,
                                where=~numpy.isnan(rmlt)).all(),
                     msg="unequal {:s}: {:} != {:}".format(mlt_str, self.out[1],
                                                           rmlt))
@@ -605,7 +605,8 @@ def eval_coords(self, hemisphere=1, tol=1.0e-7, ind=None, revert=False,
                 if not numpy.isnan(self.out[2]).all():
                     self.assertTrue(
                         numpy.less(abs(self.out[2] - self.olat[hemisphere]),
-                                   tol, where=~numpy.isnan(self.out[2])).all(),
+                                   tol, where=~numpy.isnan(self.out[2]),
+                                   out=None).all(),
                         msg="unequal OCB latitude: {:} != {:}".format(
                             self.out[2], self.olat[hemisphere]))
 
@@ -615,7 +616,8 @@ def eval_coords(self, hemisphere=1, tol=1.0e-7, ind=None, revert=False,
                 else:
                     self.assertTrue(
                         numpy.less(abs(self.out[3] - self.rcorr), tol,
-                                   where=~numpy.isnan(self.out[3])).all(),
+                                   where=~numpy.isnan(self.out[3]),
+                                   out=None).all(),
                         msg="unequal radial correction: {:} != {:}".format(
                             self.out[3], self.rcorr))
         else:
@@ -1119,7 +1121,7 @@ def test_get_current_aacgm_boundary_set(self):
 
                 # Ensure the expected values and fill values are returned
                 self.assertTrue(
-                    numpy.less(abs(abound - self.bounds[i]), 1e-7,
+                    numpy.less(abs(abound - self.bounds[i]), 1e-7, out=None,
                                where=~numpy.isnan(abound)).all(),
                     msg="unexpected boundary: {:} != {:}".format(
                         abound, self.bounds[i]))
diff --git a/ocbpy/tests/test_boundary_ocb.py b/ocbpy/tests/test_boundary_ocb.py
index 467ed8c1..993b6da3 100644
--- a/ocbpy/tests/test_boundary_ocb.py
+++ b/ocbpy/tests/test_boundary_ocb.py
@@ -590,11 +590,11 @@ def test_normal_coord_north_array(self):
         self.out = self.ocb.normal_coord(self.lat, self.mlt)
 
         self.assertTrue(numpy.all(numpy.less(abs(self.out[0] - self.ocb_lat),
-                                             1.0e-7,
+                                             1.0e-7, out=None,
                                              where=~numpy.isnan(self.out[0]))
                                   | numpy.isnan(self.out[0])))
         self.assertTrue(numpy.all(numpy.less(abs(self.out[1] - self.ocb_mlt),
-                                             1.0e-7,
+                                             1.0e-7, out=None,
                                              where=(~numpy.isnan(self.out[1])))
                                   | numpy.isnan(self.out[1])))
         self.assertTrue(numpy.where(numpy.isnan(self.out[0]))
@@ -679,11 +679,11 @@ def test_revert_coord_north_array(self):
                                          self.r_corr)
 
         self.assertTrue(numpy.all(numpy.less(abs(self.out[0] - self.lat),
-                                             1.0e-7,
+                                             1.0e-7, out=None,
                                              where=~numpy.isnan(self.out[0]))
                                   | (numpy.isnan(self.out[0]))))
         self.assertTrue(numpy.all(numpy.less(abs(self.out[1] - self.mlt),
-                                             1.0e-7,
+                                             1.0e-7, out=None,
                                              where=(~numpy.isnan(self.out[1])
                                                     & (self.lat < 90.0)))
                                   | numpy.isnan(self.out[0])
@@ -984,10 +984,10 @@ def test_normal_coord_south_array(self):
         self.out = self.ocb.normal_coord(self.lat, self.mlt)
 
         self.assertTrue(numpy.all(numpy.less(abs(self.out[0] - self.ocb_lat),
-                                             1.0e-7,
+                                             1.0e-7, out=None,
                                              where=~numpy.isnan(self.out[0]))))
         self.assertTrue(numpy.all(numpy.less(abs(self.out[1] - self.ocb_mlt),
-                                             1.0e-7,
+                                             1.0e-7, out=None,
                                              where=~numpy.isnan(self.out[1]))))
         self.assertTrue(numpy.all(numpy.where(numpy.isnan(self.out[0]))[0]
                                   == numpy.where(numpy.isnan(self.ocb_lat))[0]))
@@ -1041,11 +1041,11 @@ def test_revert_coord_south_array(self):
                                          self.r_corr)
 
         self.assertTrue(numpy.all(numpy.less(abs(self.out[0] - self.lat),
-                                             1.0e-7,
+                                             1.0e-7, out=None,
                                              where=~numpy.isnan(self.out[0]))
                                   | numpy.isnan(self.out[0])))
         self.assertTrue(numpy.all(numpy.less(abs(self.out[1] - self.mlt),
-                                             1.0e-7,
+                                             1.0e-7, out=None,
                                              where=(~numpy.isnan(self.out[1])
                                                     & (self.lat > -90.0)))
                                   | numpy.isnan(self.out[0])
diff --git a/ocbpy/tests/test_ocb_scaling.py b/ocbpy/tests/test_ocb_scaling.py
index f2936023..e326ff97 100644
--- a/ocbpy/tests/test_ocb_scaling.py
+++ b/ocbpy/tests/test_ocb_scaling.py
@@ -540,10 +540,10 @@ def test_update_loc_coords_float(self):
 
                 # Evaluate the output
                 self.assertAlmostEqual(
-                    float(self.vdata.lat), mag_out[coord][0], places=4,
+                    self.vdata.lat.item(), mag_out[coord][0], places=4,
                     msg="unexpected magnetic latitude")
                 self.assertAlmostEqual(
-                    float(self.vdata.lt), mag_out[coord][1], places=4,
+                    self.vdata.lt.item(), mag_out[coord][1], places=4,
                     msg="unexpected MLT")
                 self.assertRegex(self.vdata.loc_coord, "magnetic")
 
@@ -552,9 +552,9 @@ def test_update_loc_coords_float(self):
                                              coord=coord)
 
                 # Evaluate the output; note the loss of precision
-                self.assertAlmostEqual(float(self.vdata.lat), 75.0, places=1,
+                self.assertAlmostEqual(self.vdata.lat.item(), 75.0, places=1,
                                        msg="unexpected geographic latitude")
-                self.assertAlmostEqual(float(self.vdata.lt), 22.0, places=1,
+                self.assertAlmostEqual(self.vdata.lt.item(), 22.0, places=1,
                                        msg="unexpected SLT")
                 self.assertRegex(self.vdata.loc_coord, coord)
         return
@@ -722,15 +722,15 @@ def test_update_vect_and_loc_coords_float(self):
 
                     # Evaluate the output
                     self.assertAlmostEqual(
-                        float(self.vdata.vect_n),
+                        self.vdata.vect_n.item(),
                         mag_out[coord][loc_coord]['vect_n'],
                         places=4, msg="unexpected north component")
                     self.assertAlmostEqual(
-                        float(self.vdata.vect_e),
+                        self.vdata.vect_e.item(),
                         mag_out[coord][loc_coord]['vect_e'],
                         places=4, msg="unexpected east component")
                     self.assertAlmostEqual(
-                        float(self.vdata.vect_z),
+                        self.vdata.vect_z.item(),
                         mag_out[coord][loc_coord]['vect_z'],
                         places=4, msg="unexpected vertical component")
                     self.assertRegex(self.vdata.vect_coord, "magnetic")
@@ -738,10 +738,10 @@ def test_update_vect_and_loc_coords_float(self):
 
                     if loc_coord in mag_out.keys():
                         self.assertAlmostEqual(
-                            float(self.vdata.lat), mag_out[loc_coord]['lat'],
+                            self.vdata.lat.item(), mag_out[loc_coord]['lat'],
                             places=4, msg="unexpected magnetic latitude")
                         self.assertAlmostEqual(
-                            float(self.vdata.lt), mag_out[loc_coord]['lt'],
+                            self.vdata.lt.item(), mag_out[loc_coord]['lt'],
                             places=4, msg="unexpected MLT")
         return
 
diff --git a/ocbpy/vectors.py b/ocbpy/vectors.py
index d176eb79..83427c4d 100644
--- a/ocbpy/vectors.py
+++ b/ocbpy/vectors.py
@@ -98,10 +98,10 @@ def calc_vec_pole_angle(data_lt, data_lat, pole_lt, pole_lat):
     else:
         flat_mask = (((del_long == 0) | (abs(del_long) == np.pi))
                      & np.greater(abs(data_lat), abs(pole_lat),
-                                  where=~np.isnan(del_long)))
+                                  where=~np.isnan(del_long), out=None))
         zero_mask = (((del_long == 0) | (abs(del_long) == np.pi))
                      & np.less_equal(abs(data_lat), abs(pole_lat),
-                                     where=~np.isnan(del_long)))
+                                     where=~np.isnan(del_long), out=None))
 
         pole_angle[flat_mask] = 180.0
         pole_angle[zero_mask] = 0.0
@@ -176,39 +176,43 @@ def define_pole_quadrants(data_lt, pole_lt, pole_angle):
     pole_quad = np.zeros(shape=np.asarray(del_lt).shape)
 
     # Determine which differences need to be
-    neg_mask = np.less(del_lt, 0.0, where=~np.isnan(del_lt)) & ~np.isnan(del_lt)
+    neg_mask = np.less(del_lt, 0.0, where=~np.isnan(del_lt),
+                       out=None) & ~np.isnan(del_lt)
     while np.any(neg_mask):
         if len(del_lt.shape) == 0:
             del_lt += 24.0
-            neg_mask = np.less(del_lt, 0.0)  # Has one finite value
+            neg_mask = np.less(del_lt, 0.0, out=None)  # Has one finite value
         else:
             del_lt[neg_mask] += 24.0
-            neg_mask = np.less(del_lt, 0.0,
+            neg_mask = np.less(del_lt, 0.0, out=None,
                                where=~np.isnan(del_lt)) & ~np.isnan(del_lt)
 
-    large_mask = np.greater_equal(abs(del_lt), 24.0,
+    large_mask = np.greater_equal(abs(del_lt), 24.0, out=None,
                                   where=~np.isnan(del_lt)) & ~np.isnan(del_lt)
     while np.any(large_mask):
         if len(del_lt.shape) == 0:
             del_lt -= 24.0 * np.sign(del_lt)
-            large_mask = np.greater_equal(abs(del_lt), 24.0)  # One finite value
+            # One finite value
+            large_mask = np.greater_equal(abs(del_lt), 24.0, out=None)
         else:
             del_lt[large_mask] -= 24.0 * np.sign(del_lt[large_mask])
-            large_mask = np.greater_equal(abs(del_lt), 24.0,
+            large_mask = np.greater_equal(abs(del_lt), 24.0, out=None,
                                           where=~np.isnan(del_lt)) & ~np.isnan(
                                               del_lt)
 
     # Find the quadrant in which the OCB pole lies
     nan_mask = ~np.isnan(pole_angle) & ~np.isnan(del_lt)
-    quad1_mask = np.less(pole_angle, 90.0, where=nan_mask) & np.less(
-        del_lt, 12.0, where=nan_mask) & nan_mask
-    quad2_mask = np.less(pole_angle, 90.0, where=nan_mask) & np.greater_equal(
-        del_lt, 12.0, where=nan_mask) & nan_mask
+    quad1_mask = np.less(pole_angle, 90.0, where=nan_mask, out=None) & np.less(
+        del_lt, 12.0, where=nan_mask, out=None) & nan_mask
+    quad2_mask = np.less(
+        pole_angle, 90.0, where=nan_mask, out=None) & np.greater_equal(
+            del_lt, 12.0, where=nan_mask, out=None) & nan_mask
     quad3_mask = np.greater_equal(
-        pole_angle, 90.0, where=nan_mask) & np.greater_equal(
-            del_lt, 12.0, where=nan_mask) & nan_mask
-    quad4_mask = np.greater_equal(pole_angle, 90.0, where=nan_mask) & np.less(
-        del_lt, 12.0, where=nan_mask) & nan_mask
+        pole_angle, 90.0, where=nan_mask, out=None) & np.greater_equal(
+            del_lt, 12.0, where=nan_mask, out=None) & nan_mask
+    quad4_mask = np.greater_equal(
+        pole_angle, 90.0, where=nan_mask, out=None) & np.less(
+            del_lt, 12.0, where=nan_mask, out=None) & nan_mask
 
     if len(pole_quad.shape) == 0:
         if np.all(quad1_mask):
@@ -255,14 +259,16 @@ def define_vect_quadrants(vect_n, vect_e):
     # Get the masks for non-fill values in each quadrant
     nan_mask = ~np.isnan(vect_n) & ~np.isnan(vect_e)
     quad1_mask = np.greater_equal(
-        vect_n, 0.0, where=nan_mask) & np.greater_equal(
-            vect_e, 0.0, where=nan_mask) & nan_mask
-    quad2_mask = np.greater_equal(vect_n, 0.0, where=nan_mask) & np.less(
-        vect_e, 0.0, where=nan_mask) & nan_mask
-    quad3_mask = np.less(vect_n, 0.0, where=nan_mask) & np.less(
-        vect_e, 0.0, where=nan_mask) & nan_mask
-    quad4_mask = np.less(vect_n, 0.0, where=nan_mask) & np.greater_equal(
-        vect_e, 0.0, where=nan_mask) & nan_mask
+        vect_n, 0.0, where=nan_mask, out=None) & np.greater_equal(
+            vect_e, 0.0, where=nan_mask, out=None) & nan_mask
+    quad2_mask = np.greater_equal(
+        vect_n, 0.0, where=nan_mask, out=None) & np.less(
+            vect_e, 0.0, where=nan_mask, out=None) & nan_mask
+    quad3_mask = np.less(vect_n, 0.0, where=nan_mask, out=None) & np.less(
+        vect_e, 0.0, where=nan_mask, out=None) & nan_mask
+    quad4_mask = np.less(
+        vect_n, 0.0, where=nan_mask, out=None) & np.greater_equal(
+            vect_e, 0.0, where=nan_mask, out=None) & nan_mask
 
     # Initialize the output
     vect_quad = np.zeros(shape=nan_mask.shape)
@@ -466,13 +472,14 @@ def calc_dest_vec_sign(pole_quad, vect_quad, base_naz_angle, pole_angle,
 
         pmask = (quads[1][1] | quads[2][2] | quads[3][3] | quads[4][4]
                  | ((quads[1][4] | quads[2][3]) & np.less_equal(
-                     base_naz_angle, pole_plus, where=nan_mask))
+                     base_naz_angle, pole_plus, where=nan_mask, out=None))
                  | ((quads[1][2] | quads[2][1]) & np.less_equal(
-                     base_naz_angle, minus_pole, where=nan_mask))
+                     base_naz_angle, minus_pole, where=nan_mask, out=None))
                  | ((quads[3][4] | quads[4][3]) & np.greater_equal(
-                     base_naz_angle, 180.0 - pole_minus, where=nan_mask))
+                     base_naz_angle, 180.0 - pole_minus, where=nan_mask,
+                     out=None))
                  | ((quads[3][2] | quads[4][1]) & np.greater_equal(
-                     base_naz_angle, pole_minus, where=nan_mask)))
+                     base_naz_angle, pole_minus, where=nan_mask, out=None)))
 
         if np.any(pmask):
             if len(vsigns["north"].shape) == 0:
@@ -491,13 +498,13 @@ def calc_dest_vec_sign(pole_quad, vect_quad, base_naz_angle, pole_angle,
 
         pmask = (quads[1][4] | quads[2][1] | quads[3][2] | quads[4][3]
                  | ((quads[1][1] | quads[4][4]) & np.greater_equal(
-                     base_naz_angle, pole_angle, where=nan_mask))
+                     base_naz_angle, pole_angle, where=nan_mask, out=None))
                  | ((quads[3][1] | quads[2][4]) & np.less_equal(
-                     base_naz_angle, minus_pole, where=nan_mask))
+                     base_naz_angle, minus_pole, where=nan_mask, out=None))
                  | ((quads[4][2] | quads[1][3]) & np.greater_equal(
-                     base_naz_angle, minus_pole, where=nan_mask))
+                     base_naz_angle, minus_pole, where=nan_mask, out=None))
                  | ((quads[2][2] | quads[3][3]) & np.less_equal(
-                     base_naz_angle, pole_angle, where=nan_mask)))
+                     base_naz_angle, pole_angle, where=nan_mask, out=None)))
 
         if np.any(pmask):
             if len(vsigns["east"].shape) == 0:
@@ -618,7 +625,8 @@ def adjust_vector(vect_lt, vect_lat, vect_n, vect_e, vect_z, vect_quad,
         # Determine if the measurement is on or between the poles. This does
         # not affect the vertical direction
         sign_mask = (pole_angle == 0.0) & np.greater_equal(
-            vect_lat, pole_lat, where=~np.isnan(vect_lat)) & ~np.isnan(vect_lat)
+            vect_lat, pole_lat, where=~np.isnan(vect_lat),
+            out=None) & ~np.isnan(vect_lat)
 
         if np.any(sign_mask):
             if len(out_shape) == 0:

From 861494074d5058776f9cfe1b49a733b5be92883f Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 1 Apr 2026 16:52:34 -0400
Subject: [PATCH 18/43] DOC: added numpy changes to changelog

Added the numpy changes to the changelog.
---
 Changelog.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Changelog.rst b/Changelog.rst
index a066ae82..e01a1845 100644
--- a/Changelog.rst
+++ b/Changelog.rst
@@ -9,6 +9,7 @@ Summary of all changes made since the first stable release
 * DEP: Removed support for older versions of `zenodo_get`
 * MAINT: Added support for Python 3.14
 * MAINT: Updated GitHub CI actions
+* MAINT: Updated numpy usage
 
 0.6.0 (07-07-2025)
 ------------------

From fc6d4d1789979a8f149de2e2f57a5aae47821c7f Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 1 Apr 2026 17:03:12 -0400
Subject: [PATCH 19/43] MAINT: fixed broken links

Zenodo has disabled checks by sphinx-build, so remove the active linking in the docs.
---
 docs/citing.rst                    |  2 +-
 docs/examples/ex_dmsp.rst          | 13 ++++++-------
 ocbpy/boundaries/dmsp_ssj_files.py |  2 +-
 3 files changed, 8 insertions(+), 9 deletions(-)

diff --git a/docs/citing.rst b/docs/citing.rst
index 1e4fd674..5367259d 100644
--- a/docs/citing.rst
+++ b/docs/citing.rst
@@ -98,7 +98,7 @@ boundary method and data set are provided below.
 * **SSJ Auroral Boundaries (2010-2014)**: Kilcommons, L., et al. (2019).
   Defense Meteorology Satellite Program (DMSP) Electron Precipitation (SSJ)
   Auroral Boundaries, 2010-2014 (Version 1.0.0) [Data set]. Zenodo.
-  http://doi.org/10.5281/zenodo.3373812
+  doi:10.5281/zenodo.3373812
 
 
 .. _cite-starkov:
diff --git a/docs/examples/ex_dmsp.rst b/docs/examples/ex_dmsp.rst
index 069037ad..8ad12287 100644
--- a/docs/examples/ex_dmsp.rst
+++ b/docs/examples/ex_dmsp.rst
@@ -13,13 +13,12 @@ need the
 `ssj_auroral_boundary <https://github.com/lkilcommons/ssj_auroral_boundary>`__
 package, but now we preferentially support using the
 `zendodo_get <https://github.com/dvolgyes/zenodo_get>`__ package to obtain the
-boundary files from their
-`archive <https://zenodo.org/record/3373812>`__. Once installed,
-you can download DMSP SSJ data and obtain a boundary file for a specified time
-period (or all available times) using :py:mod:`ocbpy.boundaries.dmsp_ssj_files`.
-For this example, we'll use a single day. You can download the files into any
-directory, but this example will put them in the same directory as the other
-boundary files.
+boundary files from their archive at the website: zenodo.org/record/3373812.
+Once installed, you can download DMSP SSJ data and obtain a boundary file for a
+specified time period (or all available times) using
+:py:mod:`ocbpy.boundaries.dmsp_ssj_files`. For this example, we'll use a single
+day. You can download the files into any directory, but this example will put
+them in the same directory as the other boundary files.
 
 ::
    
diff --git a/ocbpy/boundaries/dmsp_ssj_files.py b/ocbpy/boundaries/dmsp_ssj_files.py
index b171f4e0..e556fc2a 100644
--- a/ocbpy/boundaries/dmsp_ssj_files.py
+++ b/ocbpy/boundaries/dmsp_ssj_files.py
@@ -20,7 +20,7 @@
 
 .. [7] Kilcommons, L., Redmon, R., & Knipp, D. (2019). Defense Meteorology
    Satellite Program (DMSP) Electron Precipitation (SSJ) Auroral Boundaries,
-   2010-2014 (1.0.0) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.3373812
+   2010-2014 (1.0.0) [Data set]. Zenodo. doi:10.5281/zenodo.3373812
 
 """
 

From 1d585865df89868ce72711a7d2b33cf086ba407a Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 1 Apr 2026 17:31:56 -0400
Subject: [PATCH 20/43] MAINT: update pandas access

Update the pandas index access to use `iloc`.
---
 ocbpy/instruments/pysat_instruments.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/ocbpy/instruments/pysat_instruments.py b/ocbpy/instruments/pysat_instruments.py
index 1c139f88..fdb8dd12 100644
--- a/ocbpy/instruments/pysat_instruments.py
+++ b/ocbpy/instruments/pysat_instruments.py
@@ -494,7 +494,7 @@ def add_ocb_to_data(pysat_inst, mlat_name='', mlt_name='', height_name='',
                         vname = vector_names[eattr][ikey]
                         if vname in pysat_inst.variables:
                             if time_mask is None:
-                                vector_init[ikey] = pysat_inst[vname][iout]
+                                vector_init[ikey] = pysat_inst[vname].iloc[iout]
                             else:
                                 vector_init[ikey] = reshape_pad_mask_flatten(
                                     pysat_inst[vname], time_mask)
@@ -529,7 +529,7 @@ def add_ocb_to_data(pysat_inst, mlat_name='', mlt_name='', height_name='',
                 for eattr in scale_names:
                     oattr = "{:s}_{:s}".format(eattr, bname)
                     if time_mask is None:
-                        evar = pysat_inst[eattr][iout]
+                        evar = pysat_inst[eattr].iloc[iout]
                     else:
                         evar = reshape_pad_mask_flatten(pysat_inst[eattr],
                                                         time_mask)

From fa97ed825a0fac713331ddcd3d25e9f4ed48e371 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 1 Apr 2026 17:43:04 -0400
Subject: [PATCH 21/43] Revert "MAINT: update pandas access"

This reverts commit 1d585865df89868ce72711a7d2b33cf086ba407a.
---
 ocbpy/instruments/pysat_instruments.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/ocbpy/instruments/pysat_instruments.py b/ocbpy/instruments/pysat_instruments.py
index fdb8dd12..1c139f88 100644
--- a/ocbpy/instruments/pysat_instruments.py
+++ b/ocbpy/instruments/pysat_instruments.py
@@ -494,7 +494,7 @@ def add_ocb_to_data(pysat_inst, mlat_name='', mlt_name='', height_name='',
                         vname = vector_names[eattr][ikey]
                         if vname in pysat_inst.variables:
                             if time_mask is None:
-                                vector_init[ikey] = pysat_inst[vname].iloc[iout]
+                                vector_init[ikey] = pysat_inst[vname][iout]
                             else:
                                 vector_init[ikey] = reshape_pad_mask_flatten(
                                     pysat_inst[vname], time_mask)
@@ -529,7 +529,7 @@ def add_ocb_to_data(pysat_inst, mlat_name='', mlt_name='', height_name='',
                 for eattr in scale_names:
                     oattr = "{:s}_{:s}".format(eattr, bname)
                     if time_mask is None:
-                        evar = pysat_inst[eattr].iloc[iout]
+                        evar = pysat_inst[eattr][iout]
                     else:
                         evar = reshape_pad_mask_flatten(pysat_inst[eattr],
                                                         time_mask)

From 7995ff9329de3839113b929c2d7943bf84e826a5 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Thu, 2 Apr 2026 16:11:58 -0400
Subject: [PATCH 22/43] BUG: fixed program name

Updated the program name in the pip workflow to be correct.
---
 .github/workflows/pip_rc_install.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/pip_rc_install.yml b/.github/workflows/pip_rc_install.yml
index d54da2ba..70aa2e53 100644
--- a/.github/workflows/pip_rc_install.yml
+++ b/.github/workflows/pip_rc_install.yml
@@ -27,7 +27,7 @@ jobs:
     - name: Install standard dependencies
       run: pip install -r requirements.txt
 
-    - name: Install pysat RC
+    - name: Install ocbpy RC
       run: pip install --no-deps --pre -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple/ ocbpy
 
     - name: Check that installation imports correctly

From d1732d36b891f31363ea080c74deee5ed1aee923 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Thu, 2 Apr 2026 16:12:40 -0400
Subject: [PATCH 23/43] BUG: test code with AACGMV2 release candidate

Test this code with the AACGMV2 release candidate.
---
 .github/workflows/main.yml | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 1a07e83d..03651881 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -38,6 +38,12 @@ jobs:
         mkdir pysatData
         python -c "import pysat; pysat.params['data_dirs'] = 'pysatData'"
 
+    - name: Install AACGMV2 RC
+      run: |
+        pip uninstall aacgmv2
+        pip install --no-deps --pre -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple/ aacgmv2
+        python -c "import aacgmv2; print(aacgmv2.__version__)"
+
     - name: Test PEP8 compliance
       run: flake8 . --count --show-source --statistics
 

From 2422f151f8825d3f0eaa110861ef9261ff6702ac Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Thu, 2 Apr 2026 16:24:12 -0400
Subject: [PATCH 24/43] BUG: remove uninstall line

Removed the uninstall line, just install the newer version.
---
 .github/workflows/main.yml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 03651881..f1b50d15 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -40,7 +40,6 @@ jobs:
 
     - name: Install AACGMV2 RC
       run: |
-        pip uninstall aacgmv2
         pip install --no-deps --pre -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple/ aacgmv2
         python -c "import aacgmv2; print(aacgmv2.__version__)"
 

From b36cb8d1b1683933e653ffd0b0969dd1d2fed87c Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Thu, 2 Apr 2026 16:28:47 -0400
Subject: [PATCH 25/43] BUG: fix aacgmv2 rc install

Specify the version number.
---
 .github/workflows/main.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index f1b50d15..310938e2 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -40,7 +40,7 @@ jobs:
 
     - name: Install AACGMV2 RC
       run: |
-        pip install --no-deps --pre -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple/ aacgmv2
+        pip install pip install -i https://test.pypi.org/simple/ aacgmv2==2.7.1rc1
         python -c "import aacgmv2; print(aacgmv2.__version__)"
 
     - name: Test PEP8 compliance

From dabaa5f19338ba55281dbf294adfc0a4fb4280b5 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Thu, 2 Apr 2026 16:33:59 -0400
Subject: [PATCH 26/43] BUG: fixed typo

Fixed a typo in the TestPyPi install line.
---
 .github/workflows/main.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 310938e2..465fac83 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -40,7 +40,7 @@ jobs:
 
     - name: Install AACGMV2 RC
       run: |
-        pip install pip install -i https://test.pypi.org/simple/ aacgmv2==2.7.1rc1
+        pip install --no-deps -i https://test.pypi.org/simple/ aacgmv2==2.7.1rc1
         python -c "import aacgmv2; print(aacgmv2.__version__)"
 
     - name: Test PEP8 compliance

From 82e9029a6689a446155f0ff547aeff3263d622a3 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Thu, 2 Apr 2026 16:39:35 -0400
Subject: [PATCH 27/43] BUG: re-added extra flags

Hopefully the install will go better and it won't try to install the dependencies.
---
 .github/workflows/main.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 465fac83..cbb65c50 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -40,7 +40,7 @@ jobs:
 
     - name: Install AACGMV2 RC
       run: |
-        pip install --no-deps -i https://test.pypi.org/simple/ aacgmv2==2.7.1rc1
+        pip install --no-deps --pre -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple/ aacgmv2==2.7.1rc1
         python -c "import aacgmv2; print(aacgmv2.__version__)"
 
     - name: Test PEP8 compliance

From 37f99de34888dff7b75132944b80150f22f6fe05 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Tue, 7 Apr 2026 16:39:28 -0400
Subject: [PATCH 28/43] TST: remove the aacgmv2 test installation

Now that the new version of AACGMV2 has been released, remove the test installation case.
---
 .github/workflows/main.yml | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index cbb65c50..1a07e83d 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -38,11 +38,6 @@ jobs:
         mkdir pysatData
         python -c "import pysat; pysat.params['data_dirs'] = 'pysatData'"
 
-    - name: Install AACGMV2 RC
-      run: |
-        pip install --no-deps --pre -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple/ aacgmv2==2.7.1rc1
-        python -c "import aacgmv2; print(aacgmv2.__version__)"
-
     - name: Test PEP8 compliance
       run: flake8 . --count --show-source --statistics
 

From 74c6533dd0cdc0f10b02deb4d2500680775480ce Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Tue, 7 Apr 2026 16:39:51 -0400
Subject: [PATCH 29/43] MAINT: set lower limit for aacgmv2

Set a lower limit for AACGMV2, set to the most recent version.
---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index d9ed3035..82e41de0 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -12,7 +12,7 @@ maintainers = [
 ]
 requires-python = ">=3.10"
 dependencies = [
-    "aacgmv2",
+    "aacgmv2>=2.7.1",
     "numpy",
 ]
 readme = "README.md"

From bf06e8e79713f6d15467fe53122867caa0642481 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Mon, 11 May 2026 18:02:38 -0400
Subject: [PATCH 30/43] ENH: added the CHAMP auroral boundary model

Added the CHAMP auroral boundary model.
---
 ocbpy/boundaries/models.py | 221 ++++++++++++++++++++++++++++++++++++-
 1 file changed, 220 insertions(+), 1 deletion(-)

diff --git a/ocbpy/boundaries/models.py b/ocbpy/boundaries/models.py
index fafc9457..0dec4e70 100644
--- a/ocbpy/boundaries/models.py
+++ b/ocbpy/boundaries/models.py
@@ -13,6 +13,9 @@
    Auroral Boundary, J. Geophys. Res, 88(A7), 5692-5708.
 .. [10] Umbach and Jones (2003) A Few Methods for Fitting Circles to Data,
    IEEE Transactions on Instruments and Measurements, 52(6), pp 1881-1885
+.. [11] Xiong and Luhr (2014) An empirical model of the auroral oval derived
+   from CHAMP field-aligned current signatures - Part 2, Ann. Geophys., 32,
+   pp 623-631, doi:10.5194/angeo-32-623-2014
 
 """
 
@@ -22,7 +25,7 @@
 
 
 def starkov_auroral_boundary(mlt, al=-1, bnd='ocb'):
-    """Calculate the location of the auroral boundaries.
+    """Calculate the location of the Starkov auroral boundaries.
 
     Parameters
     ----------
@@ -358,3 +361,219 @@ def circle_fit(mlt, colat):
     r_err = np.sqrt(np.mean((rad_bound - radius)**2))
 
     return phi_cent, r_cent, radius, r_err
+
+
+def ch_aurora_2014_boundary(mlt, em=0, bnd='ocb', hemi=1, obs_colat=None,
+                            obs_mlt=None):
+    """Calculate the location of the CH-Aurora-2014 auroral boundaries.
+
+    Parameters
+    ----------
+    mlt : float or array-like
+        Magnetic local time in hours
+    em : float or int
+        The time-integrated Newell Coupling Function, as described in Equation 2
+        of [11]_ (default=0)
+    bnd : str
+        Boundary to calculate, expects one of 'ocb' or 'eab' (default='ocb')
+    hemi : int
+        Sign denoting the hemisphere, 1 for North and -1 for South (default=1)
+    obs_colat : array-like
+        CHAMP, or other, boundary observation co-latitudes (default=None)
+    obs_mlt : array-like
+        MLT of the observed boundary locations (default=None)
+
+    Returns
+    -------
+    bnd_lat : float or array-like
+        Location of the boundary in degrees away from the pole in apex
+        geomagnetic coordinates for the specified magnetic local times.
+
+    References
+    ----------
+    [11]_
+
+    """
+    # Calculate each parameter for the desired IMF conditions; all are in
+    # degrees except for phi0, which is in radians.
+    semix, semiy, x0, y0, phi0 = ch_aurora_2014_coefficient_values(em, bnd,
+                                                                   hemi)
+
+    # If observations are supplied, ensure their MLT are included in the calc
+    calc_mlt = np.array(mlt)
+    iobs = 0
+    if obs_mlt is not None:
+        obs_mlt_array = np.array(obs_mlt)
+
+        if calc_mlt.ndim == 0:
+            iobs = 1
+            if obs_mlt_array.ndim == 0:
+                calc_mlt = np.array([mlt, obs_mlt])
+            else:
+                calc_mlt = list(obs_mlt)
+                calc_mlt.insert(0, mlt)
+                calc_mlt = np.array(calc_mlt)
+        else:
+            calc_mlt = list(mlt)
+            iobs = len(mlt)
+            if obs_mlt_array.ndim == 0:
+                calc_mlt.append(obs_mlt)
+            else:
+                calc_mlt.extend(list(obs_mlt))
+            calc_mlt = np.array(calc_mlt)
+
+    # Calculate the angular MLT
+    ang_lt = ocb_time.hr2rad(calc_mlt)
+
+    # First round of calculations with phi == ang_lt
+    bnd_lat = ch_aurora_2014_radius(ang_lt, semix, semiy, x0, y0, phi0)
+
+    # Now calculate the angular time from the radius
+    xprime = bnd_lat * np.cos(ang_lt + phi0) + x0
+    yprime = bnd_lat * np.sin(ang_lt + phi0) + y0
+    ang_prime = np.arctan2(yprime, xprime)
+
+    # Recalculate the radius using the new local times
+    bnd_lat = ch_aurora_2014_radius(ang_prime, semix, semiy, x0, y0, phi0)
+
+    # If desired, modify the boundary locations based on observations
+    if iobs > 0:
+        # Get the model values for the observation times
+        mod_colat = bnd_lat[iobs:]
+        del_lat = np.nanmean(obs_colat - mod_colat)
+
+        # Insert the adjustment
+        bnd_lat = ch_aurora_2014_radius(ang_lt, semix, semiy, x0, y0, phi0,
+                                        del_lat)
+
+        # Now calculate the angular time from the radius
+        xprime = bnd_lat * np.cos(ang_lt + phi0) + x0
+        yprime = bnd_lat * np.sin(ang_lt + phi0) + y0
+        ang_prime = np.arctan2(yprime, xprime)
+
+        # Recalculate the radius using the new local times
+        bnd_lat = ch_aurora_2014_radius(ang_prime, semix, semiy, x0, y0, phi0,
+                                        del_lat)
+
+        # Downselect the output to include only the desired MLT
+        bnd_lat = bnd_lat[:iobs]
+
+    # Ensure the output co-latitude is a float if the input MLT is a float
+    if calc_mlt.ndim == 0:
+        if bnd_lat.ndim == 0:
+            bnd_lat = float(bnd_lat)
+        else:
+            bnd_lat = float(bnd_lat[0])
+
+    return bnd_lat
+
+
+def ch_aurora_2014_coefficient_values(em, bnd, hemi):
+    """Retrive the CH-Aurora-2014 auroral boundary coefficients.
+
+    Parameters
+    ----------
+    em : float or int
+        The time-integrated Newell Coupling Function, as described in Equation 2
+        of [11]_ (default=0)
+    bnd : str
+        Boundary to calculate, expects one of 'ocb' or 'eab'
+    hemi : int
+        Sign denoting the hemisphere, 1 for North and -1 for South
+
+    Returns
+    -------
+    semix : float
+        Semi-axis along the x-plane
+    semiy : float
+        Semi-axis along the y-plane
+    x0 : float
+        x-coordinate of the ellipse center
+    y0 : float
+        y-coordinate of the ellipse center
+    phi0 : float
+        Orientation angle of the ellipse
+
+    References
+    ----------
+    Tables 2 and 3 in [11]_
+
+    """
+    # Set the ellipse parameters from Tables 2 and 3 from Xiong and Luhr.
+    # The list includes the p0, p1, and p2 parameters in that order, all
+    # provided in degrees
+    semix_param = {'eab': {1: [18.861, 0.95470, -0.023836],
+                           -1: [18.559, 0.81597, -0.013209]},
+                   'ocb': {1: [12.813, 0.26173, -5.9729e-4],
+                           -1: [13.251, 0.28870, -3.0559e-4]}}
+    semiy_param = {'eab': {1: [20.562, 1.1504, -0.029566],
+                           -1: [19.549, 0.10752, -0.024605]},
+                   'ocb': {1: [9.5486, 0.96759, -0.029556],
+                           -1: [11.605, 0.82006, -0.024073]}}
+    x0_param = {'eab': {1: [4.1263, 0.027827, 0.0],
+                        -1: [3.6946, 0.044667, 0.0]},
+                'ocb': {1: [4.5175, -0.025310, 0.0],
+                        -1: [4.2526, -0.021674, 0.0]}}
+    y0_param = {'eab': {1: [-0.32637, -0.028855, 1.6569e-3],
+                        -1: [-0.60436, -0.011623, 6.5985e-4]},
+                'ocb': {1: [-0.39316, -0.15732, 5.613e-3],
+                        -1: [-1.1330, -0.024479, 7.0729e-4]}}
+    phi0_param = {'eab': {1: [-3.1555, 0.31147, 0.0],
+                          -1: [-8.8836, -0.10934, 0.0]},
+                  'ocb': {1: [-8.5358, 1.2831, 0.0],
+                          -1: [3.7050, -1.5508, 0.0]}}
+
+    # Calculate each parameter for the desired IMF conditions
+    semix = semix_param[bnd][hemi][0] + semix_param[bnd][hemi][
+        1] * em + semix_param[bnd][hemi][2] * em * em
+    semiy = semiy_param[bnd][hemi][0] + semiy_param[bnd][hemi][
+        1] * em + semiy_param[bnd][hemi][2] * em * em
+    x0 = x0_param[bnd][hemi][0] + x0_param[bnd][hemi][1] * em + x0_param[bnd][
+        hemi][2] * em * em
+    y0 = y0_param[bnd][hemi][0] + y0_param[bnd][hemi][1] * em + y0_param[bnd][
+        hemi][2] * em * em
+    phi0 = np.radians(phi0_param[bnd][hemi][0] + phi0_param[bnd][hemi][1] * em
+                      + phi0_param[bnd][hemi][2] * em * em)
+
+    return semix, semiy, x0, y0, phi0
+
+
+def ch_aurora_2014_radius(ang_lt, semix, semiy, x0, y0, phi0, del_rad=0.0):
+    """Calculate the CH-Aurora-2014 ellipse radii at the desired times.
+
+    Parameters
+    ----------
+    ang_lt : float or array-like
+        Angular measure of MLT in radians
+    semix : float
+        Semi-axis along the x-plane in degrees
+    semiy : float
+        Semi-axis along the y-plane in degrees
+    x0 : float
+        x-coordinate of the ellipse center in degrees
+    y0 : float
+        y-coordinate of the ellipse center in degrees
+    phi0 : float
+        Orientation angle of the ellipse in radians
+    del_rad : float
+        Assimilative adjustment to the radius in degrees
+
+    Returns
+    -------
+    rad : float or array-like
+        Radius of the ellipse at the desired MLT in degrees
+
+    References
+    ----------
+    Equations 3 and 4 in [11]_
+
+    """
+    # Calculate the radius from the centre of the ellipse
+    r0 = del_rad + semix * semiy / np.sqrt((semix * np.sin(ang_lt + phi0))**2
+                                           + (semiy * np.cos(ang_lt + phi0))**2)
+
+    # Adjust to provide the radius from the magnetic pole
+    rad = np.sqrt((r0 * np.cos(ang_lt + phi0) + x0)**2
+                  + (r0 * np.sin(ang_lt + phi0) + y0)**2)
+
+    return rad

From dea5016ff9e1f23fda651882d559988fca56184c Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Mon, 11 May 2026 18:02:58 -0400
Subject: [PATCH 31/43] TST: added CHAMP unit tests

Added unit tests for the CHAMP auroral boundary model.
---
 ocbpy/tests/test_models.py | 148 +++++++++++++++++++++++++++++++++++++
 1 file changed, 148 insertions(+)

diff --git a/ocbpy/tests/test_models.py b/ocbpy/tests/test_models.py
index a34b9e6f..a0433a3e 100644
--- a/ocbpy/tests/test_models.py
+++ b/ocbpy/tests/test_models.py
@@ -293,3 +293,151 @@ def test_circle_fit(self):
         self.assertAlmostEqual(radius, 1.0)
         self.assertAlmostEqual(r_err, 0.0)
         return
+
+
+class TestCHAMPModel(unittest.TestCase):
+    """"Unit tests for the CH-Aurora-2014 routines."""
+
+    def setUp(self):
+        """Initialize the test case by setting some values to test against."""
+        self.mlt = np.arange(0, 24, 1)
+        self.coeff_out = {'semix': {'ocb': {1: 12.813, -1: 13.251},
+                                    'eab': {1: 18.861, -1: 18.559}},
+                          'semiy': {'ocb': {1: 9.5486, -1: 11.605},
+                                    'eab': {1: 20.562, -1: 19.549}},
+                          'x0': {'ocb': {1: 4.5175, -1: 4.2526},
+                                 'eab': {1: 4.1263, -1: 3.6946}},
+                          'y0': {'ocb': {1: -0.39316, -1: -1.1330},
+                                 'eab': {1: -0.32637, -1: -0.60436}},
+                          'phi0': {'ocb': {1: -0.1489778, -1: 0.0646644},
+                                   'eab': {1: -0.055074, -1: -0.155048}}}
+        self.em = [0, 10.5]
+        self.iobs = [0, 12]
+        self.obs_mlt = self.mlt[self.iobs]
+        self.obs_colat = {'ocb': np.array([18.0, 8.0]),
+                          'eab': np.array([25.0, 16.0])}
+        self.max_lat = {'ocb': {1: [17.31668454, 20.1346482],
+                                -1: [17.5881323, 20.741303775]},
+                        'eab': {1: [23.052935, 31.313414],
+                                -1: [22.354591, 29.69814346]}}
+        return
+
+    def tearDown(self):
+        """Clean up the test environment."""
+        del self.mlt, self.coeff_out, self.em, self.max_lat, self.obs_mlt
+        del self.obs_colat, self.iobs
+        return
+
+    def test_coeff_construction(self):
+        """Test coefficient calculation for an Em of 0."""
+
+        for coeff in self.coeff_out.keys():
+            for bnd in self.coeff_out[coeff].keys():
+                for hemi in [1, -1]:
+                    with self.subTest(coeff=coeff, bnd=bnd, hemi=hemi):
+                        # Calculate the coefficient value
+                        out = models.ch_aurora_2014_coefficient_values(
+                            self.em[0], bnd, hemi)
+
+                        # Compare the output
+                        for ic, coeff in enumerate(['semix', 'semiy', 'x0',
+                                                    'y0']):
+                            self.assertEqual(
+                                out[ic], self.coeff_out[coeff][bnd][hemi],
+                                msg="{:s} does not match".format(coeff))
+
+                        # Phi0 has been converted to radians, so equality
+                        # will not be exact. Use significance from paper
+                        self.assertAlmostEqual(
+                            out[-1], self.coeff_out['phi0'][bnd][hemi],
+                            places=5, msg="phi0 does not match")
+        return
+
+    def test_coeff_bad_hemi(self):
+        """Test a KeyError is raised for an unknown hemisphere."""
+        hemi = "north"
+        with self.assertRaisesRegex(KeyError, hemi):
+            models.ch_aurora_2014_coefficient_values(
+                self.em[0], list(self.max_lat.keys())[0], hemi)
+        return
+
+    def test_coeff_bad_bnd(self):
+        """Test a KeyError is raised for an unknown boundary name."""
+        bound = "not a boundary"
+        with self.assertRaisesRegex(KeyError, bound):
+            models.ch_aurora_2014_coefficient_values(self.em[0], bound, 1)
+        return
+
+    def test_bound_loc_array(self):
+        """Test the expected boundary location across an MLT array."""
+        # Get the boundary keys
+        bnds = list(self.max_lat.keys())
+
+        # Cycle through low and high Em values
+        for ie, in_em in enumerate(self.em):
+            for hemi in self.max_lat[bnds[0]].keys():
+                with self.subTest(em=in_em, hemi=hemi):
+                    lats = {
+                        bnd: models.ch_aurora_2014_boundary(
+                            self.mlt, em=in_em, bnd=bnd, hemi=hemi)
+                        for bnd in bnds}
+
+                    # Test the output latitude shape and values
+                    for bnd in bnds:
+                        self.assertTupleEqual(self.mlt.shape, lats[bnd].shape)
+                        self.assertAlmostEqual(max(lats[bnd]),
+                                               self.max_lat[bnd][hemi][ie],
+                                               places=5)
+
+                    # Test that the OCB is greater than zero and the EAB is
+                    # greater than the OCB
+                    self.assertGreaterEqual(min(lats['ocb']), 0)
+                    self.assertTrue(np.all(lats['eab'] > lats['ocb']))
+
+        return
+
+    def test_bound_loc_float(self):
+        """Test the expected boundary location across an MLT value."""
+        # Cycle through low and high Em values
+        for ie, in_em in enumerate(self.em):
+            for bnd in self.max_lat.keys():
+                for hemi in self.max_lat[bnd].keys():
+                    with self.subTest(em=in_em, bnd=bnd, hemi=hemi):
+                        lat = models.ch_aurora_2014_boundary(
+                            self.mlt[0], em=in_em, bnd=bnd, hemi=hemi)
+
+                        # Test the output latitude shape and values
+                        self.assertTrue(isinstance(lat, float))
+                        self.assertGreaterEqual(lat, 0)
+                        self.assertLessEqual(lat, self.max_lat[bnd][hemi][ie])
+
+        return
+
+    def test_bound_assim(self):
+        """Test the expected assimilated boundary location."""
+
+        # Cycle through low and high Em values
+        for ie, in_em in enumerate(self.em):
+            for bnd in self.max_lat.keys():
+                for hemi in self.max_lat[bnd].keys():
+                    with self.subTest(em=in_em, bnd=bnd, hemi=hemi):
+                        mlat = models.ch_aurora_2014_boundary(
+                            self.mlt, em=in_em, bnd=bnd, hemi=hemi)
+                        alat = models.ch_aurora_2014_boundary(
+                            self.mlt, em=in_em, bnd=bnd, hemi=hemi,
+                            obs_mlt=self.obs_mlt, obs_colat=self.obs_colat[bnd])
+
+                        # Test the output latitude shape
+                        self.assertTupleEqual(self.mlt.shape, alat.shape)
+
+                        # Model and assimilation should differ
+                        self.assertGreater(abs(alat - mlat).min(), 1.0e-2)
+
+                        # Assimilated output should be closer to the provided
+                        # data points than the original model in at least one
+                        # of the assimilated locations
+                        self.assertTrue(
+                            np.any(abs(alat[self.iobs] - self.obs_colat[bnd])
+                                   < abs(mlat[self.iobs]
+                                         - self.obs_colat[bnd])))
+        return

From a6548794daf26a66ad96b96e6a737f7610ed347f Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Mon, 11 May 2026 18:03:24 -0400
Subject: [PATCH 32/43] DOC: added CHAMP documentation

Added citations and a description of the CHAMP auroral boundary model.
---
 docs/citing.rst     | 20 ++++++++++++++++++++
 docs/ocb_models.rst | 13 +++++++++++++
 2 files changed, 33 insertions(+)

diff --git a/docs/citing.rst b/docs/citing.rst
index 5367259d..cb4db0b0 100644
--- a/docs/citing.rst
+++ b/docs/citing.rst
@@ -126,3 +126,23 @@ source when publishing your work.
 * Gussenhoven, M. S. et al. (1983) Systematics of the Equatorward Diffuse
   Auroral Boundary, J. Geophys. Res, 88(A7), 5692-5708.
 
+
+.. _cite-xiong:
+
+CH-Aurora-2014 Model
+--------------------
+
+The CH-Aurora-2014 auroral boundary model provides either a mathematical
+description of the poleward and equatorward auroral boundaries in each
+hemisphere or an assimilated location based on FAC observations of the auroral
+boundaries in that hemisphere.  The data that may be used to adjust the model
+are described in the first article below, while the model itself is described in
+the second article. If you use this model please cite both it and your Newell
+Couping Function data source when publishing your work.
+
+* Xiong, C., et al. (2014) Determining the boundaries of the auroral oval from
+  CHAMP field-aligned current signatures - Part 1, Ann. Geophys., 32,
+  pp 609-622, doi:10.5194/angeo-32-609-2014.
+* Xiong, C. and H. Luhr (2014) An empirical model of the auroral oval derived
+  from CHAMP field-aligned current signatures - Part 2, Ann. Geophys., 32,
+  pp 623-631, doi:10.5194/angeo-32-623-2014
diff --git a/docs/ocb_models.rst b/docs/ocb_models.rst
index e8e5be27..c89fa6fb 100644
--- a/docs/ocb_models.rst
+++ b/docs/ocb_models.rst
@@ -40,6 +40,19 @@ at the binned times ('binned'), at the nearest binned time ('closest'), or at
 the requested times using a circle fit to all of binned times ('circle').
 
 
+.. _bound-model-champ:
+
+CH-Aurora-2014
+--------------
+
+The CH-Aurora-2014 model (see :ref:`cite-xiong`) uses a mathematical
+formulation based on CHAMP field-aligned current (FAC) data and a delayed
+time-history of the Newell Coupling Function. They specify both auroral
+boundaries in each hemisphere.  Although this model may be driven entirely by
+the Newell Coupling Function, the authors recommend adjusting the model output
+with CHAMP measurements of the FAC boundaries.
+
+
 .. _bound-model-module:
 
 Boundary Models Module

From 1f0a0c71aa453880707cc09e4b30cf6cff0fcb41 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Mon, 11 May 2026 18:03:43 -0400
Subject: [PATCH 33/43] DOC: updated changelog

Added the inclusion of the new model to the changelog.
---
 Changelog.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Changelog.rst b/Changelog.rst
index e01a1845..14351ac6 100644
--- a/Changelog.rst
+++ b/Changelog.rst
@@ -6,6 +6,7 @@ Summary of all changes made since the first stable release
 0.7.0 (XX-XX-2025)
 ------------------
 * ENH: Added the Gussenhoven (1983) model for the EAB
+* ENH: Added the CH-Aurora-2014 model for the OCB and EAB
 * DEP: Removed support for older versions of `zenodo_get`
 * MAINT: Added support for Python 3.14
 * MAINT: Updated GitHub CI actions

From 4872b4750fb0eb689f704e3b151682d2b01c1993 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Tue, 12 May 2026 10:02:50 -0400
Subject: [PATCH 34/43] DOC: added AI policy

Added an AI policy rejecting the use of AI.
---
 CONTRIBUTING.rst | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst
index 5d04921d..7c3caa71 100644
--- a/CONTRIBUTING.rst
+++ b/CONTRIBUTING.rst
@@ -88,3 +88,17 @@ To run a subset of tests from the test directory for a specific environment::
 To run all the tests for a specific environment::
 
     python -m unittest discover
+
+Use of 'AI'
+-----------
+
+This project is human centered and developers do not use Artificial Inteligence
+(AI) tools such as large language models in its creation or maintenance.
+The purpose of this tool is to improve high-latitude science, and thus values
+interactions with scientists and scientific programmers. Code generated by a
+large language model or similar technology, such as Anthropic’s Claude,
+GitHub/Microsoft’s Copilot, OpenAI’s ChatGPT, Facebook/Meta’s Code Llama et al.,
+is not compliant with the covenants and representations of OCBpy’s Contributor’s
+Agreement, and is thus not acceptable as code for OCBpy.
+
+

From 4ab3a82772b160eb742206d238085f7626da5ad2 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Tue, 12 May 2026 10:03:06 -0400
Subject: [PATCH 35/43] DOC: updated pull request template

Added a check asserting AI is not used.
---
 .github/pull_request_template.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md
index 1537b560..c326d966 100644
--- a/.github/pull_request_template.md
+++ b/.github/pull_request_template.md
@@ -38,6 +38,7 @@ import ocbpy
 - [ ] Make sure you are merging into the ``develop`` (not ``main``) branch
 - [ ] My commits are formatted appropriately (following the SciPy/NumPy style) 
 - [ ] My code follows the style guidelines of this project
+- [ ] I assert that I have not used AI in the development of this pull request
 - [ ] I have performed a self-review of my own code
 - [ ] I have commented my code, particularly in hard-to-understand areas
 - [ ] I have made corresponding changes to the documentation

From 648e8982e153737b93f8045e8b93fe8ef20810d0 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Tue, 12 May 2026 10:03:20 -0400
Subject: [PATCH 36/43] DOC: updated changelog

Added the AI policy update to the changelog.
---
 Changelog.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Changelog.rst b/Changelog.rst
index 14351ac6..f81937a9 100644
--- a/Changelog.rst
+++ b/Changelog.rst
@@ -7,6 +7,7 @@ Summary of all changes made since the first stable release
 ------------------
 * ENH: Added the Gussenhoven (1983) model for the EAB
 * ENH: Added the CH-Aurora-2014 model for the OCB and EAB
+* ENH: Added an AI usage policy (no AI use allowed)
 * DEP: Removed support for older versions of `zenodo_get`
 * MAINT: Added support for Python 3.14
 * MAINT: Updated GitHub CI actions

From 62e6dc21d50f8170d151d8b739c9fb6555228622 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 13 May 2026 13:09:59 -0400
Subject: [PATCH 37/43] MAINT: updated `zenodo_get` download checks

Updated download checks to use checksum (available again).  Removed old output checks, which are now redundant, and instead use hashlib to perform the the md5 checksum evaluation.  This should be more resilient to future `zenodo_get` updates.
---
 Changelog.rst                      |  3 +-
 ocbpy/boundaries/dmsp_ssj_files.py | 62 ++++++++++++++----------------
 2 files changed, 31 insertions(+), 34 deletions(-)

diff --git a/Changelog.rst b/Changelog.rst
index f81937a9..f63db94a 100644
--- a/Changelog.rst
+++ b/Changelog.rst
@@ -8,7 +8,8 @@ Summary of all changes made since the first stable release
 * ENH: Added the Gussenhoven (1983) model for the EAB
 * ENH: Added the CH-Aurora-2014 model for the OCB and EAB
 * ENH: Added an AI usage policy (no AI use allowed)
-* DEP: Removed support for older versions of `zenodo_get`
+* DEP: Removed support for older versions of `zenodo_get`, updated download
+  checks to account for breaking changes made in the newer versions
 * MAINT: Added support for Python 3.14
 * MAINT: Updated GitHub CI actions
 * MAINT: Updated numpy usage
diff --git a/ocbpy/boundaries/dmsp_ssj_files.py b/ocbpy/boundaries/dmsp_ssj_files.py
index e556fc2a..52ca5299 100644
--- a/ocbpy/boundaries/dmsp_ssj_files.py
+++ b/ocbpy/boundaries/dmsp_ssj_files.py
@@ -25,10 +25,9 @@
 """
 
 import datetime as dt
-from io import StringIO
+import hashlib
 import numpy as np
 import os
-import sys
 import zipfile
 
 import aacgmv2
@@ -108,45 +107,42 @@ def fetch_ssj_boundary_files(stime=None, etime=None, out_dir=None,
     if not os.path.isdir(out_dir):
         raise ValueError("can't find the output directory")
 
-    # Download the zenodo archive, capturing the output
-    zenodo_io = StringIO()
-    sys.stdout = zenodo_io
-    sys.stderr = zenodo_io
-
+    # Download the zenodo archive
     zenodo_get.download(doi=doi, output_dir=out_dir)
-    zenodo_checksum = None
 
-    # Parse the output and retrieve files from the zip archive
-    sys.stdout = sys.__stdout__
-    sys.stderr = sys.__stderr__
-    zen_msg = zenodo_io.getvalue()
-    zen_split = zen_msg.split()
+    # Get the checksum file (does not download)
+    zenodo_get.download(doi=doi, output_dir=out_dir, md5=True)
+    zenodo_checksum = os.path.join(out_dir, "md5sums.txt")
 
-    if zen_msg.find('Checksum is correct') < 0 and zen_msg.find(
-            'already downloaded correctly') < 0:
-        raise IOError('Bad checksum: {:s}'.format(zen_msg))
+    # Verify the checksum and retrieve the zip archive name
+    if os.path.isfile(zenodo_checksum):
+        with open(zenodo_checksum, 'r') as zcheck:
+            csum, zip_name = zcheck.read().split()
 
-    # Remove the checksum file if the download problem wasn't found there
-    if zenodo_checksum is not None:
-        os.remove(zenodo_checksum)
+            # Set the archive name
+            archive_name = os.path.join(out_dir, zip_name)
 
-    # Get the archive name from the output
-    try:
-        link_ind = zen_split.index('Link:') + 1
+            if not os.path.isfile(archive_name):
+                raise IOError(
+                    'error downloading archive to output dir: {:}'.format(
+                        archive_name))
 
-        # If the archive is already available, message may differ
-        zip_name = os.path.split(zen_split[link_ind])
-        while zip_name[-1].find(".zip") <= 0:
-            zip_name = os.path.split(zip_name[0])
+            # Get the archive checkshum
+            with open(os.path.join(out_dir, archive_name), 'rb') as afile:
+                asum = hashlib.md5(afile.read()).hexdigest()
 
-        # Set the archive name
-        archive_name = os.path.join(out_dir, zip_name[-1])
-    except (ValueError, IndexError):
-        raise IOError('unable to identify zenodo archive: {:}'.format(zen_msg))
+        zen_msg = "Checksum is correct" if asum == csum else "Corrupted data"
+    else:
+        zen_msg = "Checksum file not created"
 
-    if not os.path.isfile(archive_name):
-        raise IOError('error downloading archive to output dir: {:}'.format(
-            archive_name))
+    # Remove the checksum file
+    if zenodo_checksum is not None:
+        os.remove(zenodo_checksum)
+
+    # Evaluate checksum output
+    if zen_msg.find('Checksum is correct') < 0 and zen_msg.find(
+            'already downloaded correctly') < 0:
+        raise IOError('Bad checksum: {:s}'.format(zen_msg))
 
     # Access the zip archive
     with zipfile.ZipFile(archive_name, 'r') as zref:

From da02ef7482f0f22ce74ff78440956b7c882f5d62 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Wed, 13 May 2026 14:34:28 -0400
Subject: [PATCH 38/43] BUG: debug statement

See if coveralls works when all tests pass.
---
 .github/workflows/main.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 1a07e83d..7b73b967 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -14,7 +14,7 @@ jobs:
       matrix:
         os: ["ubuntu-latest", "macos-latest", "windows-latest"]
         python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"]
-        install-extras: ["base", "pysat_instruments", "dmsp_ssj"]
+        install-extras: ["base", "dmsp_ssj"] # , "pysat_instruments"]
 
     name: Python ${{ matrix.python-version }} on ${{ matrix.os }} with ${{ matrix.install-extras }}
     runs-on: ${{ matrix.os }}

From 8a379aeb2ba0a1c695d7fe15c59bc050dfacbde9 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Thu, 14 May 2026 16:55:03 -0400
Subject: [PATCH 39/43] BUG: fix coveralls implementation

Debug with help from coveralls!
---
 .github/workflows/main.yml | 11 ++++-------
 1 file changed, 4 insertions(+), 7 deletions(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 7b73b967..d24f861e 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -53,7 +53,6 @@ jobs:
     - name:  Coveralls Parallel
       uses: coverallsapp/github-action@v2
       with:
-        github-token: ${{ secrets.GITHUB_TOKEN }}
         flag-name: run=${{ join(matrix.*, '-') }}
         parallel: true
         format: cobertura
@@ -65,9 +64,7 @@ jobs:
     runs-on: "ubuntu-latest"
     steps:
       - name: Coveralls Finished
-        env:
-          COVERALLS_REPO_TOKEN: ${{ secrets.COVERALLS_REPO_TOKEN }}
-          COVERALLS_PARALLEL: true
-        run: |
-          curl -sL https://coveralls.io/coveralls-linux.tar.gz | tar -xz
-          ./coveralls done --build-number ${{ github.run_number }}
+        uses: coverallsapp/github-action@v2
+        with:
+          parallel-finished: true
+          carryforward: "all"

From fc591361cd6d672b74cecfded84b640dbc1fbc28 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Thu, 14 May 2026 17:34:37 -0400
Subject: [PATCH 40/43] BUG: re-add pysat

Re-add the pysat tests, for the use case in which they are expected to pass.
---
 .github/workflows/main.yml | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index d24f861e..eaa5e41a 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -14,7 +14,12 @@ jobs:
       matrix:
         os: ["ubuntu-latest", "macos-latest", "windows-latest"]
         python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"]
-        install-extras: ["base", "dmsp_ssj"] # , "pysat_instruments"]
+        install-extras: ["base", "dmsp_ssj"]
+        include:
+          # Pysat test for only pysat-passing Python version
+          - python-version: "3.10"
+            os: "ubuntu-latest"
+            install-extras: "pysat-instrument"
 
     name: Python ${{ matrix.python-version }} on ${{ matrix.os }} with ${{ matrix.install-extras }}
     runs-on: ${{ matrix.os }}
@@ -50,7 +55,7 @@ jobs:
         coverage report
         coverage xml --rcfile=pyproject.toml
 
-    - name:  Coveralls Parallel
+    - name: Coveralls Parallel
       uses: coverallsapp/github-action@v2
       with:
         flag-name: run=${{ join(matrix.*, '-') }}

From 3f8d076a2d28f93e0a68654b5b4569a3e6055f87 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Fri, 22 May 2026 10:16:08 -0400
Subject: [PATCH 41/43] BUG: fixed `install-extras` name

Fixed the name to 'pysat_instruments' from 'pysat-instrument'
---
 .github/workflows/main.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index eaa5e41a..328f70bd 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -19,7 +19,7 @@ jobs:
           # Pysat test for only pysat-passing Python version
           - python-version: "3.10"
             os: "ubuntu-latest"
-            install-extras: "pysat-instrument"
+            install-extras: "pysat_instruments"
 
     name: Python ${{ matrix.python-version }} on ${{ matrix.os }} with ${{ matrix.install-extras }}
     runs-on: ${{ matrix.os }}

From 92da111c9e0b190e83b017806f2b96858f172a4d Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Fri, 22 May 2026 14:41:19 -0400
Subject: [PATCH 42/43] REL: updated version number

Bump the version number for release.
---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 82e41de0..f654358f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "ocbpy"
-version = "0.6.0"
+version = "0.7.0"
 license = {file = "LICENSE"}
 description = 'Convert between magnetic/geodetic and adaptive, polar boundary coordinates'
 maintainers = [

From 7c505427939020c391fd4496db077ac7f2bf3ce1 Mon Sep 17 00:00:00 2001
From: "Angeline G. Burrell" <aburrell@users.noreply.github.com>
Date: Fri, 22 May 2026 14:41:36 -0400
Subject: [PATCH 43/43] REL: set target release date

Set a target release date for next week.
---
 Changelog.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Changelog.rst b/Changelog.rst
index f63db94a..0c532ad6 100644
--- a/Changelog.rst
+++ b/Changelog.rst
@@ -3,7 +3,7 @@ Changelog
 
 Summary of all changes made since the first stable release
 
-0.7.0 (XX-XX-2025)
+0.7.0 (05-29-2026)
 ------------------
 * ENH: Added the Gussenhoven (1983) model for the EAB
 * ENH: Added the CH-Aurora-2014 model for the OCB and EAB