Fix bug with bandpasses. Add tests. Bump version number.

pyproject.toml

@@ -1,7 +1,7 @@
 [tool.poetry]
 name = "zodipy"
 homepage = "https://github.com/Cosmoglobe/zodipy"
-version = "0.7.4"
+version = "0.7.5"
 description = "Software for simulating zodiacal emission"
 authors = ["Metin San <[email protected]>"]
 readme = "README.md"

tests/_strategies.py

@@ -147,6 +147,10 @@ def random_freq(
     )
     if unit is not None:
         random_freq = random_freq.to(unit, u.spectral())
+    else:
+        micron = draw(booleans())
+        if micron:
+            random_freq = random_freq.to(u.micron, u.spectral())
 
     if bandpass:
         shape = draw(integers(min_value=2, max_value=100))

tests/test_get_emission.py

@@ -413,3 +413,68 @@ def test_weights(
         obs_time=time,
         obs=observer,
     )
+
+
+def test_custom_weights() -> None:
+    model = Zodipy()
+    time = Time("2020-01-01")
+    nside = 32
+    pix = np.arange(hp.nside2npix(nside))
+    central_freq = 25
+    sigma_freq = 3
+    freqs = (
+        np.linspace(central_freq - sigma_freq, central_freq + sigma_freq, 100)
+        * u.micron
+    )
+    weights = np.random.randn(len(freqs))
+    weights /= np.trapz(weights, freqs.value)
+
+    model.get_emission_pix(
+        freq=freqs,
+        weights=weights,
+        pixels=pix,
+        nside=nside,
+        obs_time=time,
+        obs="earth",
+    )
+
+
+def test_custom_obs_pos() -> None:
+    model = Zodipy()
+    time = Time("2020-01-01")
+    nside = 64
+    pix = np.arange(hp.nside2npix(nside))
+
+    model.get_emission_pix(
+        freq=234 * u.micron,
+        pixels=pix,
+        nside=nside,
+        obs_time=time,
+        obs_pos=[0.1, 0.2, 1] * u.AU,
+    )
+
+    model.get_emission_pix(
+        freq=234 * u.micron,
+        pixels=pix,
+        nside=nside,
+        obs_time=time,
+        obs_pos=[2, 0.1, 4] * u.AU,
+    )
+
+    with pytest.raises(TypeError):
+        model.get_emission_pix(
+            freq=234 * u.micron,
+            pixels=pix,
+            nside=nside,
+            obs_time=time,
+            obs_pos=[2, 0.1, 4],
+        )
+
+    with pytest.raises(u.UnitsError):
+        model.get_emission_pix(
+            freq=234 * u.micron,
+            pixels=pix,
+            nside=nside,
+            obs_time=time,
+            obs_pos=[2, 0.1, 4] * u.s,
+        )

tests/test_validators.py

@@ -1,6 +1,5 @@
 import astropy.units as u
 import numpy as np
-import numpy.typing as npt
 import pytest
 from hypothesis import given
 from hypothesis.strategies import DataObject, data
@@ -53,6 +52,12 @@ def test_validate_weights(freq: FrequencyOrWavelength, data: DataObject) -> None
     )
     assert np.trapz(bp_weights, freq.value) == pytest.approx(1.0)
 
+    with pytest.raises(ValueError):
+        validate_and_normalize_weights(
+            weights=None,
+            freq=BANDPASS_FREQUENCIES,
+        )
+
 
 def test_validate_weights_numbers() -> None:
     validate_and_normalize_weights(
@@ -64,10 +69,12 @@ def test_validate_weights_numbers() -> None:
             weights=np.array([1, 2, 3]),
             freq=BANDPASS_FREQUENCIES,
         )
+    with pytest.raises(ValueError):
         validate_and_normalize_weights(
             weights=BANDPASS_WEIGHTS[:10],
             freq=BANDPASS_FREQUENCIES[0],
         )
+    with pytest.raises(ValueError):
         validate_and_normalize_weights(
             weights=None,
             freq=BANDPASS_FREQUENCIES,

zodipy/__init__.py

@@ -7,7 +7,7 @@
 
 try:
     __version__ = pkg_resources.get_distribution(__name__).version
-except pkg_resources.DistributionNotFound:
+except pkg_resources.DistributionNotFound:  # pragma: no cover
     ...
 
 __all__ = (

zodipy/_interp.py

@@ -27,27 +27,30 @@ def interpolate_source_parameters(
     weights: npt.NDArray[np.float64],
 ) -> InterpolatedSourceParameters:
 
-    # Shift spectrum convention between frequency and wavelength and
-    # flip bandpass to make it strictly increasing for np.trapz.
+    # Shift spectrum convention between frequency and wavelength (if model units
+    # differes from input units) and flip bandpass to make it strictly increasing
+    # for np.trapz.
     if not freq.unit.is_equivalent(model.spectrum.unit):
-        freq_model_units = freq.to(model.spectrum.unit, u.spectral()).value
-        freq_model_units = np.flip(freq_model_units)
-        weights = np.flip(weights)
+        freq_value = freq.to(model.spectrum.unit, u.spectral()).value
+        if weights.size > 1:
+            freq_value = np.flip(freq_value)
+            weights = np.flip(weights)
+            weights = weights / np.trapz(weights, freq_value)
     else:
-        freq_model_units = freq.value
+        freq_value = freq.value
 
     interpolator = partial(interp1d, x=model.spectrum.value, fill_value="extrapolate")
     emissivities = np.asarray(
         [
-            interpolator(y=model.emissivities[comp_label])(freq_model_units)
+            interpolator(y=model.emissivities[comp_label])(freq_value)
             for comp_label in model.comps.keys()
         ]
     )
 
     if model.albedos is not None:
         albedos = np.asarray(
             [
-                interpolator(y=model.albedos[comp_label])(freq_model_units)
+                interpolator(y=model.albedos[comp_label])(freq_value)
                 for comp_label in model.comps.keys()
             ]
         )
@@ -56,18 +59,16 @@ def interpolate_source_parameters(
 
     if model.phase_coefficients is not None:
         phase_coefficients = interpolator(y=np.asarray(model.phase_coefficients))(
-            freq_model_units
+            freq_value
         )
 
     else:
         phase_coefficients = np.repeat(
-            np.zeros((3, 1)), repeats=freq_model_units.size, axis=-1
+            np.zeros((3, 1)), repeats=freq_value.size, axis=-1
         )
 
     if model.solar_irradiance is not None:
-        solar_irradiance = interpolator(y=model.solar_irradiance.value)(
-            freq_model_units
-        )
+        solar_irradiance = interpolator(y=model.solar_irradiance.value)(freq_value)
         solar_irradiance = (
             u.Quantity(solar_irradiance, model.solar_irradiance.unit)
             .to(SPECIFIC_INTENSITY_UNITS, equivalencies=u.spectral())
@@ -76,10 +77,10 @@ def interpolate_source_parameters(
     else:
         solar_irradiance = 0
 
-    if freq_model_units.size > 1:
-        bandpass_integrate = lambda quantity: partial(
-            np.trapz, x=freq_model_units, axis=-1
-        )(quantity * weights)
+    if freq_value.size > 1:
+        bandpass_integrate = lambda quantity: partial(np.trapz, x=freq_value, axis=-1)(
+            quantity * weights
+        )
 
         emissivities = bandpass_integrate(emissivities)
         albedos = bandpass_integrate(albedos)

zodipy/_line_of_sight.py

@@ -3,16 +3,15 @@
 import numpy as np
 import numpy.typing as npt
 
-EPS = float(np.finfo(float).eps)
-
 
 def get_line_of_sight_endpoints(
     cutoff: float,
     obs_pos: npt.NDArray[np.float64],
     unit_vectors: npt.NDArray[np.float64],
-) -> tuple[float, npt.NDArray[np.float64]]:
+) -> tuple[np.float64, npt.NDArray[np.float64]]:
     """Returns the start and stop positions along the line of sights."""
 
+    eps = np.finfo(float).eps
     x, y, z = obs_pos.flatten()
     r = np.sqrt(x**2 + y**2 + z**2)
     if cutoff < r:
@@ -25,7 +24,6 @@ def get_line_of_sight_endpoints(
     cos_lat = np.cos(lat)
     b = 2 * (x * cos_lat * np.cos(lon) + y * cos_lat * np.sin(lon))
     c = r**2 - cutoff**2
-
     q = -0.5 * b * (1 + np.sqrt(b**2 - 4 * c) / np.abs(b))
 
-    return EPS, np.maximum(q, c / q)
+    return eps, np.maximum(q, c / q)

zodipy/_sky_coords.py

@@ -7,7 +7,8 @@
 pointing to Earth from the Sun.
 """
 
-from __future__ import annotations
+
+from typing import Tuple, Union
 
 import astropy.units as u
 import numpy as np
@@ -19,9 +20,10 @@
 DISTANCE_TO_JUPITER = u.Quantity(5.2, u.AU)
 
 
+@u.quantity_input
 def get_obs_and_earth_positions(
-    obs: str, obs_time: Time, obs_pos: u.Quantity[u.AU] | None
-) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
+    obs: str, obs_time: Time, obs_pos: Union[u.Quantity[u.AU], None]
+) -> Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
     """Returns the position of the observer and the Earth in broadcastable shapes
     (3, `n_pointing`, `n_gauss_quad_degree`).
     """

zodipy/zodipy.py

@@ -405,7 +405,7 @@ def _compute_emission(
         normalized_weights = validate_and_normalize_weights(weights=weights, freq=freq)
 
         interpolated_source_params = interpolate_source_parameters(
-            model=self.ipd_model, freq=freq.copy(), weights=normalized_weights.copy()
+            model=self.ipd_model, freq=freq, weights=normalized_weights
         )
 
         observer_position, earth_position = get_obs_and_earth_positions(
@@ -427,15 +427,15 @@ def _compute_emission(
 
         # Prepare bandpass to be integrated in power units and in frequency convention.
         if not freq.unit.is_equivalent(u.Hz):
-            freq = freq.to(u.Hz, u.spectral())
-            freq_ndarray = np.flip(freq.value)
-            if freq_ndarray.size > 1:
+            freq_value = freq.to(u.Hz, u.spectral()).value
+            if freq_value.size > 1:
+                freq_value = np.flip(freq_value)
                 normalized_weights = np.flip(normalized_weights)
-                normalized_weights /= np.trapz(normalized_weights, freq_ndarray)
+                normalized_weights /= np.trapz(normalized_weights, freq_value)
         else:
-            freq_ndarray = freq.value
-        if freq_ndarray.size == 1:
-            freq_ndarray = np.expand_dims(freq_ndarray, axis=0)
+            freq_value = freq.value
+        if freq_value.size == 1:
+            freq_value = np.expand_dims(freq_value, axis=0)
 
         # Prepopulate ipd model and configuration specific arguments to the zodiacal emission
         # integrand function.
@@ -445,7 +445,7 @@ def _compute_emission(
             gauss_quad_degree=self.gauss_quad_degree,
             X_obs=observer_position,
             density_partials=density_partials,
-            freq=freq_ndarray,
+            freq=freq_value,
             weights=normalized_weights,
             T_0=self.ipd_model.T_0,
             delta=self.ipd_model.delta,
@@ -541,7 +541,7 @@ def __repr__(self) -> str:
 
 def _get_emission_at_step(
     r: npt.NDArray[np.float64],
-    start: float,
+    start: np.float64,
     stop: npt.NDArray[np.float64],
     gauss_quad_degree: int,
     X_obs: npt.NDArray[np.float64],
@@ -554,7 +554,7 @@ def _get_emission_at_step(
     emissivities: npt.NDArray[np.float64],
     albedos: npt.NDArray[np.float64],
     phase_coefficients: tuple[float, ...],
-    solar_irradiance: float,
+    solar_irradiance: np.float64,
 ) -> npt.NDArray[np.float64]:
     """Returns the zodiacal emission at a step along all lines of sight."""