+
+
ZodiPy is an Astropy affiliated package for simulating zodiacal light in intensity for arbitrary Solar system observers.
+
import astropy.units as u
+from astropy.time import Time
+
+from zodipy import Zodipy
+
+
+model = Zodipy("dirbe")
+
+emission = model.get_emission_ang(
+ 25 * u.micron,
+ theta=[10, 10.1, 10.2] * u.deg,
+ phi=[90, 89, 88] * u.deg,
+ obs_time=Time("2022-01-01 12:00:00"),
+ obs="earth",
+ lonlat=True,
+)
+
+print(emission)
+#> [15.53095493 15.52883577 15.53121942] MJy / sr
+What's going on here:
+Zodipy class, which is our interface, where we specify that we want to use the DIRBE interplanetary dust model.get_emission_ang method which is a method to simulate emission from angular sky coordinates (see the reference for other available simulation methods).get_emission_ang method, 25 * u.micron, specifies the wavelength of the simulated observation. Note that we use Astropy units for many of the input arguments.theta and phi represent the pointing of the observation (co-latitude and longitude, following the healpy convention). In this example we observe three sky coordinates.obs_time represents the time of observation, which we need to compute the position of the observer and all other required solar system bodies.obs represents the observer, and must be an solar system observer supported by the Astropy ephemeris used internally. If we wish to be more specific about the observer position, we can use the obs_pos keyword instead of obs, which takes in a heliocentric cartesian position in units of AU.lonlat is a boolean which converts the convention of theta and phi from co-latitude and longitude to longitude and latitude.For more information on using ZodiPy, see the usage section.
+ + + + + + +Installing ZodiPy is as simple as:
+pip install zodipy
+Note
+ZodiPy supports python versions >=3.9.
+ZodiPy has the following dependencies (these are automatically downloaded alongside ZodiPy):
+ + + + + + + +ZodiPy is an open source Python tool for simulating the zodiacal emission that a solar system observer is predicted to see given an interplanetary dust model. We attempts to make zodiacal emission simulations more accessible by providing the community with a simple interface to existing models. For other zodiacal emission tools, see Zodiacal Light Models on LAMBDA. All contributions are most welcome.
+ZodiPy supports the following interplanetary dust models:
+1.25-240 \(\boldsymbol{\mu}\)m
+100-857 GHz
+Info
+The Planck and Odegard models extend the DIRBE interplanetary dust model to CMB frequencies by fitting the blackbody emissivity of the dust in the respective DIRBE interplanetary dust components to Planck HFI data. +The distribution of the interplanetary dust is exactly the same as in the DIRBE model.
+If you see a missing model, please feel free to contact us by opening an issue on GitHub.
+For an overview of the modeling approach used in ZodiPy and other information regarding zodiacal emission and interplanetary dust modeling we refer to the scientific paper on ZodiPy:
+ + + + + + + +Zodipy
+
+
+Interface for simulating zodiacal emission.
+This class provides methods for simulating zodiacal emission given observer pointing +either in sky angles or through HEALPix pixels.
+ + + +Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
model |
+
+ str
+ |
+
+
+
+ Name of the interplanetary dust model to use in the simulations. +Defaults to DIRBE. + |
+
+ 'dirbe'
+ |
+
gauss_quad_degree |
+
+ int
+ |
+
+
+
+ Order of the Gaussian-Legendre quadrature used to evaluate +the line-of-sight integral in the simulations. Default is 50 points. + |
+
+ 50
+ |
+
interp_kind |
+
+ str
+ |
+
+
+
+ Interpolation kind used to interpolate relevant model parameters. +Defaults to 'linear'. For more information on available interpolation methods, +please visit the Scipy documentation. + |
+
+ 'linear'
+ |
+
extrapolate |
+
+ bool
+ |
+
+
+
+ If |
+
+ False
+ |
+
ephemeris |
+
+ str
+ |
+
+
+
+ Ephemeris used to compute the positions of the observer and the +Earth. Defaults to 'de432s', which requires downloading (and caching) a ~10MB +file. For more information on available ephemeridis, please visit the Astropy +documentation + |
+
+ 'de432s'
+ |
+
solar_cut |
+
+ Quantity[deg]
+ |
+
+
+
+ Cutoff angle from the sun in degrees. The emission
+for all the pointing with angular distance between the sun smaller than
+ |
+
+ None
+ |
+
solar_cut_fill_value |
+
+ float
+ |
+
+
+
+ Fill value for pixels masked with |
+
+ nan
+ |
+
n_proc |
+
+ int
+ |
+
+
+
+ Number of cores to use. If |
+
+ 1
+ |
+
zodipy/zodipy.py37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + 54 + 55 + 56 + 57 + 58 + 59 + 60 + 61 + 62 + 63 + 64 + 65 + 66 + 67 + 68 + 69 + 70 + 71 + 72 + 73 + 74 + 75 + 76 + 77 + 78 + 79 + 80 + 81 + 82 + 83 + 84 + 85 + 86 + 87 + 88 + 89 + 90 + 91 + 92 + 93 + 94 + 95 + 96 + 97 + 98 + 99 +100 +101 +102 +103 +104 +105 +106 +107 +108 +109 +110 +111 +112 +113 +114 +115 +116 +117 +118 +119 +120 +121 +122 +123 +124 +125 +126 +127 +128 +129 +130 +131 +132 +133 +134 +135 +136 +137 +138 +139 +140 +141 +142 +143 +144 +145 +146 +147 +148 +149 +150 +151 +152 +153 +154 +155 +156 +157 +158 +159 +160 +161 +162 +163 +164 +165 +166 +167 +168 +169 +170 +171 +172 +173 +174 +175 +176 +177 +178 +179 +180 +181 +182 +183 +184 +185 +186 +187 +188 +189 +190 +191 +192 +193 +194 +195 +196 +197 +198 +199 +200 +201 +202 +203 +204 +205 +206 +207 +208 +209 +210 +211 +212 +213 +214 +215 +216 +217 +218 +219 +220 +221 +222 +223 +224 +225 +226 +227 +228 +229 +230 +231 +232 +233 +234 +235 +236 +237 +238 +239 +240 +241 +242 +243 +244 +245 +246 +247 +248 +249 +250 +251 +252 +253 +254 +255 +256 +257 +258 +259 +260 +261 +262 +263 +264 +265 +266 +267 +268 +269 +270 +271 +272 +273 +274 +275 +276 +277 +278 +279 +280 +281 +282 +283 +284 +285 +286 +287 +288 +289 +290 +291 +292 +293 +294 +295 +296 +297 +298 +299 +300 +301 +302 +303 +304 +305 +306 +307 +308 +309 +310 +311 +312 +313 +314 +315 +316 +317 +318 +319 +320 +321 +322 +323 +324 +325 +326 +327 +328 +329 +330 +331 +332 +333 +334 +335 +336 +337 +338 +339 +340 +341 +342 +343 +344 +345 +346 +347 +348 +349 +350 +351 +352 +353 +354 +355 +356 +357 +358 +359 +360 +361 +362 +363 +364 +365 +366 +367 +368 +369 +370 +371 +372 +373 +374 +375 +376 +377 +378 +379 +380 +381 +382 +383 +384 +385 +386 +387 +388 +389 +390 +391 +392 +393 +394 +395 +396 +397 +398 +399 +400 +401 +402 +403 +404 +405 +406 +407 +408 +409 +410 +411 +412 +413 +414 +415 +416 +417 +418 +419 +420 +421 +422 +423 +424 +425 +426 +427 +428 +429 +430 +431 +432 +433 +434 +435 +436 +437 +438 +439 +440 +441 +442 +443 +444 +445 +446 +447 +448 +449 +450 +451 +452 +453 +454 +455 +456 +457 +458 +459 +460 +461 +462 +463 +464 +465 +466 +467 +468 +469 +470 +471 +472 +473 +474 +475 +476 +477 +478 +479 +480 +481 +482 +483 +484 +485 +486 +487 +488 +489 +490 +491 +492 +493 +494 +495 +496 +497 +498 +499 +500 +501 +502 +503 +504 +505 +506 +507 +508 +509 +510 +511 +512 +513 +514 +515 +516 +517 +518 +519 +520 +521 +522 +523 +524 +525 +526 +527 +528 +529 +530 +531 +532 +533 +534 +535 +536 +537 +538 +539 | |
supported_observers: list[str]
+
+
+ property
+
+
+Return a list of available observers given an ephemeris.
+get_binned_emission_ang(freq, theta, phi, nside, obs_time, obs='earth', obs_pos=None, weights=None, lonlat=False, return_comps=False, coord_in='E')
+
+Return the simulated binned zodiacal emission given angles on the sky.
+The pointing, for which to compute the emission, is specified in form of angles on
+the sky given by theta and phi. The emission is binned to a HEALPix map with
+resolution given by nside.
Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
freq |
+
+ FrequencyOrWavelength
+ |
+
+
+
+ Delta frequency/wavelength or a sequence of frequencies corresponding to +a bandpass over which to evaluate the zodiacal emission. The frequencies +must be strictly increasing. + |
+ + required + | +
theta |
+
+ SkyAngles
+ |
+
+
+
+ Angular co-latitude coordinate of a point, or a sequence of points, on +the celestial sphere. Must be in the range [0, π] rad. Units must be either +radians or degrees. + |
+ + required + | +
phi |
+
+ SkyAngles
+ |
+
+
+
+ Angular longitude coordinate of a point, or a sequence of points, on the +celestial sphere. Must be in the range [0, 2π] rad. Units must be either +radians or degrees. + |
+ + required + | +
nside |
+
+ int
+ |
+
+
+
+ HEALPix resolution parameter of the binned map. + |
+ + required + | +
obs_time |
+
+ Time
+ |
+
+
+
+ Time of observation. + |
+ + required + | +
obs |
+
+ str
+ |
+
+
+
+ Name of the Solar System observer. A list of all support observers (for a
+given ephemeridis) is specified in |
+
+ 'earth'
+ |
+
obs_pos |
+
+ Quantity[AU] | None
+ |
+
+
+
+ The heliocentric ecliptic cartesian position of the observer in AU.
+Overrides the |
+
+ None
+ |
+
weights |
+
+ Sequence[float] | NDArray[floating] | None
+ |
+
+
+
+ Bandpass weights corresponding the the frequencies in |
+
+ None
+ |
+
lonlat |
+
+ bool
+ |
+
+
+
+ If True, input angles |
+
+ False
+ |
+
return_comps |
+
+ bool
+ |
+
+
+
+ If True, the emission is returned component-wise. Defaults to False. + |
+
+ False
+ |
+
coord_in |
+
+ Literal['E', 'G', 'C']
+ |
+
+
+
+ Coordinate frame of the input pointing. Assumes 'E' (ecliptic +coordinates) by default. + |
+
+ 'E'
+ |
+
Returns:
+| Name | Type | +Description | +
|---|---|---|
emission |
+ Quantity[MJy / sr]
+ |
+
+
+
+ Simulated zodiacal emission in units of 'MJy/sr'. + |
+
zodipy/zodipy.py257 +258 +259 +260 +261 +262 +263 +264 +265 +266 +267 +268 +269 +270 +271 +272 +273 +274 +275 +276 +277 +278 +279 +280 +281 +282 +283 +284 +285 +286 +287 +288 +289 +290 +291 +292 +293 +294 +295 +296 +297 +298 +299 +300 +301 +302 +303 +304 +305 +306 +307 +308 +309 +310 +311 +312 +313 +314 +315 +316 +317 +318 +319 +320 +321 +322 +323 +324 +325 +326 +327 +328 +329 | |
get_binned_emission_pix(freq, pixels, nside, obs_time, obs='earth', obs_pos=None, weights=None, return_comps=False, coord_in='E')
+
+Return the simulated binned zodiacal Emission given pixel numbers.
+The pixel numbers represent the pixel indicies on a HEALPix grid with resolution
+given by nside. The emission is binned to a HEALPix map with resolution given by
+nside.
Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
freq |
+
+ FrequencyOrWavelength
+ |
+
+
+
+ Delta frequency/wavelength or a sequence of frequencies corresponding to +a bandpass over which to evaluate the zodiacal emission. The frequencies +must be strictly increasing. + |
+ + required + | +
pixels |
+
+ Pixels
+ |
+
+
+
+ HEALPix pixel indicies representing points on the celestial sphere. + |
+ + required + | +
nside |
+
+ int
+ |
+
+
+
+ HEALPix resolution parameter of the pixels and the binned map. + |
+ + required + | +
obs_time |
+
+ Time
+ |
+
+
+
+ Time of observation. + |
+ + required + | +
obs |
+
+ str
+ |
+
+
+
+ Name of the Solar System observer. A list of all support observers (for a
+given ephemeridis) is specified in |
+
+ 'earth'
+ |
+
obs_pos |
+
+ Quantity[AU] | None
+ |
+
+
+
+ The heliocentric ecliptic cartesian position of the observer in AU.
+Overrides the |
+
+ None
+ |
+
weights |
+
+ Sequence[float] | NDArray[floating] | None
+ |
+
+
+
+ Bandpass weights corresponding the the frequencies in |
+
+ None
+ |
+
return_comps |
+
+ bool
+ |
+
+
+
+ If True, the emission is returned component-wise. Defaults to False. + |
+
+ False
+ |
+
coord_in |
+
+ Literal['E', 'G', 'C']
+ |
+
+
+
+ Coordinate frame of the input pointing. Assumes 'E' (ecliptic +coordinates) by default. + |
+
+ 'E'
+ |
+
Returns:
+| Name | Type | +Description | +
|---|---|---|
emission |
+ Quantity[MJy / sr]
+ |
+
+
+
+ Simulated zodiacal emission in units of 'MJy/sr'. + |
+
zodipy/zodipy.py331 +332 +333 +334 +335 +336 +337 +338 +339 +340 +341 +342 +343 +344 +345 +346 +347 +348 +349 +350 +351 +352 +353 +354 +355 +356 +357 +358 +359 +360 +361 +362 +363 +364 +365 +366 +367 +368 +369 +370 +371 +372 +373 +374 +375 +376 +377 +378 +379 +380 +381 +382 +383 +384 +385 +386 +387 +388 +389 +390 +391 +392 | |
get_emission_ang(freq, theta, phi, obs_time, obs='earth', obs_pos=None, weights=None, lonlat=False, return_comps=False, coord_in='E')
+
+Return the simulated zodiacal emission given angles on the sky.
+The pointing, for which to compute the emission, is specified in form of angles on
+the sky given by theta and phi.
Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
freq |
+
+ FrequencyOrWavelength
+ |
+
+
+
+ Delta frequency/wavelength or a sequence of frequencies corresponding to +a bandpass over which to evaluate the zodiacal emission. The frequencies +must be strictly increasing. + |
+ + required + | +
theta |
+
+ SkyAngles
+ |
+
+
+
+ Angular co-latitude coordinate of a point, or a sequence of points, on +the celestial sphere. Must be in the range [0, π] rad. Units must be either +radians or degrees. + |
+ + required + | +
phi |
+
+ SkyAngles
+ |
+
+
+
+ Angular longitude coordinate of a point, or a sequence of points, on the +celestial sphere. Must be in the range [0, 2π] rad. Units must be either +radians or degrees. + |
+ + required + | +
obs_time |
+
+ Time
+ |
+
+
+
+ Time of observation. + |
+ + required + | +
obs |
+
+ str
+ |
+
+
+
+ Name of the Solar System observer. A list of all support observers (for a
+given ephemeridis) is specified in |
+
+ 'earth'
+ |
+
obs_pos |
+
+ Quantity[AU] | None
+ |
+
+
+
+ The heliocentric ecliptic cartesian position of the observer in AU.
+Overrides the |
+
+ None
+ |
+
weights |
+
+ Sequence[float] | NDArray[floating] | None
+ |
+
+
+
+ Bandpass weights corresponding the the frequencies in |
+
+ None
+ |
+
lonlat |
+
+ bool
+ |
+
+
+
+ If True, input angles ( |
+
+ False
+ |
+
return_comps |
+
+ bool
+ |
+
+
+
+ If True, the emission is returned component-wise. Defaults to False. + |
+
+ False
+ |
+
coord_in |
+
+ Literal['E', 'G', 'C']
+ |
+
+
+
+ Coordinate frame of the input pointing. Assumes 'E' (ecliptic +coordinates) by default. + |
+
+ 'E'
+ |
+
Returns:
+| Name | Type | +Description | +
|---|---|---|
emission |
+ Quantity[MJy / sr]
+ |
+
+
+
+ Simulated zodiacal emission in units of 'MJy/sr'. + |
+
zodipy/zodipy.py129 +130 +131 +132 +133 +134 +135 +136 +137 +138 +139 +140 +141 +142 +143 +144 +145 +146 +147 +148 +149 +150 +151 +152 +153 +154 +155 +156 +157 +158 +159 +160 +161 +162 +163 +164 +165 +166 +167 +168 +169 +170 +171 +172 +173 +174 +175 +176 +177 +178 +179 +180 +181 +182 +183 +184 +185 +186 +187 +188 +189 +190 +191 +192 +193 +194 | |
get_emission_pix(freq, pixels, nside, obs_time, obs='earth', obs_pos=None, weights=None, return_comps=False, coord_in='E')
+
+Return the simulated zodiacal emission given pixel numbers.
+The pixel numbers represent the pixel indicies on a HEALPix grid with resolution
+given by nside.
Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
freq |
+
+ FrequencyOrWavelength
+ |
+
+
+
+ Delta frequency/wavelength or a sequence of frequencies corresponding to +a bandpass over which to evaluate the zodiacal emission. The frequencies +must be strictly increasing. + |
+ + required + | +
pixels |
+
+ Pixels
+ |
+
+
+
+ HEALPix pixel indicies representing points on the celestial sphere. + |
+ + required + | +
nside |
+
+ int
+ |
+
+
+
+ HEALPix resolution parameter of the pixels and the binned map. + |
+ + required + | +
obs_time |
+
+ Time
+ |
+
+
+
+ Time of observation. + |
+ + required + | +
obs |
+
+ str
+ |
+
+
+
+ Name of the Solar System observer. A list of all support observers (for a
+given ephemeridis) is specified in |
+
+ 'earth'
+ |
+
obs_pos |
+
+ Quantity[AU] | None
+ |
+
+
+
+ The heliocentric ecliptic cartesian position of the observer in AU.
+Overrides the |
+
+ None
+ |
+
weights |
+
+ Sequence[float] | NDArray[floating] | None
+ |
+
+
+
+ Bandpass weights corresponding the the frequencies in |
+
+ None
+ |
+
return_comps |
+
+ bool
+ |
+
+
+
+ If True, the emission is returned component-wise. Defaults to False. + |
+
+ False
+ |
+
coord_in |
+
+ Literal['E', 'G', 'C']
+ |
+
+
+
+ Coordinate frame of the input pointing. Assumes 'E' (ecliptic +coordinates) by default. + |
+
+ 'E'
+ |
+
Returns:
+| Name | Type | +Description | +
|---|---|---|
emission |
+ Quantity[MJy / sr]
+ |
+
+
+
+ Simulated zodiacal emission in units of 'MJy/sr'. + |
+
zodipy/zodipy.py196 +197 +198 +199 +200 +201 +202 +203 +204 +205 +206 +207 +208 +209 +210 +211 +212 +213 +214 +215 +216 +217 +218 +219 +220 +221 +222 +223 +224 +225 +226 +227 +228 +229 +230 +231 +232 +233 +234 +235 +236 +237 +238 +239 +240 +241 +242 +243 +244 +245 +246 +247 +248 +249 +250 +251 +252 +253 +254 +255 | |
get_parameters()
+
+Return a dictionary containing the interplanetary dust model parameters.
+ +zodipy/zodipy.py103 +104 +105 | |
update_parameters(parameters)
+
+Update the interplanetary dust model parameters.
+ + + +Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
parameters |
+
+ ParameterDict
+ |
+
+
+
+ Dictionary of parameters to update. The keys must be the names
+of the parameters as defined in the model. To get the parameters dict
+of an existing model, use |
+ + required + | +
zodipy/zodipy.py107 +108 +109 +110 +111 +112 +113 +114 +115 +116 +117 +118 +119 +120 +121 +122 +123 +124 +125 +126 +127 | |
ZodiPy is an Astropy affiliated package for simulating zodiacal light in intensity for arbitrary Solar system observers.
"},{"location":"#a-simple-example","title":"A simple example","text":"import astropy.units as u\nfrom astropy.time import Time\n\nfrom zodipy import Zodipy\n\n\nmodel = Zodipy(\"dirbe\")\n\nemission = model.get_emission_ang(\n 25 * u.micron,\n theta=[10, 10.1, 10.2] * u.deg,\n phi=[90, 89, 88] * u.deg,\n obs_time=Time(\"2022-01-01 12:00:00\"),\n obs=\"earth\",\n lonlat=True,\n)\n\nprint(emission)\n#> [15.53095493 15.52883577 15.53121942] MJy / sr\nWhat's going on here:
Zodipy class, which is our interface, where we specify that we want to use the DIRBE interplanetary dust model.get_emission_ang method which is a method to simulate emission from angular sky coordinates (see the reference for other available simulation methods).get_emission_ang method, 25 * u.micron, specifies the wavelength of the simulated observation. Note that we use Astropy units for many of the input arguments.theta and phi represent the pointing of the observation (co-latitude and longitude, following the healpy convention). In this example we observe three sky coordinates.obs_time represents the time of observation, which we need to compute the position of the observer and all other required solar system bodies.obs represents the observer, and must be an solar system observer supported by the Astropy ephemeris used internally. If we wish to be more specific about the observer position, we can use the obs_pos keyword instead of obs, which takes in a heliocentric cartesian position in units of AU.lonlat is a boolean which converts the convention of theta and phi from co-latitude and longitude to longitude and latitude.For more information on using ZodiPy, see the usage section.
"},{"location":"install/","title":"Install","text":"Installing ZodiPy is as simple as:
pip install zodipy\nNote
ZodiPy supports python versions >=3.9.
"},{"location":"install/#dependencies","title":"Dependencies","text":"ZodiPy has the following dependencies (these are automatically downloaded alongside ZodiPy):
ZodiPy is an open source Python tool for simulating the zodiacal emission that a solar system observer is predicted to see given an interplanetary dust model. We attempts to make zodiacal emission simulations more accessible by providing the community with a simple interface to existing models. For other zodiacal emission tools, see Zodiacal Light Models on LAMBDA. All contributions are most welcome.
"},{"location":"introduction/#interplanetary-dust-models","title":"Interplanetary Dust Models","text":"ZodiPy supports the following interplanetary dust models:
1.25-240 \\(\\boldsymbol{\\mu}\\)m
100-857 GHz
Info
The Planck and Odegard models extend the DIRBE interplanetary dust model to CMB frequencies by fitting the blackbody emissivity of the dust in the respective DIRBE interplanetary dust components to Planck HFI data. The distribution of the interplanetary dust is exactly the same as in the DIRBE model.
If you see a missing model, please feel free to contact us by opening an issue on GitHub.
"},{"location":"introduction/#scientific-paper","title":"Scientific Paper","text":"For an overview of the modeling approach used in ZodiPy and other information regarding zodiacal emission and interplanetary dust modeling we refer to the scientific paper on ZodiPy:
Zodipy","text":"Interface for simulating zodiacal emission.
This class provides methods for simulating zodiacal emission given observer pointing either in sky angles or through HEALPix pixels.
Parameters:
Name Type Description Defaultmodel str Name of the interplanetary dust model to use in the simulations. Defaults to DIRBE.
'dirbe' gauss_quad_degree int Order of the Gaussian-Legendre quadrature used to evaluate the line-of-sight integral in the simulations. Default is 50 points.
50 interp_kind str Interpolation kind used to interpolate relevant model parameters. Defaults to 'linear'. For more information on available interpolation methods, please visit the Scipy documentation.
'linear' extrapolate bool If True all spectral quantities in the selected model are extrapolated to the requested frequencies or wavelengths. If False, an exception is raised on requested frequencies/wavelengths outside of the valid model range. Default is False.
False ephemeris str Ephemeris used to compute the positions of the observer and the Earth. Defaults to 'de432s', which requires downloading (and caching) a ~10MB file. For more information on available ephemeridis, please visit the Astropy documentation
'de432s' solar_cut Quantity[deg] Cutoff angle from the sun in degrees. The emission for all the pointing with angular distance between the sun smaller than solar_cutoff are masked. Defaults to None.
None solar_cut_fill_value float Fill value for pixels masked with solar_cut. Defaults to np.nan.
nan n_proc int Number of cores to use. If n_proc is greater than 1, the line-of-sight integrals are parallelized using the multiprocessing module. Defaults to 1.
1 Source code in zodipy/zodipy.py class Zodipy:\n\"\"\"Interface for simulating zodiacal emission.\n\n This class provides methods for simulating zodiacal emission given observer pointing\n either in sky angles or through HEALPix pixels.\n\n Args:\n model (str): Name of the interplanetary dust model to use in the simulations.\n Defaults to DIRBE.\n gauss_quad_degree (int): Order of the Gaussian-Legendre quadrature used to evaluate\n the line-of-sight integral in the simulations. Default is 50 points.\n interp_kind (str): Interpolation kind used to interpolate relevant model parameters.\n Defaults to 'linear'. For more information on available interpolation methods,\n please visit the [Scipy documentation](\n https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html).\n extrapolate (bool): If `True` all spectral quantities in the selected model are\n extrapolated to the requested frequencies or wavelengths. If `False`, an\n exception is raised on requested frequencies/wavelengths outside of the\n valid model range. Default is `False`.\n ephemeris (str): Ephemeris used to compute the positions of the observer and the\n Earth. Defaults to 'de432s', which requires downloading (and caching) a ~10MB\n file. For more information on available ephemeridis, please visit the [Astropy\n documentation](https://docs.astropy.org/en/stable/coordinates/solarsystem.html)\n solar_cut (u.Quantity[u.deg]): Cutoff angle from the sun in degrees. The emission\n for all the pointing with angular distance between the sun smaller than\n `solar_cutoff` are masked. Defaults to `None`.\n solar_cut_fill_value (float): Fill value for pixels masked with `solar_cut`.\n Defaults to `np.nan`.\n n_proc (int): Number of cores to use. If `n_proc` is greater than 1, the line-of-sight\n integrals are parallelized using the `multiprocessing` module. Defaults to 1.\n\n \"\"\"\n\n def __init__(\n self,\n model: str = \"dirbe\",\n gauss_quad_degree: int = 50,\n extrapolate: bool = False,\n interp_kind: str = \"linear\",\n ephemeris: str = \"de432s\",\n solar_cut: u.Quantity[u.deg] | None = None,\n solar_cut_fill_value: float = np.nan,\n n_proc: int = 1,\n ) -> None:\n self.model = model\n self.gauss_quad_degree = gauss_quad_degree\n self.extrapolate = extrapolate\n self.interp_kind = interp_kind\n self.ephemeris = ephemeris\n self.solar_cut = solar_cut.to(u.rad) if solar_cut is not None else solar_cut\n self.solar_cut_fill_value = solar_cut_fill_value\n self.n_proc = n_proc\n\n self._interpolator = partial(\n interp1d,\n kind=self.interp_kind,\n fill_value=\"extrapolate\" if self.extrapolate else np.nan,\n )\n self._ipd_model = model_registry.get_model(model)\n self._gauss_points_and_weights = np.polynomial.legendre.leggauss(gauss_quad_degree)\n\n @property\n def supported_observers(self) -> list[str]:\n\"\"\"Return a list of available observers given an ephemeris.\"\"\"\n return [*list(solar_system_ephemeris.bodies), \"semb-l2\"]\n\n def get_parameters(self) -> ParameterDict:\n\"\"\"Return a dictionary containing the interplanetary dust model parameters.\"\"\"\n return self._ipd_model.to_dict()\n\n def update_parameters(self, parameters: ParameterDict) -> None:\n\"\"\"Update the interplanetary dust model parameters.\n\n Args:\n parameters: Dictionary of parameters to update. The keys must be the names\n of the parameters as defined in the model. To get the parameters dict\n of an existing model, use `Zodipy(\"dirbe\").get_parameters()`.\n\n \"\"\"\n _dict = parameters.copy()\n _dict[\"comps\"] = {}\n for key, value in parameters.items():\n if key == \"comps\":\n for comp_key, comp_value in value.items():\n _dict[\"comps\"][ComponentLabel(comp_key)] = type(\n self._ipd_model.comps[ComponentLabel(comp_key)]\n )(**comp_value)\n elif isinstance(value, dict):\n _dict[key] = {ComponentLabel(k): v for k, v in value.items()}\n\n self._ipd_model = self._ipd_model.__class__(**_dict)\n\n def get_emission_ang(\n self,\n freq: FrequencyOrWavelength,\n theta: SkyAngles,\n phi: SkyAngles,\n obs_time: Time,\n obs: str = \"earth\",\n obs_pos: u.Quantity[u.AU] | None = None,\n weights: Sequence[float] | npt.NDArray[np.floating] | None = None,\n lonlat: bool = False,\n return_comps: bool = False,\n coord_in: Literal[\"E\", \"G\", \"C\"] = \"E\",\n ) -> u.Quantity[u.MJy / u.sr]:\n\"\"\"Return the simulated zodiacal emission given angles on the sky.\n\n The pointing, for which to compute the emission, is specified in form of angles on\n the sky given by `theta` and `phi`.\n\n Args:\n freq: Delta frequency/wavelength or a sequence of frequencies corresponding to\n a bandpass over which to evaluate the zodiacal emission. The frequencies\n must be strictly increasing.\n theta: Angular co-latitude coordinate of a point, or a sequence of points, on\n the celestial sphere. Must be in the range [0, \u03c0] rad. Units must be either\n radians or degrees.\n phi: Angular longitude coordinate of a point, or a sequence of points, on the\n celestial sphere. Must be in the range [0, 2\u03c0] rad. Units must be either\n radians or degrees.\n obs_time: Time of observation.\n obs: Name of the Solar System observer. A list of all support observers (for a\n given ephemeridis) is specified in `supported_observers` attribute of the\n `Zodipy` instance. Defaults to 'earth'.\n obs_pos: The heliocentric ecliptic cartesian position of the observer in AU.\n Overrides the `obs` argument. Default is None.\n weights: Bandpass weights corresponding the the frequencies in `freq`. The weights\n are assumed to be given in spectral radiance units (Jy/sr).\n lonlat: If True, input angles (`theta`, `phi`) are assumed to be longitude and\n latitude, otherwise, they are co-latitude and longitude.\n return_comps: If True, the emission is returned component-wise. Defaults to False.\n coord_in: Coordinate frame of the input pointing. Assumes 'E' (ecliptic\n coordinates) by default.\n\n Returns:\n emission: Simulated zodiacal emission in units of 'MJy/sr'.\n\n \"\"\"\n theta, phi = get_validated_ang(theta=theta, phi=phi, lonlat=lonlat)\n\n unique_angles, indicies = np.unique(np.asarray([theta, phi]), return_inverse=True, axis=1)\n unit_vectors = get_unit_vectors_from_ang(\n coord_in=coord_in,\n theta=unique_angles[0],\n phi=unique_angles[1],\n lonlat=lonlat,\n )\n\n return self._compute_emission(\n freq=freq,\n weights=weights,\n obs=obs,\n obs_time=obs_time,\n obs_pos=obs_pos,\n unit_vectors=unit_vectors,\n indicies=indicies,\n return_comps=return_comps,\n )\n\n def get_emission_pix(\n self,\n freq: FrequencyOrWavelength,\n pixels: Pixels,\n nside: int,\n obs_time: Time,\n obs: str = \"earth\",\n obs_pos: u.Quantity[u.AU] | None = None,\n weights: Sequence[float] | npt.NDArray[np.floating] | None = None,\n return_comps: bool = False,\n coord_in: Literal[\"E\", \"G\", \"C\"] = \"E\",\n ) -> u.Quantity[u.MJy / u.sr]:\n\"\"\"Return the simulated zodiacal emission given pixel numbers.\n\n The pixel numbers represent the pixel indicies on a HEALPix grid with resolution\n given by `nside`.\n\n Args:\n freq: Delta frequency/wavelength or a sequence of frequencies corresponding to\n a bandpass over which to evaluate the zodiacal emission. The frequencies\n must be strictly increasing.\n pixels: HEALPix pixel indicies representing points on the celestial sphere.\n nside: HEALPix resolution parameter of the pixels and the binned map.\n obs_time: Time of observation.\n obs: Name of the Solar System observer. A list of all support observers (for a\n given ephemeridis) is specified in `supported_observers` attribute of the\n `Zodipy` instance. Defaults to 'earth'.\n obs_pos: The heliocentric ecliptic cartesian position of the observer in AU.\n Overrides the `obs` argument. Default is None.\n weights: Bandpass weights corresponding the the frequencies in `freq`. The weights\n are assumed to be given in spectral radiance units (Jy/sr).\n return_comps: If True, the emission is returned component-wise. Defaults to False.\n coord_in: Coordinate frame of the input pointing. Assumes 'E' (ecliptic\n coordinates) by default.\n\n Returns:\n emission: Simulated zodiacal emission in units of 'MJy/sr'.\n\n \"\"\"\n pixels = get_validated_pix(pixels=pixels, nside=nside)\n\n unique_pixels, indicies = np.unique(pixels, return_inverse=True)\n unit_vectors = get_unit_vectors_from_pixels(\n coord_in=coord_in,\n pixels=unique_pixels,\n nside=nside,\n )\n\n return self._compute_emission(\n freq=freq,\n weights=weights,\n obs=obs,\n obs_time=obs_time,\n obs_pos=obs_pos,\n unit_vectors=unit_vectors,\n indicies=indicies,\n pixels=unique_pixels,\n nside=nside,\n return_comps=return_comps,\n )\n\n def get_binned_emission_ang(\n self,\n freq: FrequencyOrWavelength,\n theta: SkyAngles,\n phi: SkyAngles,\n nside: int,\n obs_time: Time,\n obs: str = \"earth\",\n obs_pos: u.Quantity[u.AU] | None = None,\n weights: Sequence[float] | npt.NDArray[np.floating] | None = None,\n lonlat: bool = False,\n return_comps: bool = False,\n coord_in: Literal[\"E\", \"G\", \"C\"] = \"E\",\n ) -> u.Quantity[u.MJy / u.sr]:\n\"\"\"Return the simulated binned zodiacal emission given angles on the sky.\n\n The pointing, for which to compute the emission, is specified in form of angles on\n the sky given by `theta` and `phi`. The emission is binned to a HEALPix map with\n resolution given by `nside`.\n\n Args:\n freq: Delta frequency/wavelength or a sequence of frequencies corresponding to\n a bandpass over which to evaluate the zodiacal emission. The frequencies\n must be strictly increasing.\n theta: Angular co-latitude coordinate of a point, or a sequence of points, on\n the celestial sphere. Must be in the range [0, \u03c0] rad. Units must be either\n radians or degrees.\n phi: Angular longitude coordinate of a point, or a sequence of points, on the\n celestial sphere. Must be in the range [0, 2\u03c0] rad. Units must be either\n radians or degrees.\n nside: HEALPix resolution parameter of the binned map.\n obs_time: Time of observation.\n obs: Name of the Solar System observer. A list of all support observers (for a\n given ephemeridis) is specified in `supported_observers` attribute of the\n `Zodipy` instance. Defaults to 'earth'.\n obs_pos: The heliocentric ecliptic cartesian position of the observer in AU.\n Overrides the `obs` argument. Default is None.\n weights: Bandpass weights corresponding the the frequencies in `freq`. The weights\n are assumed to be given in spectral radiance units (Jy/sr).\n lonlat: If True, input angles `theta`, `phi` are assumed to be longitude and\n latitude, otherwise, they are co-latitude and longitude.\n return_comps: If True, the emission is returned component-wise. Defaults to False.\n coord_in: Coordinate frame of the input pointing. Assumes 'E' (ecliptic\n coordinates) by default.\n\n Returns:\n emission: Simulated zodiacal emission in units of 'MJy/sr'.\n\n \"\"\"\n theta, phi = get_validated_ang(theta=theta, phi=phi, lonlat=lonlat)\n\n unique_angles, counts = np.unique(np.asarray([theta, phi]), return_counts=True, axis=1)\n unique_pixels = hp.ang2pix(nside, *unique_angles, lonlat=lonlat)\n unit_vectors = get_unit_vectors_from_ang(\n coord_in=coord_in,\n theta=unique_angles[0],\n phi=unique_angles[1],\n lonlat=lonlat,\n )\n\n return self._compute_emission(\n freq=freq,\n weights=weights,\n obs=obs,\n obs_time=obs_time,\n obs_pos=obs_pos,\n unit_vectors=unit_vectors,\n indicies=counts,\n binned=True,\n pixels=unique_pixels,\n nside=nside,\n return_comps=return_comps,\n )\n\n def get_binned_emission_pix(\n self,\n freq: FrequencyOrWavelength,\n pixels: Pixels,\n nside: int,\n obs_time: Time,\n obs: str = \"earth\",\n obs_pos: u.Quantity[u.AU] | None = None,\n weights: Sequence[float] | npt.NDArray[np.floating] | None = None,\n return_comps: bool = False,\n coord_in: Literal[\"E\", \"G\", \"C\"] = \"E\",\n ) -> u.Quantity[u.MJy / u.sr]:\n\"\"\"Return the simulated binned zodiacal Emission given pixel numbers.\n\n The pixel numbers represent the pixel indicies on a HEALPix grid with resolution\n given by `nside`. The emission is binned to a HEALPix map with resolution given by\n `nside`.\n\n Args:\n freq: Delta frequency/wavelength or a sequence of frequencies corresponding to\n a bandpass over which to evaluate the zodiacal emission. The frequencies\n must be strictly increasing.\n pixels: HEALPix pixel indicies representing points on the celestial sphere.\n nside: HEALPix resolution parameter of the pixels and the binned map.\n obs_time: Time of observation.\n obs: Name of the Solar System observer. A list of all support observers (for a\n given ephemeridis) is specified in `supported_observers` attribute of the\n `Zodipy` instance. Defaults to 'earth'.\n obs_pos: The heliocentric ecliptic cartesian position of the observer in AU.\n Overrides the `obs` argument. Default is None.\n weights: Bandpass weights corresponding the the frequencies in `freq`. The weights\n are assumed to be given in spectral radiance units (Jy/sr).\n return_comps: If True, the emission is returned component-wise. Defaults to False.\n coord_in: Coordinate frame of the input pointing. Assumes 'E' (ecliptic\n coordinates) by default.\n\n Returns:\n emission: Simulated zodiacal emission in units of 'MJy/sr'.\n\n \"\"\"\n pixels = get_validated_pix(pixels=pixels, nside=nside)\n\n unique_pixels, counts = np.unique(pixels, return_counts=True)\n unit_vectors = get_unit_vectors_from_pixels(\n coord_in=coord_in,\n pixels=unique_pixels,\n nside=nside,\n )\n\n return self._compute_emission(\n freq=freq,\n weights=weights,\n obs=obs,\n obs_time=obs_time,\n obs_pos=obs_pos,\n unit_vectors=unit_vectors,\n indicies=counts,\n binned=True,\n pixels=unique_pixels,\n nside=nside,\n return_comps=return_comps,\n )\n\n def _compute_emission(\n self,\n freq: FrequencyOrWavelength,\n weights: Sequence[float] | npt.NDArray[np.floating] | None,\n obs: str,\n obs_time: Time,\n unit_vectors: npt.NDArray[np.float64],\n indicies: npt.NDArray[np.int64],\n binned: bool = False,\n obs_pos: u.Quantity[u.AU] | None = None,\n pixels: npt.NDArray[np.int64] | None = None,\n nside: int | None = None,\n return_comps: bool = False,\n ) -> u.Quantity[u.MJy / u.sr]:\n\"\"\"Compute the component-wise zodiacal emission.\"\"\"\n bandpass = validate_and_get_bandpass(\n freq=freq,\n weights=weights,\n model=self._ipd_model,\n extrapolate=self.extrapolate,\n )\n\n # Get model parameters, some of which have been interpolated to the given\n # frequency or bandpass.\n source_parameters = SOURCE_PARAMS_MAPPING[type(self._ipd_model)](\n bandpass, self._ipd_model, self._interpolator\n )\n\n observer_position, earth_position = get_obs_and_earth_positions(\n obs=obs, obs_time=obs_time, obs_pos=obs_pos\n )\n\n # Get the integration limits for each zodiacal component (which may be\n # different or the same depending on the model) along all line of sights.\n start, stop = get_line_of_sight_start_and_stop_distances(\n components=self._ipd_model.comps.keys(),\n unit_vectors=unit_vectors,\n obs_pos=observer_position,\n )\n\n density_partials = construct_density_partials_comps(\n comps=self._ipd_model.comps,\n dynamic_params={\"X_earth\": earth_position},\n )\n\n # Make table of pre-computed bandpass integrated blackbody emission.\n bandpass_interpolatation_table = get_bandpass_interpolation_table(bandpass)\n\n common_integrand = partial(\n EMISSION_MAPPING[type(self._ipd_model)],\n X_obs=observer_position,\n bp_interpolation_table=bandpass_interpolatation_table,\n **source_parameters[\"common\"],\n )\n\n if self.n_proc > 1:\n unit_vector_chunks = np.array_split(unit_vectors, self.n_proc, axis=-1)\n integrated_comp_emission = np.zeros((len(self._ipd_model.comps), unit_vectors.shape[1]))\n with multiprocessing.get_context(SYS_PROC_START_METHOD).Pool(\n processes=self.n_proc\n ) as pool:\n for idx, comp_label in enumerate(self._ipd_model.comps.keys()):\n stop_chunks = np.array_split(stop[comp_label], self.n_proc, axis=-1)\n if start[comp_label].size == 1:\n start_chunks = [start[comp_label]] * self.n_proc\n else:\n start_chunks = np.array_split(start[comp_label], self.n_proc, axis=-1)\n comp_integrands = [\n partial(\n common_integrand,\n u_los=np.expand_dims(unit_vectors, axis=-1),\n start=np.expand_dims(start, axis=-1),\n stop=np.expand_dims(stop, axis=-1),\n get_density_function=density_partials[comp_label],\n **source_parameters[comp_label],\n )\n for unit_vectors, start, stop in zip(\n unit_vector_chunks, start_chunks, stop_chunks\n )\n ]\n\n proc_chunks = [\n pool.apply_async(\n _integrate_gauss_quad,\n args=(comp_integrand, *self._gauss_points_and_weights),\n )\n for comp_integrand in comp_integrands\n ]\n\n integrated_comp_emission[idx] += (\n np.concatenate([result.get() for result in proc_chunks])\n * 0.5\n * (stop[comp_label] - start[comp_label])\n )\n\n else:\n integrated_comp_emission = np.zeros((len(self._ipd_model.comps), unit_vectors.shape[1]))\n unit_vectors_expanded = np.expand_dims(unit_vectors, axis=-1)\n\n for idx, comp_label in enumerate(self._ipd_model.comps.keys()):\n comp_integrand = partial(\n common_integrand,\n u_los=unit_vectors_expanded,\n start=np.expand_dims(start[comp_label], axis=-1),\n stop=np.expand_dims(stop[comp_label], axis=-1),\n get_density_function=density_partials[comp_label],\n **source_parameters[comp_label],\n )\n\n integrated_comp_emission[idx] = (\n _integrate_gauss_quad(comp_integrand, *self._gauss_points_and_weights)\n * 0.5\n * (stop[comp_label] - start[comp_label])\n )\n\n emission = np.zeros(\n (\n len(self._ipd_model.comps),\n hp.nside2npix(nside) if binned else indicies.size,\n )\n )\n if binned:\n emission[:, pixels] = integrated_comp_emission\n else:\n emission = integrated_comp_emission[:, indicies]\n\n if self.solar_cut is not None:\n ang_dist = hp.rotator.angdist(-observer_position.flatten(), unit_vectors)\n solar_mask = ang_dist < self.solar_cut.value\n if binned and pixels is not None:\n emission[:, pixels[solar_mask]] = self.solar_cut_fill_value\n else:\n emission[:, solar_mask[indicies]] = self.solar_cut_fill_value\n\n emission = (emission << SPECIFIC_INTENSITY_UNITS).to(u.MJy / u.sr)\n\n return emission if return_comps else emission.sum(axis=0)\n\n def __repr__(self) -> str:\n repr_str = f\"{self.__class__.__name__}(\"\n for attribute_name, attribute in self.__dict__.items():\n if attribute_name.startswith(\"_\"):\n continue\n repr_str += f\"{attribute_name}={attribute!r}, \"\n\n return repr_str[:-2] + \")\"\nsupported_observers: list[str] property","text":"Return a list of available observers given an ephemeris.
"},{"location":"reference/#zodipy.zodipy.Zodipy.get_binned_emission_ang","title":"get_binned_emission_ang(freq, theta, phi, nside, obs_time, obs='earth', obs_pos=None, weights=None, lonlat=False, return_comps=False, coord_in='E')","text":"Return the simulated binned zodiacal emission given angles on the sky.
The pointing, for which to compute the emission, is specified in form of angles on the sky given by theta and phi. The emission is binned to a HEALPix map with resolution given by nside.
Parameters:
Name Type Description Defaultfreq FrequencyOrWavelength Delta frequency/wavelength or a sequence of frequencies corresponding to a bandpass over which to evaluate the zodiacal emission. The frequencies must be strictly increasing.
requiredtheta SkyAngles Angular co-latitude coordinate of a point, or a sequence of points, on the celestial sphere. Must be in the range [0, \u03c0] rad. Units must be either radians or degrees.
requiredphi SkyAngles Angular longitude coordinate of a point, or a sequence of points, on the celestial sphere. Must be in the range [0, 2\u03c0] rad. Units must be either radians or degrees.
requirednside int HEALPix resolution parameter of the binned map.
requiredobs_time Time Time of observation.
requiredobs str Name of the Solar System observer. A list of all support observers (for a given ephemeridis) is specified in supported_observers attribute of the Zodipy instance. Defaults to 'earth'.
'earth' obs_pos Quantity[AU] | None The heliocentric ecliptic cartesian position of the observer in AU. Overrides the obs argument. Default is None.
None weights Sequence[float] | NDArray[floating] | None Bandpass weights corresponding the the frequencies in freq. The weights are assumed to be given in spectral radiance units (Jy/sr).
None lonlat bool If True, input angles theta, phi are assumed to be longitude and latitude, otherwise, they are co-latitude and longitude.
False return_comps bool If True, the emission is returned component-wise. Defaults to False.
False coord_in Literal['E', 'G', 'C'] Coordinate frame of the input pointing. Assumes 'E' (ecliptic coordinates) by default.
'E' Returns:
Name Type Descriptionemission Quantity[MJy / sr] Simulated zodiacal emission in units of 'MJy/sr'.
Source code inzodipy/zodipy.py def get_binned_emission_ang(\n self,\n freq: FrequencyOrWavelength,\n theta: SkyAngles,\n phi: SkyAngles,\n nside: int,\n obs_time: Time,\n obs: str = \"earth\",\n obs_pos: u.Quantity[u.AU] | None = None,\n weights: Sequence[float] | npt.NDArray[np.floating] | None = None,\n lonlat: bool = False,\n return_comps: bool = False,\n coord_in: Literal[\"E\", \"G\", \"C\"] = \"E\",\n) -> u.Quantity[u.MJy / u.sr]:\n\"\"\"Return the simulated binned zodiacal emission given angles on the sky.\n\n The pointing, for which to compute the emission, is specified in form of angles on\n the sky given by `theta` and `phi`. The emission is binned to a HEALPix map with\n resolution given by `nside`.\n\n Args:\n freq: Delta frequency/wavelength or a sequence of frequencies corresponding to\n a bandpass over which to evaluate the zodiacal emission. The frequencies\n must be strictly increasing.\n theta: Angular co-latitude coordinate of a point, or a sequence of points, on\n the celestial sphere. Must be in the range [0, \u03c0] rad. Units must be either\n radians or degrees.\n phi: Angular longitude coordinate of a point, or a sequence of points, on the\n celestial sphere. Must be in the range [0, 2\u03c0] rad. Units must be either\n radians or degrees.\n nside: HEALPix resolution parameter of the binned map.\n obs_time: Time of observation.\n obs: Name of the Solar System observer. A list of all support observers (for a\n given ephemeridis) is specified in `supported_observers` attribute of the\n `Zodipy` instance. Defaults to 'earth'.\n obs_pos: The heliocentric ecliptic cartesian position of the observer in AU.\n Overrides the `obs` argument. Default is None.\n weights: Bandpass weights corresponding the the frequencies in `freq`. The weights\n are assumed to be given in spectral radiance units (Jy/sr).\n lonlat: If True, input angles `theta`, `phi` are assumed to be longitude and\n latitude, otherwise, they are co-latitude and longitude.\n return_comps: If True, the emission is returned component-wise. Defaults to False.\n coord_in: Coordinate frame of the input pointing. Assumes 'E' (ecliptic\n coordinates) by default.\n\n Returns:\n emission: Simulated zodiacal emission in units of 'MJy/sr'.\n\n \"\"\"\n theta, phi = get_validated_ang(theta=theta, phi=phi, lonlat=lonlat)\n\n unique_angles, counts = np.unique(np.asarray([theta, phi]), return_counts=True, axis=1)\n unique_pixels = hp.ang2pix(nside, *unique_angles, lonlat=lonlat)\n unit_vectors = get_unit_vectors_from_ang(\n coord_in=coord_in,\n theta=unique_angles[0],\n phi=unique_angles[1],\n lonlat=lonlat,\n )\n\n return self._compute_emission(\n freq=freq,\n weights=weights,\n obs=obs,\n obs_time=obs_time,\n obs_pos=obs_pos,\n unit_vectors=unit_vectors,\n indicies=counts,\n binned=True,\n pixels=unique_pixels,\n nside=nside,\n return_comps=return_comps,\n )\nget_binned_emission_pix(freq, pixels, nside, obs_time, obs='earth', obs_pos=None, weights=None, return_comps=False, coord_in='E')","text":"Return the simulated binned zodiacal Emission given pixel numbers.
The pixel numbers represent the pixel indicies on a HEALPix grid with resolution given by nside. The emission is binned to a HEALPix map with resolution given by nside.
Parameters:
Name Type Description Defaultfreq FrequencyOrWavelength Delta frequency/wavelength or a sequence of frequencies corresponding to a bandpass over which to evaluate the zodiacal emission. The frequencies must be strictly increasing.
requiredpixels Pixels HEALPix pixel indicies representing points on the celestial sphere.
requirednside int HEALPix resolution parameter of the pixels and the binned map.
requiredobs_time Time Time of observation.
requiredobs str Name of the Solar System observer. A list of all support observers (for a given ephemeridis) is specified in supported_observers attribute of the Zodipy instance. Defaults to 'earth'.
'earth' obs_pos Quantity[AU] | None The heliocentric ecliptic cartesian position of the observer in AU. Overrides the obs argument. Default is None.
None weights Sequence[float] | NDArray[floating] | None Bandpass weights corresponding the the frequencies in freq. The weights are assumed to be given in spectral radiance units (Jy/sr).
None return_comps bool If True, the emission is returned component-wise. Defaults to False.
False coord_in Literal['E', 'G', 'C'] Coordinate frame of the input pointing. Assumes 'E' (ecliptic coordinates) by default.
'E' Returns:
Name Type Descriptionemission Quantity[MJy / sr] Simulated zodiacal emission in units of 'MJy/sr'.
Source code inzodipy/zodipy.py def get_binned_emission_pix(\n self,\n freq: FrequencyOrWavelength,\n pixels: Pixels,\n nside: int,\n obs_time: Time,\n obs: str = \"earth\",\n obs_pos: u.Quantity[u.AU] | None = None,\n weights: Sequence[float] | npt.NDArray[np.floating] | None = None,\n return_comps: bool = False,\n coord_in: Literal[\"E\", \"G\", \"C\"] = \"E\",\n) -> u.Quantity[u.MJy / u.sr]:\n\"\"\"Return the simulated binned zodiacal Emission given pixel numbers.\n\n The pixel numbers represent the pixel indicies on a HEALPix grid with resolution\n given by `nside`. The emission is binned to a HEALPix map with resolution given by\n `nside`.\n\n Args:\n freq: Delta frequency/wavelength or a sequence of frequencies corresponding to\n a bandpass over which to evaluate the zodiacal emission. The frequencies\n must be strictly increasing.\n pixels: HEALPix pixel indicies representing points on the celestial sphere.\n nside: HEALPix resolution parameter of the pixels and the binned map.\n obs_time: Time of observation.\n obs: Name of the Solar System observer. A list of all support observers (for a\n given ephemeridis) is specified in `supported_observers` attribute of the\n `Zodipy` instance. Defaults to 'earth'.\n obs_pos: The heliocentric ecliptic cartesian position of the observer in AU.\n Overrides the `obs` argument. Default is None.\n weights: Bandpass weights corresponding the the frequencies in `freq`. The weights\n are assumed to be given in spectral radiance units (Jy/sr).\n return_comps: If True, the emission is returned component-wise. Defaults to False.\n coord_in: Coordinate frame of the input pointing. Assumes 'E' (ecliptic\n coordinates) by default.\n\n Returns:\n emission: Simulated zodiacal emission in units of 'MJy/sr'.\n\n \"\"\"\n pixels = get_validated_pix(pixels=pixels, nside=nside)\n\n unique_pixels, counts = np.unique(pixels, return_counts=True)\n unit_vectors = get_unit_vectors_from_pixels(\n coord_in=coord_in,\n pixels=unique_pixels,\n nside=nside,\n )\n\n return self._compute_emission(\n freq=freq,\n weights=weights,\n obs=obs,\n obs_time=obs_time,\n obs_pos=obs_pos,\n unit_vectors=unit_vectors,\n indicies=counts,\n binned=True,\n pixels=unique_pixels,\n nside=nside,\n return_comps=return_comps,\n )\nget_emission_ang(freq, theta, phi, obs_time, obs='earth', obs_pos=None, weights=None, lonlat=False, return_comps=False, coord_in='E')","text":"Return the simulated zodiacal emission given angles on the sky.
The pointing, for which to compute the emission, is specified in form of angles on the sky given by theta and phi.
Parameters:
Name Type Description Defaultfreq FrequencyOrWavelength Delta frequency/wavelength or a sequence of frequencies corresponding to a bandpass over which to evaluate the zodiacal emission. The frequencies must be strictly increasing.
requiredtheta SkyAngles Angular co-latitude coordinate of a point, or a sequence of points, on the celestial sphere. Must be in the range [0, \u03c0] rad. Units must be either radians or degrees.
requiredphi SkyAngles Angular longitude coordinate of a point, or a sequence of points, on the celestial sphere. Must be in the range [0, 2\u03c0] rad. Units must be either radians or degrees.
requiredobs_time Time Time of observation.
requiredobs str Name of the Solar System observer. A list of all support observers (for a given ephemeridis) is specified in supported_observers attribute of the Zodipy instance. Defaults to 'earth'.
'earth' obs_pos Quantity[AU] | None The heliocentric ecliptic cartesian position of the observer in AU. Overrides the obs argument. Default is None.
None weights Sequence[float] | NDArray[floating] | None Bandpass weights corresponding the the frequencies in freq. The weights are assumed to be given in spectral radiance units (Jy/sr).
None lonlat bool If True, input angles (theta, phi) are assumed to be longitude and latitude, otherwise, they are co-latitude and longitude.
False return_comps bool If True, the emission is returned component-wise. Defaults to False.
False coord_in Literal['E', 'G', 'C'] Coordinate frame of the input pointing. Assumes 'E' (ecliptic coordinates) by default.
'E' Returns:
Name Type Descriptionemission Quantity[MJy / sr] Simulated zodiacal emission in units of 'MJy/sr'.
Source code inzodipy/zodipy.py def get_emission_ang(\n self,\n freq: FrequencyOrWavelength,\n theta: SkyAngles,\n phi: SkyAngles,\n obs_time: Time,\n obs: str = \"earth\",\n obs_pos: u.Quantity[u.AU] | None = None,\n weights: Sequence[float] | npt.NDArray[np.floating] | None = None,\n lonlat: bool = False,\n return_comps: bool = False,\n coord_in: Literal[\"E\", \"G\", \"C\"] = \"E\",\n) -> u.Quantity[u.MJy / u.sr]:\n\"\"\"Return the simulated zodiacal emission given angles on the sky.\n\n The pointing, for which to compute the emission, is specified in form of angles on\n the sky given by `theta` and `phi`.\n\n Args:\n freq: Delta frequency/wavelength or a sequence of frequencies corresponding to\n a bandpass over which to evaluate the zodiacal emission. The frequencies\n must be strictly increasing.\n theta: Angular co-latitude coordinate of a point, or a sequence of points, on\n the celestial sphere. Must be in the range [0, \u03c0] rad. Units must be either\n radians or degrees.\n phi: Angular longitude coordinate of a point, or a sequence of points, on the\n celestial sphere. Must be in the range [0, 2\u03c0] rad. Units must be either\n radians or degrees.\n obs_time: Time of observation.\n obs: Name of the Solar System observer. A list of all support observers (for a\n given ephemeridis) is specified in `supported_observers` attribute of the\n `Zodipy` instance. Defaults to 'earth'.\n obs_pos: The heliocentric ecliptic cartesian position of the observer in AU.\n Overrides the `obs` argument. Default is None.\n weights: Bandpass weights corresponding the the frequencies in `freq`. The weights\n are assumed to be given in spectral radiance units (Jy/sr).\n lonlat: If True, input angles (`theta`, `phi`) are assumed to be longitude and\n latitude, otherwise, they are co-latitude and longitude.\n return_comps: If True, the emission is returned component-wise. Defaults to False.\n coord_in: Coordinate frame of the input pointing. Assumes 'E' (ecliptic\n coordinates) by default.\n\n Returns:\n emission: Simulated zodiacal emission in units of 'MJy/sr'.\n\n \"\"\"\n theta, phi = get_validated_ang(theta=theta, phi=phi, lonlat=lonlat)\n\n unique_angles, indicies = np.unique(np.asarray([theta, phi]), return_inverse=True, axis=1)\n unit_vectors = get_unit_vectors_from_ang(\n coord_in=coord_in,\n theta=unique_angles[0],\n phi=unique_angles[1],\n lonlat=lonlat,\n )\n\n return self._compute_emission(\n freq=freq,\n weights=weights,\n obs=obs,\n obs_time=obs_time,\n obs_pos=obs_pos,\n unit_vectors=unit_vectors,\n indicies=indicies,\n return_comps=return_comps,\n )\nget_emission_pix(freq, pixels, nside, obs_time, obs='earth', obs_pos=None, weights=None, return_comps=False, coord_in='E')","text":"Return the simulated zodiacal emission given pixel numbers.
The pixel numbers represent the pixel indicies on a HEALPix grid with resolution given by nside.
Parameters:
Name Type Description Defaultfreq FrequencyOrWavelength Delta frequency/wavelength or a sequence of frequencies corresponding to a bandpass over which to evaluate the zodiacal emission. The frequencies must be strictly increasing.
requiredpixels Pixels HEALPix pixel indicies representing points on the celestial sphere.
requirednside int HEALPix resolution parameter of the pixels and the binned map.
requiredobs_time Time Time of observation.
requiredobs str Name of the Solar System observer. A list of all support observers (for a given ephemeridis) is specified in supported_observers attribute of the Zodipy instance. Defaults to 'earth'.
'earth' obs_pos Quantity[AU] | None The heliocentric ecliptic cartesian position of the observer in AU. Overrides the obs argument. Default is None.
None weights Sequence[float] | NDArray[floating] | None Bandpass weights corresponding the the frequencies in freq. The weights are assumed to be given in spectral radiance units (Jy/sr).
None return_comps bool If True, the emission is returned component-wise. Defaults to False.
False coord_in Literal['E', 'G', 'C'] Coordinate frame of the input pointing. Assumes 'E' (ecliptic coordinates) by default.
'E' Returns:
Name Type Descriptionemission Quantity[MJy / sr] Simulated zodiacal emission in units of 'MJy/sr'.
Source code inzodipy/zodipy.py def get_emission_pix(\n self,\n freq: FrequencyOrWavelength,\n pixels: Pixels,\n nside: int,\n obs_time: Time,\n obs: str = \"earth\",\n obs_pos: u.Quantity[u.AU] | None = None,\n weights: Sequence[float] | npt.NDArray[np.floating] | None = None,\n return_comps: bool = False,\n coord_in: Literal[\"E\", \"G\", \"C\"] = \"E\",\n) -> u.Quantity[u.MJy / u.sr]:\n\"\"\"Return the simulated zodiacal emission given pixel numbers.\n\n The pixel numbers represent the pixel indicies on a HEALPix grid with resolution\n given by `nside`.\n\n Args:\n freq: Delta frequency/wavelength or a sequence of frequencies corresponding to\n a bandpass over which to evaluate the zodiacal emission. The frequencies\n must be strictly increasing.\n pixels: HEALPix pixel indicies representing points on the celestial sphere.\n nside: HEALPix resolution parameter of the pixels and the binned map.\n obs_time: Time of observation.\n obs: Name of the Solar System observer. A list of all support observers (for a\n given ephemeridis) is specified in `supported_observers` attribute of the\n `Zodipy` instance. Defaults to 'earth'.\n obs_pos: The heliocentric ecliptic cartesian position of the observer in AU.\n Overrides the `obs` argument. Default is None.\n weights: Bandpass weights corresponding the the frequencies in `freq`. The weights\n are assumed to be given in spectral radiance units (Jy/sr).\n return_comps: If True, the emission is returned component-wise. Defaults to False.\n coord_in: Coordinate frame of the input pointing. Assumes 'E' (ecliptic\n coordinates) by default.\n\n Returns:\n emission: Simulated zodiacal emission in units of 'MJy/sr'.\n\n \"\"\"\n pixels = get_validated_pix(pixels=pixels, nside=nside)\n\n unique_pixels, indicies = np.unique(pixels, return_inverse=True)\n unit_vectors = get_unit_vectors_from_pixels(\n coord_in=coord_in,\n pixels=unique_pixels,\n nside=nside,\n )\n\n return self._compute_emission(\n freq=freq,\n weights=weights,\n obs=obs,\n obs_time=obs_time,\n obs_pos=obs_pos,\n unit_vectors=unit_vectors,\n indicies=indicies,\n pixels=unique_pixels,\n nside=nside,\n return_comps=return_comps,\n )\nget_parameters()","text":"Return a dictionary containing the interplanetary dust model parameters.
Source code inzodipy/zodipy.py def get_parameters(self) -> ParameterDict:\n\"\"\"Return a dictionary containing the interplanetary dust model parameters.\"\"\"\n return self._ipd_model.to_dict()\nupdate_parameters(parameters)","text":"Update the interplanetary dust model parameters.
Parameters:
Name Type Description Defaultparameters ParameterDict Dictionary of parameters to update. The keys must be the names of the parameters as defined in the model. To get the parameters dict of an existing model, use Zodipy(\"dirbe\").get_parameters().
zodipy/zodipy.py def update_parameters(self, parameters: ParameterDict) -> None:\n\"\"\"Update the interplanetary dust model parameters.\n\n Args:\n parameters: Dictionary of parameters to update. The keys must be the names\n of the parameters as defined in the model. To get the parameters dict\n of an existing model, use `Zodipy(\"dirbe\").get_parameters()`.\n\n \"\"\"\n _dict = parameters.copy()\n _dict[\"comps\"] = {}\n for key, value in parameters.items():\n if key == \"comps\":\n for comp_key, comp_value in value.items():\n _dict[\"comps\"][ComponentLabel(comp_key)] = type(\n self._ipd_model.comps[ComponentLabel(comp_key)]\n )(**comp_value)\n elif isinstance(value, dict):\n _dict[key] = {ComponentLabel(k): v for k, v in value.items()}\n\n self._ipd_model = self._ipd_model.__class__(**_dict)\nBelow we illustrate how ZodiPy can be used to create timestreams of the zodiacal emission. Note that since ZodiPy assumes a constant observer position over the input pointing sequence, the output will not be real timestreams, but for small enough time intervals the error is negligible.
"},{"location":"usage/#emission-along-a-meridian","title":"Emission along a meridian","text":"In the following example we simulate what an observer on Earth is expected to see on 14 June, 2022 when looking along a meridian (line of constant longitude) at 30 microns, given the DIRBE interplanetary dust model.
import astropy.units as u\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom astropy.time import Time\n\nfrom zodipy import Zodipy\n\nmodel = Zodipy(\"dirbe\")\n\nlatitudes = np.linspace(-90, 90, 10000) * u.deg\nlongitudes = np.zeros_like(latitudes)\n\nemission = model.get_emission_ang(\n 30 * u.micron,\n theta=longitudes,\n phi=latitudes,\n lonlat=True,\n obs_time=Time(\"2022-06-14\"),\n obs=\"earth\",\n)\n\n\nplt.plot(latitudes, emission)\nplt.xlabel(\"Latitude [deg]\")\nplt.ylabel(\"Emission [MJy/sr]\")\nplt.savefig(\"../img/timestream.png\", dpi=300)\nplt.show()\nNote
ZodiPy assumes a constant observer position over an input pointing sequence. For an observer on Earth, the true zodiacal emission signal will move along the ecliptic on the sky by roughly one degree each day. To account for this effect, the full pointing sequence of an experiment should be chunked into small subsequences with timescales corresponding to at maximum a day.
"},{"location":"usage/#healpix-maps","title":"HEALPix maps","text":"Below we illustrate how ZodiPy can be used to create simulated binned HEALPix maps of the zodiacal emission.
"},{"location":"usage/#instantaneous-map-in-ecliptic-coordinates","title":"Instantaneous map in ecliptic coordinates","text":"In the following example we make an instantaneous map of of the zodiacal emission at 857 GHz as seen by an observer on earth on 14 June, 2022 given the Planck 2018 interplanetary dust model.
import astropy.units as u\nimport healpy as hp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom astropy.time import Time\n\nfrom zodipy import Zodipy\n\nmodel = Zodipy(\"planck18\")\nnside = 256\n\nbinned_emission = model.get_binned_emission_pix(\n 857 * u.GHz,\n pixels=np.arange(hp.nside2npix(nside)),\n nside=nside,\n obs_time=Time(\"2022-06-14\"),\n obs=\"earth\",\n)\n\nhp.mollview(\n binned_emission,\n title=\"Binned zodiacal emission at 857 GHz\",\n unit=\"MJy/sr\",\n min=0,\n max=1,\n cmap=\"afmhot\",\n)\nplt.savefig(\"../img/binned.png\", dpi=300)\nplt.show()\nInstruments 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.
import astropy.units as u\nimport healpy as hp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom astropy.time import Time\n\nfrom zodipy import Zodipy\n\nnside = 128\n\ncenter_freq = 800 * u.GHz\nfreqs = np.linspace(750, 850, 11) * u.GHz\nweights = np.array([2, 3, 5, 9, 11, 11.5, 11, 9, 5, 3, 2])\nplt.plot(freqs, weights)\nplt.xlabel(\"Frequency [GHz]\")\nplt.ylabel(\"Weights\")\nplt.savefig(\"../img/bandpass.png\", dpi=300)\n\nmodel = Zodipy(model=\"planck18\")\n\nemission_central_freq = model.get_binned_emission_pix(\n freq=center_freq,\n pixels=np.arange(hp.nside2npix(nside)),\n nside=nside,\n obs_time=Time(\"2022-03-10\"),\n obs=\"SEMB-L2\",\n)\n\nemission_bandpass_integrated = model.get_binned_emission_pix(\nfreq=freqs,\nweights=weights,\npixels=np.arange(hp.nside2npix(nside)),\n nside=nside,\n obs_time=Time(\"2022-03-10\"),\n obs=\"SEMB-L2\",\n)\n\nhp.mollview(\n emission_central_freq,\n title=f\"Center frequency\",\n unit=\"MJy/sr\",\n cmap=\"afmhot\",\n norm=\"log\",\n)\nplt.savefig(\"../img/center_freq.png\", dpi=300)\n\nhp.mollview(\n emission_bandpass_integrated,\n title=\"Bandpass integrated\",\n unit=\"MJy/sr\",\n cmap=\"afmhot\",\n norm=\"log\",\n)\nplt.savefig(\"../img/bandpass_integrated.png\", dpi=300)\nplt.show()\nFew 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 than the solar_cut value.
import astropy.units as u\nimport healpy as hp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom astropy.time import Time\n\nfrom zodipy import Zodipy\n\nmodel = Zodipy(\"dirbe\", solar_cut=60 * u.deg)\nnside = 256\n\nbinned_emission = model.get_binned_emission_pix(\n 25 * u.micron,\n pixels=np.arange(hp.nside2npix(nside)),\n nside=nside,\n obs_time=Time(\"2020-01-01\"),\n obs=\"earth\",\n)\n\nhp.mollview(\n binned_emission,\n title=\"Solar cutoff at 60 degrees\",\n unit=\"MJy/sr\",\n max=80,\n coord=\"E\",\n cmap=\"afmhot\",\n)\nplt.savefig(\"../img/binned_solar_cutoff.png\", dpi=300)\nplt.show()\nWe can specify the coordinate system of the input pointing with the coord_in keyword
import astropy.units as u\nimport healpy as hp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom astropy.time import Time\n\nfrom zodipy import Zodipy\n\nmodel = Zodipy(\"planck18\")\nnside = 256\n\nbinned_emission = model.get_binned_emission_pix(\n 857 * u.GHz,\n pixels=np.arange(hp.nside2npix(nside)),\n nside=nside,\n obs_time=Time(\"2022-02-20\"),\n obs=\"earth\",\ncoord_in=\"G\", # Coordinates of the input pointing\n)\n\nhp.mollview(\n binned_emission,\n title=\"Binned zodiacal emission at 857 GHz\",\n unit=\"MJy/sr\",\n cmap=\"afmhot\",\n min=0,\n max=1,\n)\nplt.savefig(\"../img/binned_gal.png\", dpi=300)\nplt.show()\nZodiPy can also return the zodiacal emission component-wise. In the following example we use the DIRBE model since the later Planck models excluded the circumsolar-ring and Earth-trailing feature components. For more information on the interplanetary dust models, please read Cosmoglobe: Simulating Zodiacal Emission with ZodiPy.
import astropy.units as u\nimport healpy as hp\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom astropy.time import Time\n\nfrom zodipy import Zodipy\n\nmodel = Zodipy(\"dirbe\")\nnside = 256\n\nbinned_emission = model.get_binned_emission_pix(\n 25 * u.micron,\n pixels=np.arange(hp.nside2npix(nside)),\n nside=nside,\n obs_time=Time(\"2022-01-01\"),\n obs=\"earth\",\nreturn_comps=True,\n)\nfig = plt.figure(figsize=(8, 6.5), constrained_layout=True)\ncomps = [\"Cloud\", \"Band1\", \"Band2\", \"Band3\", \"Ring\", \"Feature\"]\nfor idx, binned_comp_emission in enumerate(binned_emission):\n hp.mollview(\n binned_comp_emission,\n title=comps[idx],\n norm=\"log\" if idx == 0 else None,\n cmap=\"afmhot\",\n cbar=False,\n sub=(3, 2, idx + 1),\n fig=fig,\n )\n# plt.savefig(\"../img/binned_comp.png\", dpi=300)\nplt.show()\nIf you are not using ZodiPy in an already parallelized environment, you may specify the number of cores used by ZodiPy through the n_proc keyword. By default n_proc is set to 1. For values of n_proc > 1, the line-of-sight calculations are parallelized using the multiprocessing module.
import time\n\nimport astropy.units as u\nimport healpy as hp\nimport numpy as np\nfrom astropy.time import Time\n\nfrom zodipy import Zodipy\n\nnside = 256\npixels = np.arange(hp.nside2npix(nside))\nobs_time = Time(\"2020-01-01\")\nn_proc = 8\n\nmodel = Zodipy()\nmodel_parallel = Zodipy(n_proc=n_proc)\nstart = time.perf_counter()\nemission = model.get_binned_emission_pix(\n 40 * u.micron,\n pixels=pixels,\n nside=nside,\n obs_time=obs_time,\n)\nprint(\"Time spent on a single CPU:\", round(time.perf_counter() - start, 2), \"seconds\")\n# > Time spent on a single CPU: 35.23 seconds\n\nstart = time.perf_counter()\nemission_parallel = model_parallel.get_binned_emission_pix(\n 40 * u.micron,\n pixels=pixels,\n nside=nside,\n obs_time=obs_time,\n)\nprint(\n f\"Time spent on {n_proc} CPUs:\",\n round(time.perf_counter() - start, 2),\n \"seconds\",\n)\n# > Time spent on 8 CPUs: 12.85 seconds\n\nassert np.allclose(emission, emission_parallel)\nWindows users
Windows users must make sure to wrap the get_*_emission_* function calls in a if __name__ == \"__main__\" guard to avoid spawning infinite processes:
...\nif __name__ == \"__main__\":\n emission = model.get_emission_pix(\n ...\n )\nUsing ZodiPy in parallelized environments
If ZodiPy is used in a parallelized environment one may have to specifically set the environment variable OMP_NUM_THREADS=1 to avoid oversubscription. This is due automatic parallelization in third party libraries such as healpy where for instance the hp.Rotator object automatically parallelizes rotation of unit vectors. This means that when using ZodiPy with pointing in a coordinate system other than ecliptic, even if Zodipy is initialized with n_proc=1, healpy will under the hood automatically distribute the pointing to available CPU's.
It is possible to visualize the three-dimensional interplanetary dust distribution of the models used in ZodiPy by using the tabulate_density function which takes in a interplanetary dust model and a custom grid.
In the following example we tabulate the density distribution of the DIRBE interplanetary dust model and plot the cross section of the diffuse cloud components density in the yz-plane.
import matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib.colors import LogNorm\n\nfrom zodipy import tabulate_density\n\nN = 200\n\nx = np.linspace(-5, 5, N) # x-plane\ny = np.linspace(-5, 5, N) # y-plane\nz = np.linspace(-2, 2, N) # z-plane\n\ngrid = np.asarray(np.meshgrid(x, y, z))\ndensity_grid = tabulate_density(grid, model=\"dirbe\")\ndensity_grid = density_grid.sum(axis=0) # Sum over all components\n\nplt.pcolormesh(\n x,\n y,\n density_grid[N // 2].T, # cross section in the yz-plane\n cmap=\"afmhot\",\n norm=LogNorm(vmin=density_grid.min(), vmax=density_grid.max()),\n shading=\"gouraud\",\n rasterized=True,\n)\nplt.title(\"Cross section of the interplanetary dust density (yz-plane)\")\nplt.xlabel(\"x [AU]\")\nplt.ylabel(\"z [AU]\")\n# plt.savefig(\"../img/density_grid.png\", dpi=300)\nplt.show()\nBelow we illustrate how ZodiPy can be used to create timestreams of the zodiacal emission. +Note that since ZodiPy assumes a constant observer position over the input pointing sequence, the output +will not be real timestreams, but for small enough time intervals the error is negligible.
+In the following example we simulate what an observer on Earth is expected to see on 14 June, +2022 when looking along a meridian (line of constant longitude) at 30 microns, given the +DIRBE interplanetary dust model.
+import astropy.units as u
+import matplotlib.pyplot as plt
+import numpy as np
+from astropy.time import Time
+
+from zodipy import Zodipy
+
+model = Zodipy("dirbe")
+
+latitudes = np.linspace(-90, 90, 10000) * u.deg
+longitudes = np.zeros_like(latitudes)
+
+emission = model.get_emission_ang(
+ 30 * u.micron,
+ theta=longitudes,
+ phi=latitudes,
+ lonlat=True,
+ obs_time=Time("2022-06-14"),
+ obs="earth",
+)
+
+
+plt.plot(latitudes, emission)
+plt.xlabel("Latitude [deg]")
+plt.ylabel("Emission [MJy/sr]")
+plt.savefig("../img/timestream.png", dpi=300)
+plt.show()
+
Note
+ZodiPy assumes a constant observer position over an input pointing sequence. For an observer on Earth, +the true zodiacal emission signal will move along the ecliptic on the sky by roughly one degree each day. +To account for this effect, the full pointing sequence of an experiment should be chunked into small +subsequences with timescales corresponding to at maximum a day.
+Below we illustrate how ZodiPy can be used to create simulated binned HEALPix maps of the zodiacal emission.
+In the following example we make an instantaneous map of of the zodiacal emission at 857 GHz +as seen by an observer on earth on 14 June, 2022 given the Planck 2018 interplanetary dust model.
+import astropy.units as u
+import healpy as hp
+import matplotlib.pyplot as plt
+import numpy as np
+from astropy.time import Time
+
+from zodipy import Zodipy
+
+model = Zodipy("planck18")
+nside = 256
+
+binned_emission = model.get_binned_emission_pix(
+ 857 * u.GHz,
+ pixels=np.arange(hp.nside2npix(nside)),
+ nside=nside,
+ obs_time=Time("2022-06-14"),
+ obs="earth",
+)
+
+hp.mollview(
+ binned_emission,
+ title="Binned zodiacal emission at 857 GHz",
+ unit="MJy/sr",
+ min=0,
+ max=1,
+ cmap="afmhot",
+)
+plt.savefig("../img/binned.png", dpi=300)
+plt.show()
+
+Note that the color bar is logarithmic.
+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.
+
import astropy.units as u
+import healpy as hp
+import matplotlib.pyplot as plt
+import numpy as np
+from astropy.time import Time
+
+from zodipy import Zodipy
+
+nside = 128
+
+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 [GHz]")
+plt.ylabel("Weights")
+plt.savefig("../img/bandpass.png", dpi=300)
+
+model = Zodipy(model="planck18")
+
+emission_central_freq = model.get_binned_emission_pix(
+ freq=center_freq,
+ pixels=np.arange(hp.nside2npix(nside)),
+ nside=nside,
+ obs_time=Time("2022-03-10"),
+ obs="SEMB-L2",
+)
+
+emission_bandpass_integrated = model.get_binned_emission_pix(
+ freq=freqs,
+ weights=weights,
+ pixels=np.arange(hp.nside2npix(nside)),
+ nside=nside,
+ obs_time=Time("2022-03-10"),
+ obs="SEMB-L2",
+)
+
+hp.mollview(
+ emission_central_freq,
+ title=f"Center frequency",
+ unit="MJy/sr",
+ cmap="afmhot",
+ norm="log",
+)
+plt.savefig("../img/center_freq.png", dpi=300)
+
+hp.mollview(
+ emission_bandpass_integrated,
+ title="Bandpass integrated",
+ unit="MJy/sr",
+ cmap="afmhot",
+ norm="log",
+)
+plt.savefig("../img/bandpass_integrated.png", dpi=300)
+plt.show()
+
+
+
+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
+than the solar_cut value.
import astropy.units as u
+import healpy as hp
+import matplotlib.pyplot as plt
+import numpy as np
+from astropy.time import Time
+
+from zodipy import Zodipy
+
+model = Zodipy("dirbe", solar_cut=60 * u.deg)
+nside = 256
+
+binned_emission = model.get_binned_emission_pix(
+ 25 * u.micron,
+ pixels=np.arange(hp.nside2npix(nside)),
+ nside=nside,
+ obs_time=Time("2020-01-01"),
+ obs="earth",
+)
+
+hp.mollview(
+ binned_emission,
+ title="Solar cutoff at 60 degrees",
+ unit="MJy/sr",
+ max=80,
+ coord="E",
+ cmap="afmhot",
+)
+plt.savefig("../img/binned_solar_cutoff.png", dpi=300)
+plt.show()
+
+We can specify the coordinate system of the input pointing with the coord_in keyword
import astropy.units as u
+import healpy as hp
+import matplotlib.pyplot as plt
+import numpy as np
+from astropy.time import Time
+
+from zodipy import Zodipy
+
+model = Zodipy("planck18")
+nside = 256
+
+binned_emission = model.get_binned_emission_pix(
+ 857 * u.GHz,
+ pixels=np.arange(hp.nside2npix(nside)),
+ nside=nside,
+ 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",
+ unit="MJy/sr",
+ cmap="afmhot",
+ min=0,
+ max=1,
+)
+plt.savefig("../img/binned_gal.png", dpi=300)
+plt.show()
+
+ZodiPy can also return the zodiacal emission component-wise. In the following example we use +the DIRBE model since the later Planck models excluded the circumsolar-ring and Earth-trailing +feature components. For more information on the interplanetary dust models, please +read Cosmoglobe: Simulating Zodiacal Emission with ZodiPy.
+import astropy.units as u
+import healpy as hp
+import matplotlib.pyplot as plt
+import numpy as np
+from astropy.time import Time
+
+from zodipy import Zodipy
+
+model = Zodipy("dirbe")
+nside = 256
+
+binned_emission = model.get_binned_emission_pix(
+ 25 * u.micron,
+ pixels=np.arange(hp.nside2npix(nside)),
+ nside=nside,
+ obs_time=Time("2022-01-01"),
+ obs="earth",
+ return_comps=True,
+)
+fig = plt.figure(figsize=(8, 6.5), constrained_layout=True)
+comps = ["Cloud", "Band1", "Band2", "Band3", "Ring", "Feature"]
+for idx, binned_comp_emission in enumerate(binned_emission):
+ hp.mollview(
+ binned_comp_emission,
+ title=comps[idx],
+ norm="log" if idx == 0 else None,
+ cmap="afmhot",
+ cbar=False,
+ sub=(3, 2, idx + 1),
+ fig=fig,
+ )
+# plt.savefig("../img/binned_comp.png", dpi=300)
+plt.show()
+
+Note that the color for the Cloud component is logarithmic, while the others are linear.
+If you are not using ZodiPy in an already parallelized environment, you may specify the number of cores used by ZodiPy through the n_proc keyword. By default n_proc is set to 1. For values of n_proc > 1, the line-of-sight calculations are parallelized using the multiprocessing module.
import time
+
+import astropy.units as u
+import healpy as hp
+import numpy as np
+from astropy.time import Time
+
+from zodipy import Zodipy
+
+nside = 256
+pixels = np.arange(hp.nside2npix(nside))
+obs_time = Time("2020-01-01")
+n_proc = 8
+
+model = Zodipy()
+model_parallel = Zodipy(n_proc=n_proc)
+
+start = time.perf_counter()
+emission = model.get_binned_emission_pix(
+ 40 * u.micron,
+ pixels=pixels,
+ nside=nside,
+ obs_time=obs_time,
+)
+print("Time spent on a single CPU:", round(time.perf_counter() - start, 2), "seconds")
+# > Time spent on a single CPU: 35.23 seconds
+
+start = time.perf_counter()
+emission_parallel = model_parallel.get_binned_emission_pix(
+ 40 * u.micron,
+ pixels=pixels,
+ nside=nside,
+ obs_time=obs_time,
+)
+print(
+ f"Time spent on {n_proc} CPUs:",
+ round(time.perf_counter() - start, 2),
+ "seconds",
+)
+# > Time spent on 8 CPUs: 12.85 seconds
+
+assert np.allclose(emission, emission_parallel)
+Windows users
+Windows users must make sure to wrap the get_*_emission_* function calls in a if __name__ == "__main__" guard to avoid spawning infinite processes:
+
...
+if __name__ == "__main__":
+ emission = model.get_emission_pix(
+ ...
+ )
+Using ZodiPy in parallelized environments
+If ZodiPy is used in a parallelized environment one may have to specifically set the environment variable
+OMP_NUM_THREADS=1 to avoid oversubscription. This is due automatic parallelization in third party libraries such as healpy where for instance the hp.Rotator object automatically parallelizes rotation of unit vectors.
+This means that when using ZodiPy with pointing in a coordinate system other than ecliptic, even if Zodipy is initialized with n_proc=1, healpy will under the hood automatically distribute the pointing to available CPU's.
It is possible to visualize the three-dimensional interplanetary dust distribution of the models used in
+ZodiPy by using the tabulate_density function which takes in a interplanetary dust model and a custom grid.
In the following example we tabulate the density distribution of the DIRBE interplanetary dust model +and plot the cross section of the diffuse cloud components density in the yz-plane.
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.colors import LogNorm
+
+from zodipy import tabulate_density
+
+N = 200
+
+x = np.linspace(-5, 5, N) # x-plane
+y = np.linspace(-5, 5, N) # y-plane
+z = np.linspace(-2, 2, N) # z-plane
+
+grid = np.asarray(np.meshgrid(x, y, z))
+density_grid = tabulate_density(grid, model="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 interplanetary dust density (yz-plane)")
+plt.xlabel("x [AU]")
+plt.ylabel("z [AU]")
+# plt.savefig("../img/density_grid.png", dpi=300)
+plt.show()
+
+
+
+
+
+
+
+