Classes

class orbdot.star_planet.StarPlanet(settings_file, planet_num=0)

A StarPlanet class instance represents a star-planet system and acts as an interface for the core capabilities of the OrbDot package. It combines the data, methods, and attributes needed to run model fitting routines and interpret the results.

__init__(settings_file, planet_num=0)

Initializes the StarPlanet class.

Parameters:

• settings_file (str) – Path to the main settings file.

• planet_num (int, optional) – Planet number in case of multi-planet systems (default is 0).

update_default(parameter, new_value)

Updates the default (fixed) value for the specified parameter.

The default value will be used in a model fit if the parameter is not allowed to vary.

Parameters:

• parameter (str) – The parameter name.

• new_value (float) – The new parameter value.

Returns:

None – The default value for the specified parameter is updated.

update_prior(parameter, new_prior)

Updates the prior distribution for the specified parameter.

Parameters:

• parameter (str) – The parameter name.

• new_prior (list) – A list three of values specifying the prior distribution, where the first element is the type of prior ("uniform", "gaussian", or "log"), and subsequent elements define the distribution.

Returns:

None – The prior distribution for the specified parameter is updated.

---

class orbdot.analysis.Analyzer(planet, results_dic)

This class enables various computations related to the long-term variations of exoplanet orbits. It combines model fit results, star-planet system characteristics, and data to compute and summarize analyses of various physical models, such as the effects of equilibrium tides, apsidal precession, systemic proper motion, and companion objects.

__init__(planet, results_dic)

Initializes the Analyzer class.

Parameters:

• planet (object) – An instance of the StarPlanet class.

• results_dic (dict) – A dictionary containing the results of an OrbDot model fit.

apsidal_precession_fit(printout=False)

Interprets the results of an apsidal precession model fit.

This method produces a concise summary of various interpretations of the results of an apsidal precession model fit, independent of the data type(s) that are fit.

Parameters:

printout (bool, optional) – An option to print the results to the console.

Returns:

None – The results are written to a text file.

Raises:

IndexError – If the results of an apsidal precession model fit are not available.

apsidal_precession_predicted(printout=False)

Calculates various apsidal precession rates that are predicted by theory.

This method produces a concise summary of the expected rates of apsidal precession due to general relativistic effects, tides, and rotation.

Parameters:

printout (bool, optional) – An option to print the results to the console.

Returns:

None – The results are written to a text file.

model_comparison(model_2_results, printout=False)

Compares the Bayesian evidence with that of another model fit.

To compare two models, Model 1 and Model 2, this method calculates the Bayes factor, denoted as:

\[\log{B_{12}} = \log{\mathrm{Z}}_{1} - \log{\mathrm{Z}}_{2}\]

where \(\log{\mathrm{Z}}\) is the Bayesian evidence. The Bayes factor is then evaluated against the thresholds established by Kass and Raftery (1995) [1]_.

Parameters:

• model_2_results (dict) – The dictionary returned by the alternative model fit.

• printout (bool, optional) – An option to print the results to the console.

Returns:

None – The results are written to a text file.

References

[1]

Kass and Raftery [1995]. https://doi.org/10.2307/2291091

orbital_decay_fit(printout=False)

Interprets the results of an orbital decay model fit.

This method produces a concise summary of various interpretations of the results of an orbital decay model fit, independent of the data type(s) that are fit.

Parameters:

printout (bool, optional) – An option to print the results to the console.

Returns:

None – The results are written to a text file.

Raises:

IndexError – If the results of an orbital decay model fit are not available.

orbital_decay_predicted(printout=False)

Calculates various orbital decay parameters that are predicted by theory.

This method produces a concise summary of various orbital decay characteristics that are predicted by equilibrium tidal theory.

Parameters:

printout (bool, optional) – An option to print the results to the console.

Returns:

None – The results are written to a text file.

proper_motion(printout=False)

Calculates the expected TTVs and TDVs due to systemic proper motion.

This method produces a concise summary of the apparent transit timing and duration variations that are expected due to systemic proper motion.

Parameters:

printout (bool, optional) – An option to print the results to the console.

Returns:

None – The results are written to a text file.

resolved_binary(separation, secondary_mass=None, printout=False)

Calculates and summarizes observable effects and properties of a resolved companion star.

Parameters:

• separation (float) – The angular separation of the binary in arcseconds.

• secondary_mass (float, optional) – The mass of the stellar companion in solar masses.

• printout (bool, optional) – An option to print the results to the console.

Returns:

None – The results are written to a text file.

rv_trend_linear()

Plots the best-fit linear trend over the RV residuals.

Notes

The full radial velocity signal is expressed as:

\[v_r = K[\cos{(\phi\left(t\right)+\omega_p)}+e\cos{\omega_p}] + \gamma_j + \dot{\gamma} \left(t-t_0\right) + \frac{1}{2} \ddot{\gamma} \left( t-t_0\right)^2\]

where \(\dot{\gamma}\) and \(\ddot{\gamma}\) are first and second-order acceleration terms, respectively.

After subtracting the contribution from the planet and the systemic velocity \(\gamma\), the residuals are only the long-term trend [1]_:

\[RV_c(t) = \dot{\gamma} \left(t-t_0\right) + \frac{1}{2} \ddot{\gamma} \left( t-t_0\right)^2\]

In this case, \(\ddot{\gamma} = 0\).

rv_trend_quadratic()

Estimates the minimum period of an outer companion given a quadratic fit to RV residuals.

This method analyzes the best-fit quadratic radial velocity trend to estimate the minimum orbital period of an outer companion that could cause such acceleration, assuming that the orbit is circular. It is designed to work in tandem with the companion_from_quadratic_rv() method, which uses Equations 1, 3, and 4 of Kipping et al. (2011) [1]_.

Returns:

float – Lower limit for the orbital period of an outer companion in days.

Notes

The full radial velocity signal is expressed as:

\[v_r = K[\cos{(\phi\left(t\right)+\omega_p)}+e\cos{\omega_p}] + \gamma_j + \dot{\gamma} \left(t-t_0\right) + \frac{1}{2} \ddot{\gamma} \left( t-t_0\right)^2\]

where \(\dot{\gamma}\) and \(\ddot{\gamma}\) are first and second-order acceleration terms, respectively.

After subtracting the contribution from the planet and the systemic velocity \(\gamma\), the residuals are only the long term trend [1]_:

\[RV_c(t) = \dot{\gamma} \left(t-t_0\right) + \frac{1}{2} \ddot{\gamma} \left( t-t_0\right)^2\]

If both \(\dot{\gamma}\) and \(\ddot{\gamma}\) are nonzero, residuals are quadratic, but if \(\ddot{\gamma} = 0\) they are linear. For the latter case, this method is not valid.

The minimum orbital period of the companion is loosely constrained by the length of the baseline of RV observations. Assuming that the companion’s orbit is circular, for which the signal is sinusoidal, the minimum period is approximated as 4x the span of time between the vertex of the quadratic trend and furthest edge of the observed baseline:

\[P_{\mathrm{min}} = 4 \times \mathrm{max} \left[|t_{\mathrm{vertex}} - t_{ \mathrm{min}}|, |t_{\mathrm{vertex}} - t_{\mathrm{max}}|\right]\]

where,

\[t_{\mathrm{vertex}} = -\frac{\dot{\gamma}}{2 \ddot{\gamma}} + t_{\mathrm{pivot}}\]

References

[1]

Kipping et al. [2011]. https://doi.org/10.1088/0004-6256/142/3/95

unknown_companion(p2_min=0.5, p2_max=10, p2_nout=10, printout=False)

Calculates and summarizes observable effects and properties of a nonresonant companion.

Parameters:

• p2_min (float, optional) – For precession analysis, the minimum companion period in days. Default is 0.5.

• p2_max (float, optional) – For precession analysis, the minimum companion period in days. Default is 10.

• p2_nout (int, optional) – Number of companion period values to assess in precession analysis. Default is 10.

• printout (bool, optional) – An option to print the results to the console.

Returns:

None – The results are printed to the console and written to a text file.

---

class orbdot.nested_sampling.NestedSampling(fixed_values, prior)

This class contains all of the methods required to run the model fitting routines defined in the TransitTiming, RadialVelocity, TransitDuration, and JointFit classes.

The user may choose between two packages, Nestle [1]_ or PyMultiNest [2]_. PyMultiNest is generally faster and more robust, but it can be tricky to install, thus it is not a requirement to use this code. The desired sampler is specified in the settings file as "nestle" or "multinest".

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest/

__init__(fixed_values, prior)

Initializes the NestedSampling class.

Parameters:

• fixed_values (dict) – A dictionary that specifies the fixed value for each parameter.

• prior (dict) – A dictionary that specifies the prior distributions for each parameter.

clear_index()

Clears the free parameter indices to prepare for the next model fit.

This method removes instance variables that store the index (order) of the free parameters, allowing them to be redefined in a subsequent model fit.

Returns:

None

generate_random_samples(weighted_samples, num=300)

Generates a set of random samples for plotting.

This function selects random samples from a given array of weighted posterior samples and retrieves the corresponding parameter values using the get_vals method.

Parameters:

• weighted_samples (array_like) – Array of weighted posterior samples.

• num (int, optional) – Number of random samples to generate. Default is 300.

Returns:

tuple –

A tuple with the following elements:

1. A list of the orbit parameter samples.

2. A list of time-dependent parameter samples.

3. A list of radial velocity parameter samples.

get_best_fit(samples)

Retrieves the 68% credible intervals on each parameter.

This method returns the 68% credible intervals on each free parameter with a given array of weighted posterior samples. If a model parameter was not allowed to vary in the model fit, its default value is recorded for completeness.

Parameters:

samples (array_like) – The weighted posterior samples.

Returns:

dict – A dictionary containing the best-fit parameter values.

Notes

If the user has chosen to fit "ecosw" and "esinw" or "sq_ecosw" and "sq_esinw", the corresponding "e0" and "w0" values are calculated.

get_index()

Retrieves the index (order) of the free parameters.

This method iterates through the list of free parameters and determines the index (order) of each free parameter, storing them in instance variables for later use.

Returns:

None

get_vals(theta)

Retrieves and returns values for the full set of OrbDot parameters.

This function combines and returns values for the entire set of OrbDot model parameters, allowing them to be easily passed to a physical model or log-likelihood function. For any parameter that is not allowed to vary in the model fit, its default value is recorded.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

• orbital_elements (list) – A list of the orbital element parameters, in the order: "t0", "P0", "e0", "w0", "i0", "O0".

• time_dependent (list) – A list of the time-dependent parameters, in the order: "PdE", "wdE", "edE", "idE", "OdE".

• radial_velocity (list) – A list of the radial velocity parameters, in the order: "K", "v0", "jit", "dvdt", "ddvdt", K_tide.

Notes

If the user has chosen to fit "ecosw" and "esinw" or "sq_ecosw" and "sq_esinw", the corresponding "e0" and "w0" values are calculated.

print_results(dic, sampler)

Print the results of the model fit to the console.

This method prints the results of the sampler to the console, including the best-fit parameter values and Bayesian evidence.

Parameters:

• dic (dict) – Dictionary containing the model fit results.

• sampler (str) – Name of the sapling package that was used, must be "nestle" or "multinest".

Returns:

None

prior_transform(theta)

Transforms a parameter from the unit hypercube value to a physical value.

This method transforms the current state of the free parameters from the unit hypercube ( in the range 0 to 1) to their true values based on the specified prior distributions. The transformed parameters may then be passed to the log-likelihood function by the sampler.

Parameters:

theta (array_like) – Array containing the current state of the free parameters in the unit hypercube.

Returns:

trans (array_like) – Array containing the transformed values of the free parameters.

run_multinest(loglike, free_params, n_points, tol, save_dir, run_num=0, resume=False)

Runs a model fit with the nested sampling package PyMultiNest.

This method is an overhead function that performs a model fit using the PyMultiNest Python package [1]_. PyMultiNest is imported within this function, so that it does not need to be installed if Nestle is preferred.

Parameters:

• loglike (callable) – A single-argument function that computes the log-likelihood for the desired model.

• free_params (list) – A list of free parameters.

• n_points (int) – The number of live points.

• tol (int) – The evidence tolerance.

• save_dir (str) – The directory to save the PyMultiNest chains.

• run_num (int, optional) – Number for tagging the PyMultiNest chains, if the user wishes to resume a run.

• resume (bool, optional) – Resumes a run of the sampler (by run number), rather than re-starting it.

Returns:

tuple –

A tuple with the following elements:

1. A dictionary containing the results of the model fit.

2. An array of the weighted posterior samples.

3. A set of 300 random samples for plotting.

References

[1]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest/

run_nestle(loglike, free_params, method, n_points, tol)

Runs a model fit with the nested sampling package Nestle.

This method is an overhead function that performs a model fit using the Nestle Python package [1]_. Nestle is imported within this function, so that it does not need to be installed if PyMultiNest is preferred.

Parameters:

• loglike (callable) – A single-argument function that computes the log-likelihood for the desired model.

• free_params (list) – A list of free parameters.

• method (str) – The Nestle sampling type, which may be “multi” for multiple ellipsoids, “single” for single ellipsoid, or “classic” for classic MCMC exploration.

• n_points (int) – The number of live points.

• tol (int) – The evidence tolerance.

Returns:

tuple –

A tuple with the following elements:

1. A dictionary containing the results of the model fit.

2. An array of the weighted posterior samples.

3. A set of 300 random samples for plotting.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

save_random_samples(random_samples, filename)

Overhead function that saves a random set of posterior samples for plotting.

Parameters:

• random_samples (array_like) – The randomly selected samples.

• filename (str) – The output file name.

Returns:

None

save_results(random_samples, weighted_samples, res_dic, free_params, sampler_type, suffix, prefix, not_model_params)

Saves the results of the model fit by generating a set of output files.

Parameters:

• random_samples (array-like) – A set of 300 random samples for plotting.

• weighted_samples (array-like) – The weighted posterior samples.

• res_dic (dict) – A dictionary containing the results of the model fit.

• free_params (list or tuple) – The free parameters.

• sampler_type (str) – The type of sampler used ("nestle" or "multinest").

• suffix (str) – A string appended to the end of the output files.

• prefix (str) – A string added to the beginning of the output files.

• not_model_params (list or tuple) – A list of OrbDot parameters that do not belong to the model.

Returns:

None – The following output files are generated:

1. *_summary.txt: a quick visual summary of the results.

2. *_results.json: the entire model fitting results dictionary.

3. *_corner.png: a corner plot.

4. *_weighted_samples.txt: the weighted posterior samples.

5. *_random_samples.json: a random set of 300 posterior samples for plotting.

save_summary(dic, filename, sampler, not_model_params)

Saves a summary of the nested sampling results.

This method summarizes the results of the model fit in an easy-to-read text file.

Parameters:

• dic (dict) – Dictionary containing the results of the sampler.

• filename (str) – Path to the directory for the output files.

• sampler (str) – The type of sampler used ("nestle" or "multinest").

• not_model_params (list or tuple) – A list of OrbDot parameters that do not belong to the model.

Returns:

None

save_weighted_samples(weighted_samples, filename)

Overhead function that saves the weighted posterior samples.

Parameters:

• weighted_samples (array_like) – Array containing the samples generated by the model fit.

• filename (str) – Name of the output file.

Returns:

None

---

class orbdot.joint_fit.JointFit(joint_settings)

This class utilizes the capabilities of the NestedSampling class to support joint model fitting of transit/eclipse mid-times, radial velocities, and transit durations.

__init__(joint_settings)

Initializes the JointFit class.

Parameters:

joint_settings (dict) – A dictionary specifying directories and settings for the nested sampling analysis.

run_joint_fit(free_params, TTV=False, RV=False, TDV=False, model='constant', file_suffix='', make_plot=True, sigma_clip=False, clip_model='linear')

Run a model fit on multiple data types simultaneously.

This method executes a model fit for any combination of data types using one of two nested sampling packages: Nestle [1]_ or PyMultiNest [2]_.

At least two data types must be specified when calling this method. Use the TTV=True argument to include transit and/or eclipse mid-times, RV=True to include radial velocities, and TDV=True to include transit durations.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be part of the physical model.

• TTV (bool) – If True, the transit and/or eclipse timing data are included in the model fit. Default is False.

• RV (bool) – If True, the radial velocity data are included in the model fit. Default is False.

• TDV (bool) – If True, the transit duration data are included in the model fit. Default is False.

• model (str, optional) – The timing model, must be "constant", "decay", or "precession". Default is "constant".

• file_suffix (str, optional) – A string appended to the end of the output file names.

• make_plot (bool, optional) – If True, relevant plots are generated. Default is True.

• sigma_clip (bool, optional) – Option to execute a sigma-clipping routine on the transit mid-times. Default is False.

• clip_model (str, optional) – Specifies the model to be used in the sigma-clipping routine. The options are "linear" or "quadratic", with the default being "linear".

Returns:

res (dict) – A dictionary containing the model fit results and settings.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_rv_tdv_fit(free_params, model, suffix, plot)

Run a joint RV/TDV model fit.

This method executes a simultaneous model fit of the radial velocity and transit duration data. It uses one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_. The physical model is specified by the model argument, which may be equal to "constant" for an unchanging orbit, "decay" for orbital decay, or "precession" for apsidal precession.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["t0", "P0", "e0", "w0", "ecosw", "esinw", "sq_ecosw", sq_esinw", "i0", "PdE", "wdE", "K", "v0", "jit", "dvdt", "ddvdt"]

• model (string) – The desired physical model, must be "constant", "decay", or "precession". Default is "constant".

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, TDV and RV plots are generated. Default is True.

Returns:

res (dict) – A dictionary containing the results of the fit for access within a script

Note

The following output files are generated:

1. "rv_tdv_<model>_summary.txt": a quick visual summary of the results

2. "rv_tdv_<model>_results.json": the entire model fitting results dictionary.

3. "rv_tdv_<model>_corner.png": a corner plot.

4. "rv_tdv_<model>_weighted_samples.txt": the weighted posterior samples.

5. "rv_tdv_<model>_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_ttv_rv_fit(free_params, model, suffix, plot, clip, clip_method, save)

Run a joint TTV/RV model fit.

This method executes a simultaneous model fit of the transit/eclipse mid-times and radial velocity data. It uses one of two nested sampling packages, Nestle [1]_ or PyMultiNest [ 2]_. The physical model is specified by the model argument, which may be equal to "constant" for an unchanging orbit, "decay" for orbital decay, or "precession" for apsidal precession.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["t0", "P0", "e0", "w0", "ecosw", "esinw", "sq_ecosw", sq_esinw", "PdE", "wdE", "K", "v0", "jit", "dvdt", "ddvdt"].

• model (string) – The desired physical model, must be "constant", "decay", or "precession". Default is "constant".

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, TTV and RV plots are generated. Default is True.

• clip (bool) – Option to execute a sigma-clipping routine on the transit mid-times. Default is False.

• clip_method (str) – Specifies the model to be used in the sigma-clipping routine. The options are "linear" or "quadratic", with the default being "linear".

• save (bool) – If False, the output files are not saved.

Returns:

res (dict) – A dictionary containing the results of the fit for access within a script

Note

The following output files are generated:

1. "ttv_rv_<model>_summary.txt": a quick visual summary of the results

2. "ttv_rv_<model>_results.json": the entire model fitting results dictionary.

3. "ttv_rv_<model>_corner.png": a corner plot.

4. "ttv_rv_<model>_weighted_samples.txt": the weighted posterior samples.

5. "ttv_rv_<model>_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_ttv_rv_tdv_fit(free_params, model, suffix, plot, clip, clip_method, save)

Run a joint TTV/RV/TDV model fit.

This method executes a simultaneous model fit of the transit/eclipse mid-times, radial velocity, and transit duration data. It uses one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_. The physical model is specified by the model argument, which may be equal to "constant" for an unchanging orbit, "decay" for orbital decay, or "precession" for apsidal precession.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["t0", "P0", "e0", "w0", "ecosw", "esinw", "sq_ecosw", sq_esinw", "i0", "PdE", "wdE", "K", "v0", "jit", "dvdt", "ddvdt"].

• model (string) – The desired physical model, must be "constant", "decay", or "precession". Default is "constant".

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, the TTV, RV, and TDV plots are generated. Default is True.

• clip (bool) – Option to execute a sigma-clipping routine on the transit mid-times. Default is False.

• clip_method (str) – Specifies the model to be used in the sigma-clipping routine. The options are "linear" or "quadratic", with the default being "linear".

• save (bool) – If False, the output files are not saved.

Returns:

res (dict) – A dictionary containing the results of the fit for access within a script

Note

The following output files are generated:

1. "ttv_rv_tdv_<model>_summary.txt": a quick visual summary of the results

2. "ttv_rv_tdv_<model>_results.json": the entire model fitting results dictionary.

3. "ttv_rv_tdv_<model>_corner.png": a corner plot.

4. "ttv_rv_tdv_<model>_weighted_samples.txt": the weighted posterior samples.

5. "ttv_rv_tdv_<model>_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_ttv_tdv_fit(free_params, model, suffix, plot, clip, clip_method, save)

Run a joint TTV/TDV model fit.

This method executes a simultaneous model fit of the transit/eclipse mid-times and transit duration data. It uses one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_. The physical model is specified by the model argument, which may be equal to "constant" for an unchanging orbit, "decay" for orbital decay, or "precession" for apsidal precession.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["t0", "P0", "e0", "w0", "ecosw", "esinw", "sq_ecosw", sq_esinw", "i0", "PdE", "wdE"].

• model (string) – The desired physical model, must be "constant", "decay", or "precession". Default is "constant".

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, TTV and TDV plots are generated. Default is True.

• clip (bool) – Option to execute a sigma-clipping routine on the transit mid-times. Default is False.

• clip_method (str) – Specifies the model to be used in the sigma-clipping routine. The options are "linear" or "quadratic", with the default being "linear".

• save (bool) – If False, the output files are not saved.

Returns:

res (dict) – A dictionary containing the results of the fit for access within a script

Note

The following output files are generated:

1. "ttv_tdv_<model>_summary.txt": a quick visual summary of the results

2. "ttv_tdv_<model>_results.json": the entire model fitting results dictionary.

3. "ttv_tdv_<model>_corner.png": a corner plot.

4. "ttv_tdv_<model>_weighted_samples.txt": the weighted posterior samples.

5. "ttv_tdv_<model>_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

rv_tdv_loglike_constant(theta)

Calculates the joint log-likelihood for the constant-period RV/TDV model fit.

This function returns the sum of the log-likelihoods calculated by the rv_loglike_constant() and tdv_loglike_constant() methods.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The sum of the RV and TDV log-likelihood values.

rv_tdv_loglike_decay(theta)

Calculates the joint log-likelihood for the orbital decay RV/TDV model fit.

This function returns the sum of the log-likelihoods calculated by the rv_loglike_decay() and tdv_loglike_decay() methods.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The sum of the RV and TDV log-likelihood values.

rv_tdv_loglike_precession(theta)

Calculates the joint log-likelihood for the apsidal precession RV/TDV model fit.

This function returns the sum of the log-likelihoods calculated by the rv_loglike_precession() and tdv_loglike_precession() methods.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The sum of the RV and TDV log-likelihood values.

ttv_rv_loglike_constant(theta)

Calculates the joint log-likelihood for the constant-period TTV/RV model fit.

This function returns the sum of the log-likelihoods calculated by the ttv_loglike_constant() and rv_loglike_constant() methods.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The sum of the TTV and RV log-likelihood values.

ttv_rv_loglike_decay(theta)

Calculates the joint log-likelihood for the orbital decay TTV/RV model fit.

This function returns the sum of the log-likelihoods calculated by the ttv_loglike_decay() and rv_loglike_decay() methods.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The sum of the TTV and RV log-likelihood values.

ttv_rv_loglike_precession(theta)

Calculates the joint log-likelihood for the apsidal precession TTV/RV model fit.

This function returns the sum of the log-likelihoods calculated by the ttv_loglike_precession() and rv_loglike_precession() methods.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The sum of the TTV and RV log-likelihood values.

ttv_rv_tdv_loglike_constant(theta)

Calculates the joint log-likelihood for the constant-period TTV/RV/TDV model fit.

This function returns the sum of the log-likelihoods calculated by the ttv_loglike_constant(), rv_loglike_constant(), and tdv_loglike_constant() methods.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The sum of the TTV, RV, and TDV log-likelihood values.

ttv_rv_tdv_loglike_decay(theta)

Calculates the joint log-likelihood for the orbital decay TTV/RV/TDV model fit.

This function returns the sum of the log-likelihoods calculated by the ttv_loglike_decay(), rv_loglike_decay(), and tdv_loglike_decay() methods.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The sum of the TTV, RV, and TDV log-likelihood values.

ttv_rv_tdv_loglike_precession(theta)

Calculates the joint log-likelihood for the apsidal precession TTV/RV/TDV model fit.

This function returns the sum of the log-likelihoods calculated by the ttv_loglike_precession(), rv_loglike_precession(), and tdv_loglike_precession() methods.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The sum of the TTV, RV, and TDV log-likelihood values.

ttv_tdv_loglike_constant(theta)

Calculates the joint log-likelihood for the constant-period TTV/TDV model fit.

This function returns the sum of the log-likelihoods calculated by the ttv_loglike_constant() and tdv_loglike_constant() methods.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The sum of the TTV and TDV log-likelihood values.

ttv_tdv_loglike_decay(theta)

Calculates the joint log-likelihood for the orbital decay TTV/TDV model fit.

This function returns the sum of the log-likelihoods calculated by the ttv_loglike_decay() and tdv_loglike_decay() methods.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The sum of the TTV and TDV log-likelihood values.

ttv_tdv_loglike_precession(theta)

Calculates the joint log-likelihood for the apsidal precession TTV/TDV model fit.

This function returns the sum of the log-likelihoods calculated by the ttv_loglike_precession() and tdv_loglike_precession() methods.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The sum of the TTV and TDV log-likelihood values.

---

class orbdot.transit_timing.TransitTiming(ttv_settings)

This class utilizes the capabilities of the NestedSampling class to facilitate model fitting of transit and eclipse mid-times.

__init__(ttv_settings)

Initializes the TransitTiming class.

Parameters:

ttv_settings (dict) – A dictionary specifying directories and settings for the nested sampling analysis.

clip(outfile_cleaned, outfile_clipped, method='linear', max_iters=20)

Remove outliers from the transit timing data.

This method runs an iterative sigma-clipping routine on the observed transit mid-times, originally developed in Hagey et al. (2022) [1]_. At every iteration, the best-fit model (linear or quadratic) is subtracted from the observed mid-times, and any data point that lies outside of 3 standard deviations from the mean are removed. This process is repeated until no remaining data points fit this criteria, or until a maximum number of iterations has been reached.

Parameters:

• outfile_cleaned (str) – File path for saving the the cleaned data.

• outfile_clipped (str) – File path for saving the the excluded (“clipped”) data.

• method (str, optional) – The timing model to be subtracted from the data at every iteration. The options are "linear" or "quadratic", with the default being "linear".

• max_iters (int, optional) – The maximum number of iterations. Default is 20.

References

[1]

Hagey, Edwards, and Boley (2022). https://doi.org/10.3847/1538-3881/ac959a

run_ttv_constant(free_params, suffix, plot, clip, clip_method, save)

Run a fit of the constant-period timing model.

This method executes a constant-period model fit of the observed transit and/or eclipse mid-times using one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["t0", "P0", "e0", "w0", "ecosw", "esinw", "sq_ecosw", sq_esinw"].

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, a TTV (O-C) plot is generated.

• clip (bool) – If True, a sigma-clipping routine is run on the transit mid-times before model fitting.

• clip_method (str) – Specifies the model to be used in the sigma-clipping routine, either "linear" or "quadratic".

• save (bool) – If False, the output files are not saved.

Returns:

res (dict) – A dictionary containing the model fit results and settings.

Note

The following output files are generated:

1. "ttv_constant_summary.txt": a quick visual summary of the results

2. "ttv_constant_results.json": the entire model fitting results dictionary.

3. "ttv_constant_corner.png": a corner plot.

4. "ttv_constant_weighted_samples.txt": the weighted posterior samples.

5. "ttv_constant_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_ttv_decay(free_params, suffix, plot, clip, clip_method, save)

Run a fit of the orbital decay timing model.

This method executes an orbital decay model fit of the observed transit and/or eclipse mid-times using one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["t0", "P0", "PdE", "e0", "w0", "ecosw", "esinw", "sq_ecosw", sq_esinw"].

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, a TTV (O-C) plot is generated.

• save (bool) – If False, the output files are not saved.

• clip (bool) – If True, a sigma-clipping routine is run on the transit mid-times before model fitting.

• clip_method (str) – Specifies the model to be used in the sigma-clipping routine, either "linear" or "quadratic".

Returns:

res (dict) – A dictionary containing the model fit results and settings.

Note

The following output files are generated:

1. "ttv_decay_summary.txt": a quick visual summary of the results

2. "ttv_decay_results.json": the entire model fitting results dictionary.

3. "ttv_decay_corner.png": a corner plot.

4. "ttv_decay_weighted_samples.txt": the weighted posterior samples.

5. "ttv_decay_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_ttv_fit(free_params, model='constant', file_suffix='', make_plot=True, sigma_clip=False, clip_model='linear')

Run a model fit of the observed transit and/or eclipse mid-times.

This method executes a model fit of the observed transit and/or eclipse mid-times using one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must part of the physical model.

• model (str, optional) – The timing model, must be "constant", "decay", or "precession". Default is "constant".

• file_suffix (str, optional) – A string appended to the end of the output file names.

• make_plot (bool, optional) – If True, a TTV (O-C) plot is generated. Default is True.

• sigma_clip (bool, optional) – Option to execute a sigma-clipping routine on the transit mid-times. Default is False.

• clip_model (str, optional) – Specifies the model to be used in the sigma-clipping routine. The options are "linear" or "quadratic", with the default being "linear".

Returns:

res (dict) – A dictionary containing the model fit results and settings.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_ttv_precession(free_params, suffix, plot, clip, clip_method, save)

Run a fit of the apsidal precession timing model.

This method executes an apsidal precession model fit of the observed transit and/or eclipse mid-times using one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["t0", "P0", "e0", "w0", "wdE", "ecosw", "esinw", "sq_ecosw", sq_esinw"].

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, a TTV (O-C) plot is generated.

• save (bool) – If False, the output files are not saved.

• clip (bool) – If True, a sigma-clipping routine is run on the transit mid-times before model fitting.

• clip_method (str) – Specifies the model to be used in the sigma-clipping routine, either "linear" or "quadratic".

Returns:

res (dict) – A dictionary containing the model fit results and settings.

Note

The following output files are generated:

1. "ttv_precession_summary.txt": a quick visual summary of the results

2. "ttv_precession_results.json": the entire model fitting results dictionary.

3. "ttv_precession_corner.png": a corner plot.

4. "ttv_precession_weighted_samples.txt": the weighted posterior samples.

5. "ttv_precession_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

ttv_loglike_constant(theta)

Calculates the log-likelihood for the constant-period timing model.

This function returns the log-likelihood for the constant-period timing model using the ttv_constant() method.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The log-likelihood value.

ttv_loglike_decay(theta)

Calculates the log-likelihood for the orbital decay timing model.

This function returns the log-likelihood for the orbital decay timing model using the ttv_decay() method.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The log-likelihood value.

ttv_loglike_precession(theta)

Calculates the log-likelihood for the apsidal precession timing model.

This function returns the log-likelihood for the apsidal precession timing model using the ttv_precession() method.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The log-likelihood value.

---

class orbdot.radial_velocity.RadialVelocity(rv_settings)

This class utilizes the capabilities of the NestedSampling class to facilitate model fitting of radial velocity data.

__init__(rv_settings)

Initializes the RadialVelocity class.

Parameters:

rv_settings (dict) – A dictionary specifying directories and settings for the nested sampling analysis.

run_rv_constant(free_params, suffix, plot)

Run a fit of the constant-period radial velocity model.

This method executes a constant-period model fit of the observed radial velocities using one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["t0", "P0", "e0", "w0", "ecosw", "esinw", "sq_ecosw", sq_esinw", "K", "v0", "jit", "dvdt", "ddvdt"]. If the data are from multiple sources, the instrument-dependent parameters ("v0" and "jit") are split into multiple variables before the fit is performed.

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, an RV plot is generated.

Returns:

res (dict) – A dictionary containing the model fit results and settings.

Note

The following output files are generated:

1. "rv_constant_summary.txt": a quick visual summary of the results

2. "rv_constant_results.json": the entire model fitting results dictionary.

3. "rv_constant_corner.png": a corner plot.

4. "rv_constant_weighted_samples.txt": the weighted posterior samples.

5. "rv_constant_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_rv_decay(free_params, suffix, plot)

Run a fit of the orbital decay radial velocity model.

This method executes an orbital decay model fit of the observed radial velocities using one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["t0", "P0", "PdE", "e0", "w0", "ecosw", "esinw", "sq_ecosw", sq_esinw", "K", "v0", "jit", "dvdt", "ddvdt"]. If the data are from multiple sources, the instrument-dependent parameters ("v0" and "jit") are split into multiple variables before the fit is performed.

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, an RV plot is generated.

Returns:

res (dict) – A dictionary containing the model fit results and settings.

Note

The following output files are generated:

1. "rv_decay_summary.txt": a quick visual summary of the results

2. "rv_decay_results.json": the entire model fitting results dictionary.

3. "rv_decay_corner.png": a corner plot.

4. "rv_decay_weighted_samples.txt": the weighted posterior samples.

5. "rv_decay_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_rv_fit(free_params, model='constant', file_suffix='', make_plot=True)

Run a model fit of the radial velocity measurements.

This method executes a model fit of the observed radial velocities using one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be part of the physical model. If the data are from multiple sources, the instrument-dependent parameters ("v0" and "jit") are split into multiple variables before the fit is performed.

• model (str, optional) – The radial velocity model, must be "constant", "decay", or "precession". Default is "constant".

• file_suffix (str, optional) – A string appended to the end of the output file names.

• make_plot (bool, optional) – If True, an RV plot is generated. Default is True.

Returns:

res (dict) – A dictionary containing the model fit results and settings.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_rv_precession(free_params, suffix, plot)

Run a fit of the apsidal precession radial velocity model.

This method executes an apsidal precession model fit of the observed radial velocities using one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["t0", "P0", "e0", "w0", "ecosw", "esinw", "sq_ecosw", sq_esinw", "wdE", "K", "v0", "jit", "dvdt", "ddvdt"]. If the data are from multiple sources, the instrument-dependent parameters ("v0" and "jit") are split into multiple variables before the fit is performed.

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, an RV plot is generated.

Returns:

res (dict) – A dictionary containing the model fit results and settings.

Note

The following output files are generated:

1. "rv_precession_summary.txt": a quick visual summary of the results

2. "rv_precession_results.json": the entire model fitting results dictionary.

3. "rv_precession_corner.png": a corner plot.

4. "rv_precession_weighted_samples.txt": the weighted posterior samples.

5. "rv_precession_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

rv_loglike_constant(theta)

Calculates the log-likelihood for the constant-period radial velocity model.

This function returns the log-likelihood for the constant-period radial velocity model using the rv_constant() method.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The log-likelihood value.

rv_loglike_decay(theta)

Calculates the log-likelihood for the orbital decay radial velocity model.

This function returns the log-likelihood for the orbital decay radial velocity model using the rv_decay() method.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The log-likelihood value.

rv_loglike_precession(theta)

Calculates the log-likelihood for the apsidal precession radial velocity model.

This function returns the log-likelihood for the apsidal precession radial velocity model using the rv_precession() method.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The log-likelihood value.

save_rv_residuals(fit_results, outfile)

Saves the best-fit radial velocity model residuals.

This method writes to a file the difference between the observed radial velocities and the best-fit model, ie. the “residuals”.

Parameters:

• fit_results (dict) – Dictionary containing the results of the radial velocity fit.

• outfile (str) – Output file path name.

Returns:

None

---

class orbdot.transit_duration.TransitDuration(tdv_settings, system_info)

This class utilizes the capabilities of the NestedSampling class to facilitate model fitting of transit durations.

__init__(tdv_settings, system_info)

Initializes the TransitDuration class.

Parameters:

tdv_settings (dict) – A dictionary specifying directories and settings for the nested sampling analysis.

run_tdv_constant(free_params, suffix, plot)

Run a fit of the constant-period transit duration model.

This method executes a constant-period model fit of the observed transit durations using one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["P0", "e0", "w0", "ecosw", "esinw", "sq_ecosw", sq_esinw", "i0"].

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, a TDV plot is generated.

Returns:

res (dict) – A dictionary containing the model fit results and settings.

Note

The following output files are generated:

1. "tdv_constant_summary.txt": a quick visual summary of the results

2. "tdv_constant_results.json": the entire model fitting results dictionary.

3. "tdv_constant_corner.png": a corner plot.

4. "tdv_constant_weighted_samples.txt": the weighted posterior samples.

5. "tdv_constant_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_tdv_decay(free_params, suffix, plot)

Run a fit of the orbital decay transit duration model.

This method executes an orbital decay model fit of the observed transit durations using one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["P0", "e0", "w0", "ecosw", "esinw", "sq_ecosw", sq_esinw", "i0", "PdE"].

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, a TDV plot is generated.

Returns:

res (dict) – A dictionary containing the model fit results and settings.

Note

The following output files are generated:

1. "tdv_decay_summary.txt": a quick visual summary of the results

2. "tdv_decay_results.json": the entire model fitting results dictionary.

3. "tdv_decay_corner.png": a corner plot.

4. "tdv_decay_weighted_samples.txt": the weighted posterior samples.

5. "tdv_decay_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_tdv_fit(free_params, model='constant', file_suffix='', make_plot=True)

Run a model fit of the observed transit durations.

This method executes a model fit of the observed transit durations using one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be part of the physical model.

• model (str, optional) – The transit duration model, must be "constant", "decay", or "precession". Default is "constant".

• file_suffix (str, optional) – A string appended to the end of the output file names.

• make_plot (bool, optional) – If True, a TDV plot is generated. Default is True.

Returns:

res (dict) – A dictionary containing the model fit results and settings.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

run_tdv_precession(free_params, suffix, plot)

Run a fit of the apsidal precession transit duration model.

This method executes an apsidal precession model fit of the observed transit durations using one of two nested sampling packages, Nestle [1]_ or PyMultiNest [2]_.

Parameters:

• free_params (list or tuple) – The list of free parameters for the model fit, in any order. The parameter names are formatted as strings and must be in the set: ["P0", "e0", "w0", "ecosw", "esinw", "sq_ecosw", sq_esinw", "i0", "wdE"].

• suffix (str) – A string appended to the end of the output file names.

• plot (bool) – If True, a TDV plot is generated.

Returns:

res (dict) – A dictionary containing the model fit results and settings.

Note

The following output files are generated:

1. "tdv_precession_summary.txt": a quick visual summary of the results

2. "tdv_precession_results.json": the entire model fitting results dictionary.

3. "tdv_precession_corner.png": a corner plot.

4. "tdv_precession_weighted_samples.txt": the weighted posterior samples.

5. "tdv_precession_random_samples.json": a random set of 300 posterior samples.

References

[1]

Nestle by Kyle Barbary. http://kbarbary.github.io/nestle

[2]

PyMultiNest by Johannes Buchner. http://johannesbuchner.github.io/PyMultiNest

tdv_loglike_constant(theta)

Calculates the log-likelihood for the constant-period transit duration model.

This function returns the log-likelihood for the constant-period transit duration model using the tdv_constant() method.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The log-likelihood value.

tdv_loglike_decay(theta)

Calculates the log-likelihood for the orbital decay transit duration model.

This function returns the log-likelihood for the orbital decay transit duration model using the tdv_decay() method.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The log-likelihood value.

tdv_loglike_precession(theta)

Calculates the log-likelihood for the apsidal precession transit duration model.

This function returns the log-likelihood for the apsidal precession transit duration model using the tdv_precession() method.

Parameters:

theta (array_like) – An array containing parameter values, passed from the sampling algorithm.

Returns:

float – The log-likelihood value.

---

Models

Theory

This module defines analytical models for investigating the source of long-term variations of exoplanet orbits. These include equations for assessing the effects of tidal energy dissipation, apsidal precession, systemic proper motion, and more.

orbdot.models.theory.companion_doppler_pdot_from_rv_trend(P, dvdt)

Calculate the apparent variation of an orbital period due to a line-of-sight acceleration.

This method returns the time derivative of the observed orbital period of a transiting exoplanet due to the Doppler effect induced by an acceleration along the line-of-sight (i.e. a linear RV trend). It uses equation (6) from Bouma et al. (2020) [1]_, which is derived in convenient units.

Parameters:

• P (float) – The orbital period in days.

• dvdt (float) – The first-order radial acceleration term in m/s/day.

Returns:

float – The apparent time derivative of the orbital period in milliseconds per year.

Notes

When radial velocity data of a transiting exoplanet’s host star show a linear trend, and assuming this trend is a real effect, the Doppler effect can cause the observed period between transits to vary. Equation (6) from Bouma et al. (2020) [1]_ expresses the expected period derivative in convenient units:

\[\dot{P}_{\mathrm{RV}} = 105.3 \mathrm{\,ms\,yr^{-1}}\left(\frac{P}{\mathrm{ day}}\right)\left(\frac{\dot{\gamma}}{\mathrm{m\,s^{-1}\,day^{-1}}}\right)\]

where \(P\) is the orbital period of the transiting planet, \(\dot{P}_{\mathrm{ RV}}\) is the time derivative of the period, and \(\dot{\gamma}\) is the first-order RV acceleration term.

References

[1]

Bouma et al. [2020]. https://doi.org/10.3847/2041-8213/ab8563.

orbdot.models.theory.companion_doppler_rv_trend_from_pdot(P, dPdt)

Calculate the line-of-sight acceleration given the apparent variation of an orbital period.

This method returns the slope of the linear radial velocity trend (acceleration) that corresponds to a given time derivative of the observed orbital period, assuming that the latter is a due to the Doppler effect. It uses equation (6) from Bouma et al. (2020) [1]_, which is derived in convenient units.

Parameters:

• P (float) – The orbital period in days.

• dPdt (float) – The apparent time derivative of the orbital period in milliseconds per year (ms/yr).

Returns:

float – The first-order radial acceleration term in m/s/day.

Notes

For a transiting exoplanet, if a variation in the observed period between transits is measured, it may be due to the Doppler effect induced by a line-of-sight acceleration. Assuming this trend is a real effect, Equation (6) from Bouma et al. (2020) [1]_ can be used to determine the corresponding slope of the linear radial velocity trend from the period derivative:

\[\dot{\gamma} = \frac{\dot{P}_{\mathrm{RV}}}{105.3 \mathrm{\,ms\,yr^{-1}}} \left(\frac{\mathrm{day}}{P}\right)\]

where \(P\) is the orbital period of the transiting planet, \(\dot{P}_{\mathrm{ RV}}\) is the time derivative of the period, and \(\dot{\gamma}\) is the first-order RV acceleration term.

References

[1]

Bouma et al. [2020]. https://doi.org/10.3847/2041-8213/ab8563.

orbdot.models.theory.companion_from_quadratic_rv(P_min, t_pivot, dvdt, ddvdt, M_s)

Constrain properties of the orbit of an outer companion given a quadratic RV trend.

Given the first and second-order acceleration terms that parameterize a long-term, quadratic radial velocity trend, \(\dot{\gamma}\) and \(\ddot{\gamma}\), this method follows the derivations of Kipping et al. (2011) [1]_ (equations 1, 3, and 4) to constrain properties of a possible outer companion.

This method requires an estimate of minimum possible orbital period of the companion, which is constrained by the timespan of the radial velocity measurements. This is determined by the rv_trend_quadratic() method from the Analyzer class.

Parameters:

• P_min (float) – The minimum orbital period of the companion in days.

• t_pivot (float) – The “pivot” point, defined as the reference mid-time of the transiting planet.

• dvdt (float) – First-order acceleration term in m/s/day.

• ddvdt (float) – Second-order acceleration term in m/s/day^2.

• M_s (float) – The mass of the host star in solar masses.

Returns:

tuple – A tuple with the following elements: 1) the minimum companion period in days, 2) the corresponding lower limit on the RV semi-amplitude, 3) the companion mass in Earth masses, and 4) the time at which the companion RV signal is at a minimum.

References

[1]

Kipping et al. [2011]. https://doi.org/10.1088/0004-6256/142/3/95

orbdot.models.theory.companion_mass_from_precession(P, P2, dwdE, M_s)

Calculate the mass of a nonresonant companion given a precession rate and orbital period.

Given the orbital period of a nonresonant companion planet in the system, this method returns the necessary mass for it to be responsible for a measured apsidal precession rate of the studied orbit. It uses equation (8) from Heyl and Gladman (2007) [1]_, and is appropriate for a companion planet that is either interior or exterior to the observed (i.e., transiting) planet. However, it assumes that the companion mass is much less than that of the host star.

Parameters:

• P (float) – The orbital period of the transiting planet in days.

• P2 (float) – The orbital period of the outer planet in days.

• dwdE (float) – The rate of apsidal precession in radians per epoch.

• M_s (float) – The mass of the host star in solar masses.

Returns:

float – The mass of the outer companion in earth masses.

Notes

This method applies to the scenario in which a nonresonant companion in the system is driving apsidal precession of a transiting planet. It assumes that the companion object is on a circular, coplanar, and nonresonant orbit, and that its mass is far less than the host star.

In the formulation of Heyl and Gladman (2007) [1]_, the mass of the companion planet can be determined from the measured precession rate by:

\[m_2 = \dot{\omega} \frac{M_\star (\alpha+1)(\alpha-1)^2}{\alpha \left[ \left(\alpha^2+1\right) \mathcal{E}\left(\frac{2 \alpha^{1 / 2}}{\alpha+1}\right) - (\alpha-1)^2 \mathcal{K}\left(\frac{2 \alpha^{1/2}}{\alpha+1}\right)\right]}\]

where \(\dot{\omega}\) is the observed precession rate, \(\alpha = a/a_2\) is the semi major axis ratio, and \(\mathcal{K}\) and \(\mathcal{E}\) are the complete elliptic integrals of the first and second kind, respectively.

Important

The scipy implementation of the complete elliptic integrals \(\mathcal{K}\) and \(\mathcal{E}\) expect an integrand of the form \((1 - m \sin^2 t)^{-1/2}\). However, the derivation of the equation above expects an integrand of the form \((1 - m^2 \sin^2 t)^{-1/2}\) [1]_, which can be shown by expanding \(\mathcal{K}\) and \(\mathcal{E}\). To maintain consistency, the input for the scipy.special.ellipk and scipy.special.ellipe methods must but squared, ie. such that \(\mathcal{K}( 4 \alpha/( \alpha+1)^2)\) and \(\mathcal{E}(4 \alpha/( \alpha+1)^2)\).

References

[1]

Heyl and Gladman [2007]. https://doi.org/10.1111/j.1365-2966.2007.11697.x

orbdot.models.theory.companion_mass_from_rv_trend(tau, dvdt, M_s)

Calculate the minimum mass of an outer companion given the slope of a linear RV trend.

Given a measured linear trend in radial velocity data (acceleration), this method calculates the minimum mass of an unseen outer companion planet that could cause such an acceleration.

While this approach was originally described in equation (1) of Feng et al. (2015) [1]_, this method implements a version given by equation (8) in Bouma et al. (2020) [2]_, which is derived in more convenient units.

Parameters:

• tau (float) – The time span of the observations in years.

• dvdt (float) – The first-order acceleration term in in m/s/day.

• M_s (float) – The mass of the host star in solar masses.

Returns:

float – The companion planet mass in Earth masses.

Notes

Given a measured acceleration in radial velocity measurements and the time span over which there are RV data available, this method calculates the minimum mass of an outer companion that could induce such an acceleration, using equation (8) from Bouma et al. (2020) [2]_:

\[M_{c} \approx 5.99 M_{\mathrm{Jup}}\left(\frac{\tau}{\mathrm{yr}}\right)^{ 4/3}\left|\frac{\dot{\gamma}}{\mathrm{m\,s^{-1}\,day^{-1}}}\right|\left( \frac{M_{\star}}{M_{\odot}}\right)^{2 / 3}\]

where \(\tau\) is the time span of the observations in years, \(\dot{\gamma}\) is the first-order acceleration term, and \(M_\star\) is the mass of the host star.

The above equation assumes that the companion’s orbit has an eccentricity of \(e=0.5\), an argument of pericenter that is exactly \(90\) degrees, and a period that is \(1.25\times\) the time span of the observations. In this scenario, the observed linear trend manifests as a segment of the sawtooth-like RV curve of the companion planet (see Figure 1 of [1]_), for which the semi-amplitude can be approximated as half of the baseline multiplied by the acceleration, i.e. \(K_c = 0.5 \, \tau \, \dot{\gamma}\).

References

[1]

Feng et al. [2015]. https://doi.org/10.1088/0004-637X/800/1/22

[2]

Bouma et al. [2020]. https://doi.org/10.3847/2041-8213/ab8563

orbdot.models.theory.companion_precession(P, M2, P2, M_s)

Calculates the rate of apsidal precession driven by a nonresonant planetary companion.

This method returns the expected rate of apsidal precession due to a nonresonant companion planet, calculated using equation (8) from Heyl and Gladman (2007) [1]_.

It is valid for a companion planet that is either interior or exterior to the observed (i.e. transiting) planet, but it assumes that the companion mass is much less than that of the host star.

Parameters:

• P (float) – The orbital period of the transiting planet in days.

• M2 (float) – The mass of the outer planet in Earth masses.

• P2 (float) – The orbital period of the outer planet in days.

• M_s (float) – The mass of the host star in solar masses.

Returns:

float – The precession rate in radians per epoch.

Notes

This method applies to the scenario in which a nonresonant companion in the system is driving apsidal precession of a transiting planet. It assumes that the companion object is on a circular, coplanar, and nonresonant orbit, and that its mass is far less than the host star.

In the formulation of Heyl and Gladman (2007) [1]_, the induced precession rate in radians per orbit:

\[\dot{\omega} = \frac{m_2}{M_\star} \frac{\alpha}{(\alpha+1)(\alpha-1)^2}\left[ \left(\alpha^2+1\right) \mathcal{E}\left(\frac{2 \alpha^{1 / 2}}{\alpha+1}\right) -(\alpha-1)^2 \mathcal{K}\left(\frac{2 \alpha^{1/2}}{\alpha+1}\right)\right]\]

where \(m_2\) is the mass of the perturbing planet, \(\alpha = a/a_2\) is the semi major axis ratio, and \(\mathcal{K}\) and \(\mathcal{E}\) are the complete elliptic integrals of the first and second kind, respectively.

Important

The scipy implementation of the complete elliptic integrals \(\mathcal{K}\) and \(\mathcal{E}\) expect an integrand of the form \((1 - m \sin^2 t)^{-1/2}\). However, the derivation of the equation above expects an integrand of the form \((1 - m^2 \sin^2 t)^{-1/2}\) [1]_, which can be shown by expanding \(\mathcal{K}\) and \(\mathcal{E}\). To maintain consistency, the input for the scipy.special.ellipk and scipy.special.ellipe methods must but squared, ie. such that \(\mathcal{K}( 4 \alpha/( \alpha+1)^2)\) and \(\mathcal{E}(4 \alpha/( \alpha+1)^2)\).

References

[1]

Heyl and Gladman [2007]. https://doi.org/10.1111/j.1365-2966.2007.11697.x

orbdot.models.theory.companion_rv_trend_from_mass(tau, M_c, M_s)

Calculate the slope of a linear radial velocity trend given an estimate of a companion mass.

This method computes the expected slope of a linear radial velocity trend (acceleration) induced by an unseen outer companion planet with an orbital period that is far longer than the observational baseline. It requires an estimate of the companion mass.

While this approach was originally described in equation (1) of Feng et al. (2015) [1]_, this method implements a version given by equation (8) in Bouma et al. (2020) [ 2]_, which is derived in more convenient units.

Parameters:

• tau (float) – The time span of the observations in years.

• M_c (float) – The companion planet mass in earth masses.

• M_s (float) – The mass of the host star in solar masses.

Returns:

float – The first-order acceleration term in in m/s/day.

Notes

Given a companion planet’s mass, the time span over which radial velocity measurements have been taken, and the mass of the host star, this method calculates the slope of the expected linear trend in the radial velocity data using equation (8) from Bouma et al. (2020) [2]_:

\[\dot{\gamma} \approx \frac{M_c}{5.99 M_{\mathrm{Jup}}} \left(\frac{\mathrm{yr}}{ \tau}\right)^{4/3} \left(\frac{M_{\odot}}{M_{\star}}\right)^{2 / 3} \mathrm{m\, s^{-1}\,day^{-1}}\]

where \(\tau\) is the time span of the observations in years, \(M_c\) is the mass of the companion planet, and \(M_\star\) is the mass of the host star.

The above equation assumes that the companion’s orbit has an eccentricity of \(e=0.5\), an argument of pericenter that is exactly \(90\) degrees, and a period that is \(1.25\times\) the time span of the observations. In this scenario, the observed linear trend manifests as a segment of the sawtooth-like RV curve of the companion planet (see Figure 1 of [1]_), for which the semi-amplitude can be approximated as half of the baseline multiplied by the acceleration, i.e. \(K_c = 0.5 \, \tau \, \dot{\gamma}\).

References

[1]

Feng et al. [2015]. https://doi.org/10.1088/0004-637X/800/1/22

[2]

Bouma et al. [2020]. https://doi.org/10.3847/2041-8213/ab8563

orbdot.models.theory.decay_angular_momentum_loss(P, dPdE, M_s, M_p)

Calculate the rate of angular momentum loss due to tidal forces causing orbital decay.

This function uses Equation (10) from Yee et al. (2020) [1]_ to compute the rate at which the angular momentum of a planetary orbit decreases as it decays due to tidal forces.

Parameters:

• P (float) – Orbital period in days.

• dPdE (float) – Orbital decay rate in days per epoch.

• M_s (float) – Host star mass in solar masses.

• M_p (float) – Planet mass in earth masses.

Returns:

float – The rate of orbital angular momentum loss in kg m^2 / s^2.

Notes

As a planet’s orbit decays, the orbital energy and angular momentum will decrease over time. Equation (10) of Yee et al. (2020) [1]_ defines the rate of orbital angular momentum loss as:

\[\frac{dL}{dt} = \frac{M_{\mathrm{p}}}{3(2\pi)^{1/3}} \left(\frac{G M_{\star}}{ P}\right)^{2/3} \frac{dP}{dt}\]

where \(G\) is the gravitational constant, \(M_{\star}\) is the host star mass, \(M_{\mathrm{p}}\) is the planet mass, \(P\) is the orbital period, and \(\frac{dP}{dt}\) is the orbital decay rate.

References

[1]

Yee et al. [2020]. https://doi.org/10.3847/2041-8213/ab5c16

orbdot.models.theory.decay_empirical_quality_factor(P_orb, P_rot_s)

Calculate the modified stellar tidal quality factor from an empirical law.

This method calculates the tidal forcing period of a star-planet system and uses it to estimate the star’s modified tidal quality factor with an empirical law derived by Penev et al. (2018) [1]_.

Parameters:

• P_orb (float) – Orbital period in days.

• P_rot_s (float) – Period of the host star’s rotation in days.

Returns:

tuple – The star’s modified tidal quality factor and the tidal forcing period of the system in days.

Notes

The “modified” tidal quality factor is a typical parameterization of the efficiency of tidal energy dissipation within a body. It is defined as:

\[Q_\star^{'} = \frac{3}{2} \frac{Q_\star}{k_{2,\star}}\]

where \(Q_\star\) is the tidal quality factor and \(k_{2,\star}\) is the second-order potential Love number. The theoretical upper limit is \(k_{2,\star}=3/2\), that of a uniform density sphere, in which case \(Q_\star^{'} = Q_\star\).

This method estimates a value for \(Q_\star^{'}\) using an empirical law derived by Penev et al. (2018) [1]_:

\[Q^{'}_\star = \max\left(\frac{10^6}{(P_{\mathrm{tide}}/\mathrm{day})^{3.1}}, \, 10^5\right)\]

where \(P_{\mathrm{tide}}\) is the tidal forcing period of the system, defined as:

\[P_{\mathrm{tide}} = \frac{1}{2\left(P_{\mathrm{orb}}^{-1} - P_{\mathrm{rot}}^{ -1}\right)}\]

where \(P_{\mathrm{orb}}\) is the planet’s orbital period and \(P_{\mathrm{rot}}\) is the rotational period of the host star.

References

[1]

Penev et al. [2018]. https://doi.org/10.3847/1538-3881/aaaf71

orbdot.models.theory.decay_energy_loss(P, dPdE, M_s, M_p)

Calculate the rate of orbit energy loss due to tidal forces causing orbital decay.

This function uses Equation (9) from Yee et al. (2020) [1]_ to compute the rate at which the energy of a planetary orbit decreases as it decays due to tidal forces.

Parameters:

• P (float) – Orbital period in days.

• dPdE (float) – Orbital decay rate in days per epoch.

• M_s (float) – Mass of the host star in solar masses.

• M_p (float) – Mass of the planet in earth masses.

Returns:

float – The rate of orbital energy loss in watts.

Notes

As a planet’s orbit decays, the orbital energy and angular momentum will decrease over time. Equation (9) of Yee et al. (2020) [1]_ defines the rate of orbital energy loss as:

\[\frac{dE}{dt} = \frac{(2\pi)^{2/3} M_{\mathrm{p}}}{3} \left(\frac{G M_{\star}}{ P}\right)^{2/3} \frac{1}{P} \frac{dP}{dt}\]

where \(G\) is the gravitational constant, \(M_{\star}\) is the host star mass, \(M_{\mathrm{p}}\) is the planet mass, \(P\) is the orbital period, and \(\frac{dP}{dt}\) is the orbital decay rate.

References

[1]

Yee et al. [2020]. https://doi.org/10.3847/2041-8213/ab5c16

orbdot.models.theory.decay_get_equilibrium_spin_freq(P, e, epsilon_p)

Calculate the equilibrium spin frequency of a planet.

This method returns the equilibrium spin frequency of a planet with nonzero obliquity and/or orbit eccentricity.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• epsilon_p (float) – The planetary obliquity.

Returns:

float – Equilibrium spin frequency in days^{-1}.

Notes

Planetary rotation periods are typically unconstrained, and it is not always appropriate to assume that they are tidally locked (synchronous rotation). When the orbit eccentricity and/or planetary obliquity are nonzero, the rotation state of the planet may be asynchronous.

In this case, a reasonable assumption is that the planet’s spin has tidally evolved toward an equilibrium rate, given by Equation (12) of [1]_:

\[\dot{\theta}_{eq} \equiv n \frac{g(e)}{h(e)}\frac{2\cos{\varepsilon_p}}{(1 + \cos^2{\varepsilon_p})}\]

where \(\varepsilon_p\) is the planetary obliquity (defined as the angle between the planet’s spin and orbital angular momentum vectors), \(n\) is the orbital mean motion, and \(g(e)\) and \(h(e)\) are eccentricity-dependent functions that are defined in decay_get_g_e() and decay_get_h_e().

References

[1]

Leconte et al. [2010]. https://doi.org/10.1051/0004-6361/201014337

orbdot.models.theory.decay_get_f_e(e)

Computes the first eccentricity function for the equilibrium tide methods, f(e).

This method compute one of three eccentricity-dependent functions called by the equilibrium tidal theory methods, given by Equation (4) in [1]_:

\[f(e) = \frac{1+\frac{31}{2} e^2+\frac{255}{8} e^4+\frac{185}{16} e^6+\frac{25}{64} e^8}{(1-e^2)^{15/2}}\]

Parameters:

e (float) – The orbit eccentricity.

Returns:

float – The computed value of f(e).

References

[1]

Leconte et al. [2010]. https://doi.org/10.1051/0004-6361/201014337

orbdot.models.theory.decay_get_g_e(e)

Computes the second eccentricity function for the equilibrium tide methods, g(e).

This method compute one of three eccentricity-dependent functions called by the equilibrium tidal theory methods, given by Equation (3) in [1]_:

\[g(e) = \frac{1+\frac{15}{2} e^2+\frac{45}{8} e^4+\frac{5}{16} e^6}{(1-e^2)^6}\]

Parameters:

e (float) – The orbit eccentricity.

Returns:

float – The computed value of g(e).

References

[1]

Leconte et al. [2010]. https://doi.org/10.1051/0004-6361/201014337

orbdot.models.theory.decay_get_h_e(e)

Computes the third eccentricity function for the equilibrium tide methods, h(e).

This method compute one of three eccentricity-dependent functions called by the equilibrium tidal theory methods, given by Equation (11) in [1]_:

\[h(e) = \frac{1+3 e^2+\frac{3}{8} e^4}{(1-e^2)^{9/2}}\]

Parameters:

e (float) – The orbit eccentricity.

Returns:

float – The calculates value of h(e).

References

[1]

Leconte et al. [2010]. https://doi.org/10.1051/0004-6361/201014337

orbdot.models.theory.decay_planet_pdot_from_quality_factor(P, e, M_s, M_p, R_p, Q_planet, epsilon_p=0.0, P_rot_p=None)

Calculate the rate of orbital decay driven by equilibrium tides in the planet.

This method computes the orbital decay rate for a given “modified” planetary tidal quality factor and, optionally, the planet’s rotation period and obliquity. This expression accounts for both eccentricity-driven and obliquity-driven tidal dissipation, providing a unified treatment of planetary tides.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – Host star mass in solar masses.

• M_p (float) – Planet mass in earth masses.

• R_p (float) – Planet radius in earth radii.

• Q_planet (float) – The planet’s modified annual tidal quality factor.

• epsilon_p (float, optional) – Planetary obliquity in degrees. Default is zero.

• P_rot_p (float, optional) – Planet rotation period in days.

Returns:

float – Orbital decay rate in days per epoch.

Notes

If equilibrium tides are the dominant driver of orbital evolution, the rate of orbital period decay due to planetary tides is given by:

\[\frac{dP}{dt}_{\mathrm{(planet)}} = -\frac{27\pi}{Q_p^{'}} \left(\frac{M_\star}{ M_p}\right) \left(\frac{R_p}{a}\right)^5 \left[f(e) - g(e) \cos{\varepsilon_p} \frac{\dot{\theta}_p}{n}\right]\]

where \(Q_p^{'}\) is the modified annual tidal quality factor of the planet, \(\varepsilon_p\) is the planetary obliquity (defined as the angle between the planet’s spin and orbital angular momentum vectors), \(\dot{\theta}_p\) is the planet’s spin frequency, \(n\) is the orbital mean motion, \(a\) is the semi major axis, \(M_p\) is the planet mass, \(M_\star\) is the stellar mass, and \(R_p\) is the planet’s radius. The eccentricity-dependent functions \(f(e)\) and \(g(e)\) are defined in the decay_get_f_e() and decay_get_g_e() methods.

This equation is derived using Equations (2), (5), and (19) from [1]_ and is valid for any value of \(\dot{\theta}_p\), \(\varepsilon_p\), and \(e\).

If a rotation period is not provided, the planet’s spin frequency is calculated using the tidally-evolved equilibrium rate defined in decay_get_equilibrium_spin_freq(), which is dependent on the eccentricity and planetary obliquity. In this case, the above equation can be expressed as:

\[\frac{dP}{dt}_{\mathrm{(planet)}} = -\frac{27\pi}{Q_p^{'}} \left(\frac{M_\star}{ M_p}\right) \left(\frac{R_p}{a}\right)^5 \times \left[f(e) - \frac{g^2(e)}{h( e)}\frac{2 \cos^2{\varepsilon_p}}{(1 + \cos^2{\varepsilon_p})}\right]\]

We see that if the orbit is circular and aligned, meaning \(e=0\) and \(\varepsilon_p=0\), the planetary contribution to orbital decay vanishes. In this case, the functions \(f(e)\), \(g(e)\), and \(h(e)\) approach unity, and the equilibrium spin rate synchronizes with the mean motion, \(\dot{\theta}_p=n\), leading to no net tidal energy dissipation in the planet.

orbdot.models.theory.decay_planet_quality_factor_from_pdot(P, e, M_s, M_p, R_p, dPdE, epsilon_p=0.0, P_rot_p=None)

Calculate the planet’s modified planetary quality factor given an orbital decay rate.

Given an orbital decay rate and, optionally, the planet’s rotation period and obliquity, this method returns the planet’s modified tidal quality factor under the assumption that the decay is driven by equilibrium tides in the planet. See the decay_planet_pdot_from_quality_factor() docstring for a more detailed description of the relevant equations.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – Host star mass in solar masses.

• M_p (float) – Planet mass in earth masses.

• R_p (float) – Planet radius in earth radii.

• dPdE (float) – Orbital decay rate in days per epoch.

• epsilon_p (float, optional) – Planetary obliquity in degrees. Default is zero.

• P_rot_p (float, optional) – Planet rotation period in days.

Returns:

float – The planet’s modified annual tidal quality factor.

Notes

If a rotation period not provided for the planet, the planet’s spin frequency is calculated with the decay_get_equilibrium_spin_freq() method.

orbdot.models.theory.decay_star_pdot_from_quality_factor(P, e, M_s, M_p, R_s, Q_star, epsilon_s=0.0, P_rot_s=None)

Calculate the rate of orbital decay driven by equilibrium tides in the host star.

This method computes the orbital decay rate for a given “modified” stellar tidal quality factor and, optionally, the host star’s rotation period and/or obliquity.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – Host star mass in solar masses.

• M_p (float) – Planet mass in earth masses.

• R_s (float) – Host star radius in solar radii.

• Q_star (float) – The host star’s modified annual tidal quality factor.

• epsilon_s (float, optional) – Host star obliquity in degrees. Default is zero.

• P_rot_s (float, optional) – Host star rotation period in days.

Returns:

float – Orbital decay rate in days per epoch.

Notes

If equilibrium tides are the dominant driver of orbital evolution, the rate of orbital period decay due to stellar tides is given by:

\[\frac{dP}{dt}_{\mathrm{(star)}} = -\frac{27\pi}{Q_\star^{'}} \left(\frac{M_p}{ M_\star}\right) \left(\frac{R_\star}{a}\right)^5 \times \left[f(e) - g(e) \cos{ \varepsilon_\star} \frac{\dot{\theta}_\star}{n}\right]\]

where \(Q_\star^{'}\) is the modified annual tidal quality factor of the star, \(\varepsilon_\star\) is the stellar obliquity (defined as the angle between the spin and orbital angular momentum vectors), \(\dot{\theta}_\star\) is the stellar spin frequency, \(n\) is the orbital mean motion, \(a\) is the semi major axis, \(M_p\) is the planet mass, \(M_\star\) is the stellar mass, and \(R_\star\) is the stellar radius. The eccentricity-dependent functions \(f(e)\) and \(g(e)\) are defined in the decay_get_f_e() and decay_get_g_e() methods.

This equation is derived using Equations (2), (5), and (19) from [1]_ and is valid for any value of \(\dot{\theta}_\star\), \(\varepsilon_\star\), and \(e\). We see that if \(e=0\) and \(\varepsilon_\star=0\), the expression simplifies to:

\[\frac{dP}{dt}_{\mathrm{(star)}} = -\frac{27\pi}{Q_\star^{'}} \left(\frac{M_p}{ M_\star}\right) \left(\frac{R_\star}{a}\right)^5 \left[1 - \frac{\dot{ \theta}_\star}{n}\right]\]

If the rotation rate of the star is unknown, it is assumed that \(\dot{\theta}_\star \ll n\), in which case the bracketed term approaches unity, leaving:

\[\frac{dP}{dt}_{\mathrm{(star)}} = -\frac{27\pi}{Q_\star^{'}} \left(\frac{M_p}{ M_\star}\right) \left(\frac{R_\star}{a}\right)^5\]

References

[1]

Leconte et al. [2010]. https://doi.org/10.1051/0004-6361/201014337

orbdot.models.theory.decay_star_quality_factor_from_pdot(P, e, M_s, M_p, R_s, dPdE, epsilon_s=0.0, P_rot_s=None)

Calculate the host star’s modified stellar quality factor given an orbital decay rate.

Given an orbital decay rate and, optionally, the host star’s rotation period and obliquity, this method returns the host star’s modified tidal quality factor under the assumption that the decay is driven by equilibrium tides in the host star. See the decay_star_pdot_from_quality_factor() docstring for a more detailed description of the relevant equations.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – Host star mass in solar masses.

• M_p (float) – Planet mass in earth masses.

• R_s (float) – Host star radius in solar radii.

• dPdE (float) – Orbital decay rate in days per epoch.

• epsilon_s (float, optional) – Host star obliquity in degrees. Default is zero.

• P_rot_s (float, optional) – Host star rotation period in days.

Returns:

float – The host star’s modified annual tidal quality factor.

orbdot.models.theory.decay_timescale(P, dPdE)

Calculate the remaining lifetime of a planet on a decaying orbit.

Parameters:

• P (float) – Orbital period in days.

• dPdE (float) – Orbital decay rate in days per epoch.

Returns:

float – The remaining lifetime in millions of years (Myr).

Notes

For a planet on a decaying orbit, its remaining lifetime can be estimated as simply:

\[\tau = \frac{P_0}{|\dot{P}_{\mathrm{decay}}|}\]

where \(\dot{P}_{\mathrm{decay}}\) is the orbital decay rate, and \(P_0\) is the observed orbital period.

orbdot.models.theory.get_pdot_from_wdot(P, e, w, dwdE)

Calculate the apparent time derivative of the orbital period due to apsidal precession.

This method determines the apparent time derivative of the orbital period due to apsidal precession in general, independent of the physical mechanism. It uses Equation (17) from Rafikov (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• w (float) – The argument of pericenter in radians.

• dwdE (float) – The precession rate in radians per epoch.

Returns:

float – The time derivative of the orbital period in milliseconds per year (ms/yr).

Notes

For transiting exoplanets, apsidal precession causes the sidereal period, i.e. the length of time between transits, to vary over time. Because apsidal precession periods are typically far greater than observational baselines, this effect yields an apparent time derivative \(\dot{P}\).

Equation (17) of Rafikov (2009) [1]_ expresses this as:

\[\dot{P} = \frac{4\pi\left(\dot{\omega}\right)^2 e\cos\omega_p}{ n^2}\frac{\left(1-e^2\right)^{3/2}}{(1+e\sin\omega_p)^3}\]

where \(e\) is the orbit eccentricity, \(\omega_p\) is the argument of pericenter, \(\dot{\omega}\) is the apsidal precession rate, and \(n = 2\pi/P\) is the orbit mean motion.

References

[1]

Rafikov [2009]. https://doi.org/10.1088/0004-637X/700/2/965

orbdot.models.theory.get_semi_major_axis_from_period(P, M_s)

Calculate the semi major axis in meters given the orbital period and host star mass.

Parameters:

• P (float) – The orbital period in days.

• M_s (float) – The mass of the host star in solar masses.

Returns:

float – The semi major axis of the orbit in meters.

orbdot.models.theory.get_tdot_from_wdot(P, e, w, i, T, dwdE, M_s, R_s)

Calculate the expected transit duration variation (TDV) due to apsidal precession.

This method determines the expected TDV signal for apsidal precession in general, independent of the physical mechanism. It uses Equation (9) from Rafikov (2009) [1]_, assuming there is no change in the line-of-sight inclination over time.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• w (float) – The argument of pericenter in radians.

• i (float) – The line-of-sight inclination in degrees.

• T (float) – The transit duration in minutes.

• dwdE (float) – The precession rate in radians per epoch.

• M_s (float) – The host star mass in solar masses.

• R_s (float) – The host star radius in solar radii.

Returns:

float – The time derivative of the transit duration in milliseconds per year (ms/yr).

Notes

For transiting exoplanets, apsidal precession causes the duration of the transits \(T\) to vary over time. Assuming that the line-of-sight inclination \(i\) of the orbit is constant, the expected transit duration variation (TDV) is given by [1]_:

\[\dot{T} = -\frac{T}{1 + e\sin{\omega_p}} \left[e \, \dot{\omega}\cos{\omega_p} - g \, \dot{\omega}\cos{i}\frac{e\cos{\omega_p}}{1+e\sin{\omega_p}}\right]\]

where,

\[g = \frac{a}{R_\star}\frac{b}{1 - b^2}\]

and where \(e\) is the orbit eccentricity, \(\omega_p\) is the argument of pericenter, \(\dot{\omega}\) is the apsidal precession rate, \(b\) is the impact parameter, \(a\) is the semi major axis, and \(R_\star\) is the radius of the host star.

References

[1]

Rafikov [2009]. https://doi.org/10.1088/0004-637X/700/2/965

orbdot.models.theory.precession_gr(P, e, M_s)

Calculate the rate of apsidal precession predicted by general relativity.

This method returns the expected apsidal precession rate of the planet’s orbit due to general relativistic (GR) effects, using Equation (12) from Ragozzine and Wolf (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – The host star mass in solar masses.

Returns:

float – The precession rate in radians per epoch.

Notes

This method applies the lowest order of the relativistic contribution to apsidal precession, given in Equation (12) of Ragozzine and Wolf (2009) [1]_ as:

\[\dot{\omega}_{\mathrm{GR}} = \frac{3 G M_{\star} n}{a c^2(1 - e^2)}\]

where \(G\) is the gravitational constant, \(c\) is the speed of light in a vacuum, \(M_{\star}\) is the host star mass, \(a\) is the planet’s semi major axis, \(e\) is the eccentricity, and \(n = 2\pi/P\) is the mean motion.

References

[1]

Ragozzine and Wolf [2009]. https://doi.org/10.1088/0004-637X/698/2/1778

orbdot.models.theory.precession_rotational_planet(P, e, M_s, M_p, R_p, k2_p, P_rot_p)

Calculates the rate of apsidal precession driven by planetary rotation.

This method returns the expected apsidal precession rate of a planet’s orbit due to its rotational bulge, calculated using Equations (10) and (11) from Ragozzine and Wolf (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – The host star mass in solar masses.

• M_p (float) – The planet mass in earth masses.

• R_p (float) – The planet radius in earth radii.

• k2_p (float) – The planet’s Love number.

• P_rot_p (float) – The period of the planet’s rotation in days.

Returns:

float – The precession rate in radians per epoch.

Notes

The rotation of a fluid body raises an oblate distortion (i.e., flattening) that perturbs the gravitational potential, driving apsidal precession. The rate of the apsidal precession induced by a planet’s rotational bulge is given in Equation (10) of Ragozzine and Wolf (2009) [1]_ as:

\[\dot{\omega}_{\mathrm{rot,p}} = \frac{n {R_p}^5 k_{2,p} {\dot{\theta}_p}^2}{2 a^2 G M_p} g_2(e)\]

where \(\dot{\theta}_p\) is the rotational velocity of the planet, \(k_{2,p}\) is its second-order potential Love number, \(G\) is the gravitational constant, \(n = 2\pi/P\) is the orbit mean motion, \(a\) is the semi major axis, and \(M_p\) and \(R_p\) are the planet mass and radius, respectively.

The function \(g_2(e)\) represents an expansion in eccentricity \(e\), defined in Equation (11) of Ragozzine and Wolf (2009) [1]_ as:

\[g_2(e) = (1-e^2)^{-2}\]

References

[1]

Ragozzine and Wolf [2009]. https://doi.org/10.1088/0004-637X/698/2/1778

[2]

Lissauer and de Pater [2019]. https://doi.org/10.1017/9781108304061

orbdot.models.theory.precession_rotational_planet_k2(P, e, M_s, M_p, R_p, P_rot_p, dwdE)

Calculate the Love number given a rate of apsidal precession driven by planetary rotation.

Given an apsidal precession rate for a planetary orbit, this method returns the second-order potential Love number of the planet \(k_{2,p}\) under the assumption that the precession is driven entirely by the rotational bulge of the planet. See the precession_rotational_planet() docstring for a more detailed description of the relevant equations, which are adopted from Ragozzine and Wolf (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – The host star mass in solar masses.

• M_p (float) – The planet mass in earth masses.

• R_p (float) – The planet radius in earth radii.

• P_rot_p (float) – The period of the planet’s rotation in days.

• dwdE (float) – The precession rate in radians per epoch.

Returns:

float – The planet’s potential Love number of the second order.

References

[1]

Ragozzine and Wolf [2009]. https://doi.org/10.1088/0004-637X/698/2/1778

orbdot.models.theory.precession_rotational_star(P, e, M_s, R_s, k2_s, P_rot_s)

Calculate the rate of apsidal precession driven by stellar rotation.

This method returns the expected apsidal precession rate of a planet’s orbit due to the rotational bulge of the host star, calculated using Equations (10) and (11) from Ragozzine and Wolf (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – The host star mass in solar masses.

• R_s (float) – The host star radius in solar radii.

• k2_s (float) – The host star’s Love number.

• P_rot_s (float) – The period of the host star’s rotation in days.

Returns:

float – The precession rate in radians per epoch.

Notes

The rotation of a fluid body raises an oblate distortion (i.e., flattening) that perturbs the gravitational potential, driving apsidal precession. The rate of the apsidal precession induced by a star’s rotational bulge is given in Equation (10) of Ragozzine and Wolf (2009) [1]_ as:

\[\dot{\omega}_{\mathrm{rot,\star}} = \frac{n {R_\star}^5 k_{2,\star} {\dot{ \theta}_\star}^2}{2 a^2 G M_\star} g_2(e)\]

where \(\dot{\theta}_\star\) is the rotational velocity of the star, \(k_{2, \star}\) is its second-order potential Love number, \(G\) is the gravitational constant, \(a\) is the semi major axis of the orbit, \(n = 2\pi/P\) is the mean motion, and \(M_\star\) and \(R_\star\) are the stellar mass and radius, respectively.

The parameter \(g_2(e)\) represents an expansion in eccentricity \(e\), defined in Equation (11) of Ragozzine and Wolf (2009) [1]_ as:

\[g_2(e) = (1-e^2)^{-2}\]

References

[1]

Ragozzine and Wolf [2009]. https://doi.org/10.1088/0004-637X/698/2/1778

[2]

Lissauer and de Pater [2019]. https://doi.org/10.1017/9781108304061

orbdot.models.theory.precession_rotational_star_k2(P, e, M_s, R_s, P_rot_s, dwdE)

Calculate the Love number given a rate of apsidal precession driven by stellar rotation.

Given an apsidal precession rate for a planetary orbit, this method returns the second-order potential Love number of the host star \(k_{2,\star}\) under the assumption that the precession is driven entirely by the rotational bulge of the host star. See the precession_rotational_star() docstring for a more detailed description of the relevant equations, which are adopted from Ragozzine and Wolf (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – The host star mass in solar masses.

• R_s (float) – The host star radius in solar radii.

• P_rot_s (float) – The period of the host star’s rotation in days.

• dwdE (float) – The precession rate in radians per epoch.

Returns:

float – The host star’s potential Love number of the second order.

References

[1]

Ragozzine and Wolf [2009]. https://doi.org/10.1088/0004-637X/698/2/1778

orbdot.models.theory.precession_tidal_planet(P, e, M_s, M_p, R_p, k2_p)

Calculate the rate of apsidal precession driven by the planet’s tidal bulge.

This method returns the expected apsidal precession rate of the planet’s orbit due to its tidal bulge, calculated using Equations (6) and (7) from Ragozzine and Wolf (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – The host star mass in solar masses.

• M_p (float) – The planet mass in earth masses.

• R_p (float) – The planet radius in earth radii.

• k2_p (float) – The planet’s Love number.

Returns:

float – The precession rate in radians per epoch.

Notes

The host star of a massive, close-in planet may raise a significant “tidal bulge” on the planet, an ellipsoidal distortion caused by the varying gravitational force experienced across the extent of the body, driving apsidal precession. The rate of apsidal precession that is induced by the planet’s tidal bulge is given in Equation (6) of Ragozzine and Wolf (2009) [1]_ as:

\[\dot{\omega}_{\mathrm{tide,p}} = \frac{15}{2}k_{2,p} n \left(\frac{R_p}{ a}\right)^5\left(\frac{M_\star}{M_p}\right)f_2(e)\]

where \(k_{2,p}\) is the planet’s second-order potential Love number, \(n = 2\pi/P\) is the orbit mean motion, \(a\) is the semi major axis, \(M_\star\) is the host star mass, and \(M_p\) and \(R_p\) are the planet mass and radius.

The parameter \(f_2(e)\) represents an expansion in eccentricity \(e\), defined in Equation (7) of Ragozzine and Wolf (2009) [1]_ as:

\[f_2(e)= (1-e^2)^{-5} \left(1 + \frac{3}{2}e^2 + \frac{1}{8}e^4 \right).\]

References

[1]

Ragozzine and Wolf [2009]. https://doi.org/10.1088/0004-637X/698/2/1778

orbdot.models.theory.precession_tidal_planet_k2(P, e, M_s, M_p, R_p, dwdE)

Calculate a planet’s Love number given a rate of apsidal precession driven by tides.

Given an apsidal precession rate for a planetary orbit, this method returns the second-order potential Love number of the planet \(k_{2,p}\) under the assumption that the precession is driven entirely by the tidal bulge of the planet. See the precession_tidal_planet() docstring for a more detailed description of the relevant equations, which are adopted from Ragozzine and Wolf (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – The host star mass in solar masses.

• M_p (float) – The planet mass in earth masses.

• R_p (float) – The planet radius in earth radii.

• dwdE (float) – The precession rate in radians per epoch.

Returns:

float – The planet’s potential Love number of the second order.

References

[1]

Ragozzine and Wolf [2009]. https://doi.org/10.1088/0004-637X/698/2/1778

orbdot.models.theory.precession_tidal_star(P, e, M_s, M_p, R_s, k2_s)

Calculate the rate of apsidal precession driven by the star’s tidal bulge.

This method returns the expected apsidal precession rate of the planet’s orbit due to the star’s tidal bulge, calculated using Equations (6) and (7) from Ragozzine and Wolf (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – The host star mass in solar masses.

• M_p (float) – The planet mass in earth masses.

• R_s (float) – The host star radius in solar radii.

• k2_s (float) – The host star’s Love number.

Returns:

float – The precession rate in radians per epoch.

Notes

A massive, close-in planet may raise a significant “tidal bulge” on the host star, an ellipsoidal distortion caused by the varying gravitational force experienced across the extent of the body, driving apsidal precession. The rate of apsidal precession that is induced by the star’s tidal bulge is given in Equation (6) of Ragozzine and Wolf (2009) [1]_ as:

\[\dot{\omega}_{\mathrm{tide,\star}} = \frac{15}{2}k_{2,\star} n \left(\frac{ R_\star}{a}\right)^5\left(\frac{M_p}{M_\star}\right)f_2(e)\]

where \(k_{2,\star}\) is the star’s second-order potential Love number, \(n = 2\pi/P\) is the orbit mean motion, \(a\) is the semi major axis, \(M_p\) is the planet mass, and \(M_\star\) and \(R_\star\) are the stellar mass and radius.

The parameter \(f_2(e)\) represents an expansion in eccentricity \(e\), defined in Equation (7) of Ragozzine and Wolf (2009) [1]_ as:

\[f_2(e)= (1-e^2)^{-5} \left(1 + \frac{3}{2}e^2 + \frac{1}{8}e^4 \right).\]

References

[1]

Ragozzine and Wolf [2009]. https://doi.org/10.1088/0004-637X/698/2/1778

orbdot.models.theory.precession_tidal_star_k2(P, e, M_s, M_p, R_s, dwdE)

Calculate a host star’s Love number given a rate of apsidal precession driven by tides.

Given an apsidal precession rate for a planetary orbit, this method returns the second-order potential Love number of the host star \(k_{2,\star}\) under the assumption that the precession is driven entirely by the tidal bulge of the star. See the precession_tidal_star() docstring for a more detailed description of the relevant equations, which are adopted from Ragozzine and Wolf (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• M_s (float) – The host star mass in solar masses.

• M_p (float) – The planet mass in earth masses.

• R_s (float) – The host star radius in solar radii.

• dwdE (float) – The precession rate in radians per epoch.

Returns:

float – The host star’s potential Love number of the second order.

References

[1]

Ragozzine and Wolf [2009]. https://doi.org/10.1088/0004-637X/698/2/1778

orbdot.models.theory.proper_motion_idot(mu, beta)

Calculate the apparent rate of change of the inclination due to systemic proper motion.

This method returns the rate of the variation of the line-of-sight inclination of an orbit due to systemic proper motion, using Equation (3) from Rafikov (2009) [1]_.

Parameters:

• mu (float) – The proper motion of the system in mas/yr (milliarcseconds per year).

• beta (float) – The angle in radians between the proper motion vector and the angular momentum vector projected onto the sky-plane.

Returns:

float – The apparent time derivative of the line-of-sight inclination in radians per year.

Notes

The proper motion of a star-planet system alters its appearance on the sky-plane, which causes an evolution of the line-of-sight inclination \(\dot{i}_\mu\). Rafikov (2009) [1]_ derives a simple expression for this effect in terms of the magnitude of the proper motion vector, \(\mu = |\vec{\mu}|\):

\[\dot{i}_\mu = -\mu\cos{\beta}\]

where \(\beta\) is defined as the angle in the sky-plane between the proper motion vector \(\vec{\mu}\) and the projection of the orbital angular momentum vector \(\vec{L}\). Thus, the magnitude of the inclination variation \(|\dot{i}_\mu|\) is maximized when \(\beta = 0^\circ, 180^\circ\).

References

[1]

Rafikov [2009]. https://doi.org/10.1088/0004-637X/700/2/965

orbdot.models.theory.proper_motion_pdot(P, e, w, mu)

Calculate the apparent time derivative of the orbital period due to systemic proper motion.

This method returns the apparent time derivative of the orbital period that is expected as a result of the systemic proper motion-induced (apparent) apsidal precession. It uses Equation (15) from Rafikov (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• w (float) – The argument of pericenter in radians.

• mu (float) – The proper motion of the system in milliarcseconds per year (mas/yr).

Returns:

float – The apparent time derivative of the orbital period in milliseconds per year (ms/yr).

Notes

The systemic proper motion of a transiting star-planet system can have a significant influence on the observed period between transits due to the apparent apsidal precession \(\dot{\omega}_{\mu}\). This manifests as a constant change in the observed period, which Rafikov (2009) [1]_ expresses as:

\[\dot{P}_{\omega,\mu} = -\frac{2\pi}{n^2} \frac{(1-e^2)^{3/2}}{(1+e\sin{ \omega_p})^2} \times \left[\ddot{\omega}_{\mu}-2(\dot{\omega}_{\mu})^2\frac{ e\cos{\omega_p}}{1+e\sin{\omega_p}}\right]\]

where \(e\) is the orbit eccentricity, \(\omega_p\) is the argument of pericenter, \(n\) is the orbit mean motion, \(\dot{\omega}_{\mu}\) is the apparent precession rate, and \(\ddot{\omega}_{\mu}\) is the time derivative of the latter.

Rafikov [1]_ posits that for transiting systems the following approximation is valid:

\[\ddot{\omega}_{\mu} \sim \mu^2 \sim (\dot{\omega}_{\mu})^2\]

which effectively maximizes the first equation.

References

[1]

Rafikov [2009]. https://doi.org/10.1088/0004-637X/700/2/965

orbdot.models.theory.proper_motion_shklovskii(P, mu, D)

Calculate the apparent time derivative of observed period due to the Shklovskii effect.

This method returns the apparent rate of change of the period between transits due to the Shklovskii effect, using Equation (21) from Rafikov (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• mu (float) – The proper motion of the system in mas/yr (milliarcseconds per year).

• D (float) – The distance to the system in parsecs.

Returns:

float – The apparent time derivative of the orbital period in milliseconds per year (ms/yr).

Notes

As a star-planet system moves through space, the radial component of its velocity, \(v_r\), varies in time as a result of the nonzero transverse velocity component:

\[v_t = \mu D\]

where \(D\) is the distance to the system in parsecs and \(\mu = |\vec{\mu}|\) is the magnitude of the proper motion vector.

Consequently, the radial velocity component is time-dependent, which leads to a secular variation in the observed period between transits, denoted \(\dot{P}_{ \mathrm{Shk}}\). This is known as the Shklovskii Effect [1]_, expressed in the following convenient form by Rafikov (2009) [2]_:

\[\dot{P}_{\mathrm{Shk}} = 20\left(\frac{\mu}{100\,{\mathrm{mas\, yr}^{-1}}}\right)^2\frac{D}{100}\frac{P}{3\,\mathrm{day}}\,\mu {\mathrm{s\, yr}^{-1}}\]

References

[1]

Rafikov [2009]. https://doi.org/10.1088/0004-637X/700/2/965

[2]

Shklovskii [1970]. https://ui.adsabs.harvard.edu/abs/1970SvA….13..562S

orbdot.models.theory.proper_motion_tdot(P, e, w, i, T, wdot_pm, idot_pm, M_s, R_s)

Calculate the time derivative of the transit duration due to systemic proper motion.

This method returns the time derivative of the transit duration due to systemic proper motion using Equation (9) from Rafikov (2009) [1]_.

Parameters:

• P (float) – The orbital period in days.

• e (float) – The eccentricity of the orbit.

• w (float) – The argument of pericenter in radians.

• i (float) – The line-of-sight inclination of the orbit in degrees.

• T (float) – The transit duration in minutes.

• wdot_pm (float) – The apparent apsidal precession rate in radians per year.

• idot_pm (float) – The apparent time derivative of the line-of-sight inclination in radians per year.

• M_s (float) – The host star mass in solar masses.

• R_s (float) – The host star radius in solar radii.

Returns:

float – A constant time derivative of the transit duration in milliseconds per year.

Notes

The systemic proper motion of a transiting star-planet system can have a significant influence on the observed transit duration, \(T\), due to the proper motion-induced variation of the line-of-sight inclination \(\dot{i}_\mu\) and the apparent apsidal precession \(\dot{\omega}_{\mu}\). The transit duration variation (TDV) is governed by the following expression [1]_:

\[\dot{T}_{\mu} = - \frac{T}{1 + e\sin{\omega_p}} \times \left[e\dot{\omega}_{ \mu}\cos{\omega_p} - g \left(\dot{i}_{\mu}\sin{i} + \dot{\omega}_{\mu}\cos{ i}\frac{e\cos{\omega_p}}{1+e\sin{\omega_p}}\right)\right]\]

where,

\[g = \frac{a}{R_\star}\frac{b}{1 - b^2}\]

and where \(e\) is the orbit eccentricity, \(\omega_p\) is the argument of pericenter, \(b\) is the impact parameter, \(a\) is the semi major axis, and \(R_\star\) is the host star radius.

References

[1]

Rafikov [2009]. https://doi.org/10.1088/0004-637X/700/2/965

orbdot.models.theory.proper_motion_wdot(mu, i, beta)

Calculate the rate of the apparent apsidal precession due to systemic proper motion.

This method returns the rate of the apparent apsidal precession induced by systemic proper motion, using Equation (4) from Rafikov (2009) [1]_.

Parameters:

• mu (float) – The proper motion of the system in mas/yr (milliarcseconds per year).

• i (float) – The line-of-sight inclination of the orbit in degrees.

• beta (float) – The angle in radians between the proper motion vector and the angular momentum vector projected onto the sky-plane.

Returns:

float – The apparent apsidal precession rate in radians per year.

Notes

The proper motion of a star-planet system alters its appearance on the sky-plane, which causes an apparent precession of the argument of pericenter \(\dot{ \omega}_\mu\). Rafikov (2009) [1]_ derives a simple expression for this effect in terms of the magnitude of the proper motion vector, \(\mu = |\vec{\mu}|\):

\[\dot{\omega}_\mu = -\mu\frac{\sin{\beta}}{\sin{i}}\]

where \(\beta\) is defined as the angle in the sky-plane between the proper motion vector \(\vec{\mu}\) and the projection of the orbital angular momentum vector \(\vec{L}\). The magnitude of the apparent precession rate \(|\dot{\omega}_\mu|\) is maximized when \(\beta = 90^\circ, 270^\circ\).

References

[1]

Rafikov [2009]. https://doi.org/10.1088/0004-637X/700/2/965

orbdot.models.theory.resolved_binary_mass_from_rv_trend(theta, D, dvdt)

Calculate the minimum mass of a resolved secondary star given a measured radial acceleration.

This method applies to the case where an acceleration has been observed in radial velocity observations of a star, and there is a known secondary object for which an angular separation has been measured.

It uses Equation (6) from Torres (1999) [1]_ to estimate a lower limit for the mass of the secondary object, under the assumption that the radial velocity trend is due entirely to its gravitational influence. In the derivation of this equation, the author makes no assumptions about the mass or brightness of the secondary [1]_.

Parameters:

• theta (float) – The angular separation of the binary in arcseconds.

• D (float) – The distance to the system in parsecs.

• dvdt (float) – The first-order RV acceleration term in m/s/day.

Returns:

float – The mass of the secondary star in solar masses.

Notes

Given the angular separation of a known companion object, the distance to the system, and a measured acceleration in radial velocity observations, this method uses Equation (6) from Torres (1999) [1]_ to estimate the minimum mass of the companion:

\[M_B = 5.341 \times 10^{-6}(D \rho)^2\left|\dot{\gamma}\right| \Phi\]

where \(\rho\) is the angular separation of the binary in arcseconds, \(D\) is the distance to the system in parsecs, \(\dot{\gamma}\) is the first-order radial acceleration term in m/s/day, and \(M_B\) is in units of solar masses.

The parameter \(\Phi\) is a function of the eccentricity, longitude of pericenter, and inclination of the companion’s orbit. Assuming that these parameters are unconstrained, this method uses the minimum value \(\Phi = \frac{3\sqrt{3}}{2}\) to determine the minimum companion mass.

References

[1]

Torres [1999]. https://doi.org/10.1086/316313.

orbdot.models.theory.resolved_binary_rv_trend_from_mass(theta, D, M_B)

Calculate the slope of a linear RV trend given properties of a resolved stellar companion.

This method returns the minimum possible acceleration, induced by a bound stellar companion, that may be observed in radial velocity observations of the primary. The secondary star must have been resolved through imaging or astrometric measurements, such that the angular separation is known. It uses Equation (6) from Torres (1999) [1]_, which makes no assumptions about the mass or brightness of the secondary [1]_.

Parameters:

• theta (float) – The angular separation of the binary in arcseconds.

• D (float) – The distance to the system in parsecs.

• M_B (float) – The mass of the stellar companion in solar masses.

Returns:

float – The first-order radial acceleration term in m/s/day.

Notes

Given the angular separation of a known companion, the distance to the system, and the mass of the companion, this method uses Equation (6) from Torres et al. (1999) [1]_ to estimate the acceleration that could be measured in radial velocity data:

\[\left|\dot{\gamma}\right| = \frac{M_B}{5.341 \times 10^{-6}(D \rho)^2 \Phi}\]

where \(\rho\) is the angular separation of the binary in arcseconds, \(D\) is the distance to the system in parsecs, and \(M_B\) is the mass of the stellar companion in solar masses. In the derivation of this equation, the author makes no assumptions about the mass or brightness of the secondary [1]_.

The parameter \(\Phi\) is a function of the eccentricity, longitude of pericenter, and inclination of the companion’s orbit. Assuming that these parameters are unconstrained, this method uses the minimum value \(\Phi = \frac{3\sqrt{3}}{2}\) to determine the minimum companion mass.

References

[1]

Torres [1999]. https://doi.org/10.1086/316313.

---

Transit and Eclipse Timing Models

This module defines functions for modelling exoplanet transit and eclipse mid-times.

orbdot.models.ttv_models.ttv_constant(t0, P0, e0, w0, E, primary=True)

Constant-period model for transit and eclipse mid-times.

This method calculates the expected transit or eclipse mid-times for a single planet on an unchanging orbit.

Parameters:

• t0 (float) – Reference transit mid-time in \(\mathrm{BJD}_\mathrm{TDB}\).

• P0 (float) – Orbital period in days.

• e0 (float) – Eccentricity of the orbit.

• w0 (float) – Argument of pericenter of the planetary orbit in radians.

• E (array-like) – The epoch(s) at which to calculate the mid-times.

• primary (bool, optional) – If True, returns the transit mid-time. If False, returns the eclipse mid-time. Default is True.

Returns:

float or array-like – The predicted transit or eclipse times in \(\mathrm{BJD}_\mathrm{TDB}\).

Note

If the orbit is eccentric, an offset of \(\frac{2P}{\pi} e\cos{\omega}\) is added to the eclipse mid-times.

orbdot.models.ttv_models.ttv_decay(t0, P0, PdE, e0, w0, E, primary=True)

Orbital decay model for transit and eclipse mid-times.

This method calculates the expected transit or eclipse mid-times for a single planet on an orbit with a constant change in the period. Though the main application of this model is for orbital decay, a positive period derivative is allowed.

Parameters:

• t0 (float) – Reference transit time in \(\mathrm{BJD}_\mathrm{TDB}\).

• P0 (float) – Orbital period at the reference mid-time in days.

• PdE (float) – Rate of change of the orbital period in days per epoch.

• e0 (float) – Eccentricity of the orbit.

• w0 (float) – Argument of pericenter of the planetary orbit in radians.

• E (array-like) – The epoch(s) at which to calculate the mid-times.

• primary (bool, optional) – If True, returns the transit mid-time. If False, returns the eclipse mid-time. Default is True.

Returns:

float or array-like – The predicted transit or eclipse times in \(\mathrm{BJD}_\mathrm{TDB}\).

Note

If the orbit is eccentric, an offset of \(\frac{2P}{\pi} e\cos{\omega}\) is added to the eclipse mid-times.

orbdot.models.ttv_models.ttv_precession(t0, P0, e0, w0, wdE, E, primary=True)

Apsidal precession model for transit and eclipse mid-times.

This method calculates the expected transit or eclipse mid-times for an elliptical orbit undergoing apsidal precession. It uses a numerical approximation for low eccentricities that is adapted from equation (15) in Giminez and Bastero (1995) [1]_ by Patra et al. (2017) [2]_.

Parameters:

• t0 (float) – Reference transit time in \(\mathrm{BJD}_\mathrm{TDB}\).

• P0 (float) – The sidereal period in days.

• e0 (float) – Eccentricity of the orbit.

• w0 (float) – Argument of pericenter of the planetary orbit at the reference mid-time in radians.

• wdE (float) – Apsidal precession rate in radians per epoch.

• E (array-like) – The epoch(s) at which to calculate the mid-times.

• primary (bool, optional) – If True, returns the transit mid-time. If False, returns the eclipse mid-time. Default is True.

Returns:

float or array-like – The predicted transit or eclipse times in \(\mathrm{BJD}_\mathrm{TDB}\).

References

[1]

Giménez and Bastero [1995]. https://doi.org/10.1007/BF00626903

[2]

Patra et al. [2017]. https://doi.org/10.3847/1538-3881/aa6d75

---

Radial Velocity Models

This module defines functions for modelling exoplanet radial velocity observations.

orbdot.models.rv_models.rv_constant(t0, P0, e0, w0, K, v0, dvdt, ddvdt, t)

Constant-period model for radial velocities.

This method calculates the expected radial velocity signal of a star that hosts a planet on an unchanging orbit.

Parameters:

• t0 (float) – Reference transit mid-time in \(\mathrm{BJD}_\mathrm{TDB}\).

• P0 (float) – The orbital period in days.

• e0 (float) – Eccentricity of the orbit.

• w0 (float) – Argument of pericenter of the planetary orbit in radians.

• K (float) – Radial velocity semi-amplitude in m/s.

• v0 (float) – Systemic radial velocity in m/s.

• dvdt (float) – First-order acceleration term in m/s/day.

• ddvdt (float) – Second-order acceleration term in m/s/day^2.

• t (float) – Measurement time(s) in \(\mathrm{BJD}_\mathrm{TDB}\).

Returns:

float or array-like – The expected radial velocity signal at the given time(s), including long-term trends.

Notes

Tip

In this model, the orbital elements do not evolve, which means that \(\omega_p=\omega_0\), \(e=e_0\), and \(P=P_0\).

The following steps outline the implementation of this method:

1. Calculate the true anomaly, eccentric anomaly, and mean anomaly at mid-transit.

In the OrbDot coordinate system, we know that the true anomaly at the transit center is:

\[\phi_0 = \frac{\pi}{2} - \omega_p\]

where \(\omega_p\) is the argument of pericenter of the planetary orbit. The eccentric anomaly may then be calculated by:

\[\tan \left(\frac{\mathrm{E}_0}{2}\right) = \sqrt{\frac{1-e}{1+e}} \tan \left(\frac{\phi_0}{2}\right)\]

which yields the mean anomaly via Kepler’s equation:

\[\mathrm{M}_0 = \mathrm{E}_0 - e \sin\mathrm{E}_0\]

2. Determine the time of pericenter passage.

The time of pericenter passage \(t_p\) is related to the mean anomaly at \(t_0\) by:

\[\mathrm{M}_0 = n \left(t_0 - t_p\right)\]

where \(n = 2 \pi/P\) is the mean motion.

3. Calculate the true anomaly of the planet at the given time.

With the time of pericenter passage, the planet’s mean anomaly may be calculated at any time \(t\):

\[\mathrm{M} = n \left(t - t_p\right)\]

allowing for the eccentric anomaly to be determined by solving Kepler’s equation:

\[\mathrm{M} = \mathrm{E} - e \sin \mathrm{E}\]

Finally, the true anomaly at time \(t\) is determined by:

\[\tan \left(\frac{\phi}{2}\right) = \sqrt{\frac{1+e}{1-e}} \tan \left(\frac{ \mathrm{E}}{2}\right)\]

4. Return the total RV signal, including long-term trends.

Thus, the total radial velocity signal at time \(t\) is:

\[v_r = K[\cos{(\phi+\omega_p)}+e\cos{\omega_p}] + \gamma_i + \dot{\gamma} \left(t-t_0\right) + \frac{1}{2} \ddot{\gamma} \left(t-t_0\right)^2\]

where \(K\) is the semi-amplitude of the planetary signal, \(\gamma_i\) is the systemic radial velocity, and \(\dot{ \gamma}\) and \(\ddot{\gamma}\) are first and second-order acceleration terms, respectively.

orbdot.models.rv_models.rv_decay(t0, P0, e0, w0, K, v0, dvdt, ddvdt, PdE, t)

Orbital decay model for radial velocities.

This method calculates the expected radial velocity signal of a star that hosts a planet on an orbit with a constant change in the period. Though the main application of this model is for orbital decay, a positive period derivative is allowed.

Note

This function assumes that the change in the RV semi-amplitude \(K\) is negligible.

Parameters:

• t0 (float) – Reference transit mid-time in \(\mathrm{BJD}_\mathrm{TDB}\).

• P0 (float) – Orbital period at the reference mid-time, in days.

• e0 (float) – Eccentricity of the orbit.

• w0 (float) – Argument of pericenter of the planetary orbit in radians.

• K (float) – Radial velocity semi-amplitude in m/s.

• v0 (float) – Systemic radial velocity in m/s.

• dvdt (float) – First-order acceleration term in m/s/day.

• ddvdt (float) – Second-order acceleration term in m/s/day^2.

• PdE (float) – Rate of change of the orbital period in days per epoch.

• t (float) – Measurement time(s) in \(\mathrm{BJD}_\mathrm{TDB}\).

Returns:

float or array-like – The expected radial velocity signal at the given time(s), including long-term trends.

Notes

Tip

In this model, the eccentricity and argument of pericenter do not evolve, which means that \(\omega_p=\omega_0\) and \(e=e_0\).

The following steps outline the implementation of this method:

1. Calculate the orbital period at the given time.

For a decaying orbit, the orbital period is dependant on time. Assuming that the rate-of-change is constant, the period at any epoch \(E\) is:

\[P\left(E\right) = P_0 + \frac{dP}{dE}\,E\]

where \(E\) is determined by truncating the result of the following equation to an integer value:

\[E = \frac{(t - t_0)}{P_0}\]

2. Calculate the latest time of pericenter passage.

As the orbit shrinks, the relative time of pericenter passage will evolve, so it must be calculated relative to the most recent transit mid-time:

\[t_{\mathrm{I}} = t_0 + PE + \frac{1}{2}\frac{dP}{dE}\,E^2\]

In the OrbDot coordinate system, we know that the true anomaly at the transit center is:

\[\phi_{\mathrm{I}} = \frac{\pi}{2} - \omega_p\]

The eccentric anomaly may then be calculated by:

\[\tan \left(\frac{\mathrm{E}_{\mathrm{I}}}{2}\right) = \sqrt{\frac{1-e}{ 1+e}}\tan \left( \frac{\phi_{\mathrm{I}}}{2}\right)\]

which yields the mean anomaly via Kepler’s equation:

\[\mathrm{M}_{\mathrm{I}} = \mathrm{E}_{\mathrm{I}} - e \sin\mathrm{E}_{\mathrm{I}}\]

Finally, the time of pericenter passage \(t_p\) is related to the mean anomaly by:

\[\mathrm{M}_{\mathrm{I}} = n \left(t_{\mathrm{I}} - t_p\right)\]

where \(n = 2 \pi/P\) is the mean motion.

3. Calculate the true anomaly of the planet at the given time.

With the time of latest pericenter passage, the planet’s mean anomaly at time \(t\) may be calculated by:

\[\mathrm{M} = n \left(t - t_p\right)\]

allowing for the eccentric anomaly to be determined by solving Kepler’s equation:

\[\mathrm{M} = \mathrm{E} - e \sin \mathrm{E}\]

Finally, the true anomaly at time \(t\) is determined by:

\[\tan \left(\frac{\phi}{2}\right) = \sqrt{\frac{1+e}{1-e}} \tan \left(\frac{ \mathrm{E}}{2}\right)\]

4. Return the total RV signal, including long-term trends.

Thus, the total radial velocity signal at time \(t\) is:

\[v_r = K[\cos{(\phi+\omega_p)}+e\cos{\omega_p}] + \gamma_i + \dot{\gamma} \left(t-t_0\right) + \frac{1}{2} \ddot{\gamma} \left(t-t_0\right)^2\]

where \(K\) is the semi-amplitude of the planetary signal, \(\gamma_i\) is the systemic radial velocity, and \(\dot{ \gamma}\) and \(\ddot{\gamma}\) are first and second-order acceleration terms, respectively.

orbdot.models.rv_models.rv_precession(t0, P0, e0, w0, K, v0, dvdt, ddvdt, wdE, t)

Apsidal precession model for radial velocities.

This method calculates the expected radial velocity signal of a star that hosts a planet on an elliptical orbit undergoing apsidal precession.

Parameters:

• t0 (float) – Reference transit mid-time in \(\mathrm{BJD}_\mathrm{TDB}\).

• P0 (float) – The sidereal period in days.

• e0 (float) – Eccentricity of the orbit.

• w0 (float) – Argument of pericenter of the planetary orbit at the reference mid-time, in radians.

• K (float) – Radial velocity semi-amplitude in m/s.

• v0 (float) – Systemic radial velocity in m/s.

• dvdt (float) – First-order acceleration term in m/s/day.

• ddvdt (float) – Second-order acceleration term in m/s/day^2.

• wdE (float) – Apsidal precession rate in radians per epoch.

• t (float) – Measurement time(s) in \(\mathrm{BJD}_\mathrm{TDB}\).

Returns:

float or array-like – The expected radial velocity signal at the given time(s), including long-term trends.

Notes

Tip

In this model, the eccentricity and orbital period not evolve, which means that \(e=e_0\) and \(P_s=P_0\), where \(P_s\) is the sidereal period.

The following steps outline the implementation of this method:

1. Calculate the anomalistic period and argument of pericenter.

For a precessing orbit, there is a distinction between the sidereal (observed) orbital period \(P_s\) and the anomalistic orbital period \(P_a\). They are related by:

\[P_s = P_a\left(1-\frac{d\omega/{dE}}{2\pi}\right)\]

where \(d\omega/{dE}\) is the apsidal precession rate in radians per epoch. The argument of pericenter of the planet’s orbit is a function of time, such that for any epoch \(E\) it is:

\[\omega_p\left(E\right) = \omega_0 + \frac{d\omega}{dE}\,E.\]

where \(E\) is determined by truncating the result of the following equation to an integer value:

\[E = \frac{(t - t_0)}{P_s}\]

2. Calculate the latest time of pericenter passage.

As the orbit is precessing, the time of pericenter passage must be calculated relative to the most recent transit mid-time:

\[t_{\mathrm{I}} = t_0 + P_s E - \frac{e P_a}{\pi}\cos{\omega_p}\]

In the OrbDot coordinate system, we know that the true anomaly at the transit center is:

\[\phi_{\mathrm{I}} = \frac{\pi}{2} - \omega_p\left(E\right)\]

The eccentric anomaly may then be calculated by:

\[\tan\left(\frac{\mathrm{E}_{\mathrm{I}}}{2}\right) = \sqrt{\frac{1-e}{1+e}} \tan\left(\frac{\phi_{\mathrm{I}}}{2}\right)\]

which yields the mean anomaly via Kepler’s equation:

\[\mathrm{M}_{\mathrm{I}} = \mathrm{E}_{\mathrm{I}} - e \sin\mathrm{E}_{\mathrm{I}}\]

Finally, the time of pericenter passage \(t_p\) is related to the mean anomaly by:

\[\mathrm{M}_{\mathrm{I}} = n \left(t_{\mathrm{I}} - t_p\right)\]

where \(n = 2 \pi/P_a\) is the mean motion.

3. Calculate the true anomaly of the planet at the given time.

With the time of latest pericenter passage, the planet’s mean anomaly at time \(t\) may be calculated by:

\[\mathrm{M} = n \left(t - t_p\right)\]

allowing for the eccentric anomaly to be determined by solving Kepler’s equation:

\[\mathrm{M} = \mathrm{E} - e \sin \mathrm{E}\]

Finally, the true anomaly at time \(t\) is determined by:

\[\tan \left(\frac{\phi}{2}\right) = \sqrt{\frac{1+e}{1-e}} \tan \left(\frac{ \mathrm{E}}{2}\right)\]

4. Return the total RV signal, including long-term trends.

Thus, the total radial velocity signal at time \(t\) is:

\[v_r = K[\cos{(\phi+\omega_p)}+e\cos{\omega_p}] + \gamma_i + \dot{\gamma} \left(t-t_0\right) + \frac{1}{2} \ddot{\gamma} \left(t-t_0\right)^2\]

where,

\[\omega_p\left(E\right) = w_0 + \frac{d\omega}{dE}\,E\]

and where \(K\) is the semi-amplitude of the planetary signal, \(\gamma_i\) is the systemic radial velocity, and \(\dot{ \gamma}\) and \(\ddot{\gamma}\) are first and second-order acceleration terms, respectively.

orbdot.models.rv_models.solve_keplers_equation(M, e, tol=1e-08)

An iterative solver of Kepler’s equation.

Parameters:

• M (array-like) – Mean anomalies in radians.

• e (float) – Eccentricity of the orbit.

• tol (float) – Tolerance for convergence (default: 1e-8).

Returns:

array-like – Eccentric anomalies in radians.

orbdot.models.rv_models.true_anomaly(t, t_peri, nu, e)

Calculates the true anomaly at any given time.

Parameters:

• t (array-like) – Time at which to calculate the true anomaly.

• t_peri (float) – Time of known pericenter passage.

• nu (float) – The mean motion in radians per day.

• e (float) – Eccentricity of the orbit.

Returns:

array-like – The true anomaly in radians.

---

Transit Duration Models

This module defines functions for modelling exoplanet transit durations.

orbdot.models.tdv_models.tdv_constant(P0, e0, w0, i0, E, M_s, R_s)

Constant-period model for transit durations.

This method returns the expected transit duration for a single planet on an unchanging orbit using the transit_duration method defined in this module.

Parameters:

• P0 (float) – Orbital period in days.

• e0 (float) – Eccentricity of the orbit.

• w0 (float) – Argument of pericenter of the planetary orbit in radians.

• i0 (float) – Line-of-sight inclination of the orbit in degrees.

• E (array-like) – The epoch(s) at which to calculate the transit duration.

• M_s (float) – Host star mass in solar masses.

• R_s (float) – Host star radius in solar radii.

Returns:

array-like – The predicted transit durations in minutes.

orbdot.models.tdv_models.tdv_decay(P0, e0, w0, i0, PdE, E, M_s, R_s)

Orbital decay model for transit durations.

This method calculates the expected transit duration for a single planet on an orbit with a constant change in the period, using Equation 58 from Kipping (2010) [1]_. Though the main application of this model is for orbital decay, a positive period derivative is allowed.

Parameters:

• P0 (float) – Orbital period at the reference mid-time in days.

• e0 (float) – Eccentricity of the orbit.

• w0 (float) – Argument of pericenter of the planetary orbit in radians.

• i0 (float) – Line-of-sight inclination of the orbit in degrees.

• PdE (float) – Rate of change of the orbital period in days per epoch.

• E (array-like) – The epoch(s) at which to calculate the transit duration.

• M_s (float) – Host star mass in solar masses.

• R_s (float) – Host star radius in solar radii.

Returns:

array-like – The predicted transit durations in minutes.

References

[1]

Kipping [2010]. https://doi.org/10.1111/j.1365-2966.2010.16894.x

orbdot.models.tdv_models.tdv_precession(P0, e0, w0, i0, wdE, E, M_s, R_s)

Apsidal precession model for transit durations.

This method calculates the expected transit duration at the given epoch(s) for a single planet on an elliptical orbit undergoing apsidal precession, using Equation 54 from Kipping (2010) [1]_.

Parameters:

• P0 (float) – Orbital period in days.

• e0 (float) – Eccentricity of the orbit.

• w0 (float) – Argument of pericenter of the planetary orbit at the reference mid-time in radians.

• i0 (float) – Line-of-sight inclination of the orbit in degrees.

• wdE (float) – Apsidal precession rate in radians per epoch.

• E (array-like) – The epoch(s) at which to calculate the transit duration.

• M_s (float) – Host star mass in solar masses.

• R_s (float) – Host star radius in solar radii.

Returns:

array-like – The predicted transit durations in minutes.

References

[1]

Kipping [2010]. https://doi.org/10.1111/j.1365-2966.2010.16894.x

orbdot.models.tdv_models.transit_duration(P0, e0, w0, i0, M_s, R_s)

Calculates the transit duration.

This method returns the expected transit duration using Equation 15 from Kipping (2010) [1]_.

Parameters:

• P0 (float) – Orbital period in days.

• e0 (float) – Eccentricity of the orbit.

• w0 (float) – Argument of pericenter of the planet’s orbit in radians.

• i0 (float) – Line-of-sight inclination of the orbit in degrees.

• M_s (float) – Host star mass in solar masses.

• R_s (float) – Host star radius in solar radii.

Returns:

float – The transit duration in minutes.

References

[1]

Kipping [2010]. https://doi.org/10.1111/j.1365-2966.2010.16894.x

---

Tools

Priors

The methods in this module are designed for the NestedSampling class, or similar algorithms for which sampling from a unit hypercube is required to explore the parameter space.

orbdot.tools.priors.gaussian_prior(hypval, prior)

Unit hypercube transformation for a free parameter with a Gaussian (Normal) prior.

Parameters:

• hypval (float) – Value in the unit hypercube (from 0 to 1).

• prior (tuple) – Tuple containing the mean and standard deviation of the Gaussian prior distribution, ie. ["gaussian", mean, std].

Returns:

float

orbdot.tools.priors.get_prior(hypval, prior)

Retrieves the value of a parameter from the unit hypercube given a prior distribution.

Parameters:

• hypval (float) – Value in the unit hypercube (from 0 to 1).

• prior (list or tuple) –

  Tuple that defines the prior distribution for the parameter. The first element specifies the type of prior ("uniform", "gaussian", or "log"), and the subsequent elements define the bounds of the distribution. The options are:

  1. Gaussian Prior: ["gaussian", mean, std]

  2. Uniform Prior: ["uniform", min, max]

  3. Log-Uniform Prior: ["log", log10(min), log10(max)]

Returns:

float

orbdot.tools.priors.log_prior(hypval, prior)

Unit hypercube transformation for a free parameter with a Log-Uniform prior.

Parameters:

• hypval (float) – Value in the unit hypercube (from 0 to 1).

• prior (tuple) – Tuple containing the minimum and maximum values of the log-uniform prior distribution, ie. ["log", log10(min), log10(max)].

Returns:

float

orbdot.tools.priors.uniform_prior(hypval, prior)

Unit hypercube transformation for a free parameter with a Uniform prior.

Parameters:

• hypval (float) – Value in the unit hypercube (from 0 to 1).

• prior (tuple) – Tuple containing the minimum and maximum values of the uniform prior distribution, ie. ["uniform", min, max].

Returns:

float

---

Stats

This module defines methods for statistical analyses and error propagation.

orbdot.tools.stats.calc_chi2(data, model, errors)

Calculates the chi-squared value for a given model and data.

Parameters:

• data (array_like) – The data.

• model (array_like) – The model values corresponding to the data points.

• errors (array_like) – The measurement uncertainties.

Returns:

float – The chi-squared value.

orbdot.tools.stats.credible_intervals(params, samples, dic, inds, circular=False)

Calculates the 68% credible intervals from a set of weighted posterior samples.

Parameters:

• params (list) – List of the free parameters.

• samples (array_like) – Array containing the weighted posterior samples.

• dic (dict) – Dictionary to store the credible intervals.

• inds (list) – List of indices that point to the desired parameters.

• circular (bool, optional) – Flag indicating whether the parameters are circular. Default is False.

Returns:

None

orbdot.tools.stats.propagate_err_ecosw_esinw(ecosw_results, esinw_results)

Derive uncertainties on "e0" and "w0" from "ecosw" and "esinw".

This method propagates the asymmetric uncertainties on "ecosw" and "esinw" to obtain upper and lower bounds on "e0" and "w0".

Parameters:

• ecosw_results (tuple) – Tuple containing the value, upper error, and lower error of "ecosw".

• esinw_results (tuple) – Tuple containing the value, upper error, and lower error of "esinw".

Returns:

tuple – A tuple of tuples, the first being the value "e0" and its upper and lower bounds, and the second being the value "w0" and its upper and lower bounds.

orbdot.tools.stats.propagate_err_sq_ecosw_sq_esinw(sq_ecosw_results, sq_esinw_results)

Derive uncertainties on "e0" and "w0" from "sq_ecosw" and "sq_esinw".

This method propagates the asymmetric uncertainties on "sq_ecosw" and "sq_esinw" to obtain upper and lower bounds on "e0" and "w0".

Parameters:

• sq_ecosw_results (tuple) – Tuple containing the value, upper error, and lower error of "sq_ecosw".

• sq_esinw_results (tuple) – Tuple containing the value, upper error, and lower error of "sq_esinw".

Returns:

tuple – A tuple of tuples, the first being the value "e0" and its upper and lower bounds, and the second being the value "w0" and its upper and lower bounds.

orbdot.tools.stats.quantiles(samples)

Calculates quantiles for the 68% credible intervals.

Parameters:

samples (array_like) – Array containing the samples.

Returns:

tuple

---

Utilities

This module defines a set of general functions that are used throughout the OrbDot package.

orbdot.tools.utilities.assign_default_values(system_info, planet_number)

Assigns default parameter values using the system information file.

Parameters:

• system_info (dict) – Dictionary of star-planet system properties.

• planet_number (int) – Planet number/index.

Returns:

dict – A dictionary containing the default values.

Note

If the eccentricity "e0" is zero, the argument of pericenter "w0" is zeroed.

orbdot.tools.utilities.bjd_to_hjd(bjd, RA, DEC)

Converts barycentric julian days (BJD) to heliocentric julian days (HJD).

This method uses the Astropy package to convert the time of an observation from barycentric julian days (BJD_TDB) to heliocentric julian days (HJD_UTC).

Parameters:

• bjd (list) – A list of times in barycentric julian days (BJD_TDB)

• RA (str) – The right ascension of the object formatted as: 00h00m00.0s.

• DEC (str) – The declination of the object formatted as: +00d00m00.0s.

Returns:

list – A list of times converted to heliocentric julian days (HJD).

orbdot.tools.utilities.calculate_epochs(tc, P, times, primary=True)

Calculate the epoch (orbit number) of a given transit or eclipse mid-time.

Parameters:

• tc (float) – The reference transit mid-time [BJD_TDB].

• P (float) – The orbital period in days.

• times (array-like) – The mid-times at which to calculate the epochs.

• primary (bool, optional) – If True, the mid-time(s) must be transits. If False, they are treated as secondary eclipses. Default is True.

Returns:

list – The epoch of the given mid-time(s).

orbdot.tools.utilities.hjd_to_bjd(hjd, RA, DEC)

Converts heliocentric julian days (HJD) to barycentric julian days (BJD).

This method uses the Astropy package to convert the time of an observation from heliocentric julian days (HJD_UTC) to barycentric julian days (BJD_TDB).

Parameters:

• hjd (list) – A list of times in heliocentric julian days (HJD_UTC).

• RA (str) – The right ascension of the object in the form 00h00m00.0s.

• DEC (str) – The declination of the object in the form +00d00m00.0s.

Returns:

list – The list of times converted to barycentric julian days (BJD_TDB).

orbdot.tools.utilities.merge_dictionaries(default_file, user_file)

Merge user input files with defaults.

Parameters:

• default_file (str) – Path name of the default input file.

• user_file (str) – Path name of the user-defined input file.

Returns:

dict – The merged dictionary.

Note

For any given key, the user-defined setting overrides the default setting.

orbdot.tools.utilities.raise_not_valid_param_error(requested_params, allowed_params, illegal_params)

Raises an exception if a given set of free parameters are not valid.

This method checks that the given free parameters are valid, ie. that they are a part of the physical model, and that certain pairs of parameters are not fit together.

Parameters:

• requested_params (list) – The requested free parameters.

• allowed_params (list) – The parameters that are part of the model.

• illegal_params (list) – The parameters that are not relevant to the model.

Returns:

None

orbdot.tools.utilities.read_rv_data(filename, delim=' ')

Reads an RV data file with columns: [Time (BJD), Velocity (m/s), Err (m/s), Source].

Parameters:

• filename (str) – The path name of the file that contains the data.

• delim (str, optional) – The column delimiter, default is space-delimited.

Returns:

dict – A dictionary containing the RV measurements, times, errors, and sources. Each value is a list of arrays, with the different arrays corresponding to separate RV data sources.

Raises:

FileNotFoundError – If the specified data file does not exist.

Note

The data are split by the instrument/source so that the instrument-specific parameters "v0" and "jit" are fit separately.

orbdot.tools.utilities.read_tdv_data(filename, delim=' ')

Reads a transit duration data file with columns: [Epoch, Duration (min), Error (min), Source].

Parameters:

• filename (str) – The path name of the file that contains the data.

• delim (str, optional) – The column delimiter, default is space-delimited.

Returns:

dict – A dictionary containing the transit durations, errors, sources, and epoch numbers.

Raises:

FileNotFoundError – If the specified data file does not exist.

Note

OrbDot does not currently support model fitting for secondary eclipse durations.

orbdot.tools.utilities.read_ttv_data(filename, delim=' ')

Reads a timing data file with columns: [Epoch, Time (BJD), Error (BJD), Source].

Parameters:

• filename (str) – The path name of the file that contains the data.

• delim (str, optional) – The column delimiter, default is space-delimited.

Returns:

dict – A dictionary containing the epochs, mid-times, errors, and sources.

Raises:

• FileNotFoundError – If the specified data file does not exist.

• ValueError – If the epoch type in the transit data is not recognized.

Note

The epoch (orbit number) of transit observations are integers, and eclipse observations are differentiated by a half orbit. For example, the eclipse that occurs after the 100th transit has an epoch equal to 100.5.

orbdot.tools.utilities.split_rv_instrument_params(source_order, source_tags, free_params)

Splits instrument-dependent free parameters when there are multiple sources of RV data.

If "jit" and/or "v0" are given as free parameters in a model fit and there is more than one source of radial velocity data, this method removes them from the list of free parameters and appends separate variable names for each source.

Parameters:

• source_order (list) – Order of the RV sources as read from the data file.

• source_tags (list) – Tags identifying the different RV data sources.

• free_params (array-like) – The list of free parameters in the model fit.

Returns:

array-like – The updated list of free parameters.

orbdot.tools.utilities.split_rv_instrument_results(free_params, source_order, source_tags, results_dic)

Splits instrument-dependent parameter results when there are multiple sources of RV data.

If "jit" and/or "v0" are given as free parameters in a model fit and there is more than one source of radial velocity data, this method splits the relevant results dictionary entry into separate keys for each source.

Parameters:

• free_params (array-like) – The list of free parameters in the model fit.

• source_order (list) – Order of the RV sources as read from the data file.

• source_tags (list) – Tags identifying the different RV data sources.

• results_dic (dict) – The results dictionary returned by the nested sampling methods.

Returns:

dict – The updated results dictionary.

orbdot.tools.utilities.wrap(angle)

Wrap an angle to its corresponding value in the range \(0\) to \(2\pi\).

Parameters:

angle (float) – An angle in radians.

Returns:

float – The wrapped angle in radians.

---

Plots

This module defines methods for creating various plots with OrbDot model fit results.

orbdot.tools.plots.corner_plot(dic, samples, params, outfile)

Generates a corner plot from the weighted posterior samples.

This method generates a corner plot from the weighted posterior samples using the corner.py package by Daniel Foreman-Mackey [1]_.

Parameters:

• dic (dict) – Dictionary containing the results of the model fit.

• samples (array_like) – The weighted posterior samples.

• params (list) – List of the free parameters.

• outfile (str) – File name for saving the plot.

Returns:

None

References

[1]

Daniel Foreman-Mackey (2016). https://doi.org/10.21105/joss.00024

orbdot.tools.plots.make_rv_plots(plot_settings, outfile, suffix='', model='constant')

Generates a radial velocity (RV) plot.

This method creates a 3-part plot of the radial velocity model and data:

1. A plot of the radial velocity measurements vs. time, with an option to plot the best-fit model. The RV data are shifted by the systemic radial velocity and the corresponding measurement times are shifted by the reference time.

2. A plot of the residuals of subtracting the best-fit radial velocity model from the data, including the systemic radial velocity and long-term trends.

3. A plot of the phase-folded RV signal due to the planet, ie. not including the systemic velocity and long-term trends, and 300 random posterior samples from the model fit.

Parameters:

• plot_settings (dict) – A dictionary containing plot settings.

• outfile (str) – The file path for saving the plot.

• suffix (str, optional) – Optional string for matching model fit results.

• model (str, optional) – The RV model that was fit, must be "constant", "decay", or "precession". Default is "constant".

Returns:

None

orbdot.tools.plots.make_tdv_plot(plot_settings, outfile, suffix='')

Generates a transit duration variation (TDV) plot.

Parameters:

• plot_settings (dict) – A dictionary containing plot settings.

• outfile (str) – The file path for saving the plot.

• suffix (str, optional) – Optional string for matching model fit results.

Returns:

None

orbdot.tools.plots.make_ttv_plot(plot_settings, outfile, suffix='')

Generates a transit timing variation (TTV) plot.

Parameters:

• plot_settings (dict) – A dictionary containing plot settings.

• outfile (str) – The file path for saving the plot.

• suffix (str, optional) – Optional string for matching model fit results.

Returns:

None

orbdot.tools.plots.periodogram(data_file, outfile, min_period=3, max_period=100000.0, num=1000000, peak_distance=100000, peak_power=0.4)

Generates a periodogram from the residuals of an RV model fit.

Note

This plot is not currently implemented in any model fitting routines.

Parameters:

• data_file (str) – Name of the file containing the residuals.

• outfile (str) – File name of the saved plot.

• min_period (float) – The minimum period in days.

• max_period (float) – The maximum period in days.

• num (int) – The number of frequencies in between those of min_period and max_period.

• peak_distance (float) – Limit of distance between peaks to be labelled.

• peak_power (float) – Label peaks above this cutoff in power.

Returns:

None

orbdot.tools.plots.read_random_samples(data_file, delim='\t')

Reads the *_random_samples.json files.

Parameters:

• data_file (str) – Name of the file containing the random posterior samples.

• delim (str, optional) – The data file delimiter, default is tab-separated.

Returns:

None