Fix NaNs for cool atmospheres with H2+BF and H-FF opacities enabled. (#…

…219)

* fill value -> 0 for LinearNDInterpolator

* add warning if outside of interpolation range

* fix warning message for h minus FF

* slightly clearer warning

* use logger warning instead of raising a warning

* provide depth point indices in warning message

stardis/radiation_field/opacities/opacities_solvers/util.py

@@ -1,12 +1,16 @@
 import numpy as np
 import pandas as pd
+import logging
 
 from astropy import units as u, constants as const
 from tardis.util.base import species_string_to_tuple
 
 from scipy.interpolate import LinearNDInterpolator
 
 
+logger = logging.getLogger(__name__)
+
+
 def sigma_file(tracing_lambdas, temperatures, fpath, opacity_source=None):
     """
     Reads and interpolates a cross-section file.
@@ -46,12 +50,16 @@ def sigma_file(tracing_lambdas, temperatures, fpath, opacity_source=None):
         linear_interp_h2plus_bf = LinearNDInterpolator(
             np.vstack([file_waves_mesh.ravel(), file_temps_mesh.ravel()]).T,
             file_cross_sections.flatten(),
+            fill_value=0,
         )
         lambdas, temps = np.meshgrid(tracing_lambdas, temperatures)
         sigmas = (
             linear_interp_h2plus_bf(lambdas, temps) * 1e-18
         )  # Scaling from Stancil 1994 table
-
+        if np.any(sigmas == 0):
+            logger.warning(
+                f"Outside of interpolation range for H2+ BF cross-sections at depth points {np.unique(np.where(sigmas == 0)[0])}. Assuming 0 opacity from H-FF for these depth points."
+            )
     elif (
         opacity_source == "Hminus_ff"
     ):  # This section specifically ingests the Bell and Berrington 1987 h_minus_ff_B1987.dat table found in data.
@@ -68,6 +76,7 @@ def sigma_file(tracing_lambdas, temperatures, fpath, opacity_source=None):
         linear_interp_hminus_ff = LinearNDInterpolator(
             np.vstack([file_waves_mesh.ravel(), file_thetas_mesh.ravel()]).T,
             file_values.flatten(),
+            fill_value=0,
         )
         lambdas, thetas = np.meshgrid(tracing_lambdas, 5040 / temperatures)
         sigmas = (
@@ -76,6 +85,10 @@ def sigma_file(tracing_lambdas, temperatures, fpath, opacity_source=None):
             * const.k_B.cgs.value
             * temperatures[:, np.newaxis]
         )
+        if np.any(sigmas == 0):
+            logger.warning(
+                f"Outside of interpolation range for H- FF cross-sections at depth points {np.unique(np.where(sigmas == 0)[0])}. Assuming 0 opacity from H-FF for these depth points."
+            )
 
     elif (
         opacity_source == "Hminus_bf"