Fix bug where the new precomputed bandpass table resultet in unit err…

…ors for bandpasses in GHz. Update documentation

docs/examples/get_bandpass_integrated_emission.py

@@ -8,16 +8,16 @@
 
 nside = 128
 
-center_freq = 25 * u.micron
-freqs = np.linspace(20, 30, 11) * u.micron
+center_freq = 800 * u.GHz
+freqs = np.linspace(750, 850, 11) * u.GHz
 weights = np.array([2, 3, 5, 9, 11, 11.5, 11, 9, 5, 3, 2])
 
 plt.plot(freqs, weights)
-plt.xlabel("Frequency [micron]")
+plt.xlabel("Frequency [GHz]")
 plt.ylabel("Weights")
 plt.savefig("../img/bandpass.png", dpi=300)
 
-model = Zodipy()
+model = Zodipy(model="planck18")
 
 emission_central_freq = model.get_binned_emission_pix(
     freq=center_freq,
@@ -40,17 +40,17 @@
     emission_central_freq,
     title=f"Center frequency",
     unit="MJy/sr",
-    norm="log",
     cmap="afmhot",
+    norm="log",
 )
 plt.savefig("../img/center_freq.png", dpi=300)
 
 hp.mollview(
     emission_bandpass_integrated,
     title="Bandpass integrated",
     unit="MJy/sr",
-    norm="log",
     cmap="afmhot",
+    norm="log",
 )
 plt.savefig("../img/bandpass_integrated.png", dpi=300)
 plt.show()

docs/examples/get_binned_emission.py

@@ -21,11 +21,9 @@
     binned_emission,
     title="Binned zodiacal emission at 857 GHz",
     unit="MJy/sr",
+    min=0,
     max=1,
-    norm="log",
-    coord="E",
     cmap="afmhot",
 )
-hp.graticule()
-# plt.savefig("../img/binned.png", dpi=300)
+plt.savefig("../img/binned.png", dpi=300)
 plt.show()

docs/examples/get_binned_emission_solar_cutoff.py

@@ -25,6 +25,5 @@
     coord="E",
     cmap="afmhot",
 )
-hp.graticule()
-# plt.savefig("../img/binned_solar_cutoff.png", dpi=300)
+plt.savefig("../img/binned_solar_cutoff.png", dpi=300)
 plt.show()

docs/examples/get_binned_gal_emission.py

@@ -13,20 +13,18 @@
     857 * u.GHz,
     pixels=np.arange(hp.nside2npix(nside)),
     nside=nside,
-    obs_time=Time("2022-06-14"),
+    obs_time=Time("2022-02-20"),
     obs="earth",
     coord_in="G",  # Coordinates of the input pointing
 )
 
 hp.mollview(
     binned_emission,
-    title="Binned zodiacal emission  at 857 GHz",
+    title="Binned zodiacal emission at 857 GHz",
     unit="MJy/sr",
-    coord="G",
-    max=1,
-    norm="log",
     cmap="afmhot",
+    min=0,
+    max=1,
 )
-hp.graticule(coord="E")
-# plt.savefig("../img/binned_gal.png", dpi=300)
+plt.savefig("../img/binned_gal.png", dpi=300)
 plt.show()

docs/examples/get_emission_ang.py

@@ -6,20 +6,21 @@
 
 model = Zodipy("dirbe")
 
-theta = np.linspace(0, np.pi, 10000) * u.rad
-phi = np.zeros_like(theta)
+latitudes = np.linspace(-90, 90, 10000) * u.deg
+longitudes = np.zeros_like(latitudes)
 
 emission = model.get_emission_ang(
     30 * u.micron,
-    theta=theta,
-    phi=phi,
+    theta=longitudes,
+    phi=latitudes,
+    lonlat=True,
     obs_time=Time("2022-06-14"),
     obs="earth",
 )
 
 
-plt.plot(theta, emission)
-plt.xlabel("Theta [rad]")
+plt.plot(latitudes, emission)
+plt.xlabel("Latitude [deg]")
 plt.ylabel("Emission [MJy/sr]")
 plt.savefig("../img/timestream.png", dpi=300)
 plt.show()

docs/examples/get_parallel_emission.py

@@ -6,14 +6,14 @@
 import numpy as np
 from astropy.time import Time
 
-import zodipy
+from zodipy import Zodipy
 
 nside = 256
 pixels = np.arange(hp.nside2npix(nside))
 obs_time = Time("2020-01-01")
 
-model = zodipy.Zodipy(parallel=False)
-model_parallel = zodipy.Zodipy()
+model = Zodipy()
+model_parallel = Zodipy(parallel=True)
 
 start = time.perf_counter()
 emission = model.get_binned_emission_pix(
@@ -23,7 +23,7 @@
     obs_time=obs_time,
 )
 print("Time spent on a single CPU:", round(time.perf_counter() - start, 2), "seconds")
-# > Time spent on a single CPU: 91.76 seconds
+# > Time spent on a single CPU: 35.23 seconds
 
 start = time.perf_counter()
 emission_parallel = model_parallel.get_binned_emission_pix(
@@ -37,6 +37,6 @@
     round(time.perf_counter() - start, 2),
     "seconds",
 )
-# > Time spent on 8 CPUs: 26.87 seconds
+# > Time spent on 8 CPUs: 12.85 seconds
 
 assert np.allclose(emission, emission_parallel)

docs/usage.md

@@ -36,6 +36,19 @@ as seen by an observer on earth on 14 June, 2022 given the Planck 2018 interplan
 ![Zodiacal emission map](img/binned.png)
 *Note that the color bar is logarithmic.*
 
+
+### Bandpass integrated emission
+Instruments do not typically observe at delta frequencies. Usually, we are more interested in finding out
+what the emission looks like over some instrument bandpass. ZodiPy will accept a sequence of frequencies to the `freq`
+argument in addition to the corresponding bandpass weights to the `weights` argument and perform bandpass integration. 
+Note that the bandpass weights must be in spectral radiance units (Jy/sr), even though the weights them self are unitless. A top hat bandpass is assumed if a sequence of frequencies are used without providing weights.
+```python hl_lines="11 12 13 31 32"
+{!examples/get_bandpass_integrated_emission.py!}
+```
+![Generated Bandpass](img/bandpass.png)
+![Center frequency emission](img/center_freq.png)
+![Bandpass integrated emission](img/bandpass_integrated.png)
+
 ### Solar cutoff angle
 Few experiments look directly in towards the Sun. We can initialize `Zodipy` with the `solar_cut` 
 argument to mask all input pointing that looks in towards the sun with an angular distance smaller 
@@ -47,14 +60,13 @@ than the `solar_cut` value.
 ![Zodiacal emission map](img/binned_solar_cutoff.png)
 
 
-### Instantaneous map in Galactic coordinates
-We can make the same map in galactic coordinates by specifying that the input pointing is in galactic coordinates.
+### Non-ecliptic coordinates
+We can specify the coordinate system of the input pointing with the `coord_in` keyword
 
 ```python hl_lines="18"
 {!examples/get_binned_gal_emission.py!}
 ```
 ![Zodiacal emission map galactic](img/binned_gal.png)
-*Note that the color bar is logarithmic.*
 
 
 ### Component-wise maps
@@ -70,19 +82,6 @@ read [Cosmoglobe: Simulating Zodiacal Emission with ZodiPy](https://arxiv.org/ab
 *Note that the color for the Cloud component is logarithmic, while the others are linear.*
 
 
-### Bandpass integrated emission
-Instruments do not typically observe at delta frequencies. Usually, we are more interested in finding out
-what the emission looks like over some instrument bandpass. ZodiPy will accept a sequence of frequencies to the `freq`
-argument in addition to the corresponding bandpass weights to the `weights` argument and perform bandpass integration. 
-Note that the bandpass weights must be in spectral radiance units (Jy/sr), even though the weights them self are unitless. A top hat bandpass is assumed if a sequence of frequencies are used without providing weights.
-```python hl_lines="32 33"
-{!examples/get_bandpass_integrated_emission.py!}
-```
-![Generated Bandpass](img/bandpass.png)
-![Center frequency emission](img/center_freq.png)
-![Bandpass integrated emission](img/bandpass_integrated.png)
-
-
 ## Parallelization
 If you are not using ZodiPy in an already parallelized environment **and** are working with large pointing sequences, setting `parallel=True` when initializing `Zodipy` will improve the performance. ZodiPy will then automatically distribute the pointing to all available CPU's, given by `multiprocessing.cpu_count()` or to `n_proc` if this argument is provided.
 

mkdocs.yml

@@ -48,6 +48,4 @@ plugins:
     default_handler: python
     handlers:
       options:
-        merge_init_into_class: true
-    watch:
-      - zodipy
+        merge_init_into_class: true

zodipy/_bandpass.py

@@ -62,6 +62,8 @@ def get_bandpass_interpolation_table(
     # Prepare bandpass to be integrated in power units and in frequency convention.
     if not bandpass.frequencies.unit.is_equivalent(u.Hz):
         bandpass.switch_convention()
+    else:
+        bandpass.frequencies = bandpass.frequencies.to(u.Hz)
 
     integrals = np.zeros(n_points)
     temp_grid = np.linspace(min_temp, max_temp, n_points)

zodipy/_interpolate_source.py

@@ -64,12 +64,12 @@ def get_source_parameters_kelsall_comp(
         )
 
     if model.solar_irradiance is not None:
-        solar_irradiance = interpolator(y=model.solar_irradiance.value)(
+        solar_irradiance = interpolator(y=model.solar_irradiance)(
             bandpass.frequencies.value
         )
-        solar_irradiance = u.Quantity(
-            solar_irradiance, model.solar_irradiance.unit
-        ).to_value(SPECIFIC_INTENSITY_UNITS, equivalencies=u.spectral())
+        solar_irradiance = u.Quantity(solar_irradiance, "MJy /sr").to_value(
+            SPECIFIC_INTENSITY_UNITS, equivalencies=u.spectral()
+        )
     else:
         solar_irradiance = 0
 

zodipy/_ipd_model.py

@@ -4,8 +4,6 @@
 from dataclasses import dataclass, field
 from typing import Mapping, Sequence
 
-import astropy.units as u
-
 from zodipy._ipd_comps import Component, ComponentLabel
 from zodipy._types import FrequencyOrWavelength
 
@@ -40,7 +38,7 @@ class Kelsall(InterplanetaryDustModel):
     delta: float
     emissivities: Mapping[ComponentLabel, Sequence[float]]
     albedos: Mapping[ComponentLabel, Sequence[float]] | None = None
-    solar_irradiance: u.Quantity[u.MJy / u.sr] | None = None
+    solar_irradiance: Sequence[float] | None = None  # In units of MJy/sr
     phase_coefficients: Sequence[Sequence[float]] | None = None
 
 

zodipy/source_params.py

@@ -231,4 +231,4 @@
     82999.336,
     42346.605,
     14409.608,
-] * (u.MJy / u.sr)
+]

zodipy/zodipy.py

@@ -90,7 +90,7 @@ def __init__(
         self.ephemeris = ephemeris
 
         self.ipd_model = model_registry.get_model(model)
-        self.gauss_points, self.gauss_weights = np.polynomial.legendre.leggauss(
+        self.gauss_points_and_weights = np.polynomial.legendre.leggauss(
             gauss_quad_degree
         )
 
@@ -490,11 +490,7 @@ def _compute_emission(
                     proc_chunks = [
                         pool.apply_async(
                             _integrate_gauss_quad,
-                            args=(
-                                comp_integrand,
-                                self.gauss_points,
-                                self.gauss_weights,
-                            ),
+                            args=(comp_integrand, self.gauss_points_and_weights),
                         )
                         for comp_integrand in comp_integrands
                     ]
@@ -522,11 +518,7 @@ def _compute_emission(
                 )
 
                 integrated_comp_emission[idx] = (
-                    _integrate_gauss_quad(
-                        fn=comp_integrand,
-                        points=self.gauss_points,
-                        weights=self.gauss_weights,
-                    )
+                    _integrate_gauss_quad(comp_integrand, self.gauss_points_and_weights)
                     * 0.5
                     * (stop[comp_label] - start[comp_label])
                 )
@@ -572,8 +564,7 @@ def __repr__(self) -> str:
 
 def _integrate_gauss_quad(
     fn: Callable[[float], npt.NDArray[np.float64]],
-    points: npt.NDArray[np.float64],
-    weights: npt.NDArray[np.float64],
+    points_and_weights: tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]],
 ) -> npt.NDArray[np.float64]:
     """Integrate the emission from a component using Gauss-Legendre quadrature."""
-    return np.squeeze(sum(fn(point) * weight for point, weight in zip(points, weights)))
+    return np.squeeze(sum(fn(x) * w for x, w in zip(*points_and_weights)))