compute_RVs.py

Attention

This module may be moved into a separate, individual package in the future.

M_anom(t, T0, P)

Compute the mean anomaly of an orbit at a given time.

Note

The parameters must have the same units (of time), but the actual units can be freely chosen (e.g., yrs).

Parameters:

  • t (float) – The time at which the mean anomaly is evaluated.

  • T0 (float) – The reference epoch.

  • P (float) – The orbital period.

Returns:

M – The mean anomaly (radians) at time t.

Return type:

float

fzero_E_anom(E, t, T0, P, e)

Evaluate Kepler’s equation moved to one side (solve for the zeros of this function to get solutions for the eccentric anomaly, E).

Note

The parameters t, T0, and P must have the same units (of time), but the actual units can be freely chosen (e.g., yrs).

Parameters:

  • E (float) – The eccentric anomaly (radians).

  • t (float) – The time at which the anomalies are evaluated.

  • T0 (float) – The reference epoch.

  • P (float) – The orbital period.

  • e (float) – The orbital eccentricity.

Return type:

Kepler’s equation moved to one side.

E_anom(t, T0, P, e)

Solve for the eccentric anomaly at a given time.

Note

The parameters t, T0, and P must have the same units (of time), but the actual units can be freely chosen (e.g., yrs).

Parameters:

  • t (float) – The time at which the eccentric anomaly is evaluated.

  • T0 (float) – The reference epoch.

  • P (float) – The orbital period.

  • e (float) – The orbital eccentricity.

Returns:

E – The eccentric anomaly (radians).

Return type:

float

nu_anom(E, e)

Compute the true anomaly from the eccentric anomaly and orbital eccentricity.

Parameters:

  • E (float) – The eccentric anomaly (radians).

  • e (float) – The orbital eccentricity.

Returns:

nu – The true anomaly (radians).

Return type:

float

RV_true(t, K, P, T0=0.0, e=0.0, w=0.0, gamma=0.0)

Compute the true radial velocity of an orbit at a given time.

Parameters:

  • t (float) – The time (yrs) at which to compute the radial velocity.

  • K (float) – The radial velocity semi-amplitude (m/s).

  • P (float) – The orbital period (yrs).

  • T0 (float, default=0.) – The reference epoch (yrs).

  • e (float, default=0.) – The orbital eccentricity.

  • w (float, default=0.) – The argument of pericenter (radians).

  • gamma (float, default=0.) – The radial velocity reference zero-point/offset (m/s).

Returns:

V_r – The radial velocity (m/s) at time t.

Return type:

float

RV_true_sys(t, K_sys, P_sys, T0_sys, e_sys, w_sys, gamma=0.0)

Compute the true radial velocity of a planetary system at a given time.

This is the sum of the radial velocities of the individual planets in the system, each computed using syssimpyplots.compute_RVs.RV_true().

Parameters:

  • t (float) – The time (yrs) at which to compute the radial velocity.

  • K_sys (array[float]) – The radial velocity semi-amplitudes (m/s) of the planets.

  • P_sys (array[float]) – The orbital periods (yrs).

  • T0_sys (array[float]) – The reference epochs (yrs).

  • e_sys (array[float]) – The orbital eccentricities.

  • w_sys (array[float]) – The arguments of pericenter (radians).

  • gamma (float, default=0.) – The radial velocity reference zero-point/offset (m/s).

Returns:

V_r – The radial velocity (m/s) at time t.

Return type:

float

rv_K(m, P, e=None, i=None, Mstar=1.0)

Compute the radial velocity semi-amplitude of each planet in a system.

Parameters:

  • m (array[float]) – The planet masses (Earth masses).

  • P (array[float]) – The orbital periods (days).

  • e (array[float], optional) – The orbital eccentricities (assumed all zero if not provided).

  • i (array[float], optional) – The orbital inclinations (radians) relative to the sky plane (assumed to be all 90 degrees, i.e. edge-on, if not provided).

  • Mstar (float, default=1.) – The stellar mass (solar masses).

Returns:

K – The radial velocity semi-amplitudes (m/s).

Return type:

array[float]

fit_rv_K_single_planet_model_GLS(t_obs, RV_obs, covarsc, P, T0=0.0, e=0.0, w=0.0)

Fit the radial velocity (RV) semi-amplitude (K) of a single planet model given an RV time series.

Assumes fixed/known values for all other parameters of the planetary orbit, and uses generalized least squares (GLS) to solve for K.

Parameters:

  • t_obs (array[float]) – The observation times (days) corresponding to the RV observations.

  • RV_obs (array[float]) – The RV observations (m/s).

  • covarsc (array[float]) – The covariance matrix of measurement uncertainties (2-d array).

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

  • T0 (float, default=0.) – The reference epoch (days).

  • e (float, default=0.) – The orbital eccentricity.

  • w (float, default=0.) – The argument of pericenter (radians).

Returns:

  • K_hat (float) – The estimator for the RV semi-amplitude (m/s).

  • sigma_K (float) – The estimated uncertainty in the RV semi-amplitude (m/s).

fit_rv_Ks_multi_planet_model_GLS(t_obs, RV_obs, covarsc, P_all, T0_all, e_all, w_all)

Fit the radial velocity (RV) semi-amplitudes (K’s) of a multi-planet model given an RV time series.

Assumes fixed/known values for all other parameters of the planetary orbits, and uses generalized least squares (GLS) to solve for K.

Parameters:

  • t_obs (array[float]) – The observation times (days) corresponding to the RV observations.

  • RV_obs (array[float]) – The RV observations (m/s).

  • covarsc (array[float]) – The covariance matrix of measurement uncertainties (2-d array).

  • P_all (array[float]) – The orbital periods (days).

  • T0_all (array[float]) – The reference epochs (days).

  • e_all (array[float]) – The orbital eccentricities.

  • w_all (array[float]) – The arguments of pericenter (radians).

Returns:

  • K_hat_all (array[float]) – The estimators for the RV semi-amplitudes (m/s).

  • sigma_K_all (array[float]) – The estimated uncertainties in the RV semi-amplitudes (m/s).

conditionals_dict(P_cond_bounds=None, Rp_cond_bounds=None, Mp_cond_bounds=None, det=True)

Create a dictionary with conditionals for planetary systems.

Allows for setting bounds on orbital period, planet radius, planet mass, and transit detectability.

Parameters:

  • P_cond_bounds (list, optional) – The [lower, upper] bounds on orbital period (days).

  • Rp_cond_bounds (list, optional) – The [lower, upper] bounds on planet radius (Earth radii).

  • Mp_cond_bounds (list, optional) – The [lower, upper] bounds on planet mass (Earth masses).

  • det (bool, default=True) – Whether to require the planets to be detected in transits.

Returns:

conds – A dictionary of conditionals.

Return type:

dict

The dictionary contains the following fields:

  • P_lower: The lower bound on orbital period (days).

  • P_upper: The upper bound on orbital period (days).

  • Rp_lower: The lower bound on planet radius (Earth radii).

  • Rp_upper: The upper bound on planet radius (Earth radii).

  • Mp_lower: The lower bound on planet mass (Earth masses).

  • Mp_upper: The upper bound on planet mass (Earth masses).

  • det: Whether to require the planets to be detected in transits (True/False).

condition_planets_bools_per_sys(sssp_per_sys, conds)

Compute a 2-d array of booleans indicating the planets passing the conditionals in a simulated physical catalog.

The 2-d array will match the shape of the 2-d arrays of planet properties in sssp_per_sys.

Parameters:
Returns:

bools_cond_per_sys – The 2-d array of booleans indicating which planets pass all the conditionals in conds.

Return type:

array[bool]

Note

The output includes rows without any conditioned planets (i.e. rows with all False values), since it must have the same shape as the 2-d arrays of planet properties in sssp_per_sys in order to be able to index them.

condition_systems_indices(sssp_per_sys, conds)

Compute an array of indices indicating which systems in a simulated physical catalog contain at least one conditioned planet.

Parameters:
Returns:

i_cond – The array of indices indicating which systems contain at least one conditioned planet.

Return type:

array[int]

Plot a gallery of systems conditioned on a given planet, along with their radial velocity (RV) time series.

Parameters:

  • sssp_per_sys (dict) – The dictionary containing the planetary and stellar properties for each system in a physical catalog (2-d and 1-d arrays), e.g. returned by the function syssimpyplots.load_sims.compute_summary_stats_from_cat_phys().

  • sssp (dict) – A dictionary containing the planetary and stellar properties of all planets in a physical catalog (1-d arrays), e.g. returned by the function syssimpyplots.load_sims.compute_summary_stats_from_cat_phys().

  • conds (dict) – The dictionary of conditionals, e.g. returned by the function syssimpyplots.compute_RVs.conditionals_dict().

  • outputs_RVs (structured array, optional) – A table of RV simulation results (TODO: need more details).

  • mark_undet (bool, default=True) – Whether to outline the undetected planets with red.

  • fit_RVs (bool, default=False) – Whether to simulate how many RV observations are required to measure the RV semi-amplitude of the conditioned planet in each system.

  • N_obs_all (list[int], optional) – The numbers of RV observations to test, if fit_RVs=True.

  • repeat (int, default=1000) – The number of times to repeat the RV simulations for each system, if fit_RVs=True.

  • σ_1obs (float, default=1.) – The single-measurement RV precision (m/s), if fit_RVs=True.

  • t_obs_σ (float, default=0.2) – The standard deviation (days) in nightly RV observation times, if fit_RVs=True.

  • fig_size (tuple, default=(9,10)) – The figure size.

  • seed (int, optional) – A random seed, for reproducible results.

  • N_sample (int, default=20) – The number of systems with a conditioned planet to sample and plot.

  • N_per_plot (int, default=20) – The number of systems to plot per figure.

  • afs (int, default=12) – The axes fontsize.

  • tfs (int, default=12) – The text fontsize.

  • save_name_base (str, default='no_name_fig') – The start of the file names for saving the figures.

  • save_fig (bool, default=False) – Whether to save the figures. If True, will save each figure in the working directory with the file name given by save_name_base with an index appended.

fit_RVobs_systems_conditional(sssp_per_sys, sssp, conds, N_obs_all, cond_only=False, fit_sys_cond=True, fit_all_planets=False, N_sample=20, repeat=1000, σ_1obs=1.0, obs_mode='daily')

Simulate the fitting of radial velocity (RV) observations to a sample of systems with a conditioned planet from a simulated physical catalog, to measure the RV semi-amplitudes of the conditioned planets.

Parameters:

  • sssp_per_sys (dict) – The dictionary containing the planetary and stellar properties for each system in a physical catalog (2-d and 1-d arrays), e.g. returned by the function syssimpyplots.load_sims.compute_summary_stats_from_cat_phys().

  • sssp (dict) – A dictionary containing the planetary and stellar properties of all planets in a physical catalog (1-d arrays), e.g. returned by the function syssimpyplots.load_sims.compute_summary_stats_from_cat_phys().

  • conds (dict) – The dictionary of conditionals, e.g. returned by the function syssimpyplots.compute_RVs.conditionals_dict().

  • N_obs_all (list[int]) – The numbers of RV observations to test.

  • cond_only (bool, default=False) – Whether to simulate an idealized case in which there are no planets other than the conditioned planet in each system.

  • fit_sys_cond (bool, default=True) – Whether to simulate the realistic case in which we fit the RV semi-amplitude of the conditioned planet only.

  • fit_all_planets (bool, default=False) – Whether to simulate an optimistic case in which we fit the RV semi-amplitudes of all the planets in each system.

  • N_sample (int, default=20) – The number of systems with a conditioned planet to sample.

  • repeat (int, default=1000) – The number of times to repeat the RV simulations for each system.

  • σ_1obs (float, default=1.) – The single-measurement RV precision (m/s).

  • obs_mode ({'daily', 'random'}) – How the RV observation times are drawn. daily will draw times spaced by a day with some added variation to avoid aliasing; random will draw completely random times over a set baseline (150 days).

Returns:

outputs – A table with the averaged results of the RV simulations for each conditioned system.

Return type:

structured array

The table has the following default columns:

  • id_sys: The index of the system.

  • n_pl: The total number of planets in the system.

  • P_cond: The period (days) of the conditioned planet.

  • Rp_cond: The radius (Earth radii) of the conditioned planet.

  • Mp_cond: The mass (Earth masses) of the conditioned planet.

  • K_cond: The RV semi-amplitude (m/s) of the conditioned planet.

  • K_max: The maximum RV semi-amplitude (m/s) of the planets in the system.

  • K_sum: The sum of the RV semi-amplitudes (m/s) of the planets in the system.

If cond_only=True, will also have the following columns:

  • N_obs_min_{XX}p_ideal: The minimum number of RV observations required to measure K_cond to within {XX} percent error in the ideal case, where {XX} = 30, 20, 10, and 5.

  • rms_sigma_K_{XX}p_ideal: The root mean square of the uncertainties in the measured K_cond when K_cond is measured to within {XX} percent error in the ideal case, where {XX} = 30, 20, 10, and 5.

  • rmsd_best_ideal: The best (fractional) root mean square deviation of the measured K_cond from the true K_cond in the ideal case, from the number of RV observations tested.

If fit_sys_cond=True, will also have the following columns:

  • N_obs_min_{XX}p: The minimum number of RV observations required to measure K_cond to within {XX} percent error in the realistic case, where {XX} = 30, 20, 10, and 5.

  • rms_sigma_K_{XX}p: The root mean square of the uncertainties in the measured K_cond when K_cond is measured to within {XX} percent error in the realistic case, where {XX} = 30, 20, 10, and 5.

  • rmsd_best: The best (fractional) root mean square deviation of the measured K_cond from the true K_cond in the realistic case, from the number of RV observations tested.

If fit_all_planets=True, will also have the following columns:

  • N_obs_min_{XX}p_fitall: The minimum number of RV observations required to measure K_cond to within {XX} percent error in the optimistic case, where {XX} = 30, 20, 10, and 5.

  • rms_sigma_K_{XX}p_fitall: The root mean square of the uncertainties in the measured K_cond when K_cond is measured to within {XX} percent error in the optimistic case, where {XX} = 30, 20, 10, and 5.

  • rmsd_best_fitall: The best (fractional) root mean square deviation of the measured K_cond from the true K_cond in the optimistic case, from the number of RV observations tested.

fit_RVobs_single_planets_vs_K(K_array, N_obs_all, P_bounds, alpha_P=0.0, sigma_ecc=0.25, repeat=1000, σ_1obs=1.0, t_obs_σ=0.2)

Simulate the fitting of radial velocity (RV) observations of single planets as a function of their RV semi-amplitude (K).

Note

Draws single planets with periods from a power-law distribution and eccentricities from a Rayleigh distribution, while assuming a fixed array of K (thus, the mass of the planet will change).

Parameters:

  • K_array (list or array[float]) – The RV semi-amplitudes (m/s) to simulate.

  • N_obs_all (array[int]) – The numbers of RV observations to test.

  • P_bounds (list or tuple) – The (lower, upper) bounds on the orbital periods (days).

  • alpha_P (float, default=0.) – The power-law index for the period distribution.

  • sigma_ecc (float, default=0.25) – The Rayleigh scale for the eccentricity distribution.

  • repeat (int, default=1000) – The number of times to repeat the RV simulations for each system.

  • σ_1obs (float, default=1.) – The single-measurement RV precision (m/s).

  • t_obs_σ (float, default=0.2) – The standard deviation (days) in nightly RV observation times.

Returns:

outputs – A table with the averaged results of the RV simulations for each single planet.

Return type:

structured array

The table has the following columns:

  • K: The RV semi-amplitude (m/s) of the planet.

  • N_obs_min_20p: The minimum number of RV observations required to measure K to within 20 percent error.

  • rmsd_best: The best (fractional) root mean square deviation of the measured K from the true K, from the number of RV observations tested.

Warning

The results (each row) averages over the periods and eccentricities drawn for each planet, which are not independent from K and may have large effects on the minimum number of observations required!

plot_scatter_K_vs_P_conditional(sssp_per_sys, sssp, conds, log_y=False, fig_size=(8, 5), afs=20, tfs=20, lfs=16, save_name_base='no_name_fig', save_fig=False)

Plot the radial velocity semi-amplitude (K) versus orbital period (P) for all planets in systems with a conditioned planet.

Makes two figures: (1) a scatter plot of K/K_max versus P, where K_max is the maximum RV semi-amplitude of the planets in each system, and (2) a scatter plot of K versus P.

Parameters:

  • sssp_per_sys (dict) – The dictionary containing the planetary and stellar properties for each system in a physical catalog (2-d and 1-d arrays), e.g. returned by the function syssimpyplots.load_sims.compute_summary_stats_from_cat_phys().

  • sssp (dict) – A dictionary containing the planetary and stellar properties of all planets in a physical catalog (1-d arrays), e.g. returned by the function syssimpyplots.load_sims.compute_summary_stats_from_cat_phys().

  • conds (dict) – The dictionary of conditionals, e.g. returned by the function syssimpyplots.compute_RVs.conditionals_dict().

  • log_y (bool, default=False) – Whether to plot the y-axis (K/K_max or K) on a log-scale.

  • fig_size (tuple, default=(8,5)) – The figure size.

  • afs (int, default=20) – The axes fontsize.

  • tfs (int, default=20) – The text fontsize.

  • lfs (int, default=16) – The legend fontsize.

  • save_name_base (str, default='no_name_fig') – The start of the file names for saving the figures.

  • save_fig (bool, default=False) – Whether to save the figures. If True, will save each figure in the working directory with the file name given by save_name_base with either ‘v1’ (for K/K_max vs. P) or ‘v2’ (for K vs. P) appended.

generate_latex_table_RVobs_systems_conditional(outputs_RVs, nan_str='$>1000$', N_sample=40)

Make a LaTeX syntax-formatted table with the results of a simulation fitting radial velocity (RV) observations.

Parameters:

  • outputs_RVs (structured array) – A table of RV simulation results (TODO: need more details).

  • nan_str (str, default=r'$>1000$') – The string to replace NaN values.

  • N_sample (int, default=40) – The number of rows/systems from outputs_RVs to include in the table.

Returns:

table_array – An array of strings that will produce the rows of the LaTeX table.

Return type:

array[str]

fit_line_loglog_Nobs_K_single_planets(outputs_ideal, σ_1obs, p0)

Fit a line to a number of points of log(N_obs) versus log(K), where ‘N_obs’ is the minimum number of radial velocity (RV) observations required to measure the RV semi-amplitude (‘K’) to within 20 percent error.

Equivalent to fitting N_obs as a power-law function of K.

Parameters:

  • outputs_ideal (dict) – The dictionary containing the results of the RV simulations of the ideal case, including the fields K for the RV semi-amplitudes (m/s) and N_obs_min_20p for the minimum number of RV observations required to measure K to within 20 percent error.

  • σ_1obs (float) – The single measurement RV precision (m/s) to serve as the normalization point (also represented by ‘K_norm’).

  • p0 (list) – The initial guesses for the parameters of the line, [normalization, slope] where ‘normalization’ is the ‘N_obs’ at ‘K_norm’.

Returns:

  • logN (float) – The log of the normalization of the fitted line, log10(N_obs) at ‘K_norm’.

  • slope (float) – The slope of the fitted line.

linear_logNobs_logK(K, K_norm, Nobs_norm, slope, Nobs_min=5, round_to_ints=True)

Return the required number of radial velocity (RV) observations (‘N_obs’) as a function of the RV semi-amplitude (‘K’) given a linear relation for log(N_obs) versus log(K).

Parameters:

  • K (list or array[float]) – The RV semi-amplitudes (m/s) at which to predict ‘N_obs’.

  • K_norm (float) – The normalization point (m/s) corresponding to Nobs_norm.

  • Nobs_norm (int) – The required number of RV observations at the normalization point K_norm.

  • slope (float) – The slope of the linear relation.

  • Nobs_min (int, default=5) – The minimum number of RV observations possible. All predicted values less than this value will be set to this value.

  • round_to_ints (bool, default=True) – Whether to round each resulting value of ‘N_obs’ to the nearest integer.

Returns:

Nobs_K – The required number of RV observations at each RV semi-amplitude in K.

Return type:

array[int] or array[float]