Priors Module (brutus.priors)

brutus.priors implements the astrophysical prior probability distributions used by BruteForce to break degeneracies between stellar luminosity, distance, and dust reddening.

Role in Inference:

Because broadband photometry alone often permits both a nearby cool dwarf and a distant reddened giant, brutus incorporates prior distributions over:

  • Stellar Mass and Luminosity: Kroupa initial mass function (logp_imf) and Pan-STARRS \(r\)-band luminosity function (logp_ps1_luminosity_function).

  • Astrometry: Gaussian parallax likelihood (logp_parallax) and scale-factor marginalization (logp_parallax_scale).

  • Galactic Structure: Three-component thin-disk, thick-disk, and halo number density, metallicity, and age distributions (logp_galactic_structure, logn_disk, logn_halo, logp_feh, logp_age_from_feh).

  • Interstellar Extinction: Truncated Gaussian \(A_V\) prior from 3D dust maps (logp_extinction).

Typical Usage:

Priors are evaluated automatically inside BruteForce.fit, and can also be called directly:

from brutus.priors.stellar import logp_imf
from brutus.priors.galactic import logp_galactic_structure
from brutus.priors.extinction import logp_extinction
from brutus.dust import Bayestar
import numpy as np

# Evaluate IMF prior for a range of initial masses (in Msun)
masses = np.array([0.5, 1.0, 2.0, 5.0])
log_prior_imf = logp_imf(masses)

# Evaluate Galactic structure prior along a sightline
distances = np.array([0.1, 0.5, 1.0, 5.0])  # kpc
gal_l, gal_b = 45.0, 10.0  # Galactic coordinates (deg)
log_prior_gal = logp_galactic_structure(distances, (gal_l, gal_b))

# Evaluate extinction prior from a 3D dust map
dustmap = Bayestar()
av_values = np.array([0.1, 0.5, 1.0, 2.0])  # mag
log_prior_dust = logp_extinction(
    av_values, dustmap, (gal_l, gal_b), distance=distances
)

Customization:

Pass custom callables to BruteForce.fit to replace or disable specific priors. The Galactic structure prior is supplied via lngalprior(dist, coord, labels=...), where coord is a (l, b) tuple in degrees:

from brutus.analysis import BruteForce

fitter = BruteForce(grid)

# Replace the Galactic structure prior with a flat log-prior in distance
fitter.fit(
    data=flux, data_err=flux_err, data_mask=mask,
    data_labels=obj_ids, save_file="results.h5",
    data_coords=coords,
    lngalprior=lambda dist, coord, labels=None: np.zeros_like(dist, dtype=float),
)

See Also:

  • Priors - Detailed conceptual guide to all prior distributions

  • How brutus works - How priors fit into the Bayesian framework

Stellar Priors

brutus.priors.logp_imf(mgrid, alpha_low=1.3, alpha_high=2.3, mass_break=0.5, mgrid2=None, mass_min=0.08, mass_max=100.0)[source]

Compute the log-prior for a Kroupa-like broken power-law initial mass function.

Evaluates a two-part power-law IMF following Kroupa (2001), with optional independent draws for a binary companion.

Parameters:

  • mgrid (array_like) – Initial stellar masses in solar masses (\(M_\odot\)).

  • alpha_low (float, optional) – Power-law slope for mgrid <= mass_break. Default is 1.3.

  • alpha_high (float, optional) – Power-law slope for mgrid > mass_break. Default is 2.3.

  • mass_break (float, optional) – Transition mass between the low- and high-mass regimes in \(M_\odot\). Default is 0.5.

  • mgrid2 (array_like, optional) – Secondary stellar masses for binary systems in \(M_\odot\), broadcastable with mgrid.

  • mass_min (float, optional) – Minimum mass for normalization in \(M_\odot\). Default is 0.08.

  • mass_max (float, optional) – Maximum mass for normalization in \(M_\odot\). May be np.inf when alpha_high > 1. Default is 100.0.

Returns:

logp – Normalized log-prior density (per \(M_\odot\), or per \(M_\odot^2\) when mgrid2 is given) with the shape of mgrid. Masses outside [mass_min, mass_max] receive LOG_ZERO (-1e300).

Return type:

numpy.ndarray

Raises:

ValueError – If mass_min or mass_break is not finite and positive, if mass_max does not exceed mass_min, or if a slope is not finite.

See also

logp_ps1_luminosity_function

Alternative luminosity-based prior

brutus.analysis.individual.BruteForce

Uses IMF priors for fitting

Notes

The IMF follows the form:

\[\xi(M) \propto M^{-\alpha}\]

where α = α_low for M ≤ M_break and α = α_high for M > M_break.

For binary systems, assumes independent sampling from the same IMF for both components.

References

Kroupa, P. (2001), MNRAS, 322, 231

brutus.priors.logp_ps1_luminosity_function(Mr)[source]

Evaluate the Pan-STARRS 1 r-band luminosity function log-prior.

Interpolates the empirical stellar luminosity function tabulated for Pan-STARRS 1 \(M_r\) in [-1, 20] mag (Green et al. 2015, 2018) and extrapolates linearly outside that interval.

Parameters:

Mr (array_like) – Absolute \(r\)-band magnitudes in the Pan-STARRS 1 system (mag).

Returns:

logp – Unnormalized log-prior values with the shape of Mr (exp(logp) integrates to approximately 0.086 over the tabulated range).

Return type:

numpy.ndarray

See also

logp_imf

Mass-based IMF prior

brutus.analysis.individual.BruteForce

Uses luminosity priors for fitting

References

Green, G. M., et al. (2015), ApJ, 810, 25 Green, G. M., et al. (2018), MNRAS, 478, 651

Astrometric Priors

brutus.priors.logp_parallax(parallaxes, p_meas, p_err)[source]

Compute the Gaussian log-prior for parallax measurements.

Parameters:

  • parallaxes (array_like) – Model parallax values in milliarcseconds (mas).

  • p_meas (float) – Measured parallax in milliarcseconds (mas).

  • p_err (float) – Parallax measurement uncertainty in milliarcseconds (mas).

Returns:

logp – Log-prior probability density (per mas) with the shape of parallaxes. Returns zeros when p_meas or p_err is non-finite or p_err <= 0.

Return type:

numpy.ndarray

Notes

For valid measurements, the log-prior is

\[\log p(\pi \mid \pi_{\text{obs}}, \sigma_\pi) = -\frac{1}{2} \left[ \frac{(\pi - \pi_{\text{obs}})^2}{\sigma_\pi^2} + \log(2\pi\sigma_\pi^2) \right]\]
brutus.priors.logp_parallax_scale(scales, scale_errs, p_meas, p_err, snr_lim=4.0)[source]

Compute the Gaussian log-prior on flux scale factors from a parallax measurement.

Approximates the distribution of the flux scale factor \(s = 1 / d_{\text{kpc}}^2 = \pi_{\text{mas}}^2\) as a Gaussian when the parallax signal-to-noise ratio exceeds snr_lim.

Parameters:

  • scales (array_like) – Flux density scale factors (\(\text{kpc}^{-2}\), equal to \(\text{mas}^2\)).

  • scale_errs (array_like) – Uncertainties on scales (\(\text{kpc}^{-2}\)), broadcastable with scales.

  • p_meas (float) – Measured parallax in milliarcseconds (mas).

  • p_err (float) – Parallax measurement uncertainty in milliarcseconds (mas).

  • snr_lim (float, optional) – Minimum |p_meas / p_err| required to apply the Gaussian approximation. Below this threshold, zeros are returned. Default is 4.0.

Returns:

logp – Log-prior probability density with the shape of scales.

Return type:

numpy.ndarray

Notes

When |p_meas / p_err| > snr_lim, moments of \(s = \pi^2\) are computed via convert_parallax_to_scale() and combined in quadrature with scale_errs:

\[s_{\text{mean}} = \max(0, \pi_{\text{obs}})^2 + \sigma_\pi^2, \quad \sigma_s = \sqrt{2\sigma_\pi^4 + 4\max(0, \pi_{\text{obs}})^2\sigma_\pi^2}\]
brutus.priors.convert_parallax_to_scale(p_meas, p_err, snr_lim=4.0)[source]

Convert a parallax measurement to flux scale-factor moments.

Transforms a parallax measurement and uncertainty into the mean and standard deviation of the scale factor \(s = \pi^2\) via second-order moment propagation.

Parameters:

  • p_meas (float) – Measured parallax in milliarcseconds (mas).

  • p_err (float) – Parallax measurement uncertainty in milliarcseconds (mas).

  • snr_lim (float, optional) – Minimum |p_meas / p_err| required for the Gaussian conversion. Below this threshold, uninformative bounds (SCALE_MIN, SCALE_MAX) are returned. Default is 4.0.

Returns:

  • s_mean (float) – Mean of the scale factor distribution (\(\text{kpc}^{-2}\)).

  • s_std (float) – Standard deviation of the scale factor distribution (\(\text{kpc}^{-2}\)).

Notes

When |p_meas / p_err| > snr_lim, the moments are

\[ \begin{align}\begin{aligned}s_{\text{mean}} = \max(0, \pi_{\text{meas}})^2 + \sigma_\pi^2\\s_{\text{std}} = \sqrt{2\sigma_\pi^4 + 4\max(0, \pi_{\text{meas}})^2\sigma_\pi^2}\end{aligned}\end{align} \]

Examples

>>> s_mean, s_std = convert_parallax_to_scale(1.0, 0.1)
>>> print(f"Scale factor: {s_mean:.3f} +/- {s_std:.3f}")

Galactic Structure Priors

brutus.priors.logp_galactic_structure(dists, coord, labels=None, feh=None, loga=None, R_solar=8.2, Z_solar=0.025, R_thin=2.6, Z_thin=0.3, Rs_thin=2.0, R_thick=2.0, Z_thick=0.9, f_thick=0.04, Rs_thick=2.0, Rs_halo=2.0, q_halo_ctr=0.2, q_halo_inf=0.8, r_q_halo=6.0, eta_halo=4.2, f_halo=0.005, feh_thin=-0.2, feh_thin_sigma=0.3, feh_thick=-0.7, feh_thick_sigma=0.4, feh_halo=-1.6, feh_halo_sigma=0.5, max_age=13.8, min_age=0.0, feh_age_ctr=-0.5, feh_age_scale=0.5, nsigma_from_max_age=2.0, max_sigma=4.0, min_sigma=1.0, return_components=False)[source]

Compute the three-component Galactic structure log-prior.

Combines spatial number-density profiles from Bland-Hawthorn & Gerhard (2016) with volume factors and optional conditional metallicity and age priors.

Parameters:

  • dists (array_like) – Distance from observer in kpc.

  • coord (tuple of floats) – Galactic coordinates (l, b) in degrees.

  • labels (structured array, optional) – Stellar labels containing ‘feh’ and/or ‘loga’ for metallicity/age priors, with one entry per distance.

  • feh (array_like, optional) – Per-point [Fe/H] as a plain float array, an alternative to labels for the metallicity prior. Lets callers (e.g. logpost_grid) avoid tiling a structured labels array; equivalent to passing the same values via labels. Honored on both code paths (fused kernel and numpy). If both labels and feh/loga are given, the plain arrays take precedence.

  • loga (array_like, optional) – Per-point log10(age/yr) as a plain float array, an alternative to labels for the age prior (see feh).

  • R_solar (float, optional) – Solar Galactocentric radius in kpc. Default is 8.2.

  • Z_solar (float, optional) – Solar height above midplane in kpc. Default is 0.025.

  • R_thin (float, optional) – Thin disk radial scale length in kpc. Default is 2.6.

  • Z_thin (float, optional) – Thin disk vertical scale height in kpc. Default is 0.3.

  • Rs_thin (float, optional) – Thin disk smoothing radius in kpc. Default is 2.0.

  • R_thick (float, optional) – Thick disk radial scale length in kpc. Default is 2.0.

  • Z_thick (float, optional) – Thick disk vertical scale height in kpc. Default is 0.9.

  • f_thick (float, optional) – Thick disk relative normalization. Default is 0.04.

  • Rs_thick (float, optional) – Thick disk smoothing radius in kpc. Default is 2.0.

  • Rs_halo (float, optional) – Halo smoothing radius in kpc. Default is 2.0.

  • q_halo_ctr (float, optional) – Halo central oblateness. Default is 0.2.

  • q_halo_inf (float, optional) – Halo asymptotic oblateness. Default is 0.8.

  • r_q_halo (float, optional) – Halo oblateness transition radius in kpc. Default is 6.0.

  • eta_halo (float, optional) – Halo power-law index. Default is 4.2.

  • f_halo (float, optional) – Halo relative normalization. Default is 0.005.

  • feh_thin (float, optional) – Thin disk mean metallicity in dex. Default is -0.2.

  • feh_thin_sigma (float, optional) – Thin disk metallicity dispersion in dex. Default is 0.3.

  • feh_thick (float, optional) – Thick disk mean metallicity in dex. Default is -0.7.

  • feh_thick_sigma (float, optional) – Thick disk metallicity dispersion in dex. Default is 0.4.

  • feh_halo (float, optional) – Halo mean metallicity in dex. Default is -1.6.

  • feh_halo_sigma (float, optional) – Halo metallicity dispersion in dex. Default is 0.5.

  • max_age (float, optional) – Maximum stellar age in Gyr. Default is 13.8.

  • min_age (float, optional) – Minimum stellar age in Gyr. Default is 0.0.

  • feh_age_ctr (float, optional) – Central metallicity for age-metallicity relation. Default is -0.5.

  • feh_age_scale (float, optional) – Scale length for age-metallicity relation. Default is 0.5.

  • nsigma_from_max_age (float, optional) – Age dispersion parameter. Default is 2.0.

  • max_sigma (float, optional) – Maximum age dispersion in Gyr. Default is 4.0.

  • min_sigma (float, optional) – Minimum age dispersion in Gyr. Default is 1.0.

  • return_components (bool, optional) – Whether to return individual component contributions. Default is False.

Returns:

  • logp (array_like) – Total log-prior probability density.

  • components (dict, optional) – Individual component contributions (if return_components=True).

Raises:

ValueError – If labels, feh or loga does not have one entry per distance.

Notes

The Galactic model combines three stellar populations:

  1. Thin Disk: Young, metal-rich stars with small scale height

  2. Thick Disk: Intermediate-age, metal-poor stars with larger scale height

  3. Halo: Old, very metal-poor stars with flattened power-law profile

Each component has distinct spatial, metallicity, and age distributions calibrated from observations. The model accounts for:

  • Coordinate transformations from observer to Galactocentric frame

  • Volume correction factors (dV ∝ distance²)

  • Component membership probabilities

  • Conditional metallicity and age priors

When stellar labels are provided, applies population-specific priors: - Metallicity: Gaussian distributions with different means/dispersions - Age: Age-metallicity relation with truncated normal distributions

When both metallicity and age are supplied, the joint prior evaluates the mixture sum_c w_c(pos) p(feh|c) p(age|c) by weighting the age term of each component with its metallicity-updated membership probability w_c(pos) p(feh|c) / sum_c' w_c'(pos) p(feh|c'). Nearby metal-poor stars are therefore weighted toward old, halo-like ages.

References

Bland-Hawthorn & Gerhard (2016) - The Galaxy in Context

brutus.priors.logn_disk(R, Z, R_solar=8.2, Z_solar=0.025, R_scale=2.6, Z_scale=0.3, R_smooth=2.0)[source]

Compute the log-number density for an exponential Galactic disk.

Evaluates a double-exponential disk model with separate radial and vertical scale lengths, smoothed near the Galactic center to avoid singularities.

Parameters:

  • R (array_like) – Galactocentric cylindrical radius in kpc.

  • Z (array_like) – Height above the Galactic midplane in kpc.

  • R_solar (float, optional) – Solar Galactocentric radius in kpc. Default is 8.2.

  • Z_solar (float, optional) – Solar height above midplane in kpc. Default is 0.025.

  • R_scale (float, optional) – Disk radial scale length in kpc. Default is 2.6.

  • Z_scale (float, optional) – Disk vertical scale height in kpc. Default is 0.3.

  • R_smooth (float, optional) – Smoothing radius to avoid central singularity in kpc. Default is 2.0.

Returns:

logn – Normalized log-number density relative to Solar neighborhood.

Return type:

array_like

See also

logn_halo

Halo number density

logp_galactic_structure

Combined disk+halo model

Notes

The disk number density follows:

\[n_{\text{disk}}(R, Z) \propto \exp\left(-\frac{R_{\text{eff}} - R_{\text{eff},\odot}}{R_{\text{scale}}} - \frac{|Z| - |Z_\odot|}{Z_{\text{scale}}}\right)\]

where \(R_{\text{eff}} = \sqrt{R^2 + R_{\text{smooth}}^2}\) provides smoothing near the Galactic center and \(R_{\text{eff},\odot} = \sqrt{R_\odot^2 + R_{\text{smooth}}^2}\) is the smoothed solar radius, so the normalized density is exactly 1 at the solar position (consistent with logn_halo()). This keeps the mixture fractions (f_thick, f_halo) in logp_galactic_structure() exact at the solar position.

References

Bland-Hawthorn & Gerhard (2016) - The Galaxy in Context

brutus.priors.logn_halo(R, Z, R_solar=8.2, Z_solar=0.025, R_smooth=2.0, eta=4.2, q_ctr=0.2, q_inf=0.8, r_q=6.0)[source]

Compute the log-number density for a flattened power-law Galactic halo.

Evaluates a flattened power-law halo profile with radius-dependent oblateness.

Parameters:

  • R (array_like) – Galactocentric cylindrical radius in kpc.

  • Z (array_like) – Height above the Galactic midplane in kpc.

  • R_solar (float, optional) – Solar Galactocentric radius in kpc. Default is 8.2.

  • Z_solar (float, optional) – Solar height above midplane in kpc. Default is 0.025.

  • R_smooth (float, optional) – Smoothing radius to avoid central singularity in kpc. Default is 2.0.

  • eta (float, optional) – Power-law index for halo density profile. Default is 4.2.

  • q_ctr (float, optional) – Halo oblateness at Galactic center. Default is 0.2.

  • q_inf (float, optional) – Halo oblateness at large radii. Default is 0.8.

  • r_q (float, optional) – Scale radius for oblateness transition in kpc. Default is 6.0.

Returns:

logn – Normalized log-number density relative to Solar neighborhood.

Return type:

array_like

See also

logn_disk

Disk number density

logp_galactic_structure

Combined disk+halo model

Notes

The halo follows a flattened power-law profile:

\[n_{\text{halo}}(R, Z) \propto R_{\text{eff}}^{-\eta}\]

where the effective radius includes radius-dependent flattening:

\[ \begin{align}\begin{aligned}R_{\text{eff}} = \sqrt{R^2 + (Z/q)^2 + R_{\text{smooth}}^2}\\q(r) = q_\infty - (q_\infty - q_{\text{ctr}}) e^{1 - r'/r_q}\\r' = \sqrt{r^2 + r_q^2}, \quad r = \sqrt{R^2 + Z^2}\end{aligned}\end{align} \]

References

Bland-Hawthorn & Gerhard (2016) - The Galaxy in Context Bell et al. (2008) - Stellar Halo Properties from SDSS

brutus.priors.logp_feh(feh, feh_mean=-0.2, feh_sigma=0.3)[source]

Compute the Gaussian log-prior for stellar metallicity [Fe/H].

Evaluates a normal metallicity distribution for a Galactic stellar component (thin disk, thick disk, or halo).

Parameters:

  • feh (array_like) – Stellar metallicity [Fe/H] in dex.

  • feh_mean (float, optional) – Mean metallicity of the population in dex. Default is -0.2 (thin disk).

  • feh_sigma (float, optional) – Metallicity dispersion in dex. Default is 0.3.

Returns:

logp – Normalized log-probability density for the input metallicities.

Return type:

array_like

Notes

The metallicity prior follows a normal distribution:

\[\log p([\text{Fe/H}]) = -\frac{1}{2}\left[\frac{([\text{Fe/H}] - \mu_{\text{Fe/H}})^2}{\sigma_{\text{Fe/H}}^2} + \log(2\pi\sigma_{\text{Fe/H}}^2)\right]\]

Typical values for different Galactic components: - Thin disk: feh_mean = -0.2, feh_sigma = 0.3 - Thick disk: feh_mean = -0.7, feh_sigma = 0.4 - Halo: feh_mean = -1.6, feh_sigma = 0.5

References

Bland-Hawthorn & Gerhard (2016) - The Galaxy in Context

brutus.priors.logp_age_from_feh(age, feh_mean=-0.2, max_age=13.8, min_age=0.0, feh_age_ctr=-0.5, feh_age_scale=0.5, nsigma_from_max_age=2.0, max_sigma=4.0, min_sigma=1.0)[source]

Compute the truncated-normal age log-prior from a logistic age-metallicity relation.

Evaluates a truncated normal distribution over [min_age, max_age] whose mean age follows a logistic function of mean metallicity.

Parameters:

  • age (array_like) – Stellar ages in Gyr.

  • feh_mean (float, optional) – Mean metallicity of the population in dex. Default is -0.2.

  • max_age (float, optional) – Maximum allowed stellar age in Gyr. Default is 13.8 (age of Universe).

  • min_age (float, optional) – Minimum allowed stellar age in Gyr. Default is 0.0.

  • feh_age_ctr (float, optional) – Metallicity where mean age is halfway between min/max. Default is -0.5.

  • feh_age_scale (float, optional) – Scale length for metallicity-age relation in dex. Default is 0.5.

  • nsigma_from_max_age (float, optional) – Number of σ the mean age is below max_age. Default is 2.0.

  • max_sigma (float, optional) – Maximum age dispersion in Gyr. Default is 4.0.

  • min_sigma (float, optional) – Minimum age dispersion in Gyr. Default is 1.0.

Returns:

logp – Normalized log-probability density for the input ages.

Return type:

array_like

Notes

The age-metallicity relation follows a logistic function:

\[\langle t \rangle = \frac{t_{\max} - t_{\min}}{1 + \exp\left(\frac{[\text{Fe/H}] - c}{s}\right)} + t_{\min}\]

where c is the central metallicity and s is the scale length.

The age dispersion increases for younger (more metal-rich) stars:

\[\sigma_t = \min\left(\max\left(\frac{t_{\max} - \langle t \rangle}{n\sigma}, \sigma_{\min}\right), \sigma_{\max}\right)\]

The age density is a truncated normal distribution bounded by [min_age, max_age].

References

Bland-Hawthorn & Gerhard (2016) - The Galaxy in Context Nordström et al. (2004) - Age-metallicity relation in Solar neighborhood

Extinction Priors

brutus.priors.logp_extinction(avs, dustmap, coord, distance=None, return_components=False, avlim=(0.0, 20.0), scatter=0.2)[source]

Log-prior for dust extinction using 3D dust maps.

Implements Gaussian extinction priors based on dust maps with systematic uncertainty treatment. Supports both 3D dust maps (e.g., Bayestar) that return distance-resolved profiles and simpler maps that return a single mean and standard deviation.

Parameters:

  • avs (array_like) – Extinction values (A_V) in magnitudes to evaluate prior for.

  • dustmap (object) – Dust map object with a query(coord) method. For 3D dust maps (e.g., Bayestar), query returns (distances, av_mean, av_std) with distance-resolved profiles. For simpler maps, query returns (av_mean, av_std) scalars. Objects without a query method (e.g. an unloaded file path) raise TypeError.

  • coord (SkyCoord or array_like) – Sky coordinates for dust map query. Accepts astropy.coordinates.SkyCoord or [l, b] in degrees.

  • distance (float or array_like, optional) – Distance(s) in kpc at which to evaluate the extinction prior. Required for 3D dust maps to interpolate the extinction profile. If avs and distance are both arrays, they must have the same shape and the prior is evaluated element-wise.

  • return_components (bool, optional) – If True, returns tuple (logp, (av_mean, av_err)) including the dust map statistics used. Default is False.

  • avlim (tuple of 2 floats, optional) – Support (av_min, av_max) over which the Gaussian prior is normalized (truncated normal), and over which the uniform prior of uninformative map bins is normalized. Should match the avlim of the fit. Default is AVLIM, (0, 20) (the BruteForce.fit default).

  • scatter (float, optional) – Systematic scatter (in magnitudes) added in quadrature to the map’s A(V) spread, sigma = sqrt(sigma_map**2 + scatter**2). Bayestar’s per-bin spreads are typically only ~0.01-0.05 mag and do not include sub-pixel structure or map systematics; without this floor the prior pins A(V) (and, through the A(V)-distance profile, the distance) far more tightly than the map supports. Default is 0.2, matching brutus <= 0.8.3. Set to 0 to use the raw map spreads.

Returns:

  • logp (ndarray) – Log-prior probability density for the input extinction values. Where the map is uninformative (no coverage, or distance bins masked as unreliable), returns the normalized uniform density over avlim, -ln(avlim[1] - avlim[0]).

  • components (tuple, optional) – If return_components=True, returns (av_mean, av_err) containing the dust map mean and standard deviation used.

Raises:

  • ValueError – If avlim is not two finite numbers with avlim[0] < avlim[1], or if a 3D dust map is given without distance.

  • TypeError – If dustmap has no query method (e.g. a str or pathlib.Path that was never loaded into a dust map object), or if dustmap.query returns anything but a tuple, list or array of length 2 or 3.

Notes

The log-prior follows a truncated Gaussian distribution when dust map data is available:

\[\log p(A_V | A_{V,\text{map}}, \sigma_{A_V}) = -\frac{1}{2} \left[ \frac{(A_V - A_{V,\text{map}})^2}{\sigma_{A_V}^2} + \log(2\pi\sigma_{A_V}^2) \right] - \log Z\]

where \(Z = \Phi((A_{V,\max}-\mu)/\sigma) - \Phi((A_{V,\min}-\mu)/\sigma)\) normalizes the density over the physical support avlim. Because the map mean and spread vary with distance, omitting \(Z\) would under-weight distances where the profile sits near the A(V) boundary (a bias of up to ln 2).

For 3D dust maps, the expected extinction and uncertainty at the requested distance are obtained by linear interpolation of the map’s distance-resolved profiles. Distances outside the map range use the boundary values.

The spread \(\sigma_{A_V}\) is the map spread with scatter added in quadrature.

For regions without dust map coverage (NaN values, including distance bins masked by Bayestar’s reliability mask), a uniform prior normalized over avlim is returned because an unnormalized ln p = 0 would act as a density of 1 per mag and up-weight masked distances relative to reliable ones by a unit-dependent factor. If the dust map query call itself fails, the uniform prior is returned with a RuntimeWarning.

Examples

>>> from astropy.coordinates import SkyCoord
>>> coord = SkyCoord(l=90., b=0., unit='deg', frame='galactic')
>>> # With a 3D dust map:
>>> # logp = logp_extinction([0.1, 0.5], dustmap, coord, distance=1.0)
>>> # logp, (mean, err) = logp_extinction([0.1], dustmap, coord,
>>> #                                      distance=1.0,
>>> #                                      return_components=True)