feat(core): Update ReferenceAbundances to Asplund 2021 with year selection

solarwindpy/__init__.py

@@ -18,7 +18,7 @@
     spacecraft,
     alfvenic_turbulence,
 )
-from . import core, plotting, solar_activity, tools, fitfunctions, data
+from . import core, plotting, solar_activity, tools, fitfunctions
 from . import instabilities  # noqa: F401
 from . import reproducibility
 
@@ -31,6 +31,7 @@ def _configure_pandas() -> None:
 _configure_pandas()
 
 Plasma = core.plasma.Plasma
+ReferenceAbundances = core.abundances.ReferenceAbundances
 at = alfvenic_turbulence
 sc = spacecraft
 pp = plotting
@@ -41,8 +42,8 @@ def _configure_pandas() -> None:
 
 __all__ = [
     "core",
-    "data",
     "plasma",
+    "ReferenceAbundances",
     "ions",
     "tensor",
     "vector",

solarwindpy/core/__init__.py

@@ -8,7 +8,7 @@
 from .spacecraft import Spacecraft
 from .units_constants import Units, Constants
 from .alfvenic_turbulence import AlfvenicTurbulence
-from .abundances import ReferenceAbundances
+from .abundances import ReferenceAbundances, Abundance
 
 __all__ = [
     "Base",
@@ -22,4 +22,5 @@
     "Constants",
     "AlfvenicTurbulence",
     "ReferenceAbundances",
+    "Abundance",
 ]

solarwindpy/core/abundances.py

@@ -1,46 +1,131 @@
-__all__ = ["ReferenceAbundances"]
+"""Reference elemental abundances from Asplund et al. (2009, 2021).
+
+This module provides access to solar photospheric and CI chondrite
+(meteoritic) abundances from the Asplund reference papers.
+
+References
+----------
+Asplund, M., Amarsi, A. M., & Grevesse, N. (2021).
+The chemical make-up of the Sun: A 2020 vision.
+A&A, 653, A141. https://doi.org/10.1051/0004-6361/202140445
+
+Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. (2009).
+The Chemical Composition of the Sun.
+Annu. Rev. Astron. Astrophys., 47, 481-522.
+https://doi.org/10.1146/annurev.astro.46.060407.145222
+"""
+
+__all__ = ["ReferenceAbundances", "Abundance"]
 
 import numpy as np
 import pandas as pd
 from collections import namedtuple
-from pathlib import Path
+from importlib import resources
 
 Abundance = namedtuple("Abundance", "measurement,uncertainty")
 
+# Alias mapping for backward compatibility
+_KIND_ALIASES = {
+    "Meteorites": "CI_chondrites",
+}
+
 
 class ReferenceAbundances:
-    """Elemental abundances from Asplund et al. (2009).
+    """Elemental abundances from Asplund et al. (2009, 2021).
+
+    Provides photospheric and CI chondrite (meteoritic) abundances
+    in the standard dex scale: log ε_X = log(N_X/N_H) + 12.
 
-    Provides both photospheric and meteoritic abundances.
+    Parameters
+    ----------
+    year : int, default 2021
+        Reference year: 2009 or 2021. Default uses Asplund 2021.
+
+    Attributes
+    ----------
+    data : pd.DataFrame
+        MultiIndex DataFrame with abundances and uncertainties.
+    year : int
+        The reference year for the loaded data.
 
     References
     ----------
+    Asplund, M., Amarsi, A. M., & Grevesse, N. (2021).
+    The chemical make-up of the Sun: A 2020 vision.
+    A&A, 653, A141. https://doi.org/10.1051/0004-6361/202140445
+
     Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. (2009).
     The Chemical Composition of the Sun.
-    Annual Review of Astronomy and Astrophysics, 47(1), 481–522.
+    Annu. Rev. Astron. Astrophys., 47, 481-522.
     https://doi.org/10.1146/annurev.astro.46.060407.145222
+
+    Examples
+    --------
+    >>> ref = ReferenceAbundances()  # doctest: +SKIP
+    >>> fe = ref.get_element("Fe")  # doctest: +SKIP
+    >>> print(f"Fe = {fe.Ab:.2f} ± {fe.Uncert:.2f}")  # doctest: +SKIP
+    Fe = 7.46 ± 0.04
+
+    Using 2009 data:
+
+    >>> ref_2009 = ReferenceAbundances(year=2009)  # doctest: +SKIP
+    >>> fe_2009 = ref_2009.get_element("Fe")  # doctest: +SKIP
+    >>> print(f"Fe (2009) = {fe_2009.Ab:.2f}")  # doctest: +SKIP
+    Fe (2009) = 7.50
     """
 
-    def __init__(self):
-        self.load_data()
+    _VALID_YEARS = (2009, 2021)
+
+    def __init__(self, year=2021):
+        if not isinstance(year, int):
+            raise TypeError(f"year must be an integer, got {type(year).__name__}")
+        if year not in self._VALID_YEARS:
+            raise ValueError(f"year must be 2009 or 2021, got {year}")
+        self._year = year
+        self._load_data()
+
+    @property
+    def year(self):
+        """The reference year for the loaded data."""
+        return self._year
 
     @property
     def data(self):
-        r"""Elemental abundances in dex scale:
+        r"""Elemental abundances in dex scale.
 
-        log ε_X = log(N_X/N_H) + 12
+        The dex scale is defined as:
+            log ε_X = log(N_X/N_H) + 12
 
         where N_X is the number density of species X.
+
+        Returns
+        -------
+        pd.DataFrame
+            MultiIndex DataFrame with index (Z, Symbol) and columns
+            (CI_chondrites, Photosphere) × (Ab, Uncert).
         """
         return self._data
 
-    def load_data(self):
-        """Load Asplund 2009 data from package CSV."""
-        path = Path(__file__).parent / "data" / "asplund2009.csv"
-        data = pd.read_csv(path, skiprows=4, header=[0, 1], index_col=[0, 1]).astype(
-            np.float64
-        )
-        self._data = data
+    def _load_data(self):
+        """Load Asplund data from package CSV based on year."""
+        filename = f"asplund{self._year}.csv"
+        data_file = resources.files(__package__).joinpath("data", filename)
+
+        with data_file.open() as f:
+            data = pd.read_csv(f, skiprows=4, header=[0, 1], index_col=[0, 1])
+
+        # 2021 has Comment column, extract before float conversion
+        # Column is ('Comment', 'Unnamed: X_level_1') due to pandas MultiIndex parsing
+        comment_cols = [col for col in data.columns if col[0] == "Comment"]
+        if comment_cols:
+            comment_col = comment_cols[0]
+            self._comments = data[comment_col].copy()
+            data = data.drop(columns=[comment_col])
+        else:
+            self._comments = None
+
+        # Convert remaining columns to float64
+        self._data = data.astype(np.float64)
 
     def get_element(self, key, kind="Photosphere"):
         r"""Get measurements for element stored at `key`.
@@ -50,8 +135,41 @@ def get_element(self, key, kind="Photosphere"):
         key : str or int
             Element symbol ('Fe') or atomic number (26).
         kind : str, default "Photosphere"
-            Which abundance source: "Photosphere" or "Meteorites".
+            Which abundance source: "Photosphere", "CI_chondrites",
+            or "Meteorites" (alias for CI_chondrites).
+
+        Returns
+        -------
+        pd.Series
+            Series with 'Ab' (abundance in dex) and 'Uncert' (uncertainty).
+
+        Raises
+        ------
+        ValueError
+            If key is not a string or integer.
+        KeyError
+            If element not found or invalid kind.
+
+        Examples
+        --------
+        >>> ref = ReferenceAbundances()  # doctest: +SKIP
+        >>> ref.get_element("Fe")  # doctest: +SKIP
+        Ab        7.46
+        Uncert    0.04
+        Name: 26, dtype: float64
+        >>> ref.get_element(26)  # Same result using atomic number  # doctest: +SKIP
         """
+        # Handle backward compatibility alias
+        kind = _KIND_ALIASES.get(kind, kind)
+
+        # Validate kind
+        valid_kinds = ["Photosphere", "CI_chondrites"]
+        if kind not in valid_kinds:
+            raise KeyError(
+                f"Invalid kind '{kind}'. Must be one of: {valid_kinds} "
+                f"(or 'Meteorites' as alias for 'CI_chondrites')"
+            )
+
         if isinstance(key, str):
             level = "Symbol"
         elif isinstance(key, int):
@@ -63,8 +181,71 @@ def get_element(self, key, kind="Photosphere"):
         assert out.shape[0] == 1
         return out.iloc[0]
 
+    def get_comment(self, key):
+        """Get the source comment for an element (2021 data only).
+
+        The comment indicates the source methodology for elements where
+        the adopted abundance is not from photospheric spectroscopy:
+        - 'definition': H abundance is defined as 12.00
+        - 'helioseismology': Derived from helioseismology (He)
+        - 'meteorites': Adopted from CI chondrite measurements
+        - 'solar wind': Derived from solar wind measurements (Ne, Ar, Kr)
+        - 'nuclear physics': Derived from nuclear physics (Xe)
+
+        Parameters
+        ----------
+        key : str or int
+            Element symbol ('Fe') or atomic number (26).
+
+        Returns
+        -------
+        str or None
+            The comment string, or None if no comment (spectroscopic
+            measurement) or if using 2009 data.
+
+        Examples
+        --------
+        >>> ref = ReferenceAbundances()  # doctest: +SKIP
+        >>> ref.get_comment("H")  # doctest: +SKIP
+        'definition'
+        >>> print(ref.get_comment("Fe"))  # Spectroscopic, no comment  # doctest: +SKIP
+        None
+        """
+        if self._comments is None:
+            return None
+
+        if isinstance(key, str):
+            level = "Symbol"
+        elif isinstance(key, int):
+            level = "Z"
+        else:
+            raise ValueError(f"Unrecognized key type ({type(key)})")
+
+        try:
+            comment = self._comments.xs(key, axis=0, level=level)
+            if len(comment) == 1:
+                comment = comment.iloc[0]
+            # Handle empty strings and NaN
+            if pd.isna(comment) or comment == "":
+                return None
+            return comment
+        except KeyError:
+            return None
+
     @staticmethod
     def _convert_from_dex(case):
+        """Convert from dex to linear abundance ratio relative to H.
+
+        Parameters
+        ----------
+        case : pd.Series
+            Series with 'Ab' and 'Uncert' in dex.
+
+        Returns
+        -------
+        tuple
+            (measurement, uncertainty) in linear units.
+        """
         m = case.loc["Ab"]
         u = case.loc["Uncert"]
         mm = 10.0 ** (m - 12.0)
@@ -83,6 +264,21 @@ def abundance_ratio(self, numerator, denominator):
         -------
         Abundance
             namedtuple with (measurement, uncertainty).
+
+        Notes
+        -----
+        Uncertainty is propagated assuming independent uncertainties:
+            σ_ratio = ratio × ln(10) × √(σ_X² + σ_Y²)
+
+        For denominator='H', uses the special conversion from dex
+        since H is the reference element (log ε_H = 12 by definition).
+
+        Examples
+        --------
+        >>> ref = ReferenceAbundances()  # doctest: +SKIP
+        >>> fe_o = ref.abundance_ratio("Fe", "O")  # doctest: +SKIP
+        >>> print(f"Fe/O = {fe_o.measurement:.4f} ± {fe_o.uncertainty:.4f}")  # doctest: +SKIP
+        Fe/O = 0.0589 ± 0.0077
         """
         top = self.get_element(numerator)
         tu = top.Uncert