From 79aa1695423fa7e8ee0b36465d804efcaa2e39b0 Mon Sep 17 00:00:00 2001 From: Metin San Date: Mon, 29 Apr 2024 23:13:10 +0200 Subject: [PATCH] Use built in `astropy.models.physical_models.BlackBody` --- tests/_strategies.py | 2 +- tests/test_source_functions.py | 24 ------------------------ zodipy/_bandpass.py | 31 ++++++++++++++++++++----------- zodipy/_constants.py | 5 ----- zodipy/_source_funcs.py | 20 -------------------- 5 files changed, 21 insertions(+), 61 deletions(-) delete mode 100644 tests/test_source_functions.py 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]: