CH-Aurora-2014 model

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

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

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

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