cobtools.astrometry package

Submodules

cobtools.astrometry.distance module

This module provides various classes for estimating distances from parallax.

Classes

SimpleInversion: A simple model for distance estimation by inverting parallax.
FromLiterature: Fits literature values to a gamma distribution.
XRBExponentialPriorModel: A Bayesian model using an exponential prior.

The module also includes methods for sampling distances and fitting distributions.

class cobtools.astrometry.distance.FromLiterature(x_est: float, x_lo: float, x_hi: float, conf_level: float = 0.68)[source]

Bases: object

A class for taking literature distance estimates and errors. The errors can be asymmetric. The class has a method to model the asymmetric errors with a gamma distribution, and then sample distances from the fitted distribution.

x_est

The estimated distance from the literature.

Type:: float

x_lo

The lower limit of the distance estimate from the literature.

Type:: float

x_hi

The upper limit of the distance estimate from the literature.

Type:: float

conf_level

The confidence level associated with the distance estimate, by default 0.68.

Type:: float

x_loerr : float: The lower error, calculated as x_est - x_lo.

x_uperr : float: The upper error, calculated as x_hi - x_est.

sigma_0 : float: The averaged error, calculated as 0.5 * (x_loerr + x_uperr).

fit_gamma : dict: Fit the literature nominal values to a gamma distribution and return the parameters of the fitted distribution.

sample_distance : np.ndarray: Sample distance from the fitted gamma distribution.

Example

>>> from cobtools.astrometry.distance import FromLiterature
>>> model = FromLiterature(x_est=5.0, x_lo=4.0, x_hi=6.0)
>>> print(model.x_loerr)  # Lower error
1.0
>>> print(model.x_uperr)  # Upper error
1.0

fit_gamma() → dict[source]

Fit the asymmetric errors from the literature to a gamma distribution. The fitting is done by minimizing the difference between the confidence intervals of the fitted gamma distribution and the literature values.

Fitting is performed using the scipy.optimize.minimize function, with the initial guess for the scale parameter (sigma_x) set to the averaged error (self.sigma_0).

Returns:: A dictionary containing the parameters of the fitted gamma distribution, including ‘alpha’, ‘theta’, and the ‘distribution’ object itself.
Return type:: dict

sample_distance(n_samples: int = 1000) → ndarray[source]

Sample distance from the fitted gamma distribution.

Parameters:: n_samples (int, optional) – Number of random samples to generate, by default 1000
Returns:: Array of sampled random distances from the fitted gamma distribution.
Return type:: np.ndarray

Example

>>> from cobtools.astrometry.distance import FromLiterature
>>> model = FromLiterature(x_est=5.0, x_lo=4.0, x_hi=6.0)
>>> samples = model.sample_distance(n_samples=100
>>> print(samples)

property sigma_0: float: Calculate the averaged error, which could be used as an initial guess for fitting a skewed distribution.

property x_loerr: float: Calculate the lower error of the distance estimate. :returns: The lower error of the distance estimate. :rtype: float

property x_uperr: float

Calculate the upper error of the distance estimate.

Returns:: The upper error of the distance estimate.
Return type:: float

class cobtools.astrometry.distance.SimpleInversion(parallax: float, parallax_error: float)[source]

Bases: object

A class for estimating distance from parallax using simple inversion. Use this method when the measured parallax is positive and has a high signal-to-noise ratio. Commonly parallax/parallax_error >= 5 is preferred.

parallax

The measured parallax in milliarcseconds (mas). Must be positive for this method

Type:: float

parallax_error

The error in the measured parallax in milliarcseconds (mas). Must be positive.

Type:: float

d_est : float: The estimated distance in kiloparsecs (kpc) obtained by inverting the parallax, i.e., d_est = 1 / parallax.

d_est_error : float: The error in the estimated distance, calculated using error propagation from the parallax error, i.e., d_est_error = parallax_error / parallax^2.

parallax_over_error : float: The signal-to-noise ratio of the parallax measurement, calculated as parallax / parallax_error.

sample_distance : np.ndarray: A method to sample random distances based on a truncated normal distribution centred on the inversion of parallax with a standard deviation equal to parallax_error / parallax^2.

Example

>>> from cobtools.astrometry.distance import SimpleInversion
>>> model = SimpleInversion(parallax=0.5, parallax_error=0.1)
>>> print(model.d_est)  # Estimated distance
2.0
>>> print(model.d_est_error)  # Error in the estimated distance
0.4
>>> print(model.parallax_over_error)  # Signal-to-noise ratio
5.0

property d_est: float

Estimate distance by inverting the parallax.

Returns:: The estimated distance in kiloparsecs (kpc).
Return type:: float

property d_est_error: float

Estimate the error in the distance estimate.

Returns:: The error in the estimated distance.
Return type:: float

property parallax_over_error: float

Calculate the signal-to-noise ratio of the parallax measurement.

Returns:: The signal-to-noise ratio.
Return type:: float

sample_distance(n_samples: int = 1000, dist_lo: float = 0, dist_up: float = inf) → ndarray[source]

A sampler that samples random distance based on a truncated normal distribution centred on the inversion of parallax with a standard deviation equal to parallax_error / parallax^2.

Parameters:

n_samples (int, optional) – Number of random samples to generate, by default 1000
dist_lo (float, optional) – Lower bound for distance sampling, by default 0. Should keep it default.
dist_up (float, optional) – Upper bound for distance sampling, by default np.inf. Should keep it default.

Returns:

Array of sampled random distances.

Return type:

np.ndarray

Raises:

ValueError – If n_samples is not a positive integer.
TypeError – If n_samples is not an integer.

Example

>>> from cobtools.astrometry.distance import SimpleInversion
>>> model = SimpleInversion(parallax=0.5, parallax_error=0.1)
>>> samples = model.sample_distance(n_samples=1000)  # Sample distances
>>> print(samples)  # Array of sampled distances

class cobtools.astrometry.distance.XRBExponentialPriorModel(parallax: float, parallax_error: float, scale_length: float = 1.97)[source]

Bases: object

A class representing distance estimation using a Bayesian method with an exponential prior representative for X-ray binaries. The prior is formulated as:

\[p(d) \propto d^2\exp(-d / L),\]

where \(L\) is the scale_length, a parameter obtained from fitting to known X-ray binaries in the literature.

The likelihood is a Gaussian likelihood centered on 1/d with a standard deviation equal to the parallax error, i.e.,

\[p(\varpi | d) \propto \exp\left( -\frac{1}{2} \frac{(\varpi - 1/d)^2}{\sigma_\varpi^2} \right),\]

This method can be used for negative parallaxes and low signal-to-noise ratio parallaxes.

parallax

The measured parallax in milliarcseconds (mas). Could be negative.

Type:: float

parallax_error

The uncertainty of the parallax measurement in milliarcseconds (mas).

Type:: float

scale_length

The scale length of the exponential prior in kiloparsecs (kpc).

Type:: float

log_likelihood(d) → float[source]

Logarithm of the Gaussian likelihood at a given distance d. For negative distances, the log likelihood is set to negative infinity.

Parameters:: d (distance) – distance in kiloparsecs (kpc) for which to calculate the log likelihood.
Returns:: The logarithm of the Gaussian likelihood at distance d.
Return type:: float

log_posterior(d) → float[source]

Logarithm of the posterior probability at a given distance d. For negative distances, the log posterior is set to negative infinity.

Parameters:: d (float) – distance in kiloparsecs (kpc) for which to calculate the log posterior.
Returns:: Logarithm of the posterior probability at distance d.
Return type:: float

log_prior(d: float) → float[source]

Logarithm of the exponential prior probability density. For negative distances, the log prior is set to negative infinity.

Parameters:: d (float) – distance in kiloparsecs (kpc) for which to calculate the log prior.
Returns:: The logarithm of the exponential prior probability density at distance d.
Return type:: float

property parallax_over_error: float

Calculate the signal-to-noise ratio of the parallax measurement.

Returns:: The signal-to-noise ratio of the parallax measurement.
Return type:: float

sample_distance(nwalkers: int = 4, nsteps: int = 2000, burn_in: int = 500, **kwargs: Any) → ndarray[Any, dtype[float64]][source]

Sample distances from the posterior distribution using the EnsembleSampler from the emcee package.

Parameters:

nwalkers (int, optional) – Number of walkers for the MCMC sampler, by default 4
nsteps (int, optional) – Number of steps for each walker, by default 2000
burn_in (int, optional) – Number of steps to discard as burn-in, by default 500
kwargs (Any) – Additional keyword arguments to pass to the emcee.EnsembleSampler

Returns:

Array of sampled distances from the posterior distribution.

Return type:

np.ndarray[Any, np.dtype[np.float64]]

Raises:

ValueError – If nwalkers is not a positive integer.
ValueError – If nsteps is not a positive integer.
ValueError – If burn_in is not a non-negative integer.
ValueError – If burn_in is greater than or equal to nsteps.

Example

>>> from cobtools.astrometry.distance import XRBExponentialPriorModel
>>> m = XRBExponentialPriorModel(parallax=0.5, parallax_error=0.1)
>>> samples = m.sample_distance(nwalkers=4, nsteps=2000, burn_in=500)
>>> print(samples)  # Array of sampled distances

cobtools.astrometry.kinematics module

kinematics.py

This module provides a suite of functions for computing Galactic space velocities, including galactocentric 3D Cartesian velocities and peculiar velocities. The formulation is based on the work of Reid et al. 2009.

Functions

equatorial_to_galactic: Convert equatorial coordinates (RA, Dec) to Galactic coordinates (l, b).
galactic_proper_motion: Convert equatorial proper motions to Galactic proper motions.
galactocentric_cartesian_velocity: Calculate galactocentric Cartesian velocities (u, v, w) and total space velocities (square root of the quadrature sum of u, v, w) from equatorial coordinates, proper motions, distance, and radial velocity.
peculiar_velocity: Calculate the peculiar velocity and its Cartesian components by subtracting the local Galactic rotation.

Convert equatorial coordinates (ra, dec) to galactic coordinates (gal_l, gal_b).

Parameters:

ra (float or ArrayLike) – Right ascension in decimal degrees.
dec (float or ArrayLike) – Declination in decimal degrees.

Returns:

Galactic longitude (l) and latitude (b) in decimal degrees. If inputs are scalars, outputs will be numpy scalars.

Return type:

Tuple[np.ndarray, np.ndarray]

Raises:

ValueError – If ra and dec have different shapes, or if ra is not in [0, 360) degrees, or if dec is not in [-90, 90] degrees.

Example

>>> from cobtools.astrometry.kinematics import equatorial_to_galactic
>>> ra = 10.684  # degrees
>>> dec = 41.269  # degrees
>>> l, b = equatorial_to_galactic(ra, dec)
>>> print(f"Galactic coordinates: l={l:.2f}, b={b:.2f} degrees")

Calculate galactic proper motion (mu_l, mu_b) from equatorial proper motion (pmra_cosdec, pmdec) and equatorial coordinates (ra, dec).

The basic idea is to calculate differences in ra and dec due to proper motion and convert it changes in Galactic coordinate.

ra and dec should be in degrees, and their proper motion components should be in units of mas/yr; note that mu_ra_cosdec contain the cos(dec) factor. Input timestep should be in unit of years (default is 1).

Parameters:

ra (float or ArrayLike) – Right ascension in decimal degrees.
dec (float or ArrayLike) – Declination in decimal degrees.
pmra_cosdec (float or ArrayLike) – Proper motion in right ascension multiplied by cos(dec), in mas/yr.
pmdec (float or ArrayLike) – Proper motion in declination, in mas/yr.
dt (float, optional) – Time step in years for which to calculate the proper motion, by default 1.0.

Returns:

Proper motion in Galactic longitude (mu_l) and latitude (mu_b) in mas/yr. Note that mu_l does not contain the cos(b) factor, i.e., it is the proper motion in the direction of increasing Galactic longitude.

Return type:

Tuple[Union[float, np.ndarray], Union[float, np.ndarray]]

Raises:

ValueError – If the time step dt is not positive.

Example

>>> from cobtools.astrometry.kinematics import galactic_proper_motion
>>> ra = 10.684  # degrees
>>> dec = 41.269  # degrees
>>> pmra_cosdec = 0.1  # mas/yr
>>> pmdec = 0.2  # mas/yr
>>> mu_l, mu_b = galactic_proper_motion(ra, dec, pmra_cosdec, pmdec)
>>> print(f"Galactic proper motion: {mu_l:.2f}, {mu_b:.2f} mas/yr")

Calculate Galactocentric Cartesian velocity components (u1, v1, w1) and the total space velocity (vspace) from equatorial coordinates, proper motions, distance, and radial velocity.

Parameters:

ra (float or ArrayLike) – Right ascension in decimal degrees.
dec (float or ArrayLike) – Declination in decimal degrees.
pmra_cosdec (float or ArrayLike) – Proper motion in ra*cos(dec), in mas/yr.
pmdec (float or ArrayLike) – Proper motion in declination, in mas/yr.
dist (float or ArrayLike) – Distance in kpc.
rv (float or ArrayLike) – Radial velocity in km/s.
u_sun (float or ArrayLike, optional) – Solar motion u component in km/s.
v_sun (float or ArrayLike, optional) – Solar motion v component in km/s.
w_sun (float or ArrayLike, optional) – Solar motion w component in km/s.
theta_sun (float or ArrayLike, optional) – Solar rotation velocity in the Galactic plane in km/s.
dt (float, optional) – Time step in years for which to calculate the proper motion, see galactic_proper_motion, by default 1.0 (year).

Returns:

Tuple[Union[float, np.ndarray], Union[float, np.ndarray],
Union[float, np.ndarray], Union[float, np.ndarray]] – Galactocentric Cartesian velocities (u1, v1, w1, vspace) in km/s.

Raises:

ValueError – If ra is not in [0, 360) degrees, or if dec is not in [-90, 90] degrees.
ValueError – If dt is not a positive number.

Example

>>> from cobtools.astrometry.kinematics import     ... galactocentric_cartesian_velocity

>>> ra = 10.684  # degrees
>>> dec = 41.269  # degrees
>>> pmra_cosdec = 0.1  # mas/yr
>>> pmdec = 0.2  # mas/yr
>>> dist = 0.77  # kpc
>>> rv = -300  # km/s
>>> u1, v1, w1, vspace = galactocentric_cartesian_velocity(
...     ra, dec, pmra_cosdec, pmdec, dist, rv
... )

Calculate the peculiar velocity and the cartesian components given its equatorial coordinates, proper motions, distance, and radial velocity.

Parameters:

ra (float or ArrayLike) – Right ascension in decimal degrees.
dec (float or ArrayLike) – Declination in decimal degrees.
pmra_cosdec (float or ArrayLike) – Proper motion in ra*cos(dec), in mas/yr.
pmdec (float or ArrayLike) – Proper motion in declination, in mas/yr.
dist (float or ArrayLike) – Distance in kpc.
rv (float or ArrayLike) – Radial velocity in km/s.
u_sun (float or ArrayLike, optional) – Solar motion u component in km/s.
v_sun (float or ArrayLike, optional) – Solar motion v component in km/s.
w_sun (float or ArrayLike, optional) – Solar motion w component in km/s.
theta_sun (float or ArrayLike, optional) – Solar rotation velocity in the Galactic plane in km/s.
r_sun (float or ArrayLike, optional) – Distance from the Sun to the Galactic center in kpc.
dt (float, optional) – Time step in years for which to calculate the proper motion, see galactic_proper_motion, by default 1.0 (year).

Returns:

Tuple[Union[float, np.ndarray], Union[float, np.ndarray],
Union[float, np.ndarray], Union[float, np.ndarray]] – Peculiar velocity components (u_s, v_s, w_s) and the total peculiar velocity (vpec) in km/s. Note that u_s is the component toward the Galactic center, v_s is the component in the direction of Galactic rotation, and w_s is the component toward the North Galactic Pole.

Example

>>> from cobtools.astrometry.kinematics import peculiar_velocity
>>> import numpy as np
>>> # Single value example
>>> ra = 10.684  # degrees
>>> dec = 41.269  # degrees
>>> pmra_cosdec = 0.1  # mas/yr
>>> pmdec = 0.2  # mas/yr
>>> dist = 0.77  # kpc
>>> rv = -300  # km/s
>>> u_s, v_s, w_s, vpec = peculiar_velocity(
...     ra, dec, pmra_cosdec, pmdec, dist, rv
... )
>>> # Array broadcasting example
>>> ra_array = 10.684  # degrees
>>> dec_array = 41.269  # degrees
>>> pmra_cosdec_array = np.random.normal(-0.35, 0.08, 100)  # mas/yr
>>> pmdec_array = np.random.normal(0.1, 0.05, 100)  # mas/yr
>>> dist_array = np.random.normal(1.2, 0.3, 100) # kpc
>>> rv_array = 15.0  # km/s
>>> u_s_arr, v_s_arr, w_s_arr, vpec_arr = peculiar_velocity(
...     ra_array, dec_array, pmra_cosdec_array, pmdec_array,
...     dist_array, rv_array
... )

cobtools.astrometry.rotation_curve module

Load a pre-computed rotation curve, and return an interpolated function for the rotation velocity as a function of radius.

Parameters:

r_sun (float, optional) – Distance from the Sun to the Galactic center in kpc, by default con.r_sun
theta_sun (float, optional) – Circular velocity at the Sun’s position in km/s, by default con.theta_sun
data_path (Path, optional) – Path to the rotation curve data file. If not provided, the default path is resolved as resources.files(“cobtools”) / “data” / “rotcurve_mw2014.npy”.

Returns:

Interpolated rotation velocity as a function of radius. The function extrapolates beyond the provided data range, but clamps negative radii to zero, i.e., v_rot(r < 0) = v_rot(0).

Return type:

Callable[[ArrayLike], np.ndarray]

Raises:

ValueError – If r_sun or theta_sun are not positive values.
FileNotFoundError – If the rotation curve data file is not found.
ValueError – If there is an error loading the rotation curve data.

cobtools.astrometry package

Submodules

cobtools.astrometry.distance module

Classes

cobtools.astrometry.kinematics module

kinematics.py

Functions

cobtools.astrometry.rotation_curve module

Module contents