Set fvth flux to 0 for input energies above CHIANTI energy grid (#178)

* set fvth flux to 0 for input energies above CHIANTI energy grid

* continuum func: warn if energy too big, not error

* tests for zero flux above continuum grid max E

* thermal continuum: error if below lower bound

---------

Co-authored-by: William Setterberg <[email protected]>

changelog/178.feature.rst

@@ -0,0 +1 @@
+For thermal emission in `~sunkit_spex.legacy.thermal.thermal_emission`, the output flux is set to 0 for input energies above CHIANTI energy grid and raises an error if input energies are below the energy grid.

sunkit_spex/legacy/tests/test_thermal_.py

@@ -627,3 +627,34 @@ def test_energy_out_of_range_warning():
         # Trigger a warning.
         output = thermal.line_emission(np.arange(3, 28, 0.5) * u.keV, 6 * u.MK, 1e44 / u.cm**3)  # noqa
         assert issubclass(w[0].category, UserWarning)
+
+
+def test_continuum_energy_out_of_range():
+    with pytest.raises(ValueError):
+        # Use an energy range that goes out of bounds
+        # on the lower end--should error
+        _ = thermal.continuum_emission(np.arange(0.1, 28, 0.5) * u.keV, 6 * u.MK, 1e44 / u.cm**3)
+    with warnings.catch_warnings(record=True) as w:
+        warnings.simplefilter("always")
+        # The continuum emission should only warn if we go out of
+        # bounds on the upper end.
+        _ = thermal.continuum_emission(np.arange(10, 1000, 0.5) * u.keV, 6 * u.MK, 1e44 / u.cm**3)
+        assert issubclass(w[0].category, UserWarning)
+
+
+def test_empty_flux_out_of_range():
+    """The CHIANTI grid covers ~1 to 300 keV, but the values greater than
+    of the grid max energy should all be zeros."""
+    energy_edges = np.geomspace(10, 800, num=1000) << u.keV
+    midpoints = energy_edges[:-1] + np.diff(energy_edges) / 2
+
+    temperature = 20 << u.MK
+    em = 1e49 << u.cm**-3
+
+    flux = thermal.thermal_emission(energy_edges, temperature, em)
+    # the continuum is the one we need to check
+    max_e = thermal.CONTINUUM_GRID["energy range keV"][1] << u.keV
+    should_be_zeros = midpoints >= max_e
+
+    true_zero = 0 * (hopefully_zero := flux[should_be_zeros])
+    np.testing.assert_allclose(true_zero, hopefully_zero)

sunkit_spex/legacy/thermal.py

@@ -218,7 +218,10 @@ def thermal_emission(
         min(CONTINUUM_GRID["energy range keV"][0], LINE_GRID["energy range keV"][0]),
         max(CONTINUUM_GRID["energy range keV"][1], LINE_GRID["energy range keV"][1]),
     )
-    _error_if_input_outside_valid_range(energy_edges_keV, energy_range, "energy", "keV")
+    # raise an error if energy_range_min() < energy_endges_keV[1]
+    _error_if_low_energy_input_outside_valid_range(energy_edges_keV, energy_range, "energy", "keV")
+    # warning if energy_range.max() > energy_endges_keV[1], outside this range flux=0
+    _warn_if_input_outside_valid_range(energy_edges_keV, energy_range, "energy", "keV")
     temp_range = (
         min(CONTINUUM_GRID["temperature range K"][0], LINE_GRID["temperature range K"][0]),
         max(CONTINUUM_GRID["temperature range K"][1], LINE_GRID["temperature range K"][1]),
@@ -256,7 +259,10 @@ def continuum_emission(
     {doc_string_params}"""
     # Convert inputs to known units and confirm they are within range.
     energy_edges_keV, temperature_K = _sanitize_inputs(energy_edges, temperature)
-    _error_if_input_outside_valid_range(energy_edges_keV, CONTINUUM_GRID["energy range keV"], "energy", "keV")
+    _error_if_low_energy_input_outside_valid_range(
+        energy_edges_keV, CONTINUUM_GRID["energy range keV"], "energy", "keV"
+    )
+    _warn_if_input_outside_valid_range(energy_edges_keV, CONTINUUM_GRID["energy range keV"], "energy", "keV")
     _error_if_input_outside_valid_range(temperature_K, CONTINUUM_GRID["temperature range K"], "temperature", "K")
     # Calculate abundances
     abundances = _calculate_abundances(abundance_type, relative_abundances)
@@ -738,6 +744,15 @@ def _error_if_input_outside_valid_range(input_values, grid_range, param_name, pa
         raise ValueError(message)
 
 
+def _error_if_low_energy_input_outside_valid_range(input_values, grid_range, param_name, param_unit):
+    if input_values.min() < grid_range[0]:
+        message = (
+            f"Lower bound of the input {param_name} must be within the range "
+            f"{grid_range[0]}--{grid_range[1]} {param_unit}. "
+        )
+        raise ValueError(message)
+
+
 def _warn_if_input_outside_valid_range(input_values, grid_range, param_name, param_unit):
     if input_values.min() < grid_range[0] or input_values.max() > grid_range[1]:
         message = (