From a46cb20b7432b6d243979d80c07637afee2f2a1b Mon Sep 17 00:00:00 2001 From: Metin San Date: Sun, 28 Apr 2024 09:50:45 +0200 Subject: [PATCH] Add support for specifying healpix order --- tests/_strategies.py | 7 +- tests/test_get_emission.py | 190 ++++++++++++++++++++++++++----------- zodipy/_line_of_sight.py | 15 +-- zodipy/_validators.py | 41 ++++---- zodipy/zodipy.py | 26 +++-- 5 files changed, 184 insertions(+), 95 deletions(-) diff --git a/tests/_strategies.py b/tests/_strategies.py index 8bdf8e0..318b342 100644 --- a/tests/_strategies.py +++ b/tests/_strategies.py @@ -3,7 +3,7 @@ import datetime from functools import partial from math import log2 -from typing import Any, Callable, Sequence +from typing import Any, Callable, Literal, Sequence import astropy.coordinates as coords import astropy.units as u @@ -71,8 +71,9 @@ def times(draw: DrawFn) -> time.Time: @composite -def nsides(draw: Callable[[SearchStrategy[int]], int]) -> int: - return draw(integers(min_value=MIN_NSIDE_EXP, max_value=MAX_NSIDE_EXP).map(partial(pow, 2))) +def healpixes(draw: DrawFn, order: Literal["ring", "nested"] = "ring") -> hp.HEALPix: + nside = draw(integers(min_value=MIN_NSIDE_EXP, max_value=MAX_NSIDE_EXP).map(partial(pow, 2))) + return hp.HEALPix(nside=nside, order=order) @composite diff --git a/tests/test_get_emission.py b/tests/test_get_emission.py index a19b9a4..e772b29 100644 --- a/tests/test_get_emission.py +++ b/tests/test_get_emission.py @@ -4,7 +4,7 @@ import astropy.coordinates as coords import astropy.units as u -import healpy as hp +import astropy_healpix as hp import numpy as np import pytest from astropy.time import Time, TimeDelta @@ -17,7 +17,7 @@ angles, coords_in, freqs, - nsides, + healpixes, obs_positions, pixels, quantities, @@ -50,12 +50,12 @@ def test_get_emission_skycoord( assert emission.size == coordinates.size -@given(zodipy_models(), sky_coords(), nsides(), data()) +@given(zodipy_models(), sky_coords(), healpixes(), data()) @settings(deadline=None) def test_get_binned_skycoord( model: Zodipy, coordinates: coords.SkyCoord, - nside: int, + healpix: hp.HEALPix, data: DataObject, ) -> None: """Property test for get_binned_emission_pix.""" @@ -67,28 +67,28 @@ def test_get_binned_skycoord( coordinates, freq=frequency, obs_pos=observer, - nside=nside, + nside=healpix.nside, solar_cut=data.draw(quantities(20, 50, u.deg)) if cut_solar else None, ) - assert emission_binned.shape == (hp.nside2npix(nside),) + assert emission_binned.shape == (healpix.npix,) -@given(zodipy_models(), times(), nsides(), data(), coords_in()) +@given(zodipy_models(), times(), healpixes(), data(), coords_in()) @settings(deadline=None) def test_get_emission_pix( model: Zodipy, time: Time, - nside: int, + healpix: hp.HEALPix, data: DataObject, coord_in: Literal["E", "G", "C"], ) -> None: """Property test for get_emission_pix.""" observer = data.draw(obs_positions(model, time)) frequency = data.draw(freqs(model)) - pix = data.draw(pixels(nside)) + pix = data.draw(pixels(healpix.nside)) emission = model.get_emission_pix( pix, - nside=nside, + nside=healpix.nside, freq=frequency, obs_time=time, obs_pos=observer, @@ -97,30 +97,30 @@ def test_get_emission_pix( assert np.size(emission) == np.size(pix) -@given(zodipy_models(), times(), nsides(), data(), coords_in()) +@given(zodipy_models(), times(), healpixes(), data(), coords_in()) @settings(deadline=None) def test_get_binned_emission_pix( model: Zodipy, time: Time, - nside: int, + healpix: hp.HEALPix, data: DataObject, coord_in: Literal["E", "G", "C"], ) -> None: """Property test for get_binned_emission_pix.""" observer = data.draw(obs_positions(model, time)) frequency = data.draw(freqs(model)) - pix = data.draw(pixels(nside)) + pix = data.draw(pixels(healpix.nside)) cut_solar = data.draw(booleans()) emission_binned = model.get_binned_emission_pix( pix, - nside=nside, + nside=healpix.nside, freq=frequency, obs_time=time, obs_pos=observer, solar_cut=data.draw(quantities(20, 50, u.deg)) if cut_solar else None, coord_in=coord_in, ) - assert emission_binned.shape == (hp.nside2npix(nside),) + assert emission_binned.shape == (healpix.npix,) @given(zodipy_models(), times(), angles(), data(), coords_in()) @@ -149,12 +149,12 @@ def test_get_emission_ang( assert emission.size == theta.size == phi.size -@given(zodipy_models(), times(), nsides(), angles(), data(), coords_in()) +@given(zodipy_models(), times(), healpixes(), angles(), data(), coords_in()) @settings(deadline=None) def test_get_binned_emission_ang( model: Zodipy, time: Time, - nside: int, + healpix: hp.HEALPix, angles: tuple[u.Quantity[u.deg], u.Quantity[u.deg]], data: DataObject, coord_in: Literal["E", "G", "C"], @@ -169,29 +169,29 @@ def test_get_binned_emission_ang( emission_binned = model.get_binned_emission_ang( theta, phi, - nside=nside, + nside=healpix.nside, freq=frequency, obs_time=time, obs_pos=observer, coord_in=coord_in, solar_cut=data.draw(quantities(20, 50, u.deg)) if cut_solar else None, ) - assert emission_binned.shape == (hp.nside2npix(nside),) + assert emission_binned.shape == (healpix.npix,) -@given(zodipy_models(extrapolate=False), times(), angles(lonlat=True), nsides(), data()) +@given(zodipy_models(extrapolate=False), times(), angles(lonlat=True), healpixes(), data()) @settings(deadline=None) def test_invalid_freq( model: Zodipy, time: Time, angles: tuple[u.Quantity[u.deg], u.Quantity[u.deg]], - nside: int, + healpix: hp.HEALPix, data: DataObject, ) -> None: """Property test checking for unsupported spectral range.""" theta, phi = angles observer = data.draw(obs_positions(model, time)) - pix = data.draw(pixels(nside)) + pix = data.draw(pixels(healpix.nside)) freq = data.draw(random_freqs(unit=model._ipd_model.spectrum.unit)) if not (model._ipd_model.spectrum[0] <= freq <= model._ipd_model.spectrum[-1]): @@ -199,7 +199,7 @@ def test_invalid_freq( model.get_emission_pix( pix, freq=freq, - nside=nside, + nside=healpix.nside, obs_time=time, obs_pos=observer, ) @@ -217,7 +217,7 @@ def test_invalid_freq( with pytest.raises(ValueError): model.get_binned_emission_pix( pix, - nside=nside, + nside=healpix.nside, freq=freq, obs_time=time, obs_pos=observer, @@ -227,7 +227,7 @@ def test_invalid_freq( model.get_binned_emission_ang( theta, phi, - nside=nside, + nside=healpix.nside, freq=freq, obs_time=time, obs_pos=observer, @@ -269,21 +269,21 @@ def test_multiprocessing() -> None: observer = "earth" time = Time("2020-01-01") frequency = 78 * u.micron - nside = 32 - pix = np.random.randint(0, hp.nside2npix(nside), size=1000) + healpix = hp.HEALPix(32) + pix = np.random.randint(0, healpix.npix, size=1000) theta = np.random.uniform(0, 180, size=1000) * u.deg phi = np.random.uniform(0, 360, size=1000) * u.deg emission_pix = model.get_emission_pix( pix, - nside=nside, + nside=healpix.nside, freq=frequency, obs_time=time, obs_pos=observer, ) emission_pix_parallel = model_parallel.get_emission_pix( pix, - nside=nside, + nside=healpix.nside, freq=frequency, obs_time=time, obs_pos=observer, @@ -292,14 +292,14 @@ def test_multiprocessing() -> None: emission_binned_pix = model.get_binned_emission_pix( pix, - nside=nside, + nside=healpix.nside, freq=frequency, obs_time=time, obs_pos=observer, ) emission_binned_pix_parallel = model_parallel.get_binned_emission_pix( pix, - nside=nside, + nside=healpix.nside, freq=frequency, obs_time=time, obs_pos=observer, @@ -326,7 +326,7 @@ def test_multiprocessing() -> None: theta, phi, freq=frequency, - nside=nside, + nside=healpix.nside, obs_time=time, obs_pos=observer, ) @@ -335,7 +335,7 @@ def test_multiprocessing() -> None: theta, phi, freq=frequency, - nside=nside, + nside=healpix.nside, obs_time=time, obs_pos=observer, ) @@ -351,19 +351,19 @@ def test_inner_radial_cutoff_multiprocessing() -> None: observer = "earth" time = Time("2020-01-01") frequency = 78 * u.micron - nside = 32 - pix = np.random.randint(0, hp.nside2npix(nside), size=1000) + healpix = hp.HEALPix(32) + pix = np.random.randint(0, healpix.npix, size=1000) emission_pix = model.get_emission_pix( pix, - nside=nside, + nside=healpix.nside, freq=frequency, obs_time=time, obs_pos=observer, ) emission_pix_parallel = model_parallel.get_emission_pix( pix, - nside=nside, + nside=healpix.nside, freq=frequency, obs_time=time, obs_pos=observer, @@ -374,7 +374,7 @@ def test_inner_radial_cutoff_multiprocessing() -> None: @given( zodipy_models(extrapolate=True), times(), - nsides(), + healpixes(), angles(), random_freqs(bandpass=True), data(), @@ -383,7 +383,7 @@ def test_inner_radial_cutoff_multiprocessing() -> None: def test_bandpass_integration( model: Zodipy, time: Time, - nside: int, + healpix: hp.HEALPix, angles: tuple[u.Quantity[u.deg], u.Quantity[u.deg]], freqs: u.Quantity[u.Hz] | u.Quantity[u.micron], data: DataObject, @@ -397,17 +397,17 @@ def test_bandpass_integration( phi, freq=freqs, weights=bp_weights, - nside=nside, + nside=healpix.nside, obs_time=time, obs_pos=observer, ) - assert emission_binned.shape == (hp.nside2npix(nside),) + assert emission_binned.shape == (healpix.npix,) @given( zodipy_models(extrapolate=True), times(), - nsides(), + healpixes(), angles(), random_freqs(bandpass=True), data(), @@ -416,7 +416,7 @@ def test_bandpass_integration( def test_weights( model: Zodipy, time: Time, - nside: int, + healpix: hp.HEALPix, angles: tuple[u.Quantity[u.deg], u.Quantity[u.deg]], freqs: u.Quantity[u.Hz] | u.Quantity[u.micron], data: DataObject, @@ -431,7 +431,7 @@ def test_weights( phi, weights=bp_weights, freq=freqs, - nside=nside, + nside=healpix.nside, obs_time=time, obs_pos=observer, ) @@ -441,8 +441,8 @@ def test_custom_weights() -> None: """Tests bandpass integration with custom weights.""" model = Zodipy() time = Time("2020-01-01") - nside = 32 - pix = np.arange(hp.nside2npix(nside)) + healpix = hp.HEALPix(64) + pix = np.arange(healpix.npix) central_freq = 25 sigma_freq = 3 freqs = np.linspace(central_freq - sigma_freq, central_freq + sigma_freq, 100) * u.micron @@ -453,7 +453,7 @@ def test_custom_weights() -> None: pix, freq=freqs, weights=weights, - nside=nside, + nside=healpix.nside, obs_time=time, obs_pos="earth", ) @@ -463,13 +463,13 @@ def test_custom_obs_pos() -> None: """Tests a user specified observer position.""" model = Zodipy() time = Time("2020-01-01") - nside = 64 - pix = np.arange(hp.nside2npix(nside)) + healpix = hp.HEALPix(64) + pix = np.arange(healpix.npix) model.get_emission_pix( pix, freq=234 * u.micron, - nside=nside, + nside=healpix.nside, obs_time=time, obs_pos=[0.1, 0.2, 1] * u.AU, ) @@ -477,7 +477,7 @@ def test_custom_obs_pos() -> None: model.get_emission_pix( pix, freq=234 * u.micron, - nside=nside, + nside=healpix.nside, obs_time=time, obs_pos=[2, 0.1, 4] * u.AU, ) @@ -486,7 +486,7 @@ def test_custom_obs_pos() -> None: model.get_emission_pix( pix, freq=234 * u.micron, - nside=nside, + nside=healpix.nside, obs_time=time, obs_pos=[2, 0.1, 4], ) @@ -495,7 +495,91 @@ def test_custom_obs_pos() -> None: model.get_emission_pix( pix, freq=234 * u.micron, - nside=nside, + nside=healpix.nside, obs_time=time, obs_pos=[2, 0.1, 4] * u.s, ) + + +def test_pixel_ordering() -> None: + """Tests that the pixel related functions work for both healpix orderings.""" + model = Zodipy() + time = Time("2020-01-01") + + healpix_ring = hp.HEALPix(32) + healpix_nested = hp.HEALPix(32) + + pix_ring = np.arange(healpix_ring.npix) + lon, lat = healpix_ring.healpix_to_lonlat(pix_ring) + lon = u.Quantity(lon) + lat = u.Quantity(lat) + pix_nested = healpix_nested.lonlat_to_healpix(lon, lat) + + model.get_emission_pix( + pix_ring, + freq=234 * u.micron, + nside=healpix_ring.nside, + obs_time=time, + order="ring", + ) + model.get_emission_pix( + pix_nested, + freq=234 * u.micron, + nside=healpix_nested.nside, + obs_time=time, + order="nested", + ) + + model.get_binned_emission_pix( + pix_ring, + freq=234 * u.micron, + nside=healpix_ring.nside, + obs_time=time, + order="ring", + ) + model.get_binned_emission_pix( + pix_nested, + freq=234 * u.micron, + nside=healpix_nested.nside, + obs_time=time, + order="nested", + ) + + model.get_binned_emission_ang( + lon, + lat, + lonlat=True, + freq=234 * u.micron, + nside=healpix_ring.nside, + obs_time=time, + order="ring", + ) + model.get_binned_emission_ang( + lon, + lat, + lonlat=True, + freq=234 * u.micron, + nside=healpix_nested.nside, + obs_time=time, + order="nested", + ) + + sky_coord = coords.SkyCoord( + lon, + lat, + frame=coords.BarycentricMeanEcliptic, + obstime=time, + ) + + model.get_binned_emission_skycoord( + sky_coord, + freq=234 * u.micron, + nside=healpix_ring.nside, + order="ring", + ) + model.get_binned_emission_skycoord( + sky_coord, + freq=234 * u.micron, + nside=healpix_nested.nside, + order="nested", + ) diff --git a/zodipy/_line_of_sight.py b/zodipy/_line_of_sight.py index 19efb27..92a46a0 100644 --- a/zodipy/_line_of_sight.py +++ b/zodipy/_line_of_sight.py @@ -11,9 +11,7 @@ if TYPE_CHECKING: import numpy.typing as npt - from zodipy._types import Float - -DIRBE_CUTOFFS: dict[ComponentLabel, tuple[Float, Float]] = { +DIRBE_CUTOFFS: dict[ComponentLabel, tuple[float | np.floating, float | np.floating]] = { ComponentLabel.CLOUD: (R_0, R_JUPITER), ComponentLabel.BAND1: (R_0, R_JUPITER), ComponentLabel.BAND2: (R_0, R_JUPITER), @@ -22,7 +20,7 @@ ComponentLabel.FEATURE: (R_EARTH - 0.2, R_EARTH + 0.2), } -RRM_CUTOFFS: dict[ComponentLabel, tuple[Float, Float]] = { +RRM_CUTOFFS: dict[ComponentLabel, tuple[float | np.floating, float | np.floating]] = { ComponentLabel.FAN: (R_0, RRM[ComponentLabel.FAN].R_outer), # type: ignore ComponentLabel.INNER_NARROW_BAND: ( RRM[ComponentLabel.INNER_NARROW_BAND].R_inner, # type: ignore @@ -54,12 +52,7 @@ def get_sphere_intersection( unit_vectors: npt.NDArray[np.float64], cutoff: float | np.float64, ) -> npt.NDArray[np.float64]: - """Get RMAX per pixel. - - Given the observer position, return distance from the observer to the - intersection between the line-of-sights and a heliocentric sphere with radius `cutoff`. - - """ + """Returns the distance from the observer to a heliocentric sphere with radius `cutoff`.""" x, y, z = obs_pos.flatten() r_obs = np.sqrt(x**2 + y**2 + z**2) @@ -87,7 +80,7 @@ def get_line_of_sight_range( dict[ComponentLabel, npt.NDArray[np.float64]], dict[ComponentLabel, npt.NDArray[np.float64]], ]: - """Get the start and stop distances for each component.""" + """Return the component-wise start- and end-point of each line of sight.""" start = { comp: get_sphere_intersection(obs_pos, unit_vectors, COMPONENT_CUTOFFS[comp][0]) for comp in components diff --git a/zodipy/_validators.py b/zodipy/_validators.py index 9e708b9..0461226 100644 --- a/zodipy/_validators.py +++ b/zodipy/_validators.py @@ -1,33 +1,32 @@ from __future__ import annotations -import astropy.units as u import numpy as np import numpy.typing as npt +from astropy import units from ._ipd_model import InterplanetaryDustModel -from ._types import FrequencyOrWavelength, SkyAngles def get_validated_freq( - freq: FrequencyOrWavelength, model: InterplanetaryDustModel, extrapolate: bool -) -> FrequencyOrWavelength: + freq: units.Quantity, model: InterplanetaryDustModel, extrapolate: bool +) -> units.Quantity: """Validate user inputted frequency.""" - if not isinstance(freq, u.Quantity): + if not isinstance(freq, units.Quantity): msg = "Frequency must be an astropy Quantity." raise TypeError(msg) - if freq.unit.is_equivalent(u.Hz): - freq = freq.to(u.Hz) - elif freq.unit.is_equivalent(u.micron): - freq = freq.to(u.micron) + if freq.unit.is_equivalent(units.Hz): + freq = freq.to(units.Hz) + elif freq.unit.is_equivalent(units.micron): + freq = freq.to(units.micron) else: msg = "Frequency must be in units compatible with Hz or micron." - raise u.UnitsError(msg) + raise units.UnitsError(msg) if extrapolate: return freq - freq_in_spectrum_units = freq.to(model.spectrum.unit, equivalencies=u.spectral()) + freq_in_spectrum_units = freq.to(model.spectrum.unit, equivalencies=units.spectral()) lower_freq_range = model.spectrum.min() upper_freq_range = model.spectrum.max() @@ -48,7 +47,7 @@ def get_validated_freq( def get_validated_and_normalized_weights( weights: npt.ArrayLike | None, - freq: FrequencyOrWavelength, + freq: units.Quantity, ) -> npt.NDArray[np.float64]: """Validate user inputted weights.""" if weights is None and freq.size > 1: @@ -74,18 +73,22 @@ def get_validated_and_normalized_weights( def get_validated_ang( - theta: SkyAngles, phi: SkyAngles, lonlat: bool -) -> tuple[SkyAngles, SkyAngles]: + theta: units.Quantity, phi: units.Quantity, lonlat: bool +) -> tuple[units.Quantity, units.Quantity]: """Validate user inputted sky angles and make sure it adheres to the healpy convention.""" - theta = theta.to(u.deg) if lonlat else theta.to(u.rad) - phi = phi.to(u.deg) if lonlat else phi.to(u.rad) + try: + theta = theta.to(units.deg) if lonlat else theta.to(units.rad) + phi = phi.to(units.deg) if lonlat else phi.to(units.rad) + except AttributeError: + msg = "Sky angles `theta` and `phi` must be astropy Quantities." + raise TypeError(msg) from AttributeError if theta.isscalar: - theta = u.Quantity([theta]) + theta = np.expand_dims(theta, axis=0) if phi.isscalar: - phi = u.Quantity([phi]) + phi = np.expand_dims(phi, axis=0) if not lonlat: - theta, phi = phi, (np.pi / 2) * u.rad - theta + theta, phi = phi, (np.pi / 2) * units.rad - theta return theta, phi diff --git a/zodipy/zodipy.py b/zodipy/zodipy.py index 3474fa4..b36c691 100644 --- a/zodipy/zodipy.py +++ b/zodipy/zodipy.py @@ -150,18 +150,18 @@ def get_emission_skycoord( return_inverse=True, axis=1, ) + obs_time = coord.obstime coord = coords.SkyCoord( unique_lon * units.deg, unique_lat * units.deg, frame=coord.frame, - obstime=coord.obstime, ) return self._compute_emission( freq=freq, weights=weights, obs_pos=obs_pos, - obs_time=coord.obstime, + obs_time=obs_time, coordinates=coord, indicies=indicies, return_comps=return_comps, @@ -238,6 +238,7 @@ def get_emission_pix( coord_in: Literal["E", "G", "C"] = "E", weights: npt.ArrayLike | None = None, return_comps: bool = False, + order: Literal["ring", "nested"] = "ring", ) -> units.Quantity[units.MJy / units.sr]: """Return the simulated zodiacal emission given pixel numbers. @@ -259,13 +260,14 @@ def get_emission_pix( weights: Bandpass weights corresponding the the frequencies in `freq`. The weights are assumed to be given in spectral radiance units (Jy/sr). return_comps: If True, the emission is returned component-wise. Defaults to False. + order: Order of the HEALPix grid. Returns: emission: Simulated zodiacal emission in units of 'MJy/sr'. """ frame = get_frame_from_string(coord_in) - healpix = hp.HEALPix(nside=nside, order="ring", frame=frame) + healpix = hp.HEALPix(nside=nside, order=order, frame=frame) unique_pixels, indicies = np.unique(pixels, return_inverse=True) coordinates = healpix.healpix_to_skycoord(unique_pixels) @@ -288,6 +290,7 @@ def get_binned_emission_skycoord( weights: npt.ArrayLike | None = None, obs_pos: units.Quantity | str = "earth", return_comps: bool = False, + order: Literal["ring", "nested"] = "ring", solar_cut: units.Quantity | None = None, ) -> units.Quantity[units.MJy / units.sr]: """Return the simulated binned zodiacal light for all observations in a `SkyCoord` object. @@ -309,6 +312,7 @@ def get_binned_emission_skycoord( an observer in the `astropy.coordinates.solar_system_ephemeris`. This should correspond to a single position. Defaults to 'earth'. return_comps: If True, the emission is returned component-wise. Defaults to False. + order: Order of the HEALPix grid. solar_cut: Cutoff angle from the sun. The emission for all the pointing with angular distance between the sun smaller than `solar_cut` are masked. Defaults to `None`. @@ -321,18 +325,18 @@ def get_binned_emission_skycoord( return_inverse=True, axis=1, ) + obs_time = coord.obstime coord = coords.SkyCoord( unique_lon * units.deg, unique_lat * units.deg, frame=coord.frame, - obstime=coord.obstime, ) - healpix = hp.HEALPix(nside, order="ring", frame=coord.frame) + healpix = hp.HEALPix(nside, order=order, frame=coord.frame) return self._compute_emission( freq=freq, weights=weights, obs_pos=obs_pos, - obs_time=coord.obstime, + obs_time=obs_time, coordinates=coord, indicies=indicies, healpix=healpix, @@ -353,6 +357,7 @@ def get_binned_emission_ang( coord_in: Literal["E", "G", "C"] = "E", weights: npt.ArrayLike | None = None, return_comps: bool = False, + order: Literal["ring", "nested"] = "ring", solar_cut: units.Quantity | None = None, ) -> units.Quantity[units.MJy / units.sr]: """Return the simulated binned zodiacal emission given angles on the sky. @@ -383,6 +388,7 @@ def get_binned_emission_ang( lonlat: If True, input angles (`theta`, `phi`) are assumed to be longitude and latitude, otherwise, they are co-latitude and longitude. return_comps: If True, the emission is returned component-wise. Defaults to False. + order: Order of the HEALPix grid. solar_cut: Cutoff angle from the sun. The emission for all the pointing with angular distance between the sun smaller than `solar_cut` are masked. Defaults to `None`. @@ -390,9 +396,9 @@ def get_binned_emission_ang( emission: Simulated zodiacal emission in units of 'MJy/sr'. """ - theta, phi = get_validated_ang(theta=theta, phi=phi, lonlat=lonlat) + theta, phi = get_validated_ang(theta, phi, lonlat=lonlat) frame = get_frame_from_string(coord_in) - healpix = hp.HEALPix(nside, order="ring", frame=frame) + healpix = hp.HEALPix(nside, order=order, frame=frame) (theta, phi), counts = np.unique(np.vstack([theta, phi]), return_counts=True, axis=1) coordinates = coords.SkyCoord( theta, @@ -423,6 +429,7 @@ def get_binned_emission_pix( coord_in: Literal["E", "G", "C"] = "E", weights: npt.ArrayLike | None = None, return_comps: bool = False, + order: Literal["ring", "nested"] = "ring", solar_cut: units.Quantity | None = None, ) -> units.Quantity[units.MJy / units.sr]: """Return the simulated binned zodiacal Emission given pixel numbers. @@ -446,6 +453,7 @@ def get_binned_emission_pix( weights: Bandpass weights corresponding the the frequencies in `freq`. The weights are assumed to be given in spectral radiance units (Jy/sr). return_comps: If True, the emission is returned component-wise. Defaults to False. + order: Order of the HEALPix grid. solar_cut: Cutoff angle from the sun. The emission for all the pointing with angular distance between the sun smaller than `solar_cut` are masked. Defaults to `None`. @@ -454,7 +462,7 @@ def get_binned_emission_pix( """ frame = get_frame_from_string(coord_in) - healpix = hp.HEALPix(nside=nside, order="ring", frame=frame) + healpix = hp.HEALPix(nside=nside, order=order, frame=frame) unique_pixels, counts = np.unique(pixels, return_counts=True) coordinates = healpix.healpix_to_skycoord(unique_pixels)