Source code for orbitize.gaia

import numpy as np
import contextlib

with contextlib.redirect_stdout(None):
    from astroquery.gaia import Gaia
from astropy import units as u

[docs]class GaiaLogProb(object): """ Class to compute the log probability of an orbit with respect to a single astrometric position point from Gaia. Uses astroquery to look up Gaia astrometric data, and computes log-likelihood. To be used in conjunction with orbitize.hipparcos.HipLogProb; see documentation for that object for more detail. Follows Nielsen+ 2020 (studying the orbit of beta Pic b). Note that this class currently only fits for the position of the star in the Gaia epoch, not the star's proper motion. .. Note:: In orbitize, it is possible to perform a fit to just the Hipparcos IAD, but not to just the Gaia astrometric data. Args: gaia_num (int): the Gaia source ID of the object you're fitting. Note that the dr2 and edr3 source IDs are not necessarily the same. hiplogprob (orbitize.hipparcos.HipLogProb): object containing all info relevant to Hipparcos IAD fitting dr (str): either 'dr2' or 'edr3' query (bool): if True, queries the Gaia database for astrometry of the target (requires an internet connection). If False, uses user-input astrometric values (runs without internet). gaia_data (dict): see `query` keyword above. If `query` set to False, then user must supply a dictionary of Gaia astometry in the following form: gaia_data = { 'ra': 139.4 # RA in degrees 'dec': 139.4 # Dec in degrees 'ra_error': 0.004 # RA error in mas 'dec_error': 0.004 # Dec error in mas } Written: Sarah Blunt, 2021 """ def __init__(self, gaia_num, hiplogprob, dr='dr2', query=True, gaia_data=None): self.gaia_num = gaia_num self.hiplogprob = hiplogprob self.dr = dr if self.dr == 'edr3': self.gaia_epoch = 2016.0 elif self.dr == 'dr2': self.gaia_epoch = 2015.5 else: raise ValueError("`dr` must be either `dr2` or `edr3`") self.hipparcos_epoch = 1991.25 if query: query = """SELECT TOP 1 ra, dec, ra_error, dec_error FROM gaia{}.gaia_source WHERE source_id = {} """.format(self.dr, self.gaia_num) job = Gaia.launch_job_async(query) gaia_data = job.get_results() self.ra = gaia_data['ra'] self.ra_err = gaia_data['ra_error'] self.dec = gaia_data['dec'] self.dec_err = gaia_data['dec_error'] # keep this number on hand for use in lnlike computation self.mas2deg = (u.mas).to(u.degree) def _save(self, hf): """ Saves the current object to an hdf5 file Args: hf (h5py._hl.files.File): a currently open hdf5 file in which to save the object. """ hf.attrs['gaia_num'] = self.gaia_num hf.attrs['dr'] = self.dr self.hiplogprob._save(hf)
[docs] def compute_lnlike( self, raoff_model, deoff_model, samples, param_idx ): """ Computes the log likelihood of an orbit model with respect to a single Gaia astrometric point. This is added to the likelihoods calculated with respect to other data types in ``sampler._logl()``. Args: raoff_model (np.array of float): 2xM primary RA offsets from the barycenter incurred from orbital motion of companions (i.e. not from parallactic motion), where M is the number of orbits being tested, and raoff_model[0,:] are position predictions at the Hipparcos epoch, and raoff_model[1,:] are position predictions at the Gaia epoch deoff_model (np.array of float): 2xM primary decl offsets from the barycenter incurred from orbital motion of companions (i.e. not from parallactic motion), where M is the number of orbits being tested, and deoff_model[0,:] are position predictions at the Hipparcos epoch, and deoff_model[1,:] are position predictions at the Gaia epoch samples (np.array of float): R-dimensional array of fitting parameters, where R is the number of parameters being fit. Must be in the same order documented in ``System``. param_idx: a dictionary matching fitting parameter labels to their indices in an array of fitting parameters (generally set to System.basis.param_idx). Returns: np.array of float: array of length M, where M is the number of input orbits, representing the log likelihood of each orbit with respect to the Gaia position measurement. """ alpha_H0 = samples[param_idx['alpha0']] # [deg] pm_ra = samples[param_idx['pm_ra']] # [mas/yr] delta_alpha_from_pm = pm_ra * (self.gaia_epoch - self.hipparcos_epoch) # [mas] delta_H0 = samples[param_idx['delta0']] # [deg] pm_dec = samples[param_idx['pm_dec']] # [mas/yr] delta_delta_from_pm = pm_dec * (self.gaia_epoch - self.hipparcos_epoch) # [mas] # difference in position due to orbital motion between Hipparcos & Gaia epochs alpha_diff_orbit = (raoff_model[1,:] - raoff_model[0,:]) # [mas] dec_diff_orbit = (deoff_model[1,:] - deoff_model[0,:]) # RA model (not in tangent plane) alpha_model = ( # [deg] self.hiplogprob.alpha0 + self.mas2deg * ( alpha_H0 + delta_alpha_from_pm + alpha_diff_orbit # divide by cos(dec) to undo projection onto tangent plane ) / np.cos(np.radians(self.hiplogprob.delta0)) ) alpha_data = self.ra # again divide by cos(dec) to undo projection onto tangent plane alpha_unc = self.mas2deg * self.ra_err / np.cos(np.radians(self.hiplogprob.delta0)) # technically this is an angle so we should wrap it, but the precision # of Hipparcos and Gaia is so good that we'll never have to. alpha_resid = (alpha_model - alpha_data) alpha_chi2 = (alpha_resid / alpha_unc)**2 delta_model = ( # [deg] self.hiplogprob.delta0 + self.mas2deg * ( delta_H0 + delta_delta_from_pm + dec_diff_orbit ) ) dec_data = self.dec delta_unc = self.mas2deg * self.dec_err delta_chi2 = ((delta_model - dec_data) / delta_unc)**2 chi2 = alpha_chi2 + delta_chi2 lnlike = -0.5 * chi2 return lnlike