vald broadening implementation (#225)

* vald broadening implementation * apply black * apply black (again) * use vald broadening by default * slight capitalization formatting in config schema * fix inconsistent broadcasting error * fix equation bug * hwhm to fwhm * add molecules * remove redunant code * fix stark broadening mask * refactor broadening calc * documentation * use atomic number rather than ion number * use more accurate temperature scaling for vdW gamma * fix broadcasting bug * fix carsus bug * force int with loc to avoid warning * remove leftover accidental self-assignment * more comprehensive dtype object buigfix * add autoionization lines if using vald broadening * apply black (though we should use ruff instead) * fix bug with using broadening without a linelist * molecular broadening bugfix * logger messages * more logger debugging * debug logging * debug logging * more debugging * try memory efficient arrays for delta_nus * remove old logger messages * black * add comments * remove 4pi comment
tardis-sn · Dec 10, 2024 · d0b5404 · d0b5404
1 parent e3a49fe
commit d0b5404
Show file tree

Hide file tree

Showing 5 changed files with 348 additions and 41 deletions.
diff --git a/stardis/config_schema.yml b/stardis/config_schema.yml
@@ -129,6 +129,9 @@ properties:
                             shortlist:
                                 type: boolean
                                 default: false
+                            use_vald_broadening:
+                                type: boolean
+                                default: true
                         description: Options for whether to use a vald linelist. Linelist must be included the atom_data and cannot be supplied separately. 
                     include_molecules:
                         type: boolean
@@ -143,7 +146,7 @@ properties:
         description: Number of angles to use in the simulation.
     result_options:
         type: object
-        additionalProperties: False
+        additionalProperties: false
         default: {}
         properties:
             return_model:
@@ -154,7 +157,7 @@ properties:
                 default: false
             return_radiation_field:
                 type: boolean
-                default: True
+                default: true
         description: Options for what to return from the simulation.
     required:
     - stardis_config_version

diff --git a/stardis/plasma/base.py b/stardis/plasma/base.py
@@ -311,9 +311,9 @@ 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
-        )  # see 1995A&AS..112..525P for appropriate units - may be off by a factor of 4pi
+        linelist["A_ul"] = 10 ** (
+            linelist["rad"]
+        )  # see 1995A&AS..112..525P for appropriate units
 
         # 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
@@ -449,14 +449,10 @@ 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
-        )  # 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]
+        linelist["A_ul"] = 10 ** (
+            linelist["rad"]
+        )  # see 1995A&AS..112..525P for appropriate units
+        return alphas, linelist
 
 
 # Properties that haven't been used in creating stellar plasma yet,

diff --git a/stardis/plasma/molecules.py b/stardis/plasma/molecules.py
@@ -288,9 +288,9 @@ 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
-        )  # see 1995A&AS..112..525P for appropriate units - may be off by a factor of 4pi
+        linelist["A_ul"] = 10 ** (
+            linelist["rad"]
+        )  # see 1995A&AS..112..525P for appropriate units
 
         return alphas, linelist
 
@@ -414,8 +414,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
-        )  # see 1995A&AS..112..525P for appropriate units - may be off by a factor of 4pi
+        linelist["A_ul"] = 10 ** (
+            linelist["rad"]
+        )  # see 1995A&AS..112..525P for appropriate units
 
         return alphas, linelist
diff --git a/stardis/radiation_field/opacities/opacities_solvers/base.py b/stardis/radiation_field/opacities/opacities_solvers/base.py
@@ -2,6 +2,7 @@
 import numpy as np
 from pathlib import Path
 import numba
+import logging
 
 from astropy import units as u, constants as const
 
@@ -32,6 +33,8 @@
 ).value
 RYDBERG_FREQUENCY = (const.c.cgs * const.Ryd.cgs).value
 
+logger = logging.getLogger(__name__)
+
 
 # Calculate opacity from any table specified by the user
 def calc_alpha_file(stellar_plasma, stellar_model, tracing_nus, opacity_source, fpath):
@@ -406,13 +409,30 @@ 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
+
     gammas, doppler_widths = calculate_broadening(
         lines_sorted_in_range,
         stellar_model,
         stellar_plasma,
         line_opacity_config.broadening,
-    )  # This can be further improved by only calculating the broadening for the lines that are within the range.
+        use_vald_broadening=line_opacity_config.vald_linelist.use_vald_broadening
+        and line_opacity_config.vald_linelist.use_linelist,  # don't try to use vald broadening if you don't use vald linelists at all
+    )
 
+    # This line is awful for large simulations. Has to solve wavelength points times number of lines. Can be optimized.
     delta_nus = tracing_nus.value - line_nus[:, np.newaxis]
 
     # If no broadening range, compute the contribution of every line at every frequency.
@@ -484,6 +504,7 @@ def calc_molecular_alpha_line_at_nu(
     lines_sorted_in_range = lines_sorted[
         lines_sorted.nu.between(tracing_nus.min(), tracing_nus.max())
     ]
+
     line_nus = lines_sorted_in_range.nu.to_numpy()
 
     alphas_and_nu = stellar_plasma.molecule_alpha_line_from_linelist
@@ -525,7 +546,6 @@ def calc_molecular_alpha_line_at_nu(
             raise ValueError(
                 "Broadening range must be in units of length or frequency."
             )
-
     if n_threads == 1:  # Single threaded
         alpha_line_at_nu = calc_alan_entries(
             stellar_model.no_of_depth_points,