Use built in astropy.models.physical_models.BlackBody

tests/_strategies.py

@@ -30,7 +30,7 @@
 from zodipy.model_registry import model_registry
 
 MIN_FREQ = u.Quantity(10, u.GHz)
-MAX_FREQ = u.Quantity(0.1, u.micron).to(u.GHz, equivalencies=u.spectral())
+MAX_FREQ = u.Quantity(0.6, u.micron).to(u.GHz, equivalencies=u.spectral())
 N_FREQS = 1000
 FREQ_LOG_RANGE = np.geomspace(np.log(MIN_FREQ.value), np.log(MAX_FREQ.value), N_FREQS).tolist()
 

zodipy/_bandpass.py

@@ -5,13 +5,14 @@
 
 import numpy as np
 from astropy import units
+from astropy.modeling.physical_models import BlackBody
 
 from zodipy._constants import (
     MAX_INTERPOLATION_GRID_TEMPERATURE,
     MIN_INTERPOLATION_GRID_TEMPERATURE,
     N_INTERPOLATION_POINTS,
+    SPECIFIC_INTENSITY_UNITS,
 )
-from zodipy._source_funcs import get_blackbody_emission
 from zodipy._validators import get_validated_and_normalized_weights, get_validated_freq
 
 if TYPE_CHECKING:
@@ -40,6 +41,11 @@ def switch_convention(self) -> None:
             self.weights = np.flip(self.weights)
             self.weights /= np.trapz(self.weights, self.frequencies.value)
 
+    @property
+    def is_center_frequency(self) -> bool:
+        """Return True if the bandpass is centered around a single frequency."""
+        return self.frequencies.isscalar
+
 
 def validate_and_get_bandpass(
     freq: units.Quantity,
@@ -64,14 +70,17 @@ def get_bandpass_interpolation_table(
     if not bandpass.frequencies.unit.is_equivalent(units.Hz):
         bandpass.switch_convention()
 
-    bandpass_integrals = np.zeros(n_points)
-    temperature_grid = np.linspace(min_temp, max_temp, n_points)
-    for idx, temp in enumerate(temperature_grid):
-        blackbody_emission = get_blackbody_emission(bandpass.frequencies.value, temp)
-        bandpass_integrals[idx] = (
-            bandpass.integrate(blackbody_emission)
-            if bandpass.frequencies.size > 1
-            else blackbody_emission
-        )
+    temperature_grid = np.linspace(min_temp, max_temp, n_points) * units.K
+    blackbody = BlackBody(temperature_grid)
+
+    if bandpass.is_center_frequency:
+        blackbody_emission = blackbody(bandpass.frequencies)
+        return np.asarray([temperature_grid, blackbody_emission.to_value(SPECIFIC_INTENSITY_UNITS)])
 
-    return np.asarray([temperature_grid, bandpass_integrals])
+    blackbody_emission = blackbody(bandpass.frequencies[:, np.newaxis])
+    integrated_blackbody_emission = np.trapz(
+        (bandpass.weights / (1 * units.Hz)) * blackbody_emission.transpose(), bandpass.frequencies
+    )
+    return np.asarray(
+        [temperature_grid, integrated_blackbody_emission.to_value(SPECIFIC_INTENSITY_UNITS)]
+    )

zodipy/_constants.py

@@ -1,11 +1,6 @@
-import astropy.constants as const
 import astropy.units as u
 import numpy as np
 
-h: float = const.h.value
-c: float = const.c.value
-k_B: float = const.k_B.value
-
 SPECIFIC_INTENSITY_UNITS = u.W / u.Hz / u.m**2 / u.sr
 
 R_0 = np.finfo(np.float64).eps

zodipy/_source_funcs.py

@@ -5,30 +5,10 @@
 
 import numpy as np
 
-from ._constants import c, h, k_B
-
 if TYPE_CHECKING:
     import numpy.typing as npt
 
 
-def get_blackbody_emission(
-    freq: float | npt.NDArray[np.float64], T: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    """Return the blackbody emission given a sequence of frequencies and temperatures.
-
-    Args:
-        freq: Frequency [Hz].
-        T: Temperature of the blackbody [K].
-
-    Returns:
-        Blackbody emission [W / m^2 Hz sr].
-
-    """
-    term1 = (2 * h * freq**3) / c**2
-    term2 = np.expm1((h * freq) / (k_B * T))
-    return term1 / term2
-
-
 def get_dust_grain_temperature(
     R: npt.NDArray[np.float64], T_0: float, delta: float
 ) -> npt.NDArray[np.float64]: