diff --git a/examples/idl_comparisons/continuum_comparison.py b/examples/idl_comparisons/continuum_comparison.py index 219d3df9..bc05cc1e 100644 --- a/examples/idl_comparisons/continuum_comparison.py +++ b/examples/idl_comparisons/continuum_comparison.py @@ -65,7 +65,8 @@ def plot_idl_comparison(wavelength, temperature, result_fiasco, result_idl): # thermal bremsstrahlung. idl_result_freefree = read_idl_test_output('freefree_all_ions', LATEST_VERSION) ion_kwargs = {'abundance': idl_result_freefree['abundance'], - 'ionization_fraction': idl_result_freefree['ionization_fraction']} + 'ionization_fraction': idl_result_freefree['ionization_fraction'], + 'proton_electron_ratio': 0.83} all_ions = [fiasco.Ion(ion_name, idl_result_freefree['temperature'], **ion_kwargs) for ion_name in fiasco.list_ions()] all_ions = fiasco.IonCollection(*all_ions) free_free = all_ions.free_free(idl_result_freefree['wavelength']) @@ -83,7 +84,8 @@ def plot_idl_comparison(wavelength, temperature, result_fiasco, result_idl): # continuum emission. idl_result_freebound = read_idl_test_output('freebound_all_ions', LATEST_VERSION) ion_kwargs = {'abundance': idl_result_freebound['abundance'], - 'ionization_fraction': idl_result_freebound['ionization_fraction']} + 'ionization_fraction': idl_result_freebound['ionization_fraction'], + 'proton_electron_ratio': 0.83} all_ions = [fiasco.Ion(ion_name, idl_result_freebound['temperature'], **ion_kwargs) for ion_name in fiasco.list_ions()] all_ions = fiasco.IonCollection(*all_ions) free_bound = all_ions.free_bound(idl_result_freebound['wavelength']) @@ -101,7 +103,8 @@ def plot_idl_comparison(wavelength, temperature, result_fiasco, result_idl): # continuum emission. idl_result_twophoton = read_idl_test_output('twophoton_all_ions', LATEST_VERSION) ion_kwargs = {'abundance': idl_result_twophoton['abundance'], - 'ionization_fraction': idl_result_twophoton['ionization_fraction']} + 'ionization_fraction': idl_result_twophoton['ionization_fraction'], + 'proton_electron_ratio': 0.83} all_ions = [fiasco.Ion(ion_name, idl_result_twophoton['temperature'], **ion_kwargs) for ion_name in fiasco.list_ions()] all_ions = fiasco.IonCollection(*all_ions) two_photon = all_ions.two_photon(idl_result_twophoton['wavelength'], diff --git a/fiasco/base.py b/fiasco/base.py index af857f37..fa64647a 100644 --- a/fiasco/base.py +++ b/fiasco/base.py @@ -16,6 +16,7 @@ from fiasco.io.generic import GenericIonParser from fiasco.util import check_database, parse_ion_name, periodic_table_period from fiasco.util.exceptions import MissingIonError +from fiasco.util.tools import _get_atomic_symbol, _get_element_name __all__ = ['IonBase'] @@ -39,20 +40,28 @@ class IonBase: """ def __init__(self, ion_name, hdf5_dbase_root=None, **kwargs): - # base rep is a tuple of integers (atomic_number, ionization_stage) - self._base_rep = parse_ion_name(ion_name) if hdf5_dbase_root is None: self.hdf5_dbase_root = fiasco.defaults['hdf5_dbase_root'] else: self.hdf5_dbase_root = hdf5_dbase_root - check_database(self.hdf5_dbase_root, **kwargs) - if self.ion_name not in fiasco.list_ions(self.hdf5_dbase_root, sort=False): - raise MissingIonError(f'{self.ion_name} not found in {self.hdf5_dbase_root}') # Put import here to avoid circular imports from fiasco import log self.log = log - # Warn users if the database they are using is potentially stale. - self._check_dbase_fiasco_version() + # This optional kwarg is here to avoid costly checks in cases where one + # is sure the database already exists, the ion name is valid and already + # in the right forma, and the version is not potentially stale. This is + # not documented publicly as it is not meant to be used except by other + # internal fiasco functions where the inputs are already validated. + if kwargs.pop('safe_mode', True): + # base rep is a tuple of integers (atomic_number, ionization_stage) + self._base_rep = parse_ion_name(ion_name) + check_database(self.hdf5_dbase_root, **kwargs) + if self.ion_name not in fiasco.list_ions(self.hdf5_dbase_root, sort=False): + raise MissingIonError(f'{self.ion_name} not found in {self.hdf5_dbase_root}') + # Warn users if the database they are using is potentially stale. + self._check_dbase_fiasco_version() + else: + self._base_rep = ion_name def _check_dbase_fiasco_version(self): "Warn if database was generated with an earlier version of fiasco." @@ -70,17 +79,17 @@ def _check_dbase_fiasco_version(self): @property def atomic_number(self): """The atomic number of the element, :math:`Z`.""" - return plasmapy.particles.atomic_number(self._base_rep[0]) + return self._base_rep[0] @property def element_name(self): """The full name of the element, e.g. "hydrogen".""" - return plasmapy.particles.element_name(self.atomic_number) + return _get_element_name(self.atomic_number) @property def atomic_symbol(self): """The standard atomic symbol for the element, e.g. "H" for hydrogen.""" - return plasmapy.particles.atomic_symbol(self.atomic_number) + return _get_atomic_symbol(self.atomic_number) @property def ion_name(self): @@ -101,7 +110,7 @@ def charge_state(self): def isoelectronic_sequence(self): "Atomic symbol denoting to which isoelectronic sequence this ion belongs." if (Z_iso := self.atomic_number - self.charge_state) > 0: - return plasmapy.particles.atomic_symbol(Z_iso) + return _get_atomic_symbol(Z_iso) @property @u.quantity_input diff --git a/fiasco/elements.py b/fiasco/elements.py index 37c201db..7917e6aa 100644 --- a/fiasco/elements.py +++ b/fiasco/elements.py @@ -39,6 +39,9 @@ def __init__(self, element_name, temperature: u.K, **kwargs): element_name = element_name.capitalize() Z = plasmapy.particles.atomic_number(element_name) ion_list = [] + if (pe_ratio := kwargs.pop('proton_electron_ratio', None)) is None: + pe_ratio = fiasco.proton_electron_ratio(temperature, **kwargs) + kwargs['proton_electron_ratio'] = pe_ratio for i in range(Z + 1): ion = fiasco.Ion((Z, i+1), temperature, **kwargs) ion_list.append(ion) diff --git a/fiasco/fiasco.py b/fiasco/fiasco.py index 27c30978..1a50ca2b 100644 --- a/fiasco/fiasco.py +++ b/fiasco/fiasco.py @@ -48,7 +48,7 @@ def list_elements(hdf5_dbase_root=None, sort=True): return elements -def list_ions(hdf5_dbase_root=None, sort=True): +def list_ions(hdf5_dbase_root=None, sort=True, base_rep=False): """ List all available ions in the CHIANTI database @@ -58,6 +58,10 @@ def list_ions(hdf5_dbase_root=None, sort=True): If not specified, will default to that specified in ``fiasco.defaults``. sort: `bool`, optional If True, sort the list of elements by increasing atomic number. + base_rep: `bool`, optional + If True, return list of ions as tuples of ``(Z, z)``, where ``Z`` is the + atomic number and ``z`` is the ionization stage. This format may be useful + when constructing a large list of `~fiasco.Ion` objects. """ if hdf5_dbase_root is None: hdf5_dbase_root = fiasco.defaults['hdf5_dbase_root'] @@ -82,7 +86,13 @@ def list_ions(hdf5_dbase_root=None, sort=True): # NOTE: when grabbing straight from the index and not sorting, the result will be # a numpy array. Cast to a list to make sure the return type is consistent for # all possible inputs - return ions.tolist() if isinstance(ions, np.ndarray) else ions + ion_list = ions.tolist() if isinstance(ions, np.ndarray) else ions + if base_rep: + # NOTE: Explicitly not using parse_ion_name here as it can be slow. + elements = list_elements(hdf5_dbase_root=hdf5_dbase_root, sort=False) + el_map = {el: plasmapy.particles.atomic_number(el) for el in elements} + ion_list = [(el_map[el], int(ion)) for el, ion in map(lambda x: x.split(), ion_list)] + return ion_list def get_dem_model(model, hdf5_dbase_root=None): @@ -156,32 +166,40 @@ def proton_electron_ratio(temperature: u.K, **kwargs): """ # Import here to avoid circular imports from fiasco import log + + # NOTE: Set this to avoid infinite recursion. The exact value is arbitrary because it + # is not used in this calculation. + kwargs['proton_electron_ratio'] = 0.0 h_2 = fiasco.Ion('H +1', temperature, **kwargs) numerator = h_2.abundance * h_2._ion_fraction[h_2._instance_kwargs['ionization_fraction']]['ionization_fraction'] denominator = u.Quantity(np.zeros(numerator.shape)) - for el_name in list_elements(h_2.hdf5_dbase_root): - el = fiasco.Element(el_name, temperature, **h_2._instance_kwargs) + abund_file = h_2._instance_kwargs['abundance'] + ionization_file = h_2._instance_kwargs['ionization_fraction'] + for ion_name in list_ions(hdf5_dbase_root=h_2.hdf5_dbase_root, sort=False, base_rep=True): + # NOTE: Using IonBase to avoid the overhead of repeatedly constructing an Ion object. + ion = fiasco.base.IonBase(ion_name, + hdf5_dbase_root=h_2._instance_kwargs['hdf5_dbase_root'], + safe_mode=False) try: - abundance = el.abundance + abundance = ion._abund[abund_file] except KeyError: - abund_file = el[0]._instance_kwargs['abundance'] log.warning( - f'Not including {el.atomic_symbol}. Abundance not available from {abund_file}.') + f'Not including {ion.ion_name_roman}. Abundance not available from {abund_file}.') continue - for ion in el: - ionization_file = ion._instance_kwargs['ionization_fraction'] - # NOTE: We use ._ion_fraction here rather than .ionization_fraction to avoid - # doing an interpolation to the temperature array every single time and instead only - # interpolate once at the end. - # It is assumed that the ionization_fraction temperature array for each ion is the same. - try: - ionization_fraction = ion._ion_fraction[ionization_file]['ionization_fraction'] - t_ionization_fraction = ion._ion_fraction[ionization_file]['temperature'] - except KeyError: - log.warning( - f'Not including {ion.ion_name}. Ionization fraction not available from {ionization_file}.') - continue - denominator += ionization_fraction * abundance * ion.charge_state + # NOTE: We use ._ion_fraction here rather than .ionization_fraction to avoid + # doing an interpolation to the temperature array every single time and instead only + # interpolate once at the end. + # It is assumed that the ionization_fraction temperature array for each ion is the same. + try: + ion_fraction_ds = ion._ion_fraction[ionization_file] + except KeyError: + log.warning( + f'Not including {ion.ion_name}. Ionization fraction not available from {ionization_file}.') + continue + else: + ionization_fraction = ion_fraction_ds['ionization_fraction'] + t_ionization_fraction = ion_fraction_ds['temperature'] + denominator += ionization_fraction * abundance * ion.charge_state ratio = numerator / denominator f_interp = interp1d(t_ionization_fraction.to(temperature.unit).value, diff --git a/fiasco/ions.py b/fiasco/ions.py index dceede9f..5d7ebc89 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -56,6 +56,11 @@ class Ion(IonBase): ionization_potential : `str` or `~astropy.units.Quantity`, optional If a string is provided, use the appropriate "ip" dataset. If a scalar value is provided, use that value for the ionization potential. This value should be convertible to eV. + proton_electron_ratio : `float` or array-like, optional + Ratio of proton to electron densities, :math:`n_H/n_e`, as a function of temperature. Can be a scalar value or + an array with the same shape as ``temperature``. If not specified, this is calculated using + `~fiasco.proton_electron_ratio`. When instantiating many ions, it may be more efficient to precompute this + quantity and then pass it to the constructor through this keyword argument. """ @u.quantity_input @@ -65,6 +70,7 @@ def __init__(self, abundance='sun_coronal_1992_feldman_ext', ionization_fraction='chianti', ionization_potential='chianti', + proton_electron_ratio=None, *args, **kwargs): super().__init__(ion_name, *args, **kwargs) @@ -73,6 +79,7 @@ def __init__(self, self.abundance = abundance self.ionization_fraction = ionization_fraction self.ionization_potential = ionization_potential + self.proton_electron_ratio = proton_electron_ratio self.gaunt_factor = GauntFactor(hdf5_dbase_root=self.hdf5_dbase_root) def _new_instance(self, temperature=None, **kwargs): @@ -81,9 +88,14 @@ def _new_instance(self, temperature=None, **kwargs): possibly different arguments. If different arguments are not specified, this will just create a copy of itself. """ + new_kwargs = self._instance_kwargs if temperature is None: temperature = self.temperature.copy() - new_kwargs = self._instance_kwargs + else: + # If a new temperature array is specified, this could now be stale + # so we remove the old value and force it to be recomputed on the + # updated temperature array if a new value is not specified. + new_kwargs.pop('proton_electron_ratio', None) new_kwargs.update(kwargs) return type(self)(self.ion_name, temperature, **new_kwargs) @@ -146,6 +158,12 @@ def _instance_kwargs(self): kwargs['ionization_fraction'] = self.ionization_fraction if kwargs['ionization_potential'] is None: kwargs['ionization_potential'] = self.ionization_potential + # NOTE: It is possible that this property is needed prior to the proton/electron ratio + # being set. This logic guards against that. + try: + kwargs['proton_electron_ratio'] = self.proton_electron_ratio + except AttributeError: + self.log.debug('Proton/electron ratio not added to instance kwargs.') return kwargs def _has_dataset(self, dset_name): @@ -203,10 +221,19 @@ def thermal_energy(self) -> u.erg: """ return self.temperature.to('erg', equivalencies=u.equivalencies.temperature_energy()) - @cached_property + @property @u.quantity_input def proton_electron_ratio(self) -> u.dimensionless_unscaled: - return proton_electron_ratio(self.temperature, **self._instance_kwargs) + """ + Ratio of proton to electron number density as a function of temperature. + """ + return self._proton_electron_ratio + + @proton_electron_ratio.setter + def proton_electron_ratio(self, value): + if value is None: + value = proton_electron_ratio(self.temperature, **self._instance_kwargs) + self._proton_electron_ratio = u.Quantity(value*np.ones(self.temperature.shape)) def next_ion(self): """ @@ -1214,8 +1241,7 @@ def emissivity(self, density: u.cm**(-3), **kwargs) -> u.erg * u.cm**(-3) / u.s: contribution_function : Calculate contribution function, :math:`G(n,T)` """ density = np.atleast_1d(density) - pe_ratio = proton_electron_ratio(self.temperature, **self._instance_kwargs) - pe_ratio = pe_ratio[:, np.newaxis, np.newaxis] + pe_ratio = self.proton_electron_ratio[:, np.newaxis, np.newaxis] g = self.contribution_function(density, **kwargs) density_squared = density**2 couple_density_to_temperature = kwargs.get('couple_density_to_temperature', False) diff --git a/fiasco/tests/data/test_file_list.json b/fiasco/tests/data/test_file_list.json index 6f853f08..aa801929 100644 --- a/fiasco/tests/data/test_file_list.json +++ b/fiasco/tests/data/test_file_list.json @@ -151,6 +151,7 @@ "ar_18.scups", "ar_18.wgfa", "ar_19.rrparams", + "archive/sun_coronal_1992_feldman.abund", "archive/sun_coronal_1992_feldman_ext.abund", "archive/sun_photospheric_2007_grevesse.abund", "b_1.diparams", diff --git a/fiasco/tests/idl/test_idl_continuum.py b/fiasco/tests/idl/test_idl_continuum.py index 335d2905..df36ff9d 100644 --- a/fiasco/tests/idl/test_idl_continuum.py +++ b/fiasco/tests/idl/test_idl_continuum.py @@ -22,6 +22,9 @@ def ion_input_args(): return { 'abundance': 'sun_coronal_1992_feldman_ext', 'ionization_fraction': 'chianti', + # NOTE: this isn't used anywhere in the continuum calculations + # so just set it to a scalar value for efficiency. + 'proton_electron_ratio': 0.83, } diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py index 427aa3c2..f12d059f 100644 --- a/fiasco/tests/test_ion.py +++ b/fiasco/tests/test_ion.py @@ -70,7 +70,9 @@ def test_new_instance(ion): abundance = ion._instance_kwargs['abundance'] new_ion = ion._new_instance() for k in new_ion._instance_kwargs: - assert new_ion._instance_kwargs[k] == ion._instance_kwargs[k] + if k != 'proton_electron_ratio': + assert new_ion._instance_kwargs[k] == ion._instance_kwargs[k] + assert u.allclose(new_ion.proton_electron_ratio, ion.proton_electron_ratio, rtol=0) assert u.allclose(new_ion.temperature, ion.temperature, rtol=0) new_ion = ion._new_instance(temperature=ion.temperature[:1]) assert u.allclose(new_ion.temperature, ion.temperature[:1]) @@ -230,9 +232,9 @@ def test_proton_collision(fe10): def test_missing_abundance(hdf5_dbase_root): _ion = fiasco.Ion('Li 1', - temperature, - abundance='sun_coronal_1992_feldman', - hdf5_dbase_root=hdf5_dbase_root) + temperature, + abundance='sun_coronal_1992_feldman', + hdf5_dbase_root=hdf5_dbase_root) with pytest.raises(MissingDatasetException): _ion.abundance @@ -582,6 +584,28 @@ def test_ionization_fraction_setter(ion, ioneq_input, ioneq_output): assert u.allclose(ion._instance_kwargs['ionization_fraction'], ioneq_input) +@pytest.mark.parametrize('value', [ + 0.83, + 0.83 * np.ones(temperature.shape), +]) +def test_proton_electron_ratio_setter(ion, value): + ion.proton_electron_ratio = value + assert ion.proton_electron_ratio.shape == ion.temperature.shape + assert u.allclose(ion.proton_electron_ratio, value) + new_ion = ion._new_instance(proton_electron_ratio=value) + assert new_ion.proton_electron_ratio.shape == new_ion.temperature.shape + assert u.allclose(new_ion.proton_electron_ratio, value) + + +@pytest.mark.requires_dbase_version('>= 8') +def test_emissivity_uses_proton_electron_ratio(ion): + # Setting the ratio to 0 should zero the emissivity, which is only the case + # if emissivity uses the (cached) property rather than recomputing the ratio. + ion.proton_electron_ratio = 0.0 + emm = ion.emissivity(1e7 * u.cm**-3) + assert u.allclose(emm, 0 * u.erg / u.cm**3 / u.s) + + def test_ionization_fraction_setter_exception(ion): # This should fail because the input has len>1 but is not the same # shape as the temperature array diff --git a/fiasco/util/tools.py b/fiasco/util/tools.py index 15f6221b..c81b78cc 100644 --- a/fiasco/util/tools.py +++ b/fiasco/util/tools.py @@ -8,6 +8,10 @@ from functools import partial from scipy.interpolate import splev, splrep +_MAX_Z = 40 +_ATOMIC_SYMBOL_LOOKUP = {z: plasmapy.particles.atomic_symbol(z) for z in range(1, _MAX_Z+1)} +_ELEMENT_NAME_LOOKUP = {z: plasmapy.particles.element_name(z) for z in range(1, _MAX_Z+1)} + __all__ = [ 'vectorize_where', 'vectorize_where_sum', @@ -235,3 +239,25 @@ def periodic_table_period(element): if element <= r: return i+1 raise ValueError(f'No period available for {element=}.') + + +def _get_atomic_symbol(Z): + """ + Return atomic symbol for a given atomic number. + + .. note:: This function is simply a fast version of + plasmapy.particles.atomic_symbol and is only + meant for internal use. + """ + return _ATOMIC_SYMBOL_LOOKUP[Z] + + +def _get_element_name(Z): + """ + Return element name for a given atomic number. + + .. note:: This function is simply a fast version of + plasmapy.particles.element_name and is only meant for + internal use. + """ + return _ELEMENT_NAME_LOOKUP[Z]