Models ne consistency

[xandrd]

The testing and documentation were create with help of GPT. Needs to be revise for consistency and accuracy. Add example of the models as in figures in the paper.

[xandrd]

Consider the following implementation

import numpy as np

def D2004(Rmax, which="ps_D2004", MLT=None):
    r"""
    Equatorial electron density :math:`n_{e0}` vs. :math:`R_{\max}` from
    Denton et al. (2004) [#]_ — core fits used in their Figure 6.

    This routine implements the *R* dependence formulas that Denton et al.
    fitted to Polar/CRRES data and compared against earlier models. Output is
    always in :math:`\mathrm{cm^{-3}}`.

    Parameters
    ----------
    Rmax : float or ndarray
        Maximum geocentric distance of the field line (apex), in Earth radii
        (:math:`R_E`). For a dipole this is numerically equal to :math:`L`.
    which : {'ps_D2004', 'ps_CA1992', 'ps_Gal2000', 'tr_D2004', 'tr_D2004_thr14'}
        Selects the specific curve:

        - ``'ps_D2004'`` – Plasmasphere fit using eq. (5) with
          :math:`a_1=3.78`, :math:`a_2=-0.324` (their best fit; valid
          :math:`2.5 \le \bar{R}_{\max} \le 8.0`).
        - ``'ps_CA1992'`` – Carpenter & Anderson (1992) saturated
          plasmasphere density recast into eq. (5) with
          :math:`a_1^{CA}=3.90`, :math:`a_2^{CA}=-0.315`.
        - ``'ps_Gal2000'`` – Gallagher et al. (2000) recast into eq. (5) with
          :math:`a_1^{Gal}=5.3`, :math:`a_2^{Gal}=-0.79`.
        - ``'tr_D2004'`` – Plasmatrough fit using eq. (12):
          :math:`\log_{10} n_{e0} = 3.77 - 3.45 \log_{10} R_{\max}`;
          valid for :math:`2.9 \lesssim R_{\max} \le 7.5`.
        - ``'tr_D2004_thr14'`` – Alternate plasmatrough fit (using the
          Sheeley-criterion separation), eq. (15):
          :math:`\log_{10} n_{e0} = 4.13 - 4.10 \log_{10} R_{\max}`.

    MLT : float or ndarray, optional
        Magnetic local time in hours. Not used by the options implemented here
        (kept for forward-compatibility; e.g., Denton cite Sheeley et al.'s
        MLT-dependent trough form). Ignored.

    Returns
    -------
    ne0 : float or ndarray
        Equatorial electron density at the apex, :math:`\mathrm{cm^{-3}}`.
        Shape matches the broadcasted shape of inputs.

    Notes
    -----
    Denton et al. write their plasmaspheric curves in terms of (their eq. 5)

    .. math::
        \log_{10}\!\left(\frac{n_{e0}}{\mathrm{cm}^{-3}}\right)
        = a_1 + a_2 \,\bar{R}_{\max}
        + \underbrace{\left(0.00127 R_{\odot} - 0.0635\right)}_{\equiv C_{\odot}}
          \exp\!\left[-\frac{(R_{\max}-2)}{1.5}\right],
        \tag{5}

    with :math:`\bar{R}_{\max} = R_{\max}/R_E` (eq. 7), which equals
    :math:`L` for a dipole. For their data span they take the average sunspot
    number :math:`R_{\odot}\approx13`, so
    :math:`C_{\odot} \approx -0.04699`.

    The specific coefficients used here are

    .. math::
        a_1^{CA}=3.90,\quad a_2^{CA}=-0.315 , \tag{6} \\
        a_1=3.78,\quad a_2=-0.324 \quad (\text{best fit}), \tag{9}\\
        a_1^{Gal}=5.3,\quad a_2^{Gal}=-0.79. \tag{10}

    For the plasmatrough they fit power laws (their eqs. 12 and 15):

    .. math::
        \log_{10}\!\left(\frac{n_{e0}}{\mathrm{cm}^{-3}}\right)
        =(3.77) + (-3.45)\log_{10} R_{\max}, \tag{12}

    and (using the Sheeley-density criterion to separate plasmasphere/trough)

    .. math::
        \log_{10}\!\left(\frac{n_{e0}}{\mathrm{cm}^{-3}}\right)
        =(4.13) + (-4.10)\log_{10} R_{\max}. \tag{15}

    Domains (quoted/rounded from the paper):

    - Plasmasphere fits (eq. 5 with the listed :math:`a_1, a_2`):
      :math:`2.5 \le \bar{R}_{\max} \le 8.0`.
    - Trough fits (eq. 12): :math:`2.9 \lesssim R_{\max} \le 7.5`.
      Eq. (15) was obtained on a reduced range (see text around their Fig. 6).

    References
    ----------
    .. [#] Denton, R. E., Takahashi, K., Galkin, I. A., Nsumei, P. A., Huang, X.,
       Reinisch, B. W., et al. (2004). *R* :sub:`max` dependence of
       magnetospheric equatorial electron density. **J. Geophys. Res.**
       109(A09215). https://doi.org/10.1029/2003JA010201
    """
    R = np.asarray(Rmax, dtype=float)
    # Average sunspot factor used in the paper (text around eq. 5):
    C_sun = 0.00127 * 13.0 - 0.0635  # ≈ -0.04699
    exp_term = np.exp(-(R - 2.0) / 1.5)

    which = str(which).lower()

    if which == "ps_d2004":
        a1, a2 = 3.78, -0.324
        lg10 = a1 + a2 * R + C_sun * exp_term
        return 10.0 ** lg10

    if which == "ps_ca1992":
        a1, a2 = 3.90, -0.315
        lg10 = a1 + a2 * R + C_sun * exp_term
        return 10.0 ** lg10

    if which == "ps_gal2000":
        a1, a2 = 5.3, -0.79
        lg10 = a1 + a2 * R + C_sun * exp_term
        return 10.0 ** lg10

    if which == "tr_d2004":
        # eq. (12): log10(ne0) = 3.77 + (-3.45) log10 Rmax
        return 10.0 ** (3.77 - 3.45 * np.log10(R))

    if which in ("tr_d2004_thr14", "tr_d2004_alt"):
        # eq. (15): alternate trough fit when using eq. (14) separation
        return 10.0 ** (4.13 - 4.10 * np.log10(R))

    raise ValueError(
        "Unknown 'which'. Use one of: "
        "'ps_D2004', 'ps_CA1992', 'ps_Gal2000', 'tr_D2004', 'tr_D2004_thr14'."
    )