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)
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:
- Sky coordinates. Provided through Astropy's
SkyCoord
object. - Time of observeration(s). Also provided in the
SkyCoord
object. - 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"),
)
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"]),
)
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
:
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>
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>
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:
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:
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:
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
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()