import numpy as np
import astropy.units as u, astropy.constants as consts
import warnings
import abc
from orbitize import priors, kepler
from scipy.optimize import fsolve
[docs]class Basis(abc.ABC):
"""
Abstract base class for different basis sets. All new basis objects should inherit from
this class. This class is meant to control how priors are assigned to various basis sets
and how conversions are made from the basis sets to the standard keplarian set.
Author: Tirth, 2021
"""
def __init__(
self, stellar_or_system_mass, mass_err, plx, plx_err, num_secondary_bodies,
fit_secondary_mass, angle_upperlim=2*np.pi, hipparcos_IAD=None,
rv=False, rv_instruments=None
):
self.stellar_or_system_mass = stellar_or_system_mass
self.mass_err=mass_err
self.plx=plx
self.plx_err=plx_err
self.num_secondary_bodies=num_secondary_bodies
self.angle_upperlim=angle_upperlim
self.fit_secondary_mass=fit_secondary_mass
self.hipparcos_IAD = hipparcos_IAD
self.rv = rv
self.rv_instruments = rv_instruments
# Define dictionary of default priors to be updated as new basis sets are added
self.default_priors = {
'sma' : priors.LogUniformPrior(0.001, 1e4),
'per' : priors.LogUniformPrior(1e-5, 1e6),
'ecc' : priors.UniformPrior(0., 1.),
'inc' : priors.SinPrior(),
'aop' : priors.UniformPrior(0., 2.*np.pi),
'pan' : priors.UniformPrior(0., angle_upperlim),
'tau' : priors.UniformPrior(0., 1.),
'K' : priors.LogUniformPrior(1e-4, 10)
}
@abc.abstractmethod
def construct_priors(self):
pass
@abc.abstractmethod
def to_standard_basis(self, param_arr, param_idx):
pass
[docs] def verify_params(self):
'''
Displays warnings about the 'fit_secondary_mass' and 'rv' parameters for
all basis sets. Can be overriden by any basis class depending on the
necessary parameters that need to be checked.
'''
if not self.fit_secondary_mass and self.rv:
warnings.warn(
""""
Radial velocity data found in input data, but rv parameters will
not be sampled. To sample rv parameters, set 'fit_secondary_mass'
to True.
"""
)
[docs] def set_hip_iad_priors(self, priors_arr, labels_arr):
'''
Adds the necessary priors relevant to the hipparcos data to 'priors_arr'
and updates 'labels_arr' with the priors' corresponding labels.
Args:
priors_arr (list of orbitize.priors.Prior objects): holds the prior
objects for each parameter to be fitted (updated here)
labels_arr (list of strings): holds the name of all the parameters
to be fitted (updated here)
'''
priors_arr.append(priors.UniformPrior(
self.hipparcos_IAD.pm_ra0 - 10 * self.hipparcos_IAD.pm_ra0_err,
self.hipparcos_IAD.pm_ra0 + 10 * self.hipparcos_IAD.pm_ra0_err)
)
labels_arr.append('pm_ra')
priors_arr.append(priors.UniformPrior(
self.hipparcos_IAD.pm_dec0 - 10 * self.hipparcos_IAD.pm_dec0_err,
self.hipparcos_IAD.pm_dec0 + 10 * self.hipparcos_IAD.pm_dec0_err)
)
labels_arr.append('pm_dec')
priors_arr.append(priors.UniformPrior(
- 10 * self.hipparcos_IAD.alpha0_err,
10 * self.hipparcos_IAD.alpha0_err)
)
labels_arr.append('alpha0')
priors_arr.append(priors.UniformPrior(
- 10 * self.hipparcos_IAD.delta0_err,
10 * self.hipparcos_IAD.delta0_err)
)
labels_arr.append('delta0')
[docs] def set_rv_priors(self, priors_arr, labels_arr):
'''
Adds the necessary priors if radial velocity data is supplied to
'priors_arr' and updates 'labels_arr' with the priors' corresponding
labels. This function assumes that 'rv' data has been supplied and
a secondary mass is being fitted for.
Args:
priors_arr (list of orbitize.priors.Prior objects): holds the prior
objects for each parameter to be fitted (updated here)
labels_arr (list of strings): holds the name of all the parameters
to be fitted (updated here)
'''
for instrument in self.rv_instruments:
priors_arr.append(priors.UniformPrior(-5, 5)) # gamma prior in km/s
labels_arr.append('gamma_{}'.format(instrument))
priors_arr.append(priors.LogUniformPrior(1e-4, 0.05)) # jitter prior in km/s
labels_arr.append('sigma_{}'.format(instrument))
[docs] def set_default_mass_priors(self, priors_arr, labels_arr):
'''
Adds the necessary priors for the stellar and/or companion masses.
Args:
priors_arr (list of orbitize.priors.Prior objects): holds the prior
objects for each parameter to be fitted (updated here)
labels_arr (list of strings): holds the name of all the parameters
to be fitted (updated here)
'''
if self.fit_secondary_mass:
for body in np.arange(self.num_secondary_bodies)+1:
priors_arr.append(priors.LogUniformPrior(1e-6, 2)) # in Solar masses
labels_arr.append('m{}'.format(body))
labels_arr.append('m0')
else:
labels_arr.append('mtot')
if self.mass_err > 0:
priors_arr.append(priors.GaussianPrior(self.stellar_or_system_mass, self.mass_err))
else:
priors_arr.append(self.stellar_or_system_mass)
[docs]class Standard(Basis):
'''
Standard basis set based upon the 6 standard Keplarian elements: (sma, ecc, inc, aop, pan, tau).
Args:
stellar_or_system_mass (float): mass of the primary star (if fitting for
dynamical masses of both components) or total system mass (if
fitting using relative astrometry only) [M_sol]
mass_err (float): uncertainty on 'stellar_or_system_mass', in M_sol
plx (float): mean parallax of the system, in mas
plx_err (float): uncertainty on 'plx', in mas
num_secondary_bodies (int): number of secondary bodies in the system, should be at least 1
fit_secondary_mass (bool): if True, include the dynamical mass of orbitting body as fitted parameter, if False,
'stellar_or_system_mass' is taken to be total mass
angle_upperlim (float): either pi or 2pi, to restrict the prior range for 'pan' parameter (default: 2pi)
hipparcos_IAD (orbitize.HipparcosLogProb object): if not 'None', then add relevant priors to this data (default: None)
rv (bool): if True, then there is radial velocity data and assign radial velocity priors, if False, then there
is no radial velocity data and radial velocity priors are not assigned (default: False)
rv_instruments (np.array): array of unique rv instruments from the originally supplied data (default: None)
'''
def __init__(self, stellar_or_system_mass, mass_err, plx, plx_err, num_secondary_bodies, fit_secondary_mass,
angle_upperlim=2*np.pi, hipparcos_IAD=None, rv=False, rv_instruments=None):
super(Standard, self).__init__(stellar_or_system_mass, mass_err, plx, plx_err, num_secondary_bodies,
fit_secondary_mass, angle_upperlim, hipparcos_IAD, rv, rv_instruments)
[docs] def construct_priors(self):
'''
Generates the parameter label array and initializes the corresponding priors for each
parameter that's to be sampled. For the standard basis, the parameters common to each
companion are: sma, ecc, inc, aop, pan, tau. Parallax, hipparcos (optional), rv (optional),
and mass priors are added at the end.
Returns:
tuple:
list: list of strings (labels) that indicate the names of each parameter to sample
list: list of orbitize.priors.Prior objects that indicate the prior distribution of each label
'''
base_labels = ['sma', 'ecc', 'inc', 'aop', 'pan', 'tau']
basis_priors = []
basis_labels = []
# Add the priors common to each companion
for body in np.arange(self.num_secondary_bodies):
for elem in base_labels:
basis_priors.append(self.default_priors[elem])
basis_labels.append(elem + str(body+1))
# Add parallax prior
basis_labels.append('plx')
if self.plx_err > 0:
basis_priors.append(priors.GaussianPrior(self.plx, self.plx_err))
else:
basis_priors.append(self.plx)
# Add hippparcos priors if necessary
if self.hipparcos_IAD is not None:
self.set_hip_iad_priors(basis_priors, basis_labels)
# Add rv priors
if self.rv and self.fit_secondary_mass:
self.set_rv_priors(basis_priors, basis_labels)
# Add mass priors
self.set_default_mass_priors(basis_priors, basis_labels)
# Define param label dictionary in current basis & standard basis
self.param_idx = dict(zip(basis_labels, np.arange(len(basis_labels))))
self.standard_basis_idx = dict(zip(basis_labels, np.arange(len(basis_labels))))
return basis_priors, basis_labels
[docs] def to_standard_basis(self, param_arr):
'''
For standard basis, no conversion needs to be made.
Args:
param_arr (np.array of float): RxM array of fitting parameters in the standard basis,
where R is the number of parameters being fit, and M is the number of orbits. If
M=1 (for MCMC), this can be a 1d array.
Returns:
np.array of float: ``param_arr`` without any modification
'''
return param_arr
[docs]class Period(Basis):
'''
Modification of the standard basis, swapping our sma for period: (per, ecc, inc, aop, pan, tau).
Args:
stellar_or_system_mass (float): mass of the primary star (if fitting for
dynamical masses of both components) or total system mass (if
fitting using relative astrometry only) [M_sol]
mass_err (float): uncertainty on 'stellar_or_system_mass', in M_sol
plx (float): mean parallax of the system, in mas
plx_err (float): uncertainty on 'plx', in mas
num_secondary_bodies (int): number of secondary bodies in the system, should be at least 1
fit_secondary_mass (bool): if True, include the dynamical mass of orbitting body as fitted parameter, if False,
'stellar_or_system_mass' is taken to be total mass
angle_upperlim (float): either pi or 2pi, to restrict the prior range for 'pan' parameter (default: 2pi)
hipparcos_IAD (orbitize.HipparcosLogProb object): if not 'None', then add relevant priors to this data (default: None)
rv (bool): if True, then there is radial velocity data and assign radial velocity priors, if False, then there
is no radial velocity data and radial velocity priors are not assigned (default: False)
rv_instruments (np.array): array of unique rv instruments from the originally supplied data (default: None)
'''
def __init__(self, stellar_or_system_mass, mass_err, plx, plx_err, num_secondary_bodies, fit_secondary_mass,
angle_upperlim=2*np.pi, hipparcos_IAD=None, rv=False, rv_instruments=None):
super(Period, self).__init__(stellar_or_system_mass, mass_err, plx, plx_err, num_secondary_bodies,
fit_secondary_mass, angle_upperlim, hipparcos_IAD, rv, rv_instruments)
[docs] def construct_priors(self):
'''
Generates the parameter label array and initializes the corresponding priors for each
parameter that's to be sampled. For the standard basis, the parameters common to each
companion are: per, ecc, inc, aop, pan, tau. Parallax, hipparcos (optional), rv (optional),
and mass priors are added at the end.
Returns:
tuple:
list: list of strings (labels) that indicate the names of each parameter to sample
list: list of orbitize.priors.Prior objects that indicate the prior distribution of each label
'''
base_labels = ['per', 'ecc', 'inc', 'aop', 'pan', 'tau']
basis_priors = []
basis_labels = []
# Add the priors common to each companion
for body in np.arange(self.num_secondary_bodies):
for elem in base_labels:
basis_priors.append(self.default_priors[elem])
basis_labels.append(elem + str(body+1))
# Add parallax prior
basis_labels.append('plx')
if self.plx_err > 0:
basis_priors.append(priors.GaussianPrior(self.plx, self.plx_err))
else:
basis_priors.append(self.plx)
# Add hipparcos priors if necessary
if self.hipparcos_IAD is not None:
self.set_hip_iad_priors(basis_priors, basis_labels)
# Add rv priors
if self.rv and self.fit_secondary_mass:
self.set_rv_priors(basis_priors, basis_labels)
# Add mass priors
self.set_default_mass_priors(basis_priors, basis_labels)
# Define param label dictionary in current basis & standard basis
self.param_idx = dict(zip(basis_labels, np.arange(len(basis_labels))))
self.standard_basis_idx = dict(zip(basis_labels, np.arange(len(basis_labels))))
for body_num in np.arange(self.num_secondary_bodies) + 1:
self.standard_basis_idx['sma{}'.format(body_num)] = self.param_idx['per{}'.format(body_num)]
return basis_priors, basis_labels
[docs] def to_standard_basis(self, param_arr):
'''
Convert parameter array from period basis to standard basis by swapping out the period
parameter to semi-major axis for each companion.
Args:
param_arr (np.array of float): RxM array of fitting parameters in the period basis,
where R is the number of parameters being fit, and M is the number of orbits. If
M=1 (for MCMC), this can be a 1D array.
Returns:
np.array of float: modifies 'param_arr' to contain the semi-major axis for each companion
in each orbit rather than period. Shape of 'param_arr' remains the same.
'''
for body_num in np.arange(self.num_secondary_bodies)+1:
per = param_arr[self.param_idx['per{}'.format(body_num)]]
if self.fit_secondary_mass:
# Assume two-body system
m_secondary = param_arr[self.param_idx['m{}'.format(body_num)]]
m0 = param_arr[self.param_idx['m0']]
mtot = m_secondary + m0
else:
mtot = param_arr[self.param_idx['mtot']]
# Compute semi-major axis using Kepler's Third Law and replace period
sma = np.cbrt((consts.G*(mtot*u.Msun)*((per*u.year)**2))/(4*np.pi**2))
sma = sma.to(u.AU).value
param_arr[self.standard_basis_idx['per{}'.format(body_num)]] = sma
return param_arr
[docs] def to_period_basis(self, param_arr):
'''
Convert parameter array from standard basis to period basis by swapping out the semi-major
axis parameter to period for each companion. This function is used primarily for testing
purposes.
Args:
param_arr (np.array of float): RxM array of fitting parameters in the standard basis,
where R is the number of parameters being fit, and M is the number of orbits. If
M=1 (for MCMC), this can be a 1D array.
Returns:
np.array of float: modifies 'param_arr' to contain the period for each companion
in each orbit rather than semi-major axis. Shape of 'param_arr' remains the same.
'''
for body_num in np.arange(self.num_secondary_bodies)+1:
sma = param_arr[self.standard_basis_idx['sma{}'.format(body_num)]]
if self.fit_secondary_mass:
# Assume two-body system
m_secondary = param_arr[self.standard_basis_idx['m{}'.format(body_num)]]
m0 = param_arr[self.standard_basis_idx['m0']]
mtot = m_secondary + m0
else:
mtot = param_arr[self.standard_basis_idx['mtot']]
per = np.sqrt((4*(np.pi**2)*(sma*u.AU)**3) / (consts.G*(mtot*u.Msun)))
per = per.to(u.year).value
param_arr[self.param_idx['per{}'.format(body_num)]] = per
return param_arr
[docs]class SemiAmp(Basis):
'''
Modification of the standard basis, swapping our sma for period and additionally sampling in
the stellar radial velocity semi-amplitude: (per, ecc, inc, aop, pan, tau, K).
.. Note:: Ideally, 'fit_secondary_mass' is true and rv data is supplied.
Args:
stellar_or_system_mass (float): mass of the primary star (if fitting for
dynamical masses of both components) or total system mass (if
fitting using relative astrometry only) [M_sol]
mass_err (float): uncertainty on 'stellar_or_system_mass', in M_sol
plx (float): mean parallax of the system, in mas
plx_err (float): uncertainty on 'plx', in mas
num_secondary_bodies (int): number of secondary bodies in the system, should be at least 1
fit_secondary_mass (bool): if True, include the dynamical mass of orbitting body as fitted parameter, if False,
'stellar_or_system_mass' is taken to be total mass
angle_upperlim (float): either pi or 2pi, to restrict the prior range for 'pan' parameter (default: 2*pi)
hipparcos_IAD (orbitize.HipparcosLogProb object): if not 'None', then add relevant priors to this data (default: None)
rv (bool): if True, then there is radial velocity data and assign radial velocity priors, if False, then there
is no radial velocity data and radial velocity priors are not assigned (default: False)
rv_instruments (np.array): array of unique rv instruments from the originally supplied data (default: None)
'''
def __init__(self, stellar_or_system_mass, mass_err, plx, plx_err, num_secondary_bodies, fit_secondary_mass,
angle_upperlim=2*np.pi, hipparcos_IAD=None, rv=False, rv_instruments=None):
super(SemiAmp, self).__init__(stellar_or_system_mass, mass_err, plx, plx_err, num_secondary_bodies,
fit_secondary_mass, angle_upperlim, hipparcos_IAD, rv, rv_instruments)
[docs] def construct_priors(self):
'''
Generates the parameter label array and initializes the corresponding priors for each
parameter that's to be sampled. For the semi-amp basis, the parameters common to each
companion are: per, ecc, inc, aop, pan, tau, K (stellar rv semi-amplitude). Parallax,
hipparcos (optional), rv (optional), and mass priors are added at the end.
The mass parameter will always be m0.
Returns:
tuple:
list: list of strings (labels) that indicate the names of each parameter to sample
list: list of orbitize.priors.Prior objects that indicate the prior distribution of each label
'''
base_labels = ['per', 'ecc', 'inc', 'aop', 'pan', 'tau', 'K']
basis_priors = []
basis_labels = []
# Add the priors common to each companion
for body in np.arange(self.num_secondary_bodies):
for elem in base_labels:
basis_priors.append(self.default_priors[elem])
basis_labels.append(elem + str(body+1))
# Add parallax prior
basis_labels.append('plx')
if self.plx_err > 0:
basis_priors.append(priors.GaussianPrior(self.plx, self.plx_err))
else:
basis_priors.append(self.plx)
# Add hip_iad priors if necessary
if self.hipparcos_IAD is not None:
self.set_hip_iad_priors(basis_priors, basis_labels)
# Add rv priors
if self.rv and self.fit_secondary_mass:
self.set_rv_priors(basis_priors, basis_labels)
# Add star mass prior (for now, regardless of whether 'fit_secondary_mass' is true)
if self.mass_err > 0:
basis_priors.append(priors.GaussianPrior(self.stellar_or_system_mass, self.mass_err))
else:
basis_priors.append(self.stellar_or_system_mass)
basis_labels.append('m0')
# Define param label dictionary in current basis & standard basis
self.param_idx = dict(zip(basis_labels, np.arange(len(basis_labels))))
self.standard_basis_idx = dict(zip(basis_labels, np.arange(len(basis_labels))))
for body_num in np.arange(self.num_secondary_bodies) + 1:
self.standard_basis_idx['sma{}'.format(body_num)] = self.param_idx['per{}'.format(body_num)]
self.standard_basis_idx['m{}'.format(body_num)] = self.param_idx['K{}'.format(body_num)]
return basis_priors, basis_labels
[docs] def to_standard_basis(self, param_arr):
'''
Convert parameter array from semi-amp basis to standard basis by swapping out the period
parameter to semi-major axis for each companion and computing the masses of each
companion.
Args:
param_arr (np.array of float): RxM array of fitting parameters in the period basis,
where R is the number of parameters being fit, and M is the number of orbits. If
M=1 (for MCMC), this can be a 1D array.
Returns:
np.array of float: modifies 'param_arr' to contain the semi-major axis for each companion
in each orbit rather than period, removes stellar rv semi-amplitude parameters for
each companion, and appends the companion masses to 'param_arr'
'''
m0 = param_arr[self.param_idx['m0']]
# Compute each companion's mass and sma
for body_num in np.arange(self.num_secondary_bodies) + 1:
period = param_arr[self.param_idx['per{}'.format(body_num)]]
ecc = param_arr[self.param_idx['ecc{}'.format(body_num)]]
inc = param_arr[self.param_idx['inc{}'.format(body_num)]]
semi_amp = param_arr[self.param_idx['K{}'.format(body_num)]]
# Replace semi-amp with companion mass and period with sma
companion_m = self.compute_companion_mass(period, ecc, inc, semi_amp, m0)
param_arr[self.standard_basis_idx['m{}'.format(body_num)]] = companion_m
companion_sma = self.compute_companion_sma(period, m0, companion_m)
param_arr[self.standard_basis_idx['sma{}'.format(body_num)]] = companion_sma
return param_arr
[docs] def func(self, x, lhs, m0):
'''
Define function for scipy.fsolve to use when computing companion mass.
Args:
x (float): the companion mass to be calculated (Msol)
lhs (float): the left hand side of the rv semi-amplitude equation (Msol^(1/3))
m0 (float): the stellar mass (Msol)
Returns:
float: the difference between the rhs and lhs of the rv semi-amplitude equation, 'x' is a
good companion mass when this difference is very close to zero
'''
return ((x / ((x + m0)**(2/3))) - lhs)
[docs] def compute_companion_mass(self, period, ecc, inc, semi_amp, m0):
'''
Computes a single companion's mass given period, eccentricity, inclination, stellar rv semi-amplitude,
and stellar mass. Uses scipy.fsolve to compute the masses numerically.
Args:
period (np.array of float): the period values for each orbit for a single companion (can be float)
ecc (np.array of float): the eccentricity values for each orbit for a single companion (can be float)
inc (np.array of float): the inclination values for each orbit for a single companion (can be float)
semi_amp (np.array of float): the stellar rv-semi amp values for each orbit (can be float)
m0 (np.array of float): the stellar mass for each orbit (can be float)
Returns:
np.array of float: the companion mass values for each orbit (can also just be a single float)
'''
# Define LHS of equation
kms = u.km / u.s
lhs = ((semi_amp*kms)*((1-ecc**2)**(1/2))*((period*u.yr)**(1/3))*(consts.G**(-1/3))*((4*np.pi**2)**(-1/6))) / (np.sin(inc))
lhs = (lhs.to((u.solMass)**(1/3))).value
m_n = []
# Solve for companion mass numerically, making initial guess at center of uniform prior distribution (Msol)
if (not hasattr(m0, '__len__')):
comp_mass = fsolve(self.func, x0=1e-3, args=(lhs, m0))
m_n.append(comp_mass[0])
else:
for orbit in range(len(m0)):
comp_mass = fsolve(self.func, x0=1e-3, args=(lhs[orbit], m0[orbit]))
m_n.append(comp_mass[0])
# squash dimensions
if len(m_n) == 1:
m_n = m_n[0]
return m_n
[docs] def compute_companion_sma(self, period, m0, m_n):
'''
Computes a single companion's semi-major axis using Kepler's Third Law for each orbit.
Args:
period (np.array of float): the period values for each orbit for a single companion (can be float)
m0 (np.array of float): the stellar mass for each orbit (can be float)
m_n (np.array of float): the companion mass for each orbit (can be float)
Returns:
np.array of float: the semi-major axis values for each orbit
'''
sma = np.cbrt((consts.G*((m0+m_n)*u.Msun)*((period*u.yr)**2))/(4*np.pi**2))
sma = sma.to(u.AU).value
return sma
[docs] def to_semi_amp_basis(self, param_arr):
'''
Convert parameter array from standard basis to semi-amp basis by swapping out the
semi-major axis parameter to period for each companion and computing the stellar
rv semi-amplitudes for each companion.
Args:
param_arr (np.array of float): RxM array of fitting parameters in the period basis,
where R is the number of parameters being fit, and M is the number of orbits. If
M=1 (for MCMC), this can be a 1D array.
Returns:
np.array of float: modifies 'param_arr' to contain the semi-major axis for each companion
in each orbit rather than period, appends stellar rv semi-amplitude parameters, and
removes companion masses
'''
for body_num in np.arange(self.num_secondary_bodies) + 1:
# Grab necessary parameters for conversion
sma = param_arr[self.standard_basis_idx['sma{}'.format(body_num)]]
ecc = param_arr[self.standard_basis_idx['ecc{}'.format(body_num)]]
inc = param_arr[self.standard_basis_idx['inc{}'.format(body_num)]]
m_n = param_arr[self.standard_basis_idx['m{}'.format(body_num)]]
m0 = param_arr[self.standard_basis_idx['m0']]
mtot = m_n + m0
# Get stellar semi-amplitude
K_n = (np.sqrt(consts.G / (1 - ecc**2)))*(m_n*u.Msun)*(np.sin(inc))*((mtot*u.Msun)**(-1/2))*((sma*u.AU)**(-1/2))
kms = u.km / u.s
K_n = K_n.to(kms).value
# Compute Period replace in array
per = np.sqrt((4*(np.pi**2)*(sma*u.AU)**3) / (consts.G*(mtot*u.Msun)))
per = per.to(u.year).value
param_arr[self.param_idx['per{}'.format(body_num)]] = per
# Replace companion mass with semi-amplitude
param_arr[self.param_idx['K{}'.format(body_num)]] = K_n
return param_arr
[docs] def verify_params(self):
'''
Additionally warns that this basis will sample stellar mass rather than sample mass
regardless of whether 'fit_secondary_mass' is True or not.
'''
super(SemiAmp, self).verify_params()
if not self.fit_secondary_mass:
warnings.warn("This basis will not sample total mass. It will sample stellar mass.")
[docs]class XYZ(Basis):
'''
Defines an orbit using the companion's position and velocity components in XYZ space (x, y, z, xdot, ydot, zdot).
The conversion algorithms used for this basis are defined in the following paper:
http://www.dept.aoe.vt.edu/~lutze/AOE4134/9OrbitInSpace.pdf
.. Note:: Does not have support with sep,pa data yet.
.. Note:: Does not work for all multi-body data.
Args:
stellar_or_system_mass (float): mass of the primary star (if fitting for
dynamical masses of both components) or total system mass (if
fitting using relative astrometry only) [M_sol]
mass_err (float): uncertainty on 'stellar_or_system_mass', in M_sol
plx (float): mean parallax of the system, in mas
plx_err (float): uncertainty on 'plx', in mas
num_secondary_bodies (int): number of secondary bodies in the system, should be at least 1
fit_secondary_mass (bool): if True, include the dynamical mass of orbitting body as fitted parameter, if False,
'stellar_or_system_mass' is taken to be total mass
input_table (astropy.table.Table): output from 'orbitize.read_input.read_file()'
best_epoch_idx (list): indices of the epochs corresponding to the smallest uncertainties
epochs (list): all of the astrometric epochs from 'input_table'
angle_upperlim (float): either pi or 2pi, to restrict the prior range for 'pan' parameter (default: 2*pi)
hipparcos_IAD (orbitize.HipparcosLogProb object): if not 'None', then add relevant priors to this data (default: None)
rv (bool): if True, then there is radial velocity data and assign radial velocity priors, if False, then there
is no radial velocity data and radial velocity priors are not assigned (default: False)
rv_instruments (np.array): array of unique rv instruments from the originally supplied data (default: None)
Author: Rodrigo
'''
def __init__(
self, stellar_or_system_mass, mass_err, plx, plx_err, num_secondary_bodies,
fit_secondary_mass, data_table, best_epoch_idx, epochs,
angle_upperlim=2*np.pi, hipparcos_IAD=None, rv=False,
rv_instruments=None
):
super(XYZ, self).__init__(stellar_or_system_mass, mass_err, plx, plx_err,
num_secondary_bodies, fit_secondary_mass, angle_upperlim,
hipparcos_IAD, rv, rv_instruments
)
self.data_table = data_table
self.best_epoch_idx = best_epoch_idx
self.epochs = epochs
[docs] def construct_priors(self):
'''
Generates the parameter label array and initializes the corresponding priors for each
parameter that's to be sampled. For the xyz basis, the parameters common to each
companion are: x, y, z, xdot, ydot, zdot. Parallax, hipparcos (optional), rv (optional),
and mass priors are added at the end.
The xyz basis describes the position and velocity vectors with reference to the local coordinate
system (the origin of the system is star).
Returns:
tuple:
list: list of strings (labels) that indicate the names of each parameter to sample
list: list of orbitize.priors.Prior objects that indicate the prior distribution of each label
'''
basis_priors = []
basis_labels = []
# Add priors for the cartesian state vectors
for body in np.arange(self.num_secondary_bodies):
# Get the epoch with the least uncertainty for this body
# curr_idx = self.body_indices[body_num]
# radec_uncerts = self.data_table['quant1_err'][curr_idx] + self.data_table['quant2_err'][curr_idx]
# min_uncert = np.where(radec_uncerts == np.amin(radec_uncerts))[0]
# best_idx = curr_idx[0][min_uncert[0]]
datapoints_to_take = 3
best_idx = self.best_epoch_idx[body]
best_epochs = self.epochs[best_idx:(best_idx+datapoints_to_take)] # 0 is best, the others are for fitting velocity
# Get data near best epoch ASSUMING THE BEST IS NOT ONE OF THE LAST TWO EPOCHS OF A GIVEN BODY
# also assuming this is in radec
best_ras = self.data_table['quant1'][best_idx:(best_idx+datapoints_to_take)].copy()
best_ras_err = self.data_table['quant1_err'][best_idx:(best_idx+datapoints_to_take)].copy()
best_decs = self.data_table['quant2'][best_idx:(best_idx+datapoints_to_take)].copy()
best_decs_err = self.data_table['quant2_err'][best_idx:(best_idx+datapoints_to_take)].copy()
# Convert to AU for prior limits
best_xs = best_ras / self.plx
best_ys = best_decs / self.plx
best_xs_err = np.sqrt((best_ras_err / best_ras)**2 + (self.plx_err / self.plx)**2)*np.absolute(best_xs)
best_ys_err = np.sqrt((best_decs_err / best_decs)**2 + (self.plx_err / self.plx)**2)*np.absolute(best_ys)
# Least-squares fit on velocity for prior limits
A = np.vander(best_epochs, 2)
ATA_x = np.dot(A.T, A / (best_xs_err ** 2)[:, None])
cov_x = np.linalg.inv(ATA_x)
w_x = np.linalg.solve(ATA_x, np.dot(A.T, best_xs / best_xs_err ** 2))
ATA_y = np.dot(A.T, A / (best_ys_err ** 2)[:, None])
cov_y = np.linalg.inv(ATA_y)
w_y = np.linalg.solve(ATA_y, np.dot(A.T, best_ys / best_ys_err ** 2))
x_vel = w_x[0]
x_vel_err = np.sqrt(cov_x[0, 0])
y_vel = w_y[0]
y_vel_err = np.sqrt(cov_y[0, 0])
x_vel = (( x_vel* u.AU / u.day).to(u.km / u.s)).value
x_vel_err = ((x_vel_err * u.AU / u.day).to(u.km / u.s)).value
y_vel = ((y_vel * u.AU / u.day).to(u.km / u.s)).value
y_vel_err = ((y_vel_err * u.AU / u.day).to(u.km / u.s)).value
# Propose bounds on absolute Z and Z dot given the energy equation
mu = consts.G * self.stellar_or_system_mass * u.Msun
mu_vel = 2 * mu / ((x_vel**2 + y_vel**2) * (u.km / u.s * u.km / u.s))
z_bound = (np.sqrt(mu_vel**2 - (best_xs[0]**2 + best_ys[0]**2)*u.AU *u.AU)).to(u.AU)
z_bound = z_bound.value
mu_pos = 2 * mu / np.sqrt((best_xs[0]**2 + best_ys[0]**2) * (u.AU *u.AU))
z_vel_bound = (np.sqrt(mu_pos - (x_vel**2 + y_vel**2)*(u.km / u.s * u.km / u.s))).to(u.km / u.s)
z_vel_bound = z_vel_bound.value
# Add x-coordinate prior
num_uncerts_x = 5
basis_priors.append(priors.UniformPrior(best_xs[0] - num_uncerts_x*best_xs_err[0], best_xs[0] + num_uncerts_x*best_xs_err[0]))
basis_labels.append('x{}'.format(body+1))
# Add y-coordinate prior
num_uncerts_y = 5
basis_priors.append(priors.UniformPrior(best_ys[0] - num_uncerts_y*best_ys_err[0], best_ys[0] + num_uncerts_y*best_ys_err[0]))
basis_labels.append('y{}'.format(body+1))
# Add z-coordinate prior
# self.sys_priors.append(priors.UniformPrior(-z_bound,z_bound))
# self.sys_priors.append(priors.LogUniformPrior(0.0001,z_bound))
basis_priors.append(priors.GaussianPrior(0.,z_bound / 4, no_negatives=False))
basis_labels.append('z{}'.format(body+1))
# Add x-velocity prior
num_uncerts_xvel = 5
basis_priors.append(priors.UniformPrior(x_vel - num_uncerts_xvel*x_vel_err, x_vel + num_uncerts_xvel*x_vel_err))
basis_labels.append('xdot{}'.format(body+1))
# Add y-velocity prior
num_uncerts_yvel = 5
basis_priors.append(priors.UniformPrior(y_vel - num_uncerts_yvel*y_vel_err, y_vel + num_uncerts_yvel*y_vel_err))
basis_labels.append('ydot{}'.format(body+1))
# Add z-velocity prior
# self.sys_priors.append(priors.UniformPrior(-z_vel_bound,z_vel_bound))
# self.sys_priors.append(priors.LogUniformPrior(0.0001,z_vel_bound))
basis_priors.append(priors.GaussianPrior(0.,z_vel_bound / 4, no_negatives=False))
basis_labels.append('zdot{}'.format(body+1))
# Add parallax prior
basis_labels.append('plx')
if self.plx_err > 0:
basis_priors.append(priors.GaussianPrior(self.plx, self.plx_err))
else:
basis_priors.append(self.plx)
# Add hip_iad priors if necessary
if self.hipparcos_IAD is not None:
self.set_hip_iad_priors(basis_priors, basis_labels)
# Add rv priors
if self.rv and self.fit_secondary_mass:
self.set_rv_priors(basis_priors, basis_labels)
# Add mass priors
self.set_default_mass_priors(basis_priors, basis_labels)
# Define param label dictionary in current basis & standard basis
self.param_idx = dict(zip(basis_labels, np.arange(len(basis_labels))))
self.standard_basis_idx = dict(zip(basis_labels, np.arange(len(basis_labels))))
for body_num in np.arange(self.num_secondary_bodies) + 1:
self.standard_basis_idx['sma{}'.format(body_num)] = self.param_idx['x{}'.format(body_num)]
self.standard_basis_idx['ecc{}'.format(body_num)] = self.param_idx['y{}'.format(body_num)]
self.standard_basis_idx['inc{}'.format(body_num)] = self.param_idx['z{}'.format(body_num)]
self.standard_basis_idx['aop{}'.format(body_num)] = self.param_idx['xdot{}'.format(body_num)]
self.standard_basis_idx['pan{}'.format(body_num)] = self.param_idx['ydot{}'.format(body_num)]
self.standard_basis_idx['tau{}'.format(body_num)] = self.param_idx['zdot{}'.format(body_num)]
return basis_priors, basis_labels
[docs] def to_standard_basis(self, param_arr):
'''
Makes a call to 'xyz_to_standard' to convert each companion's xyz parameters
to the standard parameters an returns the updated array for conversion.
Args:
param_arr (np.array of float): RxM array of fitting parameters in the period basis,
where R is the number of parameters being fit, and M is the number of orbits. If
M=1 (for MCMC), this can be a 1D array.
Return:
np.array: Orbital elements in the standard basis for all companions.
'''
for body_num in np.arange(self.num_secondary_bodies)+1:
best_idx = self.best_epoch_idx[body_num - 1]
constrained_epoch = self.epochs[best_idx]
# Total mass is the sum of companion and stellar
if self.fit_secondary_mass:
secondary_m = param_arr[self.param_idx['m{}'.format(body_num)]]
m0 = param_arr[self.param_idx['m0']]
mtot = m0 + secondary_m
else:
mtot = param_arr[self.param_idx['mtot']]
to_convert = np.array([
param_arr[self.param_idx['x{}'.format(body_num)]],
param_arr[self.param_idx['y{}'.format(body_num)]],
param_arr[self.param_idx['z{}'.format(body_num)]],
param_arr[self.param_idx['xdot{}'.format(body_num)]],
param_arr[self.param_idx['ydot{}'.format(body_num)]],
param_arr[self.param_idx['zdot{}'.format(body_num)]],
param_arr[self.param_idx['plx']],
mtot
])
standard_params = self.xyz_to_standard(constrained_epoch, to_convert)
# Update param_arr to hold standard parameters
param_arr[self.standard_basis_idx['sma{}'.format(body_num)]] = standard_params[0]
param_arr[self.standard_basis_idx['ecc{}'.format(body_num)]] = standard_params[1]
param_arr[self.standard_basis_idx['inc{}'.format(body_num)]] = standard_params[2]
param_arr[self.standard_basis_idx['aop{}'.format(body_num)]] = standard_params[3]
param_arr[self.standard_basis_idx['pan{}'.format(body_num)]] = standard_params[4]
param_arr[self.standard_basis_idx['tau{}'.format(body_num)]] = standard_params[5]
param_arr[self.standard_basis_idx['plx']] = standard_params[6]
param_arr[self.standard_basis_idx['mtot']] = standard_params[7]
return param_arr
[docs] def xyz_to_standard(self, epoch, elems, tau_ref_epoch=58849):
"""
Converts array of orbital elements in terms of position and velocity in
xyz to the standard basis.
Args:
epoch (float): Date in MJD of observation to calculate time of
periastron passage (tau).
elems (np.array of floats): Orbital elements in xyz basis
(x-coordinate [au], y-coordinate [au], z-coordinate [au],
velocity in x [km/s], velocity in y [km/s], velocity in z [km/s],
parallax [mas], total mass of the two-body orbit (M_* + M_planet)
[Solar masses]). If more than 1 set of parameters is passed, the
first dimension must be the number of orbital parameter sets,
and the second the orbital elements.
Return:
np.array: Orbital elements in the standard basis
(sma, ecc, inc, aop, pan, tau, plx, mtot)
"""
if elems.ndim == 1:
elems = elems[:, np.newaxis]
# Velocities and positions, with units
vel = elems[3:6, :] * u.km / u.s # velocities in km / s ?
pos = elems[0:3, :] * u.AU # positions in AU ?
vel_magnitude = np.linalg.norm(vel, axis=0)
pos_magnitude = np.linalg.norm(pos, axis=0)
# Mass
mtot = elems[7, :]*u.Msun
mu = consts.G * mtot # G in m3 kg-1 s-2, mtot in msun
# Angular momentum, making sure nodal vector is not exactly zero
h = (np.cross(pos, vel, axis=0)).si
# if h[0].value == 0.0 and h[1].value == 0.0:
# pos[2] = 1e-8*u.AU
# h = (np.cross(pos, vel)).si
h_magnitude = np.linalg.norm(h, axis=0)
sma = 1 / (2.0 / pos_magnitude - (vel_magnitude**2)/mu)
sma = sma.to(u.AU)
ecc = (np.sqrt(1 - h_magnitude**2 / (sma * mu))).value
e_vector = (np.cross(vel, h, axis=0) / mu - pos / pos_magnitude).si
e_vec_magnitude = np.linalg.norm(e_vector, axis=0)
unit_k = np.array((0,0,1))[:, None]
cos_inc = (np.sum(h*unit_k, axis=0) / h_magnitude).value
inc = np.arccos(-cos_inc) # Take arccos of positive cos_inc?
#Nodal vector
n = np.cross(unit_k, h, axis=0)
n_magnitude = np.linalg.norm(n, axis=0)
# Position angle of the nodes, checking for the right quadrant
# np.arccos yields angles in [0, pi]
unit_i = np.array((1,0,0))[:, None]
unit_j = np.array((0,1,0))[:, None]
cos_pan = (np.sum(n*unit_j, axis=0) / n_magnitude).value # take dot product with i?
pan = np.arccos(cos_pan)
n_x = np.sum(n*unit_i, axis=0)
pan[n_x < 0.0] = 2*np.pi - pan[n_x < 0.0]
# Argument of periastron, checking for the right quadrant
cos_aop = (np.sum(n*e_vector, axis=0) / (n_magnitude * e_vec_magnitude)).value
aop = np.arccos(cos_aop)
e_vector_z = np.sum(e_vector*unit_k, axis=0)
aop[e_vector_z < 0.0] = 2.0*np.pi - aop[e_vector_z < 0.0]
# True anomaly, checking for the right quadrant
cos_tanom = (np.sum(pos*e_vector, axis=0) / (pos_magnitude*e_vec_magnitude)).value
tanom = np.arccos(cos_tanom)
# Check for places where tanom is nan, due to cos_tanom=1. (for some reason that was a problem)
# tanom = np.where((0.9999<cos_tanom) & (cos_tanom<1.001), 0.0, tanom)
rdotv = np.sum(pos*vel, axis=0)
tanom[rdotv < 0.0] = 2*np.pi - tanom[rdotv < 0.0]
# Eccentric anomaly to get tau, checking for the right quadrant
cos_eanom = ((1 - pos_magnitude / sma) / ecc).value
eanom = np.arccos(cos_eanom)
# Check for places where eanom is nan, due to cos_eanom = 1.(same problem as above)
# eanom = np.where((0.9999<cos_eanom ) & (cos_eanom<1.001), 0.0, eanom)
eanom[tanom > np.pi] = 2*np.pi - eanom[tanom > np.pi]
# Time of periastron passage, using Kepler's equation, in MJD:
time_tau = epoch - ((np.sqrt(sma**3 / mu)).to(u.day)).value * (eanom - ecc*np.sin(eanom))
# Calculate period from Kepler's third law, in days:
period = np.sqrt(4*np.pi**2.0*(sma)**3/mu)
period = period.to(u.day).value
# Finally, tau
tau = (time_tau - tau_ref_epoch) / period
tau = tau%1.0
mtot = mtot.value
sma = sma.value
results = np.stack([sma, ecc, inc, aop, pan, tau, elems[6, :], mtot])
return np.squeeze(results)
[docs] def to_xyz_basis(self, param_arr):
'''
Makes a call to 'standard_to_xyz' to convert each companion's standard keplerian parameters
to the xyz parameters an returns the updated array for conversion.
Args:
param_arr (np.array of float): RxM array of fitting parameters in the period basis,
where R is the number of parameters being fit, and M is the number of orbits. If
M=1 (for MCMC), this can be a 1D array.
Return:
np.array: Orbital elements in the xyz for all companions.
'''
for body_num in np.arange(self.num_secondary_bodies)+1:
best_idx = self.best_epoch_idx[body_num - 1]
constrained_epoch = self.epochs[best_idx]
# Get total mass
if self.fit_secondary_mass:
secondary_m = param_arr[self.param_idx['m{}'.format(body_num)]]
m0 = param_arr[self.standard_basis_idx['m0']]
mtot = m0 + secondary_m
else:
mtot = param_arr[self.param_idx['mtot']]
# Make conversion
to_convert = np.array([
param_arr[self.standard_basis_idx['sma{}'.format(body_num)]],
param_arr[self.standard_basis_idx['ecc{}'.format(body_num)]],
param_arr[self.standard_basis_idx['inc{}'.format(body_num)]],
param_arr[self.standard_basis_idx['aop{}'.format(body_num)]],
param_arr[self.standard_basis_idx['pan{}'.format(body_num)]],
param_arr[self.standard_basis_idx['tau{}'.format(body_num)]],
param_arr[self.standard_basis_idx['plx']],
mtot
])
xyz_params = self.standard_to_xyz(constrained_epoch, to_convert)
# Update param_arr to hold xyz parameters
param_arr[self.param_idx['x{}'.format(body_num)]] = xyz_params[0]
param_arr[self.param_idx['y{}'.format(body_num)]] = xyz_params[1]
param_arr[self.param_idx['z{}'.format(body_num)]] = xyz_params[2]
param_arr[self.param_idx['xdot{}'.format(body_num)]] = xyz_params[3]
param_arr[self.param_idx['ydot{}'.format(body_num)]] = xyz_params[4]
param_arr[self.param_idx['zdot{}'.format(body_num)]] = xyz_params[5]
param_arr[self.param_idx['plx']] = xyz_params[6]
param_arr[self.param_idx['mtot']] = xyz_params[7]
return param_arr
[docs] def standard_to_xyz(self, epoch, elems, tau_ref_epoch=58849, tolerance=1e-9, max_iter=100):
"""
Converts array of orbital elements from the regular base of Keplerian orbits to positions and velocities in xyz
Uses code from orbitize.kepler
Args:
epoch (float): Date in MJD of observation to calculate time of periastron passage (tau).
elems (np.array of floats): Orbital elements (sma, ecc, inc, aop, pan, tau, plx, mtot).
If more than 1 set of parameters is passed, the first dimension must be
the number of orbital parameter sets, and the second the orbital elements.
Return:
np.array: Orbital elements in xyz (x-coordinate [au], y-coordinate [au], z-coordinate [au],
velocity in x [km/s], velocity in y [km/s], velocity in z [km/s], parallax [mas], total mass of the two-body orbit
(M_* + M_planet) [Solar masses])
"""
# Use classical elements to obtain position and velocity in the perifocal coordinate system
# Then transform coordinates using matrix multiplication
if elems.ndim == 1:
elems = elems[:, np.newaxis]
# Define variables
sma = elems[0,:] # AU
ecc = elems[1,:] # [0.0, 1.0]
inc = elems[2,:] # rad [0, pi]
aop = elems[3,:] # rad [0, 2 pi]
pan = elems[4,:] # rad [0, 2 pi]
tau = elems[5,:] # [0.0, 1.0]
mtot = elems[7,:] # Msun
# Just in case so nothing breaks
ecc = np.where(ecc == 0.0, 1e-8, ecc)
inc = np.where(inc == 0.0, 1e-8, inc)
# Begin by calculating the eccentric anomaly
period = np.sqrt(4*np.pi**2.0*(sma*u.AU)**3/(consts.G*(mtot*u.Msun)))
period = period.to(u.day).value # Period in days
mean_motion = 2*np.pi/(period)
# Mean anomaly
manom = (mean_motion*(epoch - tau_ref_epoch) - 2*np.pi*tau) % (2.0*np.pi)
# Eccentric anomaly
eanom = kepler._calc_ecc_anom(manom, ecc, tolerance=tolerance, max_iter=max_iter)
# if eanom.ndim == 1:
# eanom = eanom[np.newaxis, :]
# Magnitude of angular momentum:
h = np.sqrt(consts.G*(mtot*u.Msun)*(sma*u.AU)*(1 - ecc**2))
# Position vector in the perifocal system in AU
pos_peri_x = (sma*(np.cos(eanom) - ecc))
pos_peri_y = (sma*np.sqrt(1 - ecc**2)*np.sin(eanom))
pos_peri_z = np.zeros(len(pos_peri_x))
pos = np.stack((pos_peri_x, pos_peri_y, pos_peri_z)).T
pos_magnitude = np.linalg.norm(pos, axis=1)
# Velocity vector in the perifocal system in km/s
vel_peri_x = - ((np.sqrt(consts.G*(mtot*u.Msun)*(sma*u.AU))*np.sin(eanom) / (pos_magnitude*u.AU)).to(u.km / u.s)).value
vel_peri_y = ((h* np.cos(eanom) / (pos_magnitude*u.AU)).to(u.km / u.s)).value
vel_peri_z = np.zeros(len(vel_peri_x))
vel = np.stack((vel_peri_x, vel_peri_y, vel_peri_z)).T
# Transformation matrix to inertial xyz system, component by component
pan = pan +np.pi / 2.0
T_11 = np.cos(pan)*np.cos(aop) - np.sin(pan)*np.sin(aop)*np.cos(inc)
T_12 = - np.cos(pan)*np.sin(aop) - np.sin(pan)*np.cos(aop)*np.cos(inc)
T_13 = np.sin(pan)*np.sin(inc)
T_21 = np.sin(pan)*np.cos(aop) + np.cos(pan)*np.sin(aop)*np.cos(inc)
T_22 = - np.sin(pan)*np.sin(aop) + np.cos(pan)*np.cos(aop)*np.cos(inc)
T_23 = - np.cos(pan)*np.sin(inc)
T_31 = np.sin(aop)*np.sin(inc)
T_32 = np.cos(aop)*np.sin(inc)
T_33 = np.cos(inc)
T = np.array([[T_11, T_12, T_13],
[T_21, T_22, T_23],
[T_31, T_32, T_33]])
pos_xyz = np.zeros((len(sma), 3))
vel_xyz = np.zeros((len(sma), 3))
for k in range(len(sma)):
pos_xyz[k,:] = np.matmul(T[:,:,k], pos[k])
vel_xyz[k,:] = np.matmul(T[:,:,k], vel[k])
result = np.stack([-pos_xyz[:,0], pos_xyz[:,1], pos_xyz[:,2], -vel_xyz[:,0], vel_xyz[:,1], vel_xyz[:,2], elems[6, :], mtot])
if len(sma) == 1:
result = result.T
return np.squeeze(result)
[docs] def verify_params(self):
'''
For now, additionally throws exceptions if data is supplied in sep/pa or if the best epoch for each
body is one of the last two (which would inevtably mess up how the priors are imposed).
'''
super(XYZ, self).verify_params()
# For now, raise error if data is in sep/pa
seppa_locs = np.where(self.data_table['quant_type'] == 'seppa')
if np.size(seppa_locs) != 0:
raise Exception("For now, the XYZ basis requires data in RA and DEC offsets.")
# For now, raise error if the best epoch for each body is one of the last two
for i in range(self.num_secondary_bodies):
body_num = i + 1
best_epoch_loc = self.best_epoch_idx[i]
body_indices = np.where(self.data_table['object'] == body_num)[0]
max_index = np.amax(body_indices)
if (max_index - best_epoch_loc < 2):
raise Exception("For now, the epoch with the lowest sepparation error should not be one of the last two entries for body{}".format(body_num))
[docs]def tau_to_tp(tau, ref_epoch, period, after_date=None):
"""
Convert tau (epoch of periastron in fractional orbital period after ref epoch) to
t_p (date in days, usually MJD, but works with whatever system ref_epoch is given in)
Args:
tau (float or np.array): value of tau to convert
ref_epoch (float or np.array): date (in days, typically MJD) that tau is defined relative to
period (float or np.array): period (in years) that tau is noralized with
after_date (float): tp will be the first periastron after this date. If None, use ref_epoch.
Returns:
float or np.array: corresponding t_p of the taus
"""
period_days = period * u.year.to(u.day)
tp = tau * (period_days) + ref_epoch
if after_date is not None:
num_periods = (after_date - tp)/period_days
num_periods = int(np.ceil(num_periods))
tp += num_periods * period_days
return tp
[docs]def tp_to_tau(tp, ref_epoch, period):
"""
Convert t_p to tau
Args:
tp (float or np.array): value to t_p to convert (days, typically MJD)
ref_epoch (float or np.array): reference epoch (in days) that tau is defined from. Same system as tp (e.g., MJD)
period (float or np.array): period (in years) that tau is defined by
Returns:
float or np.array: corresponding taus
"""
tau = (tp - ref_epoch)/(period * u.year.to(u.day))
tau %= 1
return tau
[docs]def switch_tau_epoch(old_tau, old_epoch, new_epoch, period):
"""
Convert tau to another tau that uses a different referench epoch
Args:
old_tau (float or np.array): old tau to convert
old_epoch (float or np.array): old reference epoch (days, typically MJD)
new_epoch (float or np.array): new reference epoch (days, same system as old_epoch)
period (float or np.array): orbital period (years)
Returns:
float or np.array: new taus
"""
tp = tau_to_tp(old_tau, old_epoch, period)
new_tau = tp_to_tau(tp, new_epoch, period)
return new_tau
[docs]def tau_to_manom(date, sma, mtot, tau, tau_ref_epoch):
"""
Gets the mean anomlay. Wrapper for kepler.tau_to_manom, kept here
for backwards compatibility.
Args:
date (float or np.array): MJD
sma (float): semi major axis (AU)
mtot (float): total mass (M_sun)
tau (float): epoch of periastron, in units of the orbital period
tau_ref_epoch (float): reference epoch for tau
Returns:
float or np.array: mean anomaly on that date [0, 2pi)
"""
return kepler.tau_to_manom(date, sma, mtot, tau, tau_ref_epoch)