From bf06e8e79713f6d15467fe53122867caa0642481 Mon Sep 17 00:00:00 2001 From: "Angeline G. Burrell" Date: Mon, 11 May 2026 18:02:38 -0400 Subject: [PATCH 1/4] 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" Date: Mon, 11 May 2026 18:02:58 -0400 Subject: [PATCH 2/4] 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" Date: Mon, 11 May 2026 18:03:24 -0400 Subject: [PATCH 3/4] 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" Date: Mon, 11 May 2026 18:03:43 -0400 Subject: [PATCH 4/4] 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