diff --git a/tests/_strategies.py b/tests/_strategies.py
index 3039908..f993b5a 100644
--- a/tests/_strategies.py
+++ b/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()
 
diff --git a/tests/test_source_functions.py b/tests/test_source_functions.py
deleted file mode 100644
index d0315f2..0000000
--- a/tests/test_source_functions.py
+++ /dev/null
@@ -1,24 +0,0 @@
-import numpy as np
-import pytest
-
-from zodipy._source_funcs import get_blackbody_emission
-
-TEMPERATURE = np.array([30])
-TEMPERATURE_ARRAY = np.array([31, 45, 53])
-R = 3
-R_ARRAY = np.array([4, 5.3, 6])
-DELTA = 0.324
-FREQUENCY = 549e9
-
-
-def test_blackbody_emission_value() -> None:
-    """Tests that return value."""
-    emission = get_blackbody_emission(T=TEMPERATURE, freq=FREQUENCY)
-    assert emission == pytest.approx(1.73442848898e-15, abs=1e-20)
-
-
-def test_blackbody_emission_value_array() -> None:
-    """Tests the return value given a temperature array."""
-    emission = get_blackbody_emission(T=TEMPERATURE_ARRAY, freq=FREQUENCY)
-    true_values = np.array([1.82147825366e-15, 3.06550295038e-15, 3.78860400626e-15])
-    assert emission == pytest.approx(true_values, abs=1e-20)
diff --git a/zodipy/_bandpass.py b/zodipy/_bandpass.py
index fba3428..696cb48 100644
--- a/zodipy/_bandpass.py
+++ b/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)]
+    )
diff --git a/zodipy/_constants.py b/zodipy/_constants.py
index d513e49..d01d265 100644
--- a/zodipy/_constants.py
+++ b/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
diff --git a/zodipy/_source_funcs.py b/zodipy/_source_funcs.py
index e6ebb5d..9b04809 100644
--- a/zodipy/_source_funcs.py
+++ b/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]: