Reproducing real TTVs with ttvfast-python #20
Description
Hello,
I am trying to reproduce the TTVs from (Holman, M., J, et al 2010): Kepler-9 - A system of Multiple Planets Transiting a Sun-like Star, Confirmed bu Timing Variations. A google image of the TTVs can be seen hopefully here (top figure you see Kepler9-b in blue and kepler9-c in red) https://www.google.com/imgres?imgurl=https%3A%2F%2Fscience.sciencemag.org%2Fcontent%2Fsci%2F330%2F6000%2F51%2FF3.large.jpg&imgrefurl=https%3A%2F%2Fscience.sciencemag.org%2Fcontent%2F330%2F6000%2F51&tbnid=mZIzIzB4zCNX6M&vet=12ahUKEwjwtZrwzbrxAhUXkZUCHbPLAkoQMygDegUIARCWAQ..i&docid=Xlb5WQ8q-INRpM&w=755&h=1280&q=Kepler-9%3A%20A%20System%20of%20Multiple%20Planets%20Transiting%20a%20Sun-Like%20Star%2C%20Confirmed%20by%20Timing%20Variations&ved=2ahUKEwjwtZrwzbrxAhUXkZUCHbPLAkoQMygDegUIARCWAQ. The idea is to use TTVfast-python for a different target for a paper, however I am struggling a bit to understand the code behavior. The piece of code I am using follows:
import ttvfast
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook
# gravity = 0.000295994511 # AU^3/day^2/M_sun
stellar_mass = 1.0#0.95573417954 # M_sun
#planet parameters come mostly from exofop and posterior dist. from previous modeling. except for longnode, argument, mean_anomaly
planet1 = ttvfast.models.Planet(
mass=0.0001303490487, # M_sun
period=19.23891, # days
eccentricity=0.0609,
inclination=88.982, # degrees
longnode=0., # degrees
argument=90., # degrees
mean_anomaly=0., # degrees
)
planet2 = ttvfast.models.Planet(
mass=8.98076785e-5,
period=38.9853,
eccentricity=0.06691,
inclination=89.188,
longnode=0.,
argument=90.,
mean_anomaly=0.,
)
planets = [planet1, planet2]
Time = 0 # days
dt = 1e-3 # days
Total = 1000 # days
results = ttvfast.ttvfast(planets, stellar_mass, Time, dt, Total)#, input_flag=0)
#Planet1
idx_planet1 = np.array(results['positions'][0]) == 0
planet1_epochs = np.array(results['positions'][1])[idx_planet1]
planet1_times = np.array(results['positions'][2])[idx_planet1]
#make epoch matrix
A = np.vstack([np.ones(planet1_epochs.size), planet1_epochs]).T
#clean -2 from arrays
cond = planet1_times == -2.
planet1_epochs = planet1_epochs[~cond]
planet1_times = planet1_times[~cond] #np.where(~cond, planet1_times, 1)#
A = np.copy(A[~cond])
c, m = np.linalg.lstsq(A, planet1_times, rcond=-1)[0] #leastsq to get best t0 and p
#Planet2
idx_planet2 = np.array(results['positions'][0]) == 1
planet2_epochs = np.array(results['positions'][1])[idx_planet2]
planet2_times = np.array(results['positions'][2])[idx_planet2]
A2 = np.vstack([np.ones(planet2_epochs.size), planet2_epochs]).T
#clean same as before
cond2 = planet2_times == -2.
planet2_epochs = planet2_epochs[~cond2]
planet2_times = planet2_times[~cond2]
A2 = np.copy(A2[~cond2])
c2, m2 = np.linalg.lstsq(A2, planet2_times, rcond=-1)[0]
fig, ax = plt.subplots(1,1)
ax.plot(A.T[1], (planet1_times-m*A.T[1]-c)*(60.*24.), '.', label='Kepler-9b')
ax.plot(A2.T[1], (planet2_times-m2*A2.T[1]-c2)*(60.*24.), '.', label='Kepler-9c')
ax.set_xlabel("Transit number")
ax.set_ylabel("TTV [min]")
ax.legend()
Questions
(0) Why do I constantly get -2, is there a way to simulate the TTVs where I do not get it. From the c version README "BAD_TRANSIT = -1 /* If the bisection method is passed a window that does not bracket the root uniquely, or if the bisection method does not converge within MAX_ITER iteractions, the code returns -1, as discussed in the paper", However I get -2 not -1. Not sure if I should worry about that. In the code above I remove these -2 making sure I do not lose track of the transit times.
(1) If I run the code for different Total = 250, 500, 1000, days I get a different TTV amplitude and shape. I was expecting that a larger Total time would only give me more TTV points, therefore for a zoom at x = [0,12] I should have the same TTV amplitudes more or less. Why not? The figures from top to bottom are for total 250, 500, 1000. I did some tests changing dt to shorter integration steps, and also initial times to negative values.
(2) Usually we do not have the following parameters longnode, argument, mean_anomaly. If I want to reproduce observed TTVs, how exactly may I set these parameters by hand. In the example above, I believe I set both planets to start integration at periapse passage. Moreover, as in question (1) If I change the initial orbital configuration, the TTV shape changes drastically (depending on the values) do you know why? I cannot get a TTV simulation close to the amplitudes from the paper as seen in the image.
I really appreciate any comment/help, thank you!