Usage

Breaking API changes in v1.0.0

In version v1.0.0 the get_emission_* and get_binned_emission_* methods were deprecated and removed. Users wanting to simulate zodiacal light directly from HEALPix pixel indices should see the HEALPix section under examples.

As an Astropy-affiliated package, ZodiPy is highly integrated with the Astropy ecosystem, particularly with the astropy.units, astropy.coordinates, and astropy.time modules.

Initializing a zodiacal light model

To make zodiacal light simulations we must first import and initialize a zodiacal light model

import astropy.units as u
import zodipy

model = zodipy.Model(25 * u.micron)

The zodipy.Model object has one required positional argument, x, which can either represent a center wavelength/frequency or the points of an empirical bandpass. If x represents the points of an instrument bandpass, the weights argument must also be provided

import astropy.units as u
import zodipy

points = [3, 4, 5, 6] * u.micron
weights = [0.2, 0.4, 0.3, 0.1]

model = zodipy.Model(points, weights=weights)

ZodiPy supports several zodiacal light models valid at different wavelength/frequency ranges. By default, a Model will initialize on the DIRBE model. To select another model, for instance the Planck 2013 model , we need to specify the name keyword in the model constructor

import astropy.units as u
import zodipy

model = zodipy.Model(25 * u.micron, name="planck13")

Other possible keyword arguments to the zodipy.Model are gauss_quad_degree, which determines the number of discrete points evaluated along each line-of-sight, extrapolate, which, if set to True, uses linear interpolation to extrapolate the frequency/wavelength dependent model parameters allowing the model to be evaluated outside of it's original bounds, and finally, ephemeris, which one can use to specify the solar system ephemeris. used to compute the position of the Earth and optionally the observer.

Evaluating a zodiacal light model

To make zodiacal light simulations ZodiPy needs some inputs data from the user:

1. Sky coordinates. Provided through Astropy's SkyCoord object.
2. Time of observeration(s). Also provided in the SkyCoord object.
3. Position(s) of the observer. Provided in a separate argument when evaluating the model.

The SkyCoord object

Users unfamiliar with Astropy's SkyCoord should first visit the official Astropy docs to learn the basics.

For a single observation in galactic coordinates, the SkyCoord object may look something like the following:

import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.time import Time

skycoord = SkyCoord(
    40 * u.deg, 
    60 * u.deg, 
    obstime=Time("2020-01-01"), 
)

where the coordinates here are specified as longitude and latitude values. Note that the SkyCoord object is very flexible and supports many different formats for the coordinates.

The obstime keyword is mandatory and is given by the Astropy Time object, which can represent time in many formats, including regular and modified Julian dates (see the Time docs for more information).

For several observations, each with their own obstime input, the SkyCoord may instead look like the following:

import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.time import Time

skycoord = SkyCoord(
    [40, 41, 42] * u.deg, 
    [60, 59, 58] * u.deg, 
    obstime=Time(["2020-01-01", "2020-01-02", "2020-01-03"]), 
)

If a single value is given for obstime, all coordinates are assumed to be viewed instantaneously at that time from a single position in the solar system. Otherwise, each coordinate must have its own obstime value.

The sky coordinates should represent observer-centric coordinates. The observer-position is therefore required to compute the line-of-sight integrals, but this is provided not in the SkyCoord object, but rather in the evaluate method, which we will see soon.

The coordinate frame in the SkyCoord object defaults to the ICRS frame, but can be changed by providing the frame keyword argument:

import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.time import Time

skycoord = SkyCoord(
    [40, 41, 42] * u.deg, 
    [60, 59, 58] * u.deg, 
    obstime=Time(["2020-01-01", "2020-01-02", "2020-01-03"]), 
    frame="galactic",
)

Common coordinate frames are the Ecliptic, Galactic, and Celestial frames (these are represented as "E", "G", "C" in healpy), which can be specified either through string representations:

• "barycentricmeanecliptic" (Ecliptic)
• "galactic" (Galactic)
• "icrs" (Celestial)

or through frame objects imported from astropy.coordinates:

• BarycentricMeanEcliptic
• Galactic
• ICRS

Notes on coordinate frames

Note that these built-in Astropy frames do not inherently represent observer-centric coordinate frames. However this is fine, since we only need to know the rotation of the coordinates with respect to the ecliptic plane (internally, the coordinates are manually shifted to heliocentric coordinates using the obspos value)-

The evaluate method

The zodiacal light is evaluated by providing the SkyCoord object to the evaluate method.

import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.time import Time
import zodipy

model = zodipy.Model(25 * u.micron)

skycoord = SkyCoord(
    40 * u.deg, 
    60 * u.deg, 
    obstime=Time("2020-01-01"), 
    frame="galactic"
)

emission = model.evaluate(skycoord)
print(emission)
# <Quantity [25.08189292] MJy / sr>

By default, the observer is assumed to be the center of the Earth. In this case the Earth position is computed internally using Astropy's solar system ephemeris. The position of the observer can be explicitly provided through the keyword argument obspos in the evaluate method

import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.time import Time
import zodipy

model = zodipy.Model(25 * u.micron)

skycoord = SkyCoord(
    40 * u.deg, 
    60 * u.deg, 
    obstime=Time("2020-01-01"), 
    frame="galactic"
)

emission = model.evaluate(skycoord, obspos="mars")
print(emission)
# <Quantity [8.36985535] MJy / sr>

emission = model.evaluate(skycoord, obspos=[0.87, -0.53, 0.001] * u.AU)
print(emission)
# <Quantity [20.37750965] MJy / sr>

This argument accepts both a string representing a body recognized by the solar system ephemeris, or a heliocentric ecliptic cartesian position.

Similar to with the obstime attribute in the SkyCoord object, the value provided to the obspos keyword must have shape (3,) or (3, ncoords), where ncoords is the number of coordinates in the SkyCoord object.

Multiprocessing

ZodiPy will distribute the input coordinates to available cores using Python's multiprocessing module if the keyword argument nprocesses in the evaluate method is greater or equal to 1.

import multiprocessing

nprocesses = multiprocessing.cpu_count() # 8 cores

emission = model.evaluate(skycoord, nprocesses=nprocesses)

Examples

Zodiacal light along an Ecliptic scan

In the following, we simulate a scan across the Ecliptic plane:

ecliptic_scan.py

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.coordinates import BarycentricMeanEcliptic, SkyCoord
from astropy.time import Time, TimeDelta

from zodipy import Model

model = Model(30 * u.micron)

lats = np.linspace(-90, 90, 100) * u.deg
lons = np.zeros_like(lats)

t0 = Time("2022-06-14")
dt = TimeDelta(1, format="sec")
obstimes = t0 + np.arange(lats.size) * dt

coords = SkyCoord(
    lons,
    lats,
    frame=BarycentricMeanEcliptic,
    obstime=obstimes,
)

emission = model.evaluate(coords)

plt.plot(lats, emission)
plt.xlabel("Latitude [deg]")
plt.ylabel("Emission [MJy/sr]")
plt.savefig("../img/ecliptic_scan.png", dpi=300, bbox_inches="tight")
plt.show()

HEALPix

We can use healpy or Astropy-healpix to create a SkyCoord object directly from a HEALPix pixelization. In the following two examples, we produce an instantaneous map of the sky in HEALPix representation:

healpyastropy-healpix

import multiprocessing

import astropy.units as u
import healpy as hp
import matplotlib.pyplot as plt
import numpy as np
from astropy.coordinates import SkyCoord
from astropy.time import Time

import zodipy

model = zodipy.Model(30 * u.micron)

nside = 256
pixels = np.arange(hp.nside2npix(nside))
lon, lat = hp.pix2ang(nside, pixels, lonlat=True)

skycoord = SkyCoord(lon, lat, unit=u.deg, frame="galactic", obstime=Time("2022-01-14"))

emission = model.evaluate(skycoord, nprocesses=multiprocessing.cpu_count())

hp.mollview(
    emission,
    unit="MJy/sr",
    cmap="afmhot",
    min=0,
    max=80,
    title="Zodiacal light at 30 µm (2022-01-14)",
)
plt.savefig("../img/healpix_map.png", dpi=300, bbox_inches="tight")
plt.show()

import multiprocessing

import astropy.units as u
import astropy_healpix as ahp
import healpy as hp
import matplotlib.pyplot as plt
import numpy as np
from astropy.time import Time

import zodipy

model = zodipy.Model(30 * u.micron)

healpix = ahp.HEALPix(nside=256, frame="galactic")
pixels = np.arange(healpix.npix)
skycoord = healpix.healpix_to_skycoord(pixels)

# Note that we manually set the obstime attribute
skycoord.obstime = Time("2022-01-14")

emission = model.evaluate(skycoord, nprocesses=multiprocessing.cpu_count())

# Plot with healpy
hp.mollview(
    emission,
    unit="MJy/sr",
    cmap="afmhot",
    min=0,
    max=80,
    title="Zodiacal light at 30 µm (2022-01-14)",
)
# plt.savefig("../img/healpix_map.png", dpi=300, bbox_inches="tight")
plt.show()

Component-wise zodiacal light

We can return the zodiacal light for each component by using setting the keyword argument return_comps to True in the evaluate method:

component_maps.py

import multiprocessing

import astropy.units as u
import healpy as hp
import matplotlib.pyplot as plt
import numpy as np
from astropy.coordinates import SkyCoord
from astropy.time import Time

import zodipy

COMP_NAMES = [
    "Smooth cloud",
    "Dust band 1",
    "Dust band 2",
    "Dust band 3",
    "Circum-solar Ring",
    "Earth-trailing Feature",
]

model = zodipy.Model(24 * u.micron)

nside = 128
pixels = np.arange(hp.nside2npix(nside))
lon, lat = hp.pix2ang(nside, pixels, lonlat=True)

skycoord = SkyCoord(
    lon,
    lat,
    unit=u.deg,
    frame="barycentricmeanecliptic",
    obstime=Time("2022-01-14"),
)

emission = model.evaluate(skycoord, return_comps=True, nprocesses=multiprocessing.cpu_count())

fig = plt.figure(figsize=(8, 7))
for idx, comp_emission in enumerate(emission):
    hp.mollview(
        comp_emission,
        title=COMP_NAMES[idx],
        norm="log" if idx == 0 else None,
        cmap="afmhot",
        cbar=False,
        sub=(3, 2, idx + 1),
        fig=fig,
    )
plt.savefig("../img/component_maps.png", dpi=250, bbox_inches="tight")
plt.show()

Visualizing the interplanetary dust distribution

We can visualize the number density of a supported zodiacal light model by using the grid_number_density function

number_density.py

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.time import Time
from matplotlib.colors import LogNorm

from zodipy import grid_number_density

N = 200

x = np.linspace(-5, 5, N) * u.AU  # x-plane
y = np.linspace(-5, 5, N) * u.AU  # y-plane
z = np.linspace(-2, 2, N) * u.AU  # z-plane

density_grid = grid_number_density(
    x,
    y,
    z,
    obstime=Time("2021-01-01T00:00:00", scale="utc"),
    name="DIRBE",
)
density_grid = density_grid.sum(axis=0)  # Sum over all components

plt.pcolormesh(
    x,
    y,
    density_grid[N // 2].T,  # cross section in the yz-plane
    cmap="afmhot",
    norm=LogNorm(vmin=density_grid.min(), vmax=density_grid.max()),
    shading="gouraud",
    rasterized=True,
)
plt.title("Cross section of the number density in the DIRBE model")
plt.xlabel("x [AU]")
plt.ylabel("z [AU]")
plt.savefig("../img/number_density.png", dpi=300, bbox_inches="tight")
plt.show()