From 3974df77e90e1892f74d03e624f4aba32fceb577 Mon Sep 17 00:00:00 2001
From: Joshua Shields <54691495+jvshields@users.noreply.github.com>
Date: Mon, 14 Oct 2024 15:00:02 -0400
Subject: [PATCH] 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
---
 .../opacities/opacities_solvers/util.py           | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/stardis/radiation_field/opacities/opacities_solvers/util.py b/stardis/radiation_field/opacities/opacities_solvers/util.py
index 7e297e80..0d83b2cb 100644
--- a/stardis/radiation_field/opacities/opacities_solvers/util.py
+++ b/stardis/radiation_field/opacities/opacities_solvers/util.py
@@ -1,5 +1,6 @@
 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
@@ -7,6 +8,9 @@
 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"