add autoionization lines if using vald broadening

stardis/plasma/base.py

@@ -311,8 +311,8 @@ def calculate(
         linelist["level_energy_upper"] = ((linelist["e_up"].values * u.eV).cgs).value
 
         # Radiation broadening parameter is approximated as the einstein A coefficient. Vald parameters are in log scale.
-        linelist["A_ul"] = 10 ** (
-            linelist["rad"]
+        linelist["A_ul"] = (
+            10 ** (linelist["rad"])
         )  # see 1995A&AS..112..525P for appropriate units - may be off by a factor of 4pi
 
         # Need to remove autoionization lines - can't handle with current broadening treatment because can't calculate effective principal quantum number
@@ -449,15 +449,11 @@ def calculate(
         linelist["level_energy_upper"] = ((linelist["e_up"].values * u.eV).cgs).value
 
         # Radiation broadening parameter is approximated as the einstein A coefficient. Vald parameters are in log scale.
-        linelist["A_ul"] = 10 ** (
-            linelist["rad"]
+        linelist["A_ul"] = (
+            10 ** (linelist["rad"])
         )  # see 1995A&AS..112..525P for appropriate units - may be off by a factor of 4pi
 
-        # Need to remove autoionization lines - can't handle with current broadening treatment because can't calculate effective principal quantum number
-        valid_indices = linelist.level_energy_upper < linelist.ionization_energy
-
-        return alphas[valid_indices], linelist[valid_indices]
-
+        return alphas, linelist
 
 # Properties that haven't been used in creating stellar plasma yet,
 # might be useful in future ----------------------------------------------------

stardis/plasma/molecules.py

@@ -54,9 +54,7 @@ def calculate(self, ion_number_density, t_electrons, atomic_data):
         for (
             molecule,
             molecule_row,
-        ) in (
-            molecules_df.iterrows()
-        ):  # Loop over all molecules, calculate number densities using Barklem and Collet 2016 equilibrium constants - if a component ion does not exist in the plasma or is negative, assume no molecule
+        ) in molecules_df.iterrows():  # Loop over all molecules, calculate number densities using Barklem and Collet 2016 equilibrium constants - if a component ion does not exist in the plasma or is negative, assume no molecule
             if (molecule_row.Ion1_charge == -1) or (molecule_row.Ion2_charge == -1):
                 logger.warning(
                     f"Negative ionic molecules not currently supported. Assuming no {molecule}."
@@ -288,8 +286,8 @@ def calculate(
         linelist["level_energy_upper"] = ((linelist["e_up"].values * u.eV).cgs).value
 
         # Radiation broadening parameter is approximated as the einstein A coefficient. Vald parameters are in log scale.
-        linelist["A_ul"] = 10 ** (linelist["rad"]) / (
-            4 * np.pi
+        linelist["A_ul"] = (
+            10 ** (linelist["rad"])
         )  # see 1995A&AS..112..525P for appropriate units - may be off by a factor of 4pi
 
         return alphas, linelist
@@ -414,8 +412,8 @@ def calculate(
         linelist["level_energy_upper"] = ((linelist["e_up"].values * u.eV).cgs).value
 
         # Radiation broadening parameter is approximated as the einstein A coefficient. Vald parameters are in log scale.
-        linelist["A_ul"] = 10 ** (linelist["rad"]) / (
-            4 * np.pi
+        linelist["A_ul"] = (
+            10 ** (linelist["rad"])
         )  # see 1995A&AS..112..525P for appropriate units - may be off by a factor of 4pi
 
         return alphas, linelist

stardis/radiation_field/opacities/opacities_solvers/base.py

@@ -406,6 +406,16 @@ def calc_alpha_line_at_nu(
         .to_numpy()
     )
 
+    if not line_opacity_config.vald_linelist.use_vald_broadening:
+        autoionization_lines = (
+            lines_sorted_in_range.level_energy_upper
+            > lines_sorted_in_range.ionization_energy
+        ).values
+
+        lines_sorted_in_range = lines_sorted_in_range[~autoionization_lines].copy()
+        alphas_array = alphas_array[~autoionization_lines].copy()
+        line_nus = line_nus[~autoionization_lines].copy()
+
     lines_sorted_in_range = lines_sorted_in_range.apply(
         pd.to_numeric
     )  # weird bug cropped up with ion_number being an object instead of an int