diff --git a/pyproject.toml b/pyproject.toml index 66b70ab4..2a4b2e0a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -117,6 +117,9 @@ analysis = [ "Bug Tracker" = "https://github.com/blalterman/SolarWindPy/issues" "Source" = "https://github.com/blalterman/SolarWindPy" +[tool.setuptools.package-data] +solarwindpy = ["core/data/*.csv"] + [tool.pip-tools] # pip-compile configuration for lockfile generation generate-hashes = false # Set to true for security-critical deployments diff --git a/setup.cfg b/setup.cfg index 9a3d1227..0cbe0c2d 100644 --- a/setup.cfg +++ b/setup.cfg @@ -10,7 +10,7 @@ tests_require = [flake8] extend-select = D402, D413, D205, D406 -ignore = E501, W503, D100, D101, D102, D103, D104, D105, D200, D202, D209, D214, D215, D300, D302, D400, D401, D403, D404, D405, D409, D412, D414 +ignore = E231, E501, W503, D100, D101, D102, D103, D104, D105, D200, D202, D209, D214, D215, D300, D302, D400, D401, D403, D404, D405, D409, D412, D414 enable = W605 docstring-convention = numpy max-line-length = 88 diff --git a/solarwindpy/core/__init__.py b/solarwindpy/core/__init__.py index b4e4bc06..db86118f 100644 --- a/solarwindpy/core/__init__.py +++ b/solarwindpy/core/__init__.py @@ -8,6 +8,7 @@ from .spacecraft import Spacecraft from .units_constants import Units, Constants from .alfvenic_turbulence import AlfvenicTurbulence +from .abundances import ReferenceAbundances __all__ = [ "Base", @@ -20,4 +21,5 @@ "Units", "Constants", "AlfvenicTurbulence", + "ReferenceAbundances", ] diff --git a/solarwindpy/core/abundances.py b/solarwindpy/core/abundances.py new file mode 100644 index 00000000..9cec4d69 --- /dev/null +++ b/solarwindpy/core/abundances.py @@ -0,0 +1,103 @@ +__all__ = ["ReferenceAbundances"] + +import numpy as np +import pandas as pd +from collections import namedtuple +from pathlib import Path + +Abundance = namedtuple("Abundance", "measurement,uncertainty") + + +class ReferenceAbundances: + """Elemental abundances from Asplund et al. (2009). + + Provides both photospheric and meteoritic abundances. + + References + ---------- + 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. + https://doi.org/10.1146/annurev.astro.46.060407.145222 + """ + + def __init__(self): + self.load_data() + + @property + def data(self): + r"""Elemental abundances in dex scale: + + log ε_X = log(N_X/N_H) + 12 + + where N_X is the number density of species X. + """ + 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 get_element(self, key, kind="Photosphere"): + r"""Get measurements for element stored at `key`. + + Parameters + ---------- + key : str or int + Element symbol ('Fe') or atomic number (26). + kind : str, default "Photosphere" + Which abundance source: "Photosphere" or "Meteorites". + """ + if isinstance(key, str): + level = "Symbol" + elif isinstance(key, int): + level = "Z" + else: + raise ValueError(f"Unrecognized key type ({type(key)})") + + out = self.data.loc[:, kind].xs(key, axis=0, level=level) + assert out.shape[0] == 1 + return out.iloc[0] + + @staticmethod + def _convert_from_dex(case): + m = case.loc["Ab"] + u = case.loc["Uncert"] + mm = 10.0 ** (m - 12.0) + uu = mm * np.log(10) * u + return mm, uu + + def abundance_ratio(self, numerator, denominator): + r"""Calculate abundance ratio N_X/N_Y with uncertainty. + + Parameters + ---------- + numerator, denominator : str or int + Element symbols ('Fe', 'O') or atomic numbers. + + Returns + ------- + Abundance + namedtuple with (measurement, uncertainty). + """ + top = self.get_element(numerator) + tu = top.Uncert + if np.isnan(tu): + tu = 0 + + if denominator != "H": + bottom = self.get_element(denominator) + bu = bottom.Uncert + if np.isnan(bu): + bu = 0 + + rat = 10.0 ** (top.Ab - bottom.Ab) + uncert = rat * np.log(10) * np.sqrt((tu**2) + (bu**2)) + else: + rat, uncert = self._convert_from_dex(top) + + return Abundance(rat, uncert) diff --git a/solarwindpy/core/data/asplund2009.csv b/solarwindpy/core/data/asplund2009.csv new file mode 100644 index 00000000..32d1ea3a --- /dev/null +++ b/solarwindpy/core/data/asplund2009.csv @@ -0,0 +1,90 @@ +Chemical composition of the Sun from Table 1 in [1]. + +[1] 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. https://doi.org/10.1146/annurev.astro.46.060407.145222 + +Kind,,Meteorites,Meteorites,Photosphere,Photosphere +,,Ab,Uncert,Ab,Uncert +Z,Symbol,,,, +1,H,8.22 , 0.04,12.00, +2,He,1.29,,10.93 , 0.01 +3,Li,3.26 , 0.05,1.05 , 0.10 +4,Be,1.30 , 0.03,1.38 , 0.09 +5,B,2.79 , 0.04,2.70 , 0.20 +6,C,7.39 , 0.04,8.43 , 0.05 +7,N,6.26 , 0.06,7.83 , 0.05 +8,O,8.40 , 0.04,8.69 , 0.05 +9,F,4.42 , 0.06,4.56 , 0.30 +10,Ne,-1.12,,7.93 , 0.10 +11,Na,6.27 , 0.02,6.24 , 0.04 +12,Mg,7.53 , 0.01,7.60 , 0.04 +13,Al,6.43 , 0.01,6.45 , 0.03 +14,Si,7.51 , 0.01,7.51 , 0.03 +15,P,5.43 , 0.04,5.41 , 0.03 +16,S,7.15 , 0.02,7.12 , 0.03 +17,Cl,5.23 , 0.06,5.50 , 0.30 +18,Ar,-0.05,,6.40 , 0.13 +19,K,5.08 , 0.02,5.03 , 0.09 +20,Ca,6.29 , 0.02,6.34 , 0.04 +21,Sc,3.05 , 0.02,3.15 , 0.04 +22,Ti,4.91 , 0.03,4.95 , 0.05 +23,V,3.96 , 0.02,3.93 , 0.08 +24,Cr,5.64 , 0.01,5.64 , 0.04 +25,Mn,5.48 , 0.01,5.43 , 0.04 +26,Fe,7.45 , 0.01,7.50 , 0.04 +27,Co,4.87 , 0.01,4.99 , 0.07 +28,Ni,6.20 , 0.01,6.22 , 0.04 +29,Cu,4.25 , 0.04,4.19 , 0.04 +30,Zn,4.63 , 0.04,4.56 , 0.05 +31,Ga,3.08 , 0.02,3.04 , 0.09 +32,Ge,3.58 , 0.04,3.65 , 0.10 +33,As,2.30 , 0.04,, +34,Se,3.34 , 0.03,, +35,Br,2.54 , 0.06,, +36,Kr,-2.27,,3.25 , 0.06 +37,Rb,2.36 , 0.03,2.52 , 0.10 +38,Sr,2.88 , 0.03,2.87 , 0.07 +39,Y,2.17 , 0.04,2.21 , 0.05 +40,Zr,2.53 , 0.04,2.58 , 0.04 +41,Nb,1.41 , 0.04,1.46 , 0.04 +42,Mo,1.94 , 0.04,1.88 , 0.08 +44,Ru,1.76 , 0.03,1.75 , 0.08 +45,Rh,1.06 , 0.04,0.91 , 0.10 +46,Pd,1.65 , 0.02,1.57 , 0.10 +47,Ag,1.20 , 0.02,0.94 , 0.10 +48,Cd,1.71 , 0.03,, +49,In,0.76 , 0.03,0.80 , 0.20 +50,Sn,2.07 , 0.06,2.04 , 0.10 +51,Sb,1.01 , 0.06,, +52,Te,2.18 , 0.03,, +53,I,1.55 , 0.08,, +54,Xe,-1.95,,2.24 , 0.06 +55,Cs,1.08 , 0.02,, +56,Ba,2.18 , 0.03,2.18 , 0.09 +57,La,1.17 , 0.02,1.10 , 0.04 +58,Ce,1.58 , 0.02,1.58 , 0.04 +59,Pr,0.76 , 0.03,0.72 , 0.04 +60,Nd,1.45 , 0.02,1.42 , 0.04 +62,Sm,0.94 , 0.02,0.96 , 0.04 +63,Eu,0.51 , 0.02,0.52 , 0.04 +64,Gd,1.05 , 0.02,1.07 , 0.04 +65,Tb,0.32 , 0.03,0.30 , 0.10 +66,Dy,1.13 , 0.02,1.10 , 0.04 +67,Ho,0.47 , 0.03,0.48 , 0.11 +68,Er,0.92 , 0.02,0.92 , 0.05 +69,Tm,0.12 , 0.03,0.10 , 0.04 +70,Yb,0.92 , 0.02,0.84 , 0.11 +71,Lu,0.09 , 0.02,0.10 , 0.09 +72,Hf,0.71 , 0.02,0.85 , 0.04 +73,Ta,-0.12 , 0.04,, +74,W,0.65 , 0.04,0.85 , 0.12 +75,Re,0.26 , 0.04,, +76,Os,1.35 , 0.03,1.40 , 0.08 +77,Ir,1.32 , 0.02,1.38 , 0.07 +78,Pt,1.62 , 0.03,, +79,Au,0.80 , 0.04,0.92 , 0.10 +80,Hg,1.17 , 0.08,, +81,Tl,0.77 , 0.03,0.90 , 0.20 +82,Pb,2.04 , 0.03,1.75 , 0.10 +83,Bi,0.65 , 0.04,, +90,Th,0.06 , 0.03,0.02 , 0.10 +92,U,-0.54 , 0.03,, diff --git a/tests/core/test_abundances.py b/tests/core/test_abundances.py new file mode 100644 index 00000000..a045add1 --- /dev/null +++ b/tests/core/test_abundances.py @@ -0,0 +1,213 @@ +"""Tests for ReferenceAbundances class. + +Tests verify: +1. Data structure matches expected CSV format +2. Values match published Asplund 2009 Table 1 +3. Uncertainty propagation formula is correct +4. Edge cases (NaN, H denominator) handled properly + +Run: pytest tests/core/test_abundances.py -v +""" + +import numpy as np +import pandas as pd +import pytest + +from solarwindpy.core.abundances import ReferenceAbundances, Abundance + + +class TestDataStructure: + """Verify CSV loads with correct structure.""" + + @pytest.fixture + def ref(self): + return ReferenceAbundances() + + def test_data_is_dataframe(self, ref): + # NOT: assert ref.data is not None (trivial) + # GOOD: Verify specific type + assert isinstance( + ref.data, pd.DataFrame + ), f"Expected DataFrame, got {type(ref.data)}" + + def test_data_has_83_elements(self, ref): + # Verify row count matches Asplund Table 1 + assert ( + ref.data.shape[0] == 83 + ), f"Expected 83 elements (Asplund Table 1), got {ref.data.shape[0]}" + + def test_index_is_multiindex_with_z_symbol(self, ref): + assert isinstance( + ref.data.index, pd.MultiIndex + ), f"Expected MultiIndex, got {type(ref.data.index)}" + assert list(ref.data.index.names) == [ + "Z", + "Symbol", + ], f"Expected index levels ['Z', 'Symbol'], got {ref.data.index.names}" + + def test_columns_have_photosphere_and_meteorites(self, ref): + top_level = ref.data.columns.get_level_values(0).unique().tolist() + assert "Photosphere" in top_level, "Missing 'Photosphere' column group" + assert "Meteorites" in top_level, "Missing 'Meteorites' column group" + + def test_data_dtype_is_float64(self, ref): + # All values should be float64 after .astype(np.float64) + for col in ref.data.columns: + assert ( + ref.data[col].dtype == np.float64 + ), f"Column {col} has dtype {ref.data[col].dtype}, expected float64" + + def test_h_has_nan_photosphere_uncertainty(self, ref): + # H photosphere uncertainty is NaN (by definition, H is the reference) + h = ref.get_element("H") + assert np.isnan(h.Uncert), f"H uncertainty should be NaN, got {h.Uncert}" + + def test_arsenic_photosphere_is_nan(self, ref): + # As (Z=33) has no photospheric measurement (only meteoritic) + arsenic = ref.get_element("As", kind="Photosphere") + assert np.isnan( + arsenic.Ab + ), f"As photosphere Ab should be NaN, got {arsenic.Ab}" + + +class TestGetElement: + """Verify element lookup by symbol and Z.""" + + @pytest.fixture + def ref(self): + return ReferenceAbundances() + + def test_get_element_by_symbol_returns_series(self, ref): + fe = ref.get_element("Fe") + assert isinstance(fe, pd.Series), f"Expected Series, got {type(fe)}" + + def test_iron_photosphere_matches_asplund(self, ref): + # Asplund 2009 Table 1: Fe = 7.50 +/- 0.04 + fe = ref.get_element("Fe") + assert np.isclose( + fe.Ab, 7.50, atol=0.01 + ), f"Fe photosphere Ab: expected 7.50, got {fe.Ab}" + assert np.isclose( + fe.Uncert, 0.04, atol=0.01 + ), f"Fe photosphere Uncert: expected 0.04, got {fe.Uncert}" + + def test_get_element_by_z_matches_symbol(self, ref): + # Z=26 is Fe, should return identical data values + # Note: Series names differ (26 vs 'Fe') but values are identical + by_symbol = ref.get_element("Fe") + by_z = ref.get_element(26) + pd.testing.assert_series_equal(by_symbol, by_z, check_names=False) + + def test_get_element_meteorites_differs_from_photosphere(self, ref): + # Fe meteorites: 7.45 vs photosphere: 7.50 + photo = ref.get_element("Fe", kind="Photosphere") + meteor = ref.get_element("Fe", kind="Meteorites") + assert ( + photo.Ab != meteor.Ab + ), "Photosphere and Meteorites should have different values" + assert np.isclose( + meteor.Ab, 7.45, atol=0.01 + ), f"Fe meteorites Ab: expected 7.45, got {meteor.Ab}" + + def test_invalid_key_type_raises_valueerror(self, ref): + with pytest.raises(ValueError, match="Unrecognized key type"): + ref.get_element(3.14) # float is invalid + + def test_unknown_element_raises_keyerror(self, ref): + with pytest.raises(KeyError, match="Xx"): + ref.get_element("Xx") # No element Xx + + def test_invalid_kind_raises_keyerror(self, ref): + with pytest.raises(KeyError, match="Invalid"): + ref.get_element("Fe", kind="Invalid") + + +class TestAbundanceRatio: + """Verify ratio calculation with uncertainty propagation.""" + + @pytest.fixture + def ref(self): + return ReferenceAbundances() + + def test_returns_abundance_namedtuple(self, ref): + result = ref.abundance_ratio("Fe", "O") + assert isinstance( + result, Abundance + ), f"Expected Abundance namedtuple, got {type(result)}" + assert hasattr(result, "measurement"), "Missing 'measurement' attribute" + assert hasattr(result, "uncertainty"), "Missing 'uncertainty' attribute" + + def test_fe_o_ratio_matches_computed_value(self, ref): + # Fe/O = 10^(7.50 - 8.69) = 0.06457 + result = ref.abundance_ratio("Fe", "O") + expected = 10.0 ** (7.50 - 8.69) + assert np.isclose( + result.measurement, expected, rtol=0.01 + ), f"Fe/O ratio: expected {expected:.5f}, got {result.measurement:.5f}" + + def test_fe_o_uncertainty_matches_formula(self, ref): + # sigma = ratio * ln(10) * sqrt(sigma_Fe^2 + sigma_O^2) + # sigma = 0.06457 * 2.303 * sqrt(0.04^2 + 0.05^2) = 0.00951 + result = ref.abundance_ratio("Fe", "O") + expected_ratio = 10.0 ** (7.50 - 8.69) + expected_uncert = expected_ratio * np.log(10) * np.sqrt(0.04**2 + 0.05**2) + assert np.isclose( + result.uncertainty, expected_uncert, rtol=0.01 + ), f"Fe/O uncertainty: expected {expected_uncert:.5f}, got {result.uncertainty:.5f}" + + def test_c_o_ratio_matches_computed_value(self, ref): + # C/O = 10^(8.43 - 8.69) = 0.5495 + result = ref.abundance_ratio("C", "O") + expected = 10.0 ** (8.43 - 8.69) + assert np.isclose( + result.measurement, expected, rtol=0.01 + ), f"C/O ratio: expected {expected:.4f}, got {result.measurement:.4f}" + + def test_ratio_destructuring_works(self, ref): + # Verify namedtuple can be destructured + measurement, uncertainty = ref.abundance_ratio("Fe", "O") + assert isinstance(measurement, float), "measurement should be float" + assert isinstance(uncertainty, float), "uncertainty should be float" + + +class TestHydrogenDenominator: + """Verify special case when denominator is H.""" + + @pytest.fixture + def ref(self): + return ReferenceAbundances() + + def test_fe_h_uses_convert_from_dex(self, ref): + # Fe/H = 10^(7.50 - 12) = 3.162e-5 + result = ref.abundance_ratio("Fe", "H") + expected = 10.0 ** (7.50 - 12.0) + assert np.isclose( + result.measurement, expected, rtol=0.01 + ), f"Fe/H ratio: expected {expected:.3e}, got {result.measurement:.3e}" + + def test_fe_h_uncertainty_from_numerator_only(self, ref): + # H has no uncertainty, so sigma = Fe_linear * ln(10) * sigma_Fe + result = ref.abundance_ratio("Fe", "H") + fe_linear = 10.0 ** (7.50 - 12.0) + expected_uncert = fe_linear * np.log(10) * 0.04 + assert np.isclose( + result.uncertainty, expected_uncert, rtol=0.01 + ), f"Fe/H uncertainty: expected {expected_uncert:.3e}, got {result.uncertainty:.3e}" + + +class TestNaNHandling: + """Verify NaN uncertainties are replaced with 0 in ratio calculations.""" + + @pytest.fixture + def ref(self): + return ReferenceAbundances() + + def test_ratio_with_nan_uncertainty_uses_zero(self, ref): + # H/O should use 0 for H's uncertainty + # sigma = ratio * ln(10) * sqrt(0^2 + sigma_O^2) = ratio * ln(10) * sigma_O + result = ref.abundance_ratio("H", "O") + expected_ratio = 10.0 ** (12.00 - 8.69) + expected_uncert = expected_ratio * np.log(10) * 0.05 # Only O contributes + assert np.isclose( + result.uncertainty, expected_uncert, rtol=0.01 + ), f"H/O uncertainty: expected {expected_uncert:.2f}, got {result.uncertainty:.2f}"