From cdaaa95fcc467d9f4c2570b754fdb0e55d8c01ff Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Wed, 7 Jan 2026 17:47:22 -0700 Subject: [PATCH 01/19] Adding logic to ion.n_levels to account for cases where more energy levels are described in the database than are populated in the model --- fiasco/ions.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/fiasco/ions.py b/fiasco/ions.py index 9dfd35e7..156494e0 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -162,7 +162,23 @@ def n_levels(self): Number of energy levels in the CHIANTI model """ try: - return len(self.levels) + n_levels_elvlc= len(self.levels) + + # some ions don't have .wgfa + try: + n_levels_wgfa = self._wgfa['upper_level'].max() + except KeyError: + return 0 + + n_levels_scups = self._scups['upper_level'].max() + + if n_levels_elvlc > n_levels_wgfa or n_levels_elvlc > n_levels_scups: + levels_list = [n_levels_elvlc, n_levels_wgfa, n_levels_scups] + n_levels = np.max([n_levels_scups, n_levels_wgfa]) + else: + n_levels = n_levels_elvlc + + return n_levels except MissingDatasetException: return 0 From 0bdc69d8503a62f31d95b1e6a900cabc05cda766 Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Thu, 8 Jan 2026 10:48:07 -0700 Subject: [PATCH 02/19] Implementing suggested changes, cleaning up print statements, adding test for N4. --- fiasco/ions.py | 18 ++++++++++-------- fiasco/tests/test_ion.py | 11 +++++++++++ 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/fiasco/ions.py b/fiasco/ions.py index 156494e0..c2d655cf 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -159,21 +159,23 @@ def _has_dataset(self, dset_name): @property def n_levels(self): """ - Number of energy levels in the CHIANTI model + Number of energy levels in the CHIANTI model. It is possible this number + will not match the number of levels in the elvlc file. """ try: n_levels_elvlc= len(self.levels) - # some ions don't have .wgfa - try: - n_levels_wgfa = self._wgfa['upper_level'].max() - except KeyError: - return 0 + if self._has_dataset('wgfa'): + n_levels_wgfa = self.transitions.upper_level.max() + else: + n_levels_wgfa = 0 - n_levels_scups = self._scups['upper_level'].max() + if self._has_dataset('scups'): + n_levels_scups = self._scups['upper_level'].max() + else: + n_levels_scups = 0 if n_levels_elvlc > n_levels_wgfa or n_levels_elvlc > n_levels_scups: - levels_list = [n_levels_elvlc, n_levels_wgfa, n_levels_scups] n_levels = np.max([n_levels_scups, n_levels_wgfa]) else: n_levels = n_levels_elvlc diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py index b6465d77..3c9eb073 100644 --- a/fiasco/tests/test_ion.py +++ b/fiasco/tests/test_ion.py @@ -50,6 +50,12 @@ def c5(hdf5_dbase_root): def c6(hdf5_dbase_root): return fiasco.Ion('C VI', temperature, hdf5_dbase_root=hdf5_dbase_root) +@pytest.fixture +def n4(hdf5_dbase_root): + # NOTE: This ion was added because there are more levels listed in elvlc than in wgfa. + # It is used in test_n_levels. + return fiasco.Ion('N 4', temperature, hdf5_dbase_root=hdf5_dbase_root) + @pytest.fixture def fe20(hdf5_dbase_root): @@ -93,6 +99,11 @@ def test_level(ion): assert level.orbital_angular_momentum == 2 assert level.is_observed +# @pytest.mark.requires_dbase_version('>= 10') +def test_n_levels(n4): + n_levels = n4.n_levels + assert n_levels != len(n4.levels) + assert n_levels == np.max([n4.transitions.upper_level.max(), n4._scups['upper_level'].max()]) def test_level_energy_parsing(fe10): # Test falling back to theoretical energies if energies are not observed From f37c8b16acc549bb53cc9394b83821a814a9b310 Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Thu, 8 Jan 2026 14:38:42 -0700 Subject: [PATCH 03/19] fixing typo docstring describing location of test_file_list.json --- tools/generate_test_file_list.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/generate_test_file_list.py b/tools/generate_test_file_list.py index 046a3f3b..8e76db92 100644 --- a/tools/generate_test_file_list.py +++ b/tools/generate_test_file_list.py @@ -17,7 +17,7 @@ def write_test_file_list(files): """ Write a unique, sorted list of files to use in the test database. - These are written to fiasco/util/data/test_files.json + These are written to fiasco/tests/data/test_files.json """ current_files = get_test_file_list() # Read current files current_files += files From 5ed63a56d0d8d16a4b1d02cbefffa26a4333c903 Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Thu, 8 Jan 2026 14:46:38 -0700 Subject: [PATCH 04/19] add files for test_n_levels on N 4 ion --- fiasco/tests/data/test_file_list.json | 3 +++ 1 file changed, 3 insertions(+) diff --git a/fiasco/tests/data/test_file_list.json b/fiasco/tests/data/test_file_list.json index c823b321..57e60df4 100644 --- a/fiasco/tests/data/test_file_list.json +++ b/fiasco/tests/data/test_file_list.json @@ -1128,11 +1128,14 @@ "n_3.drparams", "n_3.fblvl", "n_3.rrparams", + "n_5.elvlc", "n_4.diparams", "n_4.drparams", "n_4.easplups", "n_4.fblvl", "n_4.rrparams", + "n_4.scups", + "n_4.wgfa", "n_5.diparams", "n_5.drparams", "n_5.easplups", From 63d414013f47974d8229a1852b8720929ccadec7 Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Thu, 8 Jan 2026 14:47:16 -0700 Subject: [PATCH 05/19] typo --- fiasco/tests/data/test_file_list.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fiasco/tests/data/test_file_list.json b/fiasco/tests/data/test_file_list.json index 57e60df4..61441cfb 100644 --- a/fiasco/tests/data/test_file_list.json +++ b/fiasco/tests/data/test_file_list.json @@ -1129,7 +1129,7 @@ "n_3.fblvl", "n_3.rrparams", "n_5.elvlc", - "n_4.diparams", + "n_4.elvlc", "n_4.drparams", "n_4.easplups", "n_4.fblvl", From 67305000531f5544dca5590ab2d3c26af02bbe6d Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Thu, 8 Jan 2026 14:48:17 -0700 Subject: [PATCH 06/19] typo again --- fiasco/tests/data/test_file_list.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fiasco/tests/data/test_file_list.json b/fiasco/tests/data/test_file_list.json index 61441cfb..65646d15 100644 --- a/fiasco/tests/data/test_file_list.json +++ b/fiasco/tests/data/test_file_list.json @@ -1128,8 +1128,8 @@ "n_3.drparams", "n_3.fblvl", "n_3.rrparams", - "n_5.elvlc", "n_4.elvlc", + "n_4.diparams", "n_4.drparams", "n_4.easplups", "n_4.fblvl", From fe4b3079a7ce36e25a559c0d8834ae12c7e43915 Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Mon, 12 Jan 2026 17:59:08 -0700 Subject: [PATCH 07/19] Applying suggested changes --- fiasco/ions.py | 18 ++++-------------- fiasco/tests/test_ion.py | 2 -- 2 files changed, 4 insertions(+), 16 deletions(-) diff --git a/fiasco/ions.py b/fiasco/ions.py index c2d655cf..4f685e9e 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -160,25 +160,15 @@ def _has_dataset(self, dset_name): def n_levels(self): """ Number of energy levels in the CHIANTI model. It is possible this number - will not match the number of levels in the elvlc file. + will not match the number of levels in the energy level file. """ try: n_levels_elvlc= len(self.levels) - if self._has_dataset('wgfa'): - n_levels_wgfa = self.transitions.upper_level.max() - else: - n_levels_wgfa = 0 - - if self._has_dataset('scups'): - n_levels_scups = self._scups['upper_level'].max() - else: - n_levels_scups = 0 + n_levels_wgfa = self.transitions.upper_level.max() if self._has_dataset('wgfa') else 0 + n_levels_scups = self._scups['upper_level'].max() if self._has_dataset('scups') else 0 - if n_levels_elvlc > n_levels_wgfa or n_levels_elvlc > n_levels_scups: - n_levels = np.max([n_levels_scups, n_levels_wgfa]) - else: - n_levels = n_levels_elvlc + n_levels = min(n_levels_elvlc, np.max([n_levels_scups, n_levels_wgfa])) return n_levels except MissingDatasetException: diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py index 3c9eb073..9fcbff5f 100644 --- a/fiasco/tests/test_ion.py +++ b/fiasco/tests/test_ion.py @@ -99,10 +99,8 @@ def test_level(ion): assert level.orbital_angular_momentum == 2 assert level.is_observed -# @pytest.mark.requires_dbase_version('>= 10') def test_n_levels(n4): n_levels = n4.n_levels - assert n_levels != len(n4.levels) assert n_levels == np.max([n4.transitions.upper_level.max(), n4._scups['upper_level'].max()]) def test_level_energy_parsing(fe10): From 04cf01d6f267f43913f560dbc33b8ffc2db2c688 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Mon, 26 Jan 2026 14:14:57 -0500 Subject: [PATCH 08/19] account for more n_levels cases --- fiasco/ions.py | 30 ++++++++++++++++----------- fiasco/tests/data/test_file_list.json | 4 +++- fiasco/tests/test_ion.py | 20 +++++++++++++++--- 3 files changed, 38 insertions(+), 16 deletions(-) diff --git a/fiasco/ions.py b/fiasco/ions.py index 4f685e9e..6c7c9b95 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -159,20 +159,26 @@ def _has_dataset(self, dset_name): @property def n_levels(self): """ - Number of energy levels in the CHIANTI model. It is possible this number - will not match the number of levels in the energy level file. - """ - try: - n_levels_elvlc= len(self.levels) - - n_levels_wgfa = self.transitions.upper_level.max() if self._has_dataset('wgfa') else 0 - n_levels_scups = self._scups['upper_level'].max() if self._has_dataset('scups') else 0 - - n_levels = min(n_levels_elvlc, np.max([n_levels_scups, n_levels_wgfa])) + Number of energy levels in the atomic model. - return n_levels - except MissingDatasetException: + .. note:: It is possible this number will not match the number of levels + in the energy level file. + """ + # There is no atomic model if there are no energies, collisions, and decays + # NOTE: This logic is here rather than using a decorator so that it returns + # zero rather than throwing an exception. + if not all([self._has_dataset('elvlc'), + self._has_dataset('wgfa'), + self._has_dataset('scups')]): return 0 + n_elvlc = len(self.levels) + n_wgfa = max(self.transitions.lower_level.max(), self.transitions.upper_level.max()) + n_scups = max(self._scups['upper_level'].max(), self._scups['lower_level'].max()) + # If there is autoionization data associated with this ion, ensure that the model + # has enough levels to include these rates. + if self._has_dataset('auto'): + n_scups = max(n_scups, self._auto['upper_level'].max()) + return np.min([n_elvlc, n_wgfa, n_scups]) @property def n_transitions(self): diff --git a/fiasco/tests/data/test_file_list.json b/fiasco/tests/data/test_file_list.json index 65646d15..9d9d085d 100644 --- a/fiasco/tests/data/test_file_list.json +++ b/fiasco/tests/data/test_file_list.json @@ -701,6 +701,7 @@ "fe_4.fblvl", "fe_4.rrparams", "fe_4.trparams", + "fe_5.auto", "fe_5.diparams", "fe_5.drparams", "fe_5.easplups", @@ -1128,10 +1129,10 @@ "n_3.drparams", "n_3.fblvl", "n_3.rrparams", - "n_4.elvlc", "n_4.diparams", "n_4.drparams", "n_4.easplups", + "n_4.elvlc", "n_4.fblvl", "n_4.rrparams", "n_4.scups", @@ -1771,6 +1772,7 @@ "ti_21.rrparams", "ti_22.diparams", "ti_22.drparams", + "ti_22.elvlc", "ti_22.rrparams", "ti_23.rrparams", "v_1.diparams", diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py index 9fcbff5f..14ece480 100644 --- a/fiasco/tests/test_ion.py +++ b/fiasco/tests/test_ion.py @@ -99,9 +99,23 @@ def test_level(ion): assert level.orbital_angular_momentum == 2 assert level.is_observed -def test_n_levels(n4): - n_levels = n4.n_levels - assert n_levels == np.max([n4.transitions.upper_level.max(), n4._scups['upper_level'].max()]) + +@pytest.mark.parametrize(('ion_name', 'n_levels'), [ + # NOTE: The expected number of levels may change depending on the database version + # When they do, use a marker per parameterization to mark which ones apply to which + # version. + # NOTE: These ions are chosen because they cover the different cases of how to choose + # the number of atomic levels in a model based upon the available atomic data. + ('Fe V', 34), + ('N IV', 136), + ('Ar XVIII', 24), + ('Ti XXII', 0), + ('C IV', 923), +]) +def test_n_levels(ion_name, n_levels, hdf5_dbase_root): + test_ion = fiasco.Ion(ion_name, temperature, hdf5_dbase_root=hdf5_dbase_root) + test_ion.n_levels == n_levels + def test_level_energy_parsing(fe10): # Test falling back to theoretical energies if energies are not observed From f2bfd6545aa37f9037f7939aa67fb4837d98aece Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Mon, 26 Jan 2026 23:26:19 -0500 Subject: [PATCH 09/19] add all auto files to test database --- fiasco/tests/data/test_file_list.json | 46 +++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/fiasco/tests/data/test_file_list.json b/fiasco/tests/data/test_file_list.json index 9d9d085d..9836d0b4 100644 --- a/fiasco/tests/data/test_file_list.json +++ b/fiasco/tests/data/test_file_list.json @@ -5,6 +5,7 @@ "al_2.diparams", "al_2.drparams", "al_2.rrparams", + "al_3.auto", "al_3.diparams", "al_3.drparams", "al_3.easplups", @@ -39,11 +40,13 @@ "al_10.easplups", "al_10.fblvl", "al_10.rrparams", + "al_11.auto", "al_11.diparams", "al_11.drparams", "al_11.easplups", "al_11.fblvl", "al_11.rrparams", + "al_12.auto", "al_12.diparams", "al_12.drparams", "al_12.elvlc", @@ -88,6 +91,7 @@ "ar_7.easplups", "ar_7.fblvl", "ar_7.rrparams", + "ar_8.auto", "ar_8.diparams", "ar_8.drparams", "ar_8.easplups", @@ -122,11 +126,13 @@ "ar_15.easplups", "ar_15.fblvl", "ar_15.rrparams", + "ar_16.auto", "ar_16.diparams", "ar_16.drparams", "ar_16.easplups", "ar_16.fblvl", "ar_16.rrparams", + "ar_17.auto", "ar_17.diparams", "ar_17.drparams", "ar_17.elvlc", @@ -184,6 +190,7 @@ "c_3.easplups", "c_3.fblvl", "c_3.rrparams", + "c_4.auto", "c_4.diparams", "c_4.drparams", "c_4.easplups", @@ -192,6 +199,7 @@ "c_4.rrparams", "c_4.scups", "c_4.wgfa", + "c_5.auto", "c_5.diparams", "c_5.drparams", "c_5.elvlc", @@ -245,6 +253,7 @@ "ca_9.easplups", "ca_9.fblvl", "ca_9.rrparams", + "ca_10.auto", "ca_10.diparams", "ca_10.drparams", "ca_10.easplups", @@ -283,11 +292,13 @@ "ca_17.easplups", "ca_17.fblvl", "ca_17.rrparams", + "ca_18.auto", "ca_18.diparams", "ca_18.drparams", "ca_18.easplups", "ca_18.fblvl", "ca_18.rrparams", + "ca_19.auto", "ca_19.diparams", "ca_19.drparams", "ca_19.elvlc", @@ -429,6 +440,7 @@ "co_16.drparams", "co_16.easplups", "co_16.rrparams", + "co_17.auto", "co_17.diparams", "co_17.drparams", "co_17.easplups", @@ -514,6 +526,7 @@ "cr_13.drparams", "cr_13.easplups", "cr_13.rrparams", + "cr_14.auto", "cr_14.diparams", "cr_14.drparams", "cr_14.easplups", @@ -545,6 +558,7 @@ "cr_22.easplups", "cr_22.fblvl", "cr_22.rrparams", + "cr_23.auto", "cr_23.diparams", "cr_23.drparams", "cr_23.fblvl", @@ -794,6 +808,7 @@ "fe_17.rrparams", "fe_17.scups", "fe_17.wgfa", + "fe_18.auto", "fe_18.diparams", "fe_18.drparams", "fe_18.elvlc", @@ -801,6 +816,7 @@ "fe_18.rrparams", "fe_18.scups", "fe_18.wgfa", + "fe_19.auto", "fe_19.diparams", "fe_19.drparams", "fe_19.fblvl", @@ -827,20 +843,24 @@ "fe_21.rrparams", "fe_21.scups", "fe_21.wgfa", + "fe_22.auto", "fe_22.diparams", "fe_22.drparams", "fe_22.fblvl", "fe_22.rrparams", + "fe_23.auto", "fe_23.diparams", "fe_23.drparams", "fe_23.easplups", "fe_23.fblvl", "fe_23.rrparams", + "fe_24.auto", "fe_24.diparams", "fe_24.drparams", "fe_24.easplups", "fe_24.fblvl", "fe_24.rrparams", + "fe_25.auto", "fe_25.diparams", "fe_25.drparams", "fe_25.elvlc", @@ -914,6 +934,7 @@ "k_8.drparams", "k_8.easplups", "k_8.rrparams", + "k_9.auto", "k_9.diparams", "k_9.drparams", "k_9.easplups", @@ -974,6 +995,7 @@ "li_4.rrparams", "mg_1.diparams", "mg_1.fblvl", + "mg_2.auto", "mg_2.diparams", "mg_2.drparams", "mg_2.fblvl", @@ -1007,11 +1029,13 @@ "mg_9.easplups", "mg_9.fblvl", "mg_9.rrparams", + "mg_10.auto", "mg_10.diparams", "mg_10.drparams", "mg_10.easplups", "mg_10.fblvl", "mg_10.rrparams", + "mg_11.auto", "mg_11.diparams", "mg_11.drparams", "mg_11.elvlc", @@ -1082,6 +1106,7 @@ "mn_14.drparams", "mn_14.easplups", "mn_14.rrparams", + "mn_15.auto", "mn_15.diparams", "mn_15.drparams", "mn_15.easplups", @@ -1137,11 +1162,13 @@ "n_4.rrparams", "n_4.scups", "n_4.wgfa", + "n_5.auto", "n_5.diparams", "n_5.drparams", "n_5.easplups", "n_5.fblvl", "n_5.rrparams", + "n_6.auto", "n_6.diparams", "n_6.drparams", "n_6.elvlc", @@ -1228,11 +1255,13 @@ "ne_7.easplups", "ne_7.fblvl", "ne_7.rrparams", + "ne_8.auto", "ne_8.diparams", "ne_8.drparams", "ne_8.easplups", "ne_8.fblvl", "ne_8.rrparams", + "ne_9.auto", "ne_9.diparams", "ne_9.drparams", "ne_9.elvlc", @@ -1319,6 +1348,7 @@ "ni_17.easplups", "ni_17.fblvl", "ni_17.rrparams", + "ni_18.auto", "ni_18.diparams", "ni_18.drparams", "ni_18.easplups", @@ -1353,11 +1383,13 @@ "ni_25.easplups", "ni_25.fblvl", "ni_25.rrparams", + "ni_26.auto", "ni_26.diparams", "ni_26.drparams", "ni_26.easplups", "ni_26.fblvl", "ni_26.rrparams", + "ni_27.auto", "ni_27.diparams", "ni_27.drparams", "ni_27.elvlc", @@ -1394,6 +1426,7 @@ "o_5.easplups", "o_5.fblvl", "o_5.rrparams", + "o_6.auto", "o_6.diparams", "o_6.drparams", "o_6.easplups", @@ -1402,6 +1435,7 @@ "o_6.rrparams", "o_6.scups", "o_6.wgfa", + "o_7.auto", "o_7.diparams", "o_7.drparams", "o_7.elvlc", @@ -1431,6 +1465,7 @@ "p_4.drparams", "p_4.easplups", "p_4.rrparams", + "p_5.auto", "p_5.diparams", "p_5.drparams", "p_5.easplups", @@ -1496,6 +1531,7 @@ "s_5.easplups", "s_5.fblvl", "s_5.rrparams", + "s_6.auto", "s_6.diparams", "s_6.drparams", "s_6.easplups", @@ -1530,11 +1566,13 @@ "s_13.easplups", "s_13.fblvl", "s_13.rrparams", + "s_14.auto", "s_14.diparams", "s_14.drparams", "s_14.easplups", "s_14.fblvl", "s_14.rrparams", + "s_15.auto", "s_15.diparams", "s_15.drparams", "s_15.elvlc", @@ -1636,6 +1674,7 @@ "si_3.easplups", "si_3.fblvl", "si_3.rrparams", + "si_4.auto", "si_4.diparams", "si_4.drparams", "si_4.easplups", @@ -1670,11 +1709,13 @@ "si_11.easplups", "si_11.fblvl", "si_11.rrparams", + "si_12.auto", "si_12.diparams", "si_12.drparams", "si_12.easplups", "si_12.fblvl", "si_12.rrparams", + "si_13.auto", "si_13.diparams", "si_13.drparams", "si_13.elvlc", @@ -1735,6 +1776,7 @@ "ti_11.drparams", "ti_11.easplups", "ti_11.rrparams", + "ti_12.auto", "ti_12.diparams", "ti_12.drparams", "ti_12.easplups", @@ -1766,6 +1808,7 @@ "ti_20.easplups", "ti_20.fblvl", "ti_20.rrparams", + "ti_21.auto", "ti_21.diparams", "ti_21.drparams", "ti_21.fblvl", @@ -1926,6 +1969,7 @@ "zn_19.drparams", "zn_19.easplups", "zn_19.rrparams", + "zn_20.auto", "zn_20.diparams", "zn_20.drparams", "zn_20.easplups", @@ -1952,10 +1996,12 @@ "zn_27.drparams", "zn_27.easplups", "zn_27.rrparams", + "zn_28.auto", "zn_28.diparams", "zn_28.drparams", "zn_28.easplups", "zn_28.rrparams", + "zn_29.auto", "zn_29.diparams", "zn_29.drparams", "zn_29.elvlc", From ab202584406dd7142e640f08995b371871657f61 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Mon, 2 Feb 2026 11:30:33 -0500 Subject: [PATCH 10/19] limit wgfa data to number of levels in the model --- .pre-commit-config.yaml | 2 +- fiasco/ions.py | 34 +++++++++++++++++++++------------- fiasco/levels.py | 32 ++++++++++++++++++++++---------- 3 files changed, 44 insertions(+), 24 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index c748bc36..020aa5a9 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -4,7 +4,7 @@ ci: repos: - repo: https://github.com/astral-sh/ruff-pre-commit - rev: "v0.12.7" + rev: "v0.14.4" hooks: - id: ruff args: ["--fix"] diff --git a/fiasco/ions.py b/fiasco/ions.py index 6c7c9b95..97412f7d 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -172,7 +172,7 @@ def n_levels(self): self._has_dataset('scups')]): return 0 n_elvlc = len(self.levels) - n_wgfa = max(self.transitions.lower_level.max(), self.transitions.upper_level.max()) + n_wgfa = max(self._wgfa['lower_level'].max(), self._wgfa['upper_level'].max()) n_scups = max(self._scups['upper_level'].max(), self._scups['lower_level'].max()) # If there is autoionization data associated with this ion, ensure that the model # has enough levels to include these rates. @@ -229,7 +229,7 @@ def previous_ion(self): @needs_dataset('elvlc', 'wgfa') def transitions(self): "A `~fiasco.Transitions` object holding the information about transitions for this ion." - return Transitions(self.levels, self._wgfa) + return Transitions(self.levels, self._wgfa, n_levels=self.n_levels) @property @u.quantity_input @@ -709,12 +709,16 @@ def _rate_matrix_radiative_decay(self) -> u.Unit('s-1'): @u.quantity_input def _rate_matrix_collisional_electron(self) -> u.Unit('cm3 s-1'): rate_matrix = u.Quantity(np.zeros(self.temperature.shape + (self.n_levels, self.n_levels,)), 'cm3 s-1') - lower_index = self._scups['lower_level'] - 1 - upper_index = self._scups['upper_level'] - 1 + # NOTE: For some ions, there may be more rate data available than + # there are levels in the model. + idx = np.where(np.logical_and(self._scups['lower_level']<=self.n_levels, + self._scups['upper_level']<=self.n_levels)) + lower_index = self._scups['lower_level'][idx] - 1 + upper_index = self._scups['upper_level'][idx] - 1 # De-excitation from upper states - rate_matrix[:, lower_index, upper_index] += self.electron_collision_deexcitation_rate + rate_matrix[:, lower_index, upper_index] += self.electron_collision_deexcitation_rate[..., *idx] # Excitation from lower states - rate_matrix[:, upper_index, lower_index] += self.electron_collision_excitation_rate + rate_matrix[:, upper_index, lower_index] += self.electron_collision_excitation_rate[..., *idx] return rate_matrix @cached_property @@ -722,10 +726,14 @@ def _rate_matrix_collisional_electron(self) -> u.Unit('cm3 s-1'): @u.quantity_input def _rate_matrix_collisional_proton(self) -> u.Unit('cm3 s-1'): rate_matrix = u.Quantity(np.zeros(self.temperature.shape + (self.n_levels, self.n_levels,)), 'cm3 s-1') - lower_index = self._psplups['lower_level'] - 1 - upper_index = self._psplups['upper_level'] - 1 - rate_matrix[:, lower_index, upper_index] += self.proton_collision_deexcitation_rate - rate_matrix[:, upper_index, lower_index] += self.proton_collision_excitation_rate + # NOTE: For some ions, there may be more rate data available than + # there are levels in the model. + idx = np.where(np.logical_and(self._psplups['lower_level']<=self.n_levels, + self._psplups['upper_level']<=self.n_levels)) + lower_index = self._psplups['lower_level'][idx] - 1 + upper_index = self._psplups['upper_level'][idx] - 1 + rate_matrix[:, lower_index, upper_index] += self.proton_collision_deexcitation_rate[..., *idx] + rate_matrix[:, upper_index, lower_index] += self.proton_collision_excitation_rate[..., *idx] return rate_matrix def _empty_rate_matrix(self, temperature_dependent=True, unit='cm3 s-1'): @@ -849,10 +857,10 @@ def _rate_matrix_dielectronic_recombination(self) -> u.Unit('cm3 s-1'): self._auto['upper_level'][idx_ground], self._auto['autoionization_rate']) # Sum radiative decay rates between upper levels and lower bound levels - idx_bound = self._wgfa['lower_level'] < self._auto['upper_level'].min() - A_rad_sum = vectorize_where_sum(self._wgfa['upper_level'][idx_bound], + idx_bound = self.transitions.lower_level < self._auto['upper_level'].min() + A_rad_sum = vectorize_where_sum(self.transitions.upper_level[idx_bound], self._auto['upper_level'][idx_ground], - self._wgfa['A'][idx_bound],) + self.transitions.A[idx_bound],) branching_ratio = A_rad_sum / (A_rad_sum + A_auto_sum) # Get needed levels for recombined and recombining ions dc_rate = self._dielectronic_capture_rate(self._auto['lower_level'][idx_ground], diff --git a/fiasco/levels.py b/fiasco/levels.py index 265a0de9..6bb6fad4 100644 --- a/fiasco/levels.py +++ b/fiasco/levels.py @@ -133,13 +133,26 @@ class Transitions: Data structure holding information about all of the energy levels for a given ion in the CHIANTI database. wgfa: `~fiasco.io.datalayer.DataIndexer` - Pointer to the transition information for a given ion in + Table of transition information for a given ion in the CHIANTI database. + n_levels: `int`, optional + Maximum number of levels in the CHIANTI atomic model. This is used to + appropriately limit the transition information to the size of the atomic + model. Typically, this is calculated by `fiasco.Ion.n_levels`. """ - def __init__(self, levels, wgfa): + def __init__(self, levels, wgfa, n_levels=None): self._levels = levels self._wgfa = wgfa + if n_levels is None: + self._idx = np.s_[:] + else: + # NOTE: For some ions, there may be more rate data available than + # there are levels in the model. + self._idx = np.where(np.logical_and( + self._wgfa['lower_level']<=n_levels, + self._wgfa['upper_level']<=n_levels, + )) def __len__(self): return self.wavelength.shape[0] @@ -149,8 +162,7 @@ def is_twophoton(self): """ True if the transition is a two-photon decay """ - return np.logical_and(self.wavelength == 0, - self.upper_level<=10) + return np.logical_and(self.wavelength == 0, self.upper_level<=10) @property def is_autoionization(self): @@ -164,14 +176,14 @@ def is_bound_bound(self): """ True for bound-bound transitions. """ - return self._wgfa['wavelength'] != 0 + return self._wgfa['wavelength'][self._idx] != 0 @property def is_observed(self): """ True for transitions that connect two observed energy levels """ - return self._wgfa['wavelength'] > 0 + return self._wgfa['wavelength'][self._idx] > 0 @property @u.quantity_input @@ -179,23 +191,23 @@ def A(self) -> u.s**(-1): """ Spontaneous transition probability due to radiative decay """ - return self._wgfa['A'] + return self._wgfa['A'][self._idx] @property @u.quantity_input def wavelength(self) -> u.angstrom: "Wavelength of each transition." - return np.fabs(self._wgfa['wavelength']) + return np.fabs(self._wgfa['wavelength'][self._idx]) @property def upper_level(self): "Index of the upper level of the transition." - return self._wgfa['upper_level'] + return self._wgfa['upper_level'][self._idx] @property def lower_level(self): "Index of the lower level of the transition." - return self._wgfa['lower_level'] + return self._wgfa['lower_level'][self._idx] @property @u.quantity_input From d56923cf9a13d0d1fed7871cf9aafb2f20ab09ea Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Mon, 2 Feb 2026 11:30:57 -0500 Subject: [PATCH 11/19] move ruff config to separate file --- .ruff.toml | 49 +++++++++++++++++++++++++++++++++++++++++++++++++ pyproject.toml | 50 -------------------------------------------------- 2 files changed, 49 insertions(+), 50 deletions(-) create mode 100644 .ruff.toml diff --git a/.ruff.toml b/.ruff.toml new file mode 100644 index 00000000..c52100ac --- /dev/null +++ b/.ruff.toml @@ -0,0 +1,49 @@ +target-version = "py312" +line-length = 110 +exclude=[ + ".git,", + "__pycache__", + "build", + "fiasco/version.py", +] +show-fixes = true +output-format = "full" + +[lint] +select = [ + "E", + "F", + "W", + "UP", + "PT", + #"RET", + #"TID", +] +extend-ignore = [ + # pycodestyle (E, W) + "E501", # LineTooLong # TODO! fix + "E741", # Ambiguous variable name + + # pytest (PT) + "PT001", # Always use pytest.fixture() + "PT004", # Fixtures which don't return anything should have leading _ + "PT007", # Parametrize should be lists of tuples # TODO! fix + "PT011", # Too broad exception assert # TODO! fix + "PT023", # Always use () on pytest decorators +] + +[lint.per-file-ignores] +# Part of configuration, not a package. +"setup.py" = ["INP001"] +"conftest.py" = ["INP001"] +# implicit-namespace-package. The examples are not a package. +"docs/*.py" = ["INP001"] +# Module level imports do not need to be at the top of a file here +"docs/conf.py" = ["E402"] +"__init__.py" = ["E402", "F401", "F403"] +"test_*.py" = ["B011", "D", "E402", "S101"] +# Allow star import for the factory +"fiasco/io/factory.py" = ["F403"] + +[lint.pydocstyle] +convention = "numpy" diff --git a/pyproject.toml b/pyproject.toml index f3180ace..29c92063 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -143,53 +143,3 @@ lines_between_types = 1 [tool.codespell] skip = "*.fts,*.fits,venv,*.pro,*.bib,*.asdf,*.json" ignore-words-list = "te,emiss" - -[tool.ruff] -target-version = "py310" -line-length = 110 -exclude=[ - ".git,", - "__pycache__", - "build", - "fiasco/version.py", -] -show-fixes = true -output-format = "full" - -lint.select = [ - "E", - "F", - "W", - "UP", - "PT", - #"RET", - #"TID", -] -lint.extend-ignore = [ - # pycodestyle (E, W) - "E501", # LineTooLong # TODO! fix - "E741", # Ambiguous variable name - - # pytest (PT) - "PT001", # Always use pytest.fixture() - "PT004", # Fixtures which don't return anything should have leading _ - "PT007", # Parametrize should be lists of tuples # TODO! fix - "PT011", # Too broad exception assert # TODO! fix - "PT023", # Always use () on pytest decorators -] - -[tool.ruff.lint.per-file-ignores] -# Part of configuration, not a package. -"setup.py" = ["INP001"] -"conftest.py" = ["INP001"] -# implicit-namespace-package. The examples are not a package. -"docs/*.py" = ["INP001"] -# Module level imports do not need to be at the top of a file here -"docs/conf.py" = ["E402"] -"__init__.py" = ["E402", "F401", "F403"] -"test_*.py" = ["B011", "D", "E402", "S101"] -# Allow star import for the factory -"fiasco/io/factory.py" = ["F403"] - -[tool.ruff.lint.pydocstyle] -convention = "numpy" From 1c3c201057eac0d325bea1d5f28ba3cfb056b141 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Mon, 2 Feb 2026 11:31:39 -0500 Subject: [PATCH 12/19] use .transitions everywhere --- fiasco/tests/test_ion.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/fiasco/tests/test_ion.py b/fiasco/tests/test_ion.py index 14ece480..669dda2a 100644 --- a/fiasco/tests/test_ion.py +++ b/fiasco/tests/test_ion.py @@ -241,7 +241,8 @@ def test_level_populations_proton_data_toggle(ion): @pytest.mark.requires_dbase_version('>= 8') def test_contribution_function(ion): cont_func = ion.contribution_function(1e7 * u.cm**-3) - assert cont_func.shape == ion.temperature.shape + (1, ) + ion._wgfa['wavelength'].shape + wvl_bb = ion.transitions.wavelength[ion.transitions.is_bound_bound] + assert cont_func.shape == ion.temperature.shape + (1, ) + wvl_bb.shape # This value has not been tested for correctness assert u.allclose(cont_func[0, 0, 0], 2.51408088e-30 * u.cm**3 * u.erg / u.s) @@ -326,7 +327,8 @@ def test_two_ion_model_couple_density_temperature(fe20): @pytest.mark.requires_dbase_version('>= 8') def test_emissivity(ion): emm = ion.emissivity(1e7 * u.cm**-3) - assert emm.shape == ion.temperature.shape + (1, ) + ion._wgfa['wavelength'].shape + wvl_bb = ion.transitions.wavelength[ion.transitions.is_bound_bound] + assert emm.shape == ion.temperature.shape + (1, ) + wvl_bb.shape # This value has not been tested for correctness assert u.allclose(emm[0, 0, 0], 2.18000422e-16 * u.erg / u.cm**3 / u.s) @@ -338,7 +340,7 @@ def test_emissivity(ion): 1e29 * np.ones(temperature.shape) * u.cm**-5, ]) def test_intensity(ion, em): - wave_shape = ion._wgfa['wavelength'].shape + wave_shape = ion.transitions.wavelength[ion.transitions.is_bound_bound].shape intens = ion.intensity(1e7 * u.cm**-3, em) assert intens.shape == ion.temperature.shape + (1, ) + wave_shape # Test density varying along independent axis From d0e17a5224ea23b5b3c5b7c6a388fe03c99e8df7 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Mon, 2 Feb 2026 14:13:28 -0500 Subject: [PATCH 13/19] try to reduce memory footprint of goft example --- examples/user_guide/plot_contribution_function.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/user_guide/plot_contribution_function.py b/examples/user_guide/plot_contribution_function.py index 386c849c..634fe577 100644 --- a/examples/user_guide/plot_contribution_function.py +++ b/examples/user_guide/plot_contribution_function.py @@ -19,7 +19,7 @@ # Specify the plasma properties; note that an `~fiasco.Ion` has to be created with a # range of temperatures, but the density is only used later in the contribution # function calculation. -Te = np.geomspace(0.1, 100, 51) * u.MK +Te = np.geomspace(0.1, 10, 26) * u.MK ne = 1e8 * u.cm**-3 ############################################################################### From 89ef767049e93739fc3b2e57229fe99df174812e Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Mon, 2 Feb 2026 14:13:59 -0500 Subject: [PATCH 14/19] minor bugfix in goft indexing; docs clarifications --- fiasco/ions.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/fiasco/ions.py b/fiasco/ions.py index 97412f7d..778e466d 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -108,6 +108,12 @@ def __repr__(self): @cached_property @needs_dataset('elvlc') def levels(self): + """ + Information for each energy level of the ion. + + .. note:: For some ions, there may be more levels listed here than + given by `~fiasco.Ion.n_levels`. + """ return Levels(self._elvlc) def __getitem__(self, key): @@ -582,7 +588,9 @@ def level_populations(self, `~astropy.units.Quantity` A ``(l, m, n)`` shaped quantity, where ``l`` is the number of temperatures, ``m`` is the number of densities, and ``n`` is the number of energy - levels. If ``couple_density_to_temperature=True``, then ``m=1`` and ``l`` + levels in the ion model. Note that ``n`` will always be the same as the `~fiasco.Ion.n_levels`, + but may be different than the number of levels returned by `~fiasco.Levels`. + If ``couple_density_to_temperature=True``, then ``m=1`` and ``l`` represents the number of temperatures and densities. """ if use_two_ion_model: @@ -1134,8 +1142,13 @@ def contribution_function(self, density: u.cm**(-3), **kwargs) -> u.cm**3 * u.er upper_level = self.transitions.upper_level[self.transitions.is_bound_bound] wavelength = self.transitions.wavelength[self.transitions.is_bound_bound] A = self.transitions.A[self.transitions.is_bound_bound] - energy = const.h * const.c / wavelength - i_upper = vectorize_where(self._elvlc['level'], upper_level) + energy = wavelength.to('erg', equivalencies=u.equivalencies.spectral()) + # NOTE: The first array below provides the correspondence between the last + # dimension of the level populations array and the energy level index of + # the model. The upper level of the transition is used to make this selection + # because the contribution function is proportional to the population of the + # level from which the transition is happening. + i_upper = vectorize_where(np.arange(1, self.n_levels+1), upper_level) g = term * populations[:, :, i_upper] * (A * energy) return g From e7044e3aa0686930b76272aab2b305faa46c451e Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Mon, 2 Feb 2026 14:32:24 -0500 Subject: [PATCH 15/19] more memory footprint reduction --- examples/user_guide/aia_response.py | 2 +- examples/user_guide/ion_collection.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/user_guide/aia_response.py b/examples/user_guide/aia_response.py index d391292e..e8265b49 100644 --- a/examples/user_guide/aia_response.py +++ b/examples/user_guide/aia_response.py @@ -42,7 +42,7 @@ ############################################################ # Next, we construct for the `~fiasco.Ion` object for Fe 18. Note # that by default we use the coronal abundances of :cite:t:`feldman_potential_1992`. -temperature = 10.**(np.arange(4.5, 8, 0.05)) * u.K +temperature = 10.**(np.arange(4.5, 8, 0.1)) * u.K fe18 = fiasco.Ion('Fe 18', temperature) ############################################################ diff --git a/examples/user_guide/ion_collection.py b/examples/user_guide/ion_collection.py index 538ecde1..b1e3479e 100644 --- a/examples/user_guide/ion_collection.py +++ b/examples/user_guide/ion_collection.py @@ -25,7 +25,7 @@ # # We will create `~fiasco.IonCollection` containing both O VI # and Fe XVIII -temperature = np.geomspace(1e4,1e8,100) * u.K +temperature = np.geomspace(10**5, 10**7.5, 30) * u.K fe18 = fiasco.Ion('Fe XVIII', temperature) o6 = fiasco.Ion('O VI', temperature) col = fiasco.IonCollection(fe18, o6) From 1f7625ab4f8c32112ef5115e2aef1c894684361f Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Tue, 10 Feb 2026 14:35:41 -0500 Subject: [PATCH 16/19] ignore two-ion effects in comparison calculation to reduce memory footprint on RtD --- examples/idl_comparisons/continuum_comparison.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/idl_comparisons/continuum_comparison.py b/examples/idl_comparisons/continuum_comparison.py index e7c02050..219d3df9 100644 --- a/examples/idl_comparisons/continuum_comparison.py +++ b/examples/idl_comparisons/continuum_comparison.py @@ -104,7 +104,9 @@ def plot_idl_comparison(wavelength, temperature, result_fiasco, result_idl): 'ionization_fraction': idl_result_twophoton['ionization_fraction']} 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'], idl_result_twophoton['density']).squeeze() +two_photon = all_ions.two_photon(idl_result_twophoton['wavelength'], + idl_result_twophoton['density'], + use_two_ion_model=False).squeeze() plot_idl_comparison(idl_result_twophoton['wavelength'], idl_result_twophoton['temperature'], two_photon, From 642acb3e9fe8c0666379c1c9bc70fa880952a5ce Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Tue, 10 Feb 2026 14:37:17 -0500 Subject: [PATCH 17/19] throw exception if there are missing entries in the array being searched --- fiasco/util/tools.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/fiasco/util/tools.py b/fiasco/util/tools.py index 0d446a1a..140b2ddc 100644 --- a/fiasco/util/tools.py +++ b/fiasco/util/tools.py @@ -29,6 +29,8 @@ def vectorize_where(x_1, x_2): """ x_1 = np.atleast_1d(x_1) x_2 = np.atleast_1d(x_2) + if not np.isin(x_2, x_1).all(): + raise ValueError('All elements of second input array must be found in first input array.') return np.array([np.where(x_1==x)[0] for x in x_2]).squeeze() From 2e16261e28e2afd755ee7a48f686221888d54e60 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Tue, 10 Feb 2026 14:37:47 -0500 Subject: [PATCH 18/19] be more careful about indexing energy level arrays --- fiasco/ions.py | 45 ++++++++++++++++++++------------------------- fiasco/levels.py | 42 +++++++++++++++++++++++++++++------------- 2 files changed, 49 insertions(+), 38 deletions(-) diff --git a/fiasco/ions.py b/fiasco/ions.py index 778e466d..810078d7 100644 --- a/fiasco/ions.py +++ b/fiasco/ions.py @@ -110,9 +110,6 @@ def __repr__(self): def levels(self): """ Information for each energy level of the ion. - - .. note:: For some ions, there may be more levels listed here than - given by `~fiasco.Ion.n_levels`. """ return Levels(self._elvlc) @@ -168,7 +165,9 @@ def n_levels(self): Number of energy levels in the atomic model. .. note:: It is possible this number will not match the number of levels - in the energy level file. + in `~fiasco.Ion.levels`. The number of levels in a model is + determined by the number of energy levels as well as the level + information available for radiative decays and collisions. """ # There is no atomic model if there are no energies, collisions, and decays # NOTE: This logic is here rather than using a decorator so that it returns @@ -177,7 +176,7 @@ def n_levels(self): self._has_dataset('wgfa'), self._has_dataset('scups')]): return 0 - n_elvlc = len(self.levels) + n_elvlc = self._elvlc['level'].max() n_wgfa = max(self._wgfa['lower_level'].max(), self._wgfa['upper_level'].max()) n_scups = max(self._scups['upper_level'].max(), self._scups['lower_level'].max()) # If there is autoionization data associated with this ion, ensure that the model @@ -429,7 +428,7 @@ def effective_collision_strength(self) -> u.dimensionless_unscaled: return upsilon.T @cached_property - @needs_dataset('elvlc', 'scups') + @needs_dataset('scups') @u.quantity_input def electron_collision_deexcitation_rate(self) -> u.cm**3 / u.s: r""" @@ -454,11 +453,12 @@ def electron_collision_deexcitation_rate(self) -> u.cm**3 / u.s: """ c = const.h**2 / (2. * np.pi * const.m_e)**(1.5) upsilon = self.effective_collision_strength - omega_upper = 2. * self._elvlc['J'][self._scups['upper_level'] - 1] + 1. + idx = vectorize_where(self.levels.level, self._scups['upper_level']) + omega_upper = 2. * self.levels.total_angular_momentum[idx] + 1. return c * upsilon / np.sqrt(self.thermal_energy[:, np.newaxis]) / omega_upper @cached_property - @needs_dataset('elvlc', 'scups') + @needs_dataset('scups') @u.quantity_input def electron_collision_excitation_rate(self) -> u.cm**3 / u.s: r""" @@ -483,8 +483,9 @@ def electron_collision_excitation_rate(self) -> u.cm**3 / u.s: -------- electron_collision_deexcitation_rate : De-excitation rate due to collisions """ - omega_upper = 2. * self._elvlc['J'][self._scups['upper_level'] - 1] + 1. - omega_lower = 2. * self._elvlc['J'][self._scups['lower_level'] - 1] + 1. + J = self.levels.total_angular_momentum + omega_upper = 2. * J[vectorize_where(self.levels.level, self._scups['upper_level'])] + 1. + omega_lower = 2. * J[vectorize_where(self.levels.level, self._scups['lower_level'])] + 1. kBTE = np.outer(1./self.thermal_energy, self._scups['delta_energy']) return omega_upper / omega_lower * self.electron_collision_deexcitation_rate * np.exp(-kBTE) @@ -517,7 +518,7 @@ def proton_collision_excitation_rate(self) -> u.cm**3 / u.s: return u.Quantity(np.where(ex_rate > 0., ex_rate, 0.), u.cm**3/u.s).T @cached_property - @needs_dataset('elvlc', 'psplups') + @needs_dataset('psplups') @u.quantity_input def proton_collision_deexcitation_rate(self) -> u.cm**3 / u.s: r""" @@ -541,8 +542,9 @@ def proton_collision_deexcitation_rate(self) -> u.cm**3 / u.s: proton_collision_excitation_rate : Excitation rate due to collisions with protons """ kBTE = np.outer(self.thermal_energy, 1.0 / self._psplups['delta_energy']) - omega_upper = 2. * self._elvlc['J'][self._psplups['upper_level'] - 1] + 1. - omega_lower = 2. * self._elvlc['J'][self._psplups['lower_level'] - 1] + 1. + J = self.levels.total_angular_momentum + omega_upper = 2. * J[vectorize_where(self.levels.level, self._psplups['upper_level'])] + 1. + omega_lower = 2. * J[vectorize_where(self.levels.level, self._psplups['lower_level'])] + 1. dex_rate = (omega_lower / omega_upper) * self.proton_collision_excitation_rate * np.exp(1. / kBTE) return dex_rate @@ -1969,7 +1971,6 @@ def _karzas_cross_section(self, photon_energy: u.erg, ionization_energy: u.erg, cross_section[np.where(photon_energy < ionization_energy)] = 0.*cross_section.unit return cross_section - @needs_dataset('elvlc') @u.quantity_input def two_photon(self, wavelength: u.angstrom, @@ -2008,22 +2009,16 @@ def two_photon(self, A_ji = self._hseq['A'] psi_norm = self._hseq['psi_norm'] x_interp, y_interp = self._hseq['y'], self._hseq['psi'] - config = '2s' # Get the index of the 2S1/2 state for H-like - J = 0.5 + label = '2s 2S1/2' elif self.helium_like: A_ji = self._heseq['A'] psi_norm = 1.0 * u.dimensionless_unscaled x_interp, y_interp = self._heseq['y'], self._heseq['psi'] - config = '1s.2s' # Get the index of the 1s2s 1S0 state for He-like: - J = 0 + label = '1s 2s 1S0' else: return u.Quantity(np.zeros(final_shape), 'erg cm^3 s^-1 Angstrom^-1') - level_index = np.where((self._elvlc['config'] == config) & (np.isclose(self._elvlc['J'], J)) )[0][0] - - E_obs = self._elvlc['E_obs'][level_index] - E_th = self._elvlc['E_th'][level_index] - E_2p = E_obs if E_obs > 0.0 else E_th - rest_wavelength = 1 / E_2p + level_index = np.where(self.levels.label==label) + rest_wavelength = self.levels.energy[level_index].to('AA', equivalencies=u.equivalencies.spectral()) # NOTE: Explicitly setting the boundary condition type here to match the behavior of the # IDL spline interpolation functions. See https://github.com/wtbarnes/fiasco/pull/297 for @@ -2039,7 +2034,7 @@ def two_photon(self, # for more details. kwargs.setdefault('include_protons', False) level_population = self.level_populations(electron_density, **kwargs) - level_population = level_population[..., level_index] + level_population = level_population[..., self.levels[level_index].level-1] if couple_density_to_temperature: electron_density = electron_density[:, np.newaxis] diff --git a/fiasco/levels.py b/fiasco/levels.py index 13be9436..a90499d0 100644 --- a/fiasco/levels.py +++ b/fiasco/levels.py @@ -3,6 +3,7 @@ """ import astropy.units as u import numpy as np + from fractions import Fraction from fiasco.util import vectorize_where @@ -37,7 +38,7 @@ class Levels: def __init__(self, elvlc, index=None): self._elvlc = elvlc - self._index = np.s_[:] if index is None else index + self._idx = np.s_[:] if index is None else index def __len__(self): try: @@ -48,7 +49,7 @@ def __len__(self): def __getitem__(self, index): # NOTE: This throws an IndexError to stop iteration _ = self.level[index] - return type(self)(self._elvlc, index=index) + return type(self)(self._elvlc, n_levels=self.level.max(), index=index) def __repr__(self): return f"""Level: {self.level} @@ -61,17 +62,26 @@ def __repr__(self): @property def level(self): "Index of each level." - return self._elvlc['level'][self._index] + return self._elvlc['level'][self._idx] @property def configuration(self): "Label denoting the electronic configuration." - return np.char.replace(self._elvlc['config'][self._index], ".", " ") + configuration = self._elvlc['config'][self._idx] + # NOTE: This conditional is necessary because np.char.replace does not + # handle empty arrays. + if configuration.size == 0: + return configuration + return np.char.replace( + configuration, + ".", + " " + ) @property def multiplicity(self): "Multiplicity, :math:`2S+1`" - return self._elvlc['multiplicity'][self._index] + return self._elvlc['multiplicity'][self._idx] @property @u.quantity_input @@ -82,13 +92,19 @@ def spin(self) -> u.dimensionless_unscaled: @property def total_angular_momentum(self): "Total angular momentum number :math:`J`." - return self._elvlc['J'][self._index] + return self._elvlc['J'][self._idx] @property def label(self): - "Label denoting level configuration, multiplicity, angular momentum label, and total angular momentum." - zipped = zip(self.configuration, self.multiplicity, self.orbital_angular_momentum_label, self.total_angular_momentum) - return np.array([f"{i} {j}{k}{str(Fraction(l.value))}" for i,j,k,l in zipped]) + """ + Label denoting level configuration, multiplicity, angular momentum label, + and total angular momentum. + """ + zipped = zip(self.configuration, + self.multiplicity, + self.orbital_angular_momentum_label, + self.total_angular_momentum) + return np.array([f"{i} {j}{k}{Fraction(l.value)}" for i,j,k,l in zipped]) @property def weight(self): @@ -98,7 +114,7 @@ def weight(self): @property def orbital_angular_momentum_label(self): "Orbital angular momentum label." - return self._elvlc['L_label'][self._index] + return self._elvlc['L_label'][self._idx] @property @u.quantity_input @@ -112,7 +128,7 @@ def orbital_angular_momentum(self) -> u.dimensionless_unscaled: @property def is_observed(self): "True if the energy of the level is from laboratory measurements." - return self._elvlc['E_obs'][self._index].to_value('cm-1') != -1 + return self._elvlc['E_obs'][self._idx].to_value('cm-1') != -1 @property @u.quantity_input @@ -122,8 +138,8 @@ def energy(self) -> u.eV: theoretical energy if no measured energy is available. """ energy = np.where(self.is_observed, - self._elvlc['E_obs'][self._index], - self._elvlc['E_th'][self._index]) + self._elvlc['E_obs'][self._idx], + self._elvlc['E_th'][self._idx]) return energy.to('eV', equivalencies=u.equivalencies.spectral()) From c5aabf525fb4cfd8e2b3c80813e89e3cdb7ef0f6 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Tue, 10 Feb 2026 14:55:21 -0500 Subject: [PATCH 19/19] levels indexing bugfix --- fiasco/levels.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fiasco/levels.py b/fiasco/levels.py index a90499d0..5250f868 100644 --- a/fiasco/levels.py +++ b/fiasco/levels.py @@ -49,7 +49,7 @@ def __len__(self): def __getitem__(self, index): # NOTE: This throws an IndexError to stop iteration _ = self.level[index] - return type(self)(self._elvlc, n_levels=self.level.max(), index=index) + return type(self)(self._elvlc, index=index) def __repr__(self): return f"""Level: {self.level}