<jupyter_start><jupyter_text>TOTEMS - Tidal Orbital decay Timing Extrapolation & Modelling Software**<NAME>, last update 2020-12-24 | [Github](https://github.com/Ross-Dobson/TOTEMS)** IntroductionThis code analyses Exoplanet mid-transit times $t_\text{mid}$ for signs of orbital decay, specifically those caused by tidal interaction. It also attempts to ascertain between these, or the similar effects of apsidal precession of the orbit. In this notebook, I provide an example of doing this with WASP-12b. This code was originally written as my contribution to a group project during my undergrad at University College London; it has been considerably rewritten and improved upon since then. Future plans are to write this as standalone Python code, to be incorporated into a PyPI package. TheoryConsider that a stable, circular orbit with no decay can be modelled as:$$ t(E)=t_0+P\times{E}, $$ where $E$ is the epoch number of periods/transits since the reference mid-time $t_0$, and $P$ is the orbital period.A decaying orbit due to tidal interaction can be modelled as$$ t(E)=t_0+P\times{E} + \frac{1}{2}\frac{dP}{dE}E^2, $$and apsidal precession can be modelled as$$ t(E)=t_0+P_s\times{E} - \frac{eP_a}{\pi}\cos\left({\omega(E)}\right)$$ where orbital period is now $P_s$, the sidereal period, to distinguish from anomalistic period $P_a=P_s\left(1-\frac{d\omega/dE}{2\pi}\right)^{-1} $, and the argument of pericenter evolves as $ \omega(E) = \omega_0 + \frac{d\omega}{dE}E $. Reference: [<NAME> *et al* 2017 AJ 154 4](https://doi.org/10.3847/1538-3881/aa6d75) Data PreparationPrepare a CSV with three columns: the mid-time in Julian Date, the uncertainty on this, and the type of JD: either HJD_UTC, HJD_TDB, BJD_UTC, BJD_TDB. All times will be converted if necessary to BJD_TDB automatically.Alternatively, or in addition, mid-time data can be imported from the [Exoplanet Transit Database](http://var2.astro.cz/ETD/). However, this comes from a variety of sources of differing quality, so it in the authors experience is best to avoid using solely this data.If you have values for the (initial) orbital period $P_0$ and the mid-time for epoch 0, $t_0$, you can input these; else the value for $P_0$ will be obtained from ETD, the value for $t_0$ will be set at the earliest $t_\text{mid}$. To convert mid-times into BJD$_\text{TDB}$, the [`pylightcurve`](https://github.com/ucl-exoplanets/pylightcurve) module is utilised. This was initially provided to me by Dr. <NAME> of University College London.<jupyter_code>!pip install pylightcurve
import numpy as np
import pandas as pd
from tqdm import tqdm
import matplotlib.pyplot as plt

from scipy.optimize import curve_fit
from scipy.stats import sigmaclip

# Used for BJD HJD conversions - thanks to Angelos Tsiaras
import pylightcurve as plc
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Planet propertiesUsing the [Open Exoplanet Catalogue](https://github.com/OpenExoplanetCatalogue/open_exoplanet_catalogue), we can automatically obtain a large number of parameters for most exoplanets. A lot of this information is loaded from NASA Exoplanet Archive, see the linked Github page for more information.<jupyter_code>import xml.etree.ElementTree as ET, urllib.request, gzip, io
oec_url = "https://github.com/OpenExoplanetCatalogue/oec_gzip/raw/master/systems.xml.gz"
oec = ET.parse(gzip.GzipFile(fileobj=io.BytesIO(urllib.request.urlopen(oec_url).read())))<jupyter_output><empty_output><jupyter_text>A list of all planets in the database (to check if your target is there) can be generated with:```for planet in oec.findall(".//planet"):    print(planet.findtext("name"))```Beware this produces a long list. The below function returns the Right Ascension (hh mm ss.ss), Declination (dd mm ss.ss), Eccentricity $e$, orbital period $P$ (days), and Argument of Periastron $\omega$ (deg), if available.<jupyter_code>def get_oec_planet_data(planet_name):
    """Obtains RA, Dec, e, P, and Periastron from Open Exoplanet Catalogue.

    Also prints the full system and planetary properties. Planet name needs to 
    match one on the list. To see the full list, run::

        for planet in oec.findall(".//planet"):
            print(planet.findtext("name"))
    
    Parameters
    ----------
    planet_name : string
        The name of the exoplanet to get the OEC data of.

    Returns
    -------
    ra : string
        Right Ascension of the system, format hh mm ss.ss
    dec : string
        Declination of the system, format dd mm ss.ss
    period : float
        Period of the planetary orbit in days
    eccentricity : float
        Eccentricity of the planetary orbit
    periastron : float
        Argument of periastron of the orbit (degrees)

    Notes
    -----
    If a value is not present in the catalogue, a None will be returned.
    """
    system_name = planet_name.split()[0]

    print("System Properties:")
    for properties in oec.findall(".//system[name='{0}']/*".format(system_name)):
        print("\t" + properties.tag + ":", properties.text)

    print("Planet Properties:")
    for properties in oec.findall(".//planet[name='{0}']/*".format(planet_name)):
        print("\t" + properties.tag + ":", properties.text)

    for system in oec.findall(".//system[name='{0}']".format(system_name)):
        ra = system.findtext("rightascension")
        dec = system.findtext("declination")

    for planet in oec.findall(".//planet[name='{0}']".format(planet_name)):
        period = planet.findtext("period")
        eccentricity = planet.findtext("eccentricity")
        periastron = planet.findtext("periastron")

    return ra, dec, period, eccentricity, periastron<jupyter_output><empty_output><jupyter_text>For the example planet, WASP-12b:<jupyter_code>planet_name='WASP-12 b'
wasp12_ra, wasp12_dec, wasp12_period, wasp12_eccentricity, wasp12_periastron = get_oec_planet_data(planet_name)
wasp12_period = float(wasp12_period)
wasp12_eccentricity = float(wasp12_eccentricity)
wasp12_periastron = float(wasp12_eccentricity)<jupyter_output>System Properties:
	name: WASP-12
	rightascension: 06 30 33
	declination: +29 40 20
	distance: 250
	binary: 
		
Planet Properties:
	name: WASP-12 b
	name: WASP-12 A b
	name: 2MASS J06303279+2940202 b
	list: Confirmed planets
	mass: 1.404
	radius: 1.736
	period: 1.0914222
	semimajoraxis: 0.02293
	eccentricity: 0.0
	inclination: 86
	spinorbitalignment: 59
	transittime: 2454508.98074
	discoverymethod: transit
	lastupdate: 13/07/31
	description: The host star of WASP-12 b is the primary of a hierarchical triple system. The distant binary consists of two M dwarfs.
	discoveryyear: 2008
	temperature: 2319.7
	list: Planets in binary systems, S-type
	istransiting: 1<jupyter_text>Not every quantity is always present, or given in a suitable level of precision. We can see that for WASP-12, the $\frac{dP}{dE}$ and the argument of periastron are not given, so we need to manually override it with the value from another source. The eccentricity is also not given in enough precision, so we use a value from the same source. 

Source: [<NAME> *et al* 2017 AJ 154 4](https://doi.org/10.3847/1538-3881/aa6d75)<jupyter_code>wasp12_eccentricity = .00208 # unitless
wasp12_periastron = 2.92  # degrees
wasp12_dPdE = -1.02e-9  # days/epoch
wasp12_dwdE = 0.00133  # rad / epoch<jupyter_output><empty_output><jupyter_text>Transit times One option is to get data from the [Exoplanet Transit Database](http://var2.astro.cz/ETD). This data tends to be from multiple sources of differing techniques and quality, however. In the author's experience, it should be used only as a last resort, or to bolster data from any recent published paper if few mid-times are available. If using solely ETD data, it usually makes sense to use their values of $P_0$ and $t_0$ as to inherit their epoch system.<jupyter_code>def get_etd_p0_t0(etd_url):
    """Obtains the values of P_0 (period) and t_0 (tmid at epoch 0) from ETD.
    
    Parameters
    ----------
    etd_url : string
        The URL to the .csv file for the planet on ETD.

    Returns
    -------
    p0
    """
    # Load in the row, split in two via delimeter choice.
    p_str, t_str = np.loadtxt(etd_url, dtype=str, encoding='unicode_escape', delimiter=', ', skiprows=2, max_rows=1)
    
    # Now get just the number, lose everything else.
    for token in p_str.split():
        try:
            p0 = float(token)
        except ValueError:
            pass # if its not the number, move on
    for token in t_str.split():
        try:
            t0 = float(token)
        except ValueError:
            pass # if its not the number, move on

    return p0, t0
wasp12_etd_url = 'http://var2.astro.cz/ETD/ascii-etd.php?id=246&STARNAME=WASP-12&PLANET=b&PER=1.0914222&EPOCH=2454508.97605'
wasp12_etd_p0, wasp12_etd_t0 = get_etd_p0_t0(wasp12_etd_url)
print("P0:", wasp12_etd_p0, "t0:",wasp12_etd_t0)<jupyter_output>P0: 1.0914222 t0: 2454508.97605<jupyter_text>This function gets the transit data for the exoplanet from the [Exoplanet Transit Database](http://var2.astro.cz/ETD).Requires RA and Dec to convert from HJD_UTC to BJD_TDB.<jupyter_code>def get_etd_data(right_ascension, declination, etd_url):
    """Imports the data from ETD. Uses RA & Dec to convert HJD_UTC to BJD_TDB.
    
    Parameters
    ----------
    right_ascension : string
        The right ascension in hh mm ss.ss
    declination : string
        The declination in dd mm ss.ss
    etd_url : string
        The URL of the ETD .csv datafile

    Returns
    -------
    data : Dataframe
        A Pandas DataFrame with bjd, mid_err, dq. These are the mid-time 
        in BJD_TDB, the uncertainty of the mid-time, and the ETD Data Quality 
        factor, where 1 is best and 5 is worst. The index is the row numbers 
        from ETD, we keep this as an arbitrary index for now.

    Notes
    -----
    It might be worth checking that each individual data point is actually in 
    HJD, rather than e.g. BJD already, as ETD doesn't specify. Check the source
    paper for the datapoint.
    """

    # Note that we keep ETD's index numbers for now, for convenience. Having
    # tmid as the index makes manipulation of it a lot harder.

    # Columns that we want from the CSV.
    csv_headers = ['#', 'HJDmid', 'HJDmid Error', 'DQ']

    # What we want to name the columns in our DataFrame.
    cols = ['tmid', 'err', 'DQ']

    # Now to import the actual transit data.
    data = pd.read_csv(etd_url, delimiter=';', header=0, index_col=0,
                       usecols=csv_headers, skiprows=4)
    
    # Set the columns to the names we want.
    data.index.rename('number', inplace=True)  # row header
    data.columns = cols  # columns

    # Convert from the truncated HJD to the full HJD.
    data['tmid'] = data['tmid'].apply(lambda x: x + 2400000)

    # Assuming (fairly certain) ETD uses HJD_UTC, let's convert to BJD_TDB.
    ra_dec_string = right_ascension+" "+declination
    ra, dec = plc.ra_dec_string_to_deg(ra_dec_string)

    # Convert to BJD_TDB; we assume uncertainty remains the same for HJD or BJD
    bjd = np.array([])
    for (col, value) in data['tmid'].iteritems():
        bjd = np.append(bjd, plc.hjd_utc_to_bjd_tdb(ra, dec, value))

    return data
etd_data = get_etd_data(wasp12_ra, wasp12_dec, wasp12_etd_url)
etd_data<jupyter_output>/usr/local/lib/python3.6/dist-packages/scipy/optimize/minpack.py:808: OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)<jupyter_text>Next, filter the ETD data. A maximum DQ (1 best, 5 worst) can be specified. It also removes mid-times that do not have a published uncertainty.<jupyter_code>def filter_etd_data(data, filter):
    """Filters ETD data to DQ<=filter. Checks uncertainty exists else deletes.

    ETD has it so that the best DQ is 1, the worst DQ is 5. By checking whether
    the transit mid-time has an associated uncertainty, we can discard the most
    egregious datapoints that are scientifically useless for us.

    Parameters
    ----------
    data : DataFrame
        The dataframe of ETD data to be filtered.
    filter : int
        The value of DQ for which transits with DQ greater than this should be 
        discarded.
    
    Returns
    -------
    filtered_data : DataFrame
        DataFrame of ETD data with DQ <= filter
    """

    # Use only values where an uncertainty exists, and isn't somehow negative 
    # (e.g. debug values of -9999 making it into the dataset!)
    filtered_data = data[data.err.notnull()]
    filtered_data = filtered_data[filtered_data.err > 0]
    
    # filter by DQ
    print("Before filtering:", len(data.index), "transits.")
    print("Filtering values with DQ filter {0}.".format(filter))
    filtered_data = filtered_data[filtered_data.DQ <= filter]
    print("After filtering:", len(filtered_data.index), "transits.")

    filtered_data.drop('DQ', axis='columns', inplace=True)

    return filtered_data
# Setting to lower filter values is better. Set to 2 as an example.
filtered_etd_data = filter_etd_data(etd_data, 2)<jupyter_output>Before filtering: 266 transits.
Filtering values with DQ filter 2.
After filtering: 103 transits.<jupyter_text>To add datapoints from a CSV, use this function. It can also convert from HJD to BJD (either UTC or TDB) if the type is specified as a third column in the data, and if parameter `convert` is set to `True`.<jupyter_code>def add_own_data(right_ascension, declination, fpath, convert=False):
    """Function to add your own data from csv, with format tmid, err, JD type.

    JD type needs to be one of "HJD_UTC", "HJD_TDB", "BJD_UTC", or "BJD_TDB". If
    necessary, these will be converted to BJD_TDB.

    Parameters
    ----------
    right_ascension : string
        The right ascension in hh mm ss.ss
    declination : string
        The declination in dd mm ss.ss
    fpath : string
        The path to your .csv datafile, local path or URL

    Returns
    -------
    data : pandas DataFrame 
        Pandas Dataframe with bjd, mid_err. These are the mid-time in BJD_TDB,
        and the uncertainty of said mid-time. The index is an arbitrary index
        for now, automatically generated.
    """

    # import all the data as NumPy arrays. We'll make a df at the end.
    # genfromtxt has a bug/"feature" at not unpacking if dtype=None (where it 
    # guesses the dtype), so we have to import all of itas strings, then change 
    # mid and err to floats afterwards.
    if convert==True:
        mid, err, jd_type = np.genfromtxt(fpath, unpack=True, dtype=str, skip_header=1, encoding='unicode_escape', delimiter=',', missing_values='')
        mid, err = mid.astype(float), err.astype(float)
        
        # Prepare for BJD conversion
        print("Mid-times will be converted to BJD_TDB.")
        bjd=np.array([])
        ra_dec_string = right_ascension+" "+declination
        ra, dec = plc.ra_dec_string_to_deg(ra_dec_string)

        # Convert the mid-time
        for i, jd in enumerate(mid):
            if jd_type[i] == "BJD_TDB":
                bjd = np.append(bjd, jd)
            elif jd_type[i] == "HJD_UTC":
                bjd = np.append(bjd, plc.hjd_utc_to_bjd_tdb(ra, dec, jd))
            elif jd_type[i] == "HJD_TDB":
                bjd = np.append(bjd, plc.hjd_tdb_to_bjd_tdb(ra, dec, jd))
            elif jd_type[i] == "BJD_UTC":
                bjd = np.append(bjd, plc.bjd_utc_to_bjd_tdb(ra, dec, jd))
            else:
                # i+2 as csv will count from 1, and +1 for the header
                raise ValueError("Type needs to be \'BJD_TDB\', \'HJD_UTC\', \'HJD_TDB\' or \`BJD_UTC\`. Check csv row",i+2)

    elif convert==False:
        bjd, err = np.genfromtxt(fpath, unpack=True, dtype=float, skip_header=1, usecols=(0,1), encoding='unicode_escape', delimiter=',', missing_values='')
        print("Mid-times assumed to already be BJD_TDB.")
    
    else:
        raise TypeError('Convert needs to be True or False')

    # Make a DataFrame for our NumPy arrays, set the columns
    own_data = pd.DataFrame(data=(np.array([bjd,err]).T))
    own_data.columns = ['tmid', 'err']
    return own_data
own_data = add_own_data(wasp12_ra, wasp12_dec, 'example_data_WASP12b.csv', convert=True)
own_data<jupyter_output>Mid-times will be converted to BJD_TDB.<jupyter_text>This function calculates the epochs for your mid-times. It makes the assumption of rounding any mid-time to the nearest epoch. It is extremely unlikely to have enough data recorded to have a situation where any TTV has caused it to "skip" an epoch or have an extra epoch, so this is a valid assumption to make.<jupyter_code>def epoch_calculator(df, p0, t0):
    """Calculates epochs for a df. Needs a P0 & T0 (or use ETD & earliest tmid). 

    Parameters
    ----------
    df : pandas DataFrame
        DataFrame of tmids and uncertainties
    p0 : float
        Period of exoplanet orbit in days
    t0 : float
        Reference mid-time for epoch 0, in BJD_TDB
    
    Returns
    -------
    sorted_df : pandas DataFrame 
        DataFrame with new epoch column, sorted by tmid.

    Raises
    ------
    ValueError
        If P0 or T0 are negative.

    Notes
    -----
    Supply a P0 and T0. If you have your own values, e.g. from a paper, pass
    these in. Else, consider using either ETD's values with get_etd_p0_t0(),
    or use the period provided by EOC, and use one of your tmid values; 
    earliest or median are good choices to see trends clearly, use 
    `df.tmid.min()` or `df.tmid.median()` to obtain.

    """
    if (p0<=0):
        raise ValueError("Period cannot be negative")

    if (t0<=0):
        raise ValueError("Mid-times cannot be negative")

    # Calculate the epoch system, assume closest rounding is correct via np.rint
    # Epoch is calculated assuming a constant-period stable circular orbit.
    df['epoch'] = data['tmid'].apply(lambda x: np.rint((x-t0)/p0) )
    df['epoch'] = pd.Series(df['epoch'], dtype=int)  # make sure it's integer

    # sort by tmid, and replace the index with sequential numbers
    sorted_df = df.sort_values('tmid', ignore_index=True)
    return sorted_df<jupyter_output><empty_output><jupyter_text>For example, the ETD and the data from the .csv will be combined and run through the function, using ETD's epoch system.<jupyter_code>data = pd.concat([filtered_etd_data, own_data])
data = epoch_calculator(data, wasp12_etd_p0, wasp12_etd_t0)
data<jupyter_output><empty_output><jupyter_text>Fitting 

In each of the three scenarios discussed earlier, the Levenberg-Marquadt algorithm is used to find the model parameters that best fit the observed data. An O-C graph of predicted data using this best-fit model, and the real observed data will be produced. Outlier datapoints will be sigma-clipped.

Note that an initial guess for each model parameter has to be specified, because LM only finds local minima, not global minima. This will be adressed in a future update.  1. Linear: Circular orbit, constant periodIf this is the case, and the exoplanet has no orbital decay, then we expect that using the first formula $t = t_0 + PE$ should have no difference from our recorded mid-times. We fit the linear model to the data using `curve_fit`, and then we generate an O-C graph of this scenario:<jupyter_code>def linear_model(E,t0,P):
    '''Formula for circular orbit with constant period.

    Parameters
    ----------
    E : int
        The epoch for this mid-time
    t0 : float
        Reference mid-time for epoch 0, in BJD_TDB
    P : float
        Orbital period in days

    Returns
    -------
    tmid : float
        the tmid for a constant period circular orbit, at that epoch
    '''
    return t0+(P*E)<jupyter_output><empty_output><jupyter_text>To get the parameters that best fit the data to the function, we use the Levenberg-Marquadt algorithm in `scipy.curve_fit`:<jupyter_code>param_guess = [wasp12_etd_t0, wasp12_etd_p0]  # Initial guess for LM to find.
pars_linear, cov_linear = curve_fit(linear_model, data.epoch, data.tmid, param_guess, method='lm')

# so this should be the t_0 and P_guess that minimise chi-sq
t0_linear, p0_linear = pars_linear

print("Original linear t0:", param_guess[0])
print("Best-fit linear t0:", t0_linear)

print("\nOriginal linear p0:", param_guess[1])
print("Best-fit linear p0:", p0_linear)
plt.figure()

# Use the formula to generate the expected tmid at each epoch we have observations at
tmid_linear = data['epoch'].apply(lambda E: linear_model(E, t0_linear, p0_linear))

# O-C between our observed tmids and these expected tmids
OC_linear = pd.Series(data=(data['tmid'] - tmid_linear))

plt.errorbar(data['epoch'], OC_linear, xerr=0, yerr=data['err'], linestyle='None', marker='.', mfc="w", mec="b")
plt.axhline(y=0, xmin=0.0, xmax=1.0, color='black', linestyle='dashed', linewidth=0.5)
plt.ylabel('Deviation of BJD$_{TDB} (days)$')
plt.xlabel('Epoch')
plt.title('O-C between observed $t_{mid}$ and predicted constant-P $t_{mid}$')
plt.show()<jupyter_output><empty_output><jupyter_text>We can see some outliers, so we will sigma-clip the O-C, and use the index of the resulting Series to clip our DataFrame.<jupyter_code># 3 Sigma-clip the OC Series
OC_linear_clipped, low, upp = sigmaclip(OC_linear, low=3.0, high=3.0)

OC_linear_clipped = pd.Series(data=OC_linear_clipped)

# Apply the index to the observed mid-times dataframe
data_linear = data.loc[OC_linear_clipped.index]

# Apply the index to the calculated mid-times
tmid_linear = tmid_linear.loc[OC_linear_clipped.index]

plt.errorbar(data_linear['epoch'], OC_linear_clipped, xerr=0, yerr=data_linear['err'], linestyle='None', marker='.', mfc="w", mec="b")
plt.axhline(y=0, xmin=0.0, xmax=1.0, color='black', linestyle='dashed', linewidth=0.5)
plt.ylabel('Deviation of BJD$_{TDB} (days)$')
plt.xlabel('Epoch')
plt.title('O-C between observed $t_{mid}$ and predicted constant-P $t_{mid}$')
plt.show()<jupyter_output><empty_output><jupyter_text>In this case, we can see that it isn't a straight line, it follows a trend, that could be descibed as possibly quadratic (orbital decay) or the beginning of a sinusoidal curve (apsidal precession).  2. Quadratic: Orbital decayIf the planet is undergoing orbital decay due to tidal interactions, then we expect that there should be no difference between the observed $t_\text{mid}$ and $t=t_0+PE+\frac{1}{2}\frac{dP}{dE}E^2$<jupyter_code>def decay_model(E,t0,P,dPdE):
    '''Formula for orbital decay due to tidal interaction. Quadratic decay term.

    Parameters
    ----------
    E : int
        The epoch for this mid-time
    t0 : float
        Reference mid-time for epoch 0, in BJD_TDB
    P : float
        Orbital period in days
    dPdE : float
        The value of dP/dE, the rate of decay. Unitless (days/days)

    Returns
    -------
    tmid : float
        the tmid for a an orbit undergoing tidal orbital decay, at that epoch
    '''
    return t0+(P*E)+(0.5*dPdE*(E**2))  #defining the equation to be modelled
param_guess = [wasp12_etd_t0, wasp12_etd_p0, wasp12_dPdE]
pars_decay, cov_decay = curve_fit(decay_model, data.epoch, data.tmid, param_guess, method='lm')

# The t_0 and P_guess that best fit decay_model to the observed mid-times
t0_decay, p0_decay, dPdE_decay = pars_decay

print("Original decay model t0:", param_guess[0])
print("Best-fit decay model decay t0:", t0_linear)

print("\nOriginal decay model p0:", param_guess[1])
print("Best-fit decay model p0:", p0_linear)

print("\nOriginal decay model dP/dE:", param_guess[2])
print("Best-fit decay model dP/dE:", dPdE_decay)
plt.figure()

# Use the formula to generate the expected tmid at each epoch we have observations at
tmid_decay = data['epoch'].apply(lambda E: decay_model(E, t0_decay, p0_decay, dPdE_decay))

# O-C between our observed tmids and these expected tmids
OC_decay = pd.Series(data=(data['tmid'] - tmid_decay))

plt.errorbar(data['epoch'], OC_decay, xerr=0, yerr=data['err'], linestyle='None', marker='.', mfc="w", mec="b")
plt.axhline(y=0, xmin=0.0, xmax=1.0, color='black', linestyle='dashed', linewidth=0.5)
plt.ylabel('Deviation of BJD$_{TDB} (days)$')
plt.xlabel('Epoch')
plt.title('O-C: observed $t_{mid}$ minus predicted decay model $t_{mid}$')
plt.show()<jupyter_output><empty_output><jupyter_text>The data seems to follow this trend much more closely (supporting the conclusion we drew earlier that the data looked more quadratic/sinusoidal, than linear). Sigma clip again to remove those outliers.<jupyter_code># 3 Sigma-clip the OC Series
OC_decay_clipped, low, upp = sigmaclip(OC_decay, low=3.0, high=3.0)

OC_decay_clipped = pd.Series(data=OC_decay_clipped)

# Apply the index to the observed mid-times dataframe
data_decay = data.loc[OC_decay_clipped.index]

# Apply the index to the calculated mid-times
tmid_decay = tmid_decay.loc[OC_decay_clipped.index]

plt.errorbar(data_decay['epoch'], OC_decay_clipped, xerr=0, yerr=data_decay['err'], linestyle='None', marker='.', mfc="w", mec="b")
plt.axhline(y=0, xmin=0.0, xmax=1.0, color='black', linestyle='dashed', linewidth=0.5)
plt.ylabel('Deviation of BJD$_{TDB} (days)$')
plt.xlabel('Epoch')
plt.title('O-C: observed $t_{mid}$ minus predicted decay-model $t_{mid}$')
plt.show()<jupyter_output><empty_output><jupyter_text>This looks more promising, but to be thorough, we should repeat the process for apsidal precession, and then quantify all three fits to see which is best, and if we can draw a conclusion with any confidence.  3. Sinusoidal: Apsidal precession<jupyter_code>def precession_model(E,t0,Ps,ecc,w0,dwdE):
    '''Formula for orbit with apsidal precession (sinusoidal term).

    Parameters
    ----------
    E : int
        The epoch for this mid-time
    t0 : float
        Reference mid-time for epoch 0, in BJD_TDB
    Ps : float
        Sideral (=orbital) period in days
    ecc: float
        Eccentricity of orbital path, unitless
    w0 : float
        The argument of periastron, in radians
    dwdE : float
        The value of dw/dE, the rate the argument of periastron is precessing,
        in radians/orbit

    Returns
    -------
    tmid : float
        the tmid for an orbit undergoing apsidal precession, at that epoch
    '''

    Pa = Ps/(1-(dwdE/(2*np.pi)))  # anomalistic period (time for one precession)
    return t0 + (Ps*E) - ( ((ecc*Pa)/np.pi) * np.cos(w0+(dwdE*E)) )<jupyter_output><empty_output><jupyter_text>Running this with the values given in the paper earlier doesn't enable `curve_fit` to find an acceptable fit. However, as an uncertainty were given, we can try use values within the uncertainty range to see if there is a set of numbers that allows a fit. This can be judged by looking at the covariance matrix and seeing whether it is actually populated (i.e. not $\infty$ or `NaN`). Later on, it can be attempted to fine-tune this to achieve the best fit/lowest $\chi^2$.<jupyter_code># Move within bound of published uncertainties
wasp12_eccentricity = 0.00208 + 0.00047
wasp12_periastron = 2.92 + 0.19
wasp12_dwdE = 0.00133 + 0.00018

param_guess = [wasp12_etd_t0, wasp12_etd_p0, wasp12_eccentricity, wasp12_periastron, wasp12_dwdE]
pars_precession, cov_precession = curve_fit(precession_model, data.epoch, data.tmid, param_guess, method='lm')

# If this is non-NaN, curve_fit has found an acceptable fit
print("Covariance matrix:")
print(repr(cov_precession))

t0_precession, p_precession, e_precession, w_precession, dwdE_precession = pars_precession

print("\nOriginal apsidal precession t0:", param_guess[0])
print("Best-fit apsidal precession t0:", t0_precession)

print("\nOriginal apsidal precesssion p0:", param_guess[1])
print("Best-fit apsidal precession pS:", p_precession)

print("\nOriginal apsidal precession eccentricity:", param_guess[2])
print("Best-fit apsidal precession:", e_precession)

print("\nOriginal apsidal precession argument of periastron:", param_guess[3])
print("Best-fit apsidal precession:", w_precession)

print("\nOriginal apsidal precession dw/dE:", param_guess[4])
print("Best-fit apsidal precession:", dwdE_precession)
plt.figure()

# Use the formula to generate the expected tmid at each epoch we have observations at
tmid_precession = data['epoch'].apply(lambda E: precession_model(E, t0_precession, p_precession, e_precession, w_precession, dwdE_precession))

# O-C between our observed tmids and these expected tmids
OC_precession = pd.Series(data=(data['tmid'] - tmid_precession))

plt.errorbar(data['epoch'], OC_precession, xerr=0, yerr=data['err'], linestyle='None', marker='.', mfc="w", mec="b")
plt.axhline(y=0, xmin=0.0, xmax=1.0, color='black', linestyle='dashed', linewidth=0.5)
plt.ylabel('Deviation of BJD$_{TDB} (days)$')
plt.xlabel('Epoch')
plt.title('O-C between observed $t_{mid}$ and predicted apsidal-precession $t_{mid}$')
plt.show()<jupyter_output><empty_output><jupyter_text>Again, sigma-clip the outliers.<jupyter_code># 3 Sigma-clip the OC Series
OC_precession_clipped, low, upp = sigmaclip(OC_precession, low=3.0, high=3.0)

OC_precession_clipped = pd.Series(data=OC_precession_clipped)

# Apply the index to the main DataFrame
data_precession = data.loc[OC_precession_clipped.index]

plt.errorbar(data_precession['epoch'], OC_precession_clipped, xerr=0, yerr=data_precession['err'], linestyle='None', marker='.', mfc="w", mec="b")
plt.axhline(y=0, xmin=0.0, xmax=1.0, color='black', linestyle='dashed', linewidth=0.5)
plt.ylabel('Deviation of BJD$_{TDB} (days)$')
plt.xlabel('Epoch')
plt.title('O-C: observed $t_{mid}$ minus predicted precession model $t_{mid}$')
plt.show()<jupyter_output><empty_output><jupyter_text>This doesn't look as promising as the decay model at first glance - a perfect fit would be a straight line along y=0, this is showing a curve - but let us quantify the fits.  Statistical Analysis  Judging models: $\chi^2$ calculationThese functions are for the $\chi^2$ statistic, which can measure how well a model $y(x_i)$ is fitted to the observed data $y_i$ with uncertainty $\alpha_i$. $$ \chi^2 = \sum_i{\frac{y_i-y(x_i)}{\alpha_i^2}}$$We divide $\chi^2$ by $\nu$, the degrees of freedom (the number of datapoints minus the number of fitted parameters) to obtain the reduced chi-squared:$$ \chi^2_\text{reduced} = \frac{\chi^2}{\nu} $$For each model, the null hypothesis is that the model is the best representation (thus the closest fit) to the observed mid-time/epochs. As the three models have different numbers of parameters, we used the reduced chi-squared. If a model is a good-fit, then $\chi^2\sim\nu$, so $\chi{^2_\nu}\sim1$.  (Models could also be deemed a good fit if $\chi^2$ is within two standard deviations $\pm2\sqrt{2\nu}$ of the mean, $\nu$. The probability of being outside this is approximately 0.04, similar to the 5% significance level: so if over 1 s.d. larger than the mean, it suggests the null hypothesis should be rejected; if smaller by more than one 1 s.d., the model might show over-fitting.)Reference: <NAME> and <NAME> (2010). *Measurements and their Uncertainties: A practical guide to modern error analysis.* Oxford University Press.<jupyter_code>def chi_squared(y_i,yx,alpha):
    """Calculates the unreduced chi squared.
    
    Parameters
    ----------
    y_i : float
        Observed y value (the y_i in the formula)
    yx : float
        y values of the fitted line (the y(x) in the formula )
    alpha : float
        Error in observed y, in the same units as y
    
    Returns
    -------
    chi2 : float
        the unreduced chi squared value
    """

    chi2=0
    for i,y in enumerate(y_i):
        chi2 += ((y-yx[i])**2) / (alpha[i]**2)
    return chi2

def reduced_chi_squared(y_i,yx,alpha,m):
    """Calculates the reduced chi squared from the raw data.
    
    Parameters
    ----------
    y_i : float
        Observed y value (the y_i in the formula)
    yx : float
        y values of the fitted line (the y(x) in the formula )
    alpha : float
        Error in observed y, in the same units as y
    m : int
        number of fitted parameters

    Returns
    -------
    reduced_chi2 : float
        Unreduced chi squared value
    
    Notes
    -----

    Just chisq divided by degrees of freedom. Degree of freedom is just number 
    of observations n minus number of fitted parameters m, where n is just the 
    number of y (or y(x) or x) values.
    """
    n = len(y_i) # number of observations/datapoints
    dof = n-m # degree of freedom, m is fitted params
    chi2 = chi_squared(y_i,yx,alpha) # get the unreduced chi squared
    reduced_chi2 = chi2/dof # reduce it
    return reduced_chi2
print("Constant period:")
print("Chi-squared", chi_squared(data_linear['tmid'], tmid_linear, data_linear['err']))
print("Reduced chi-squared", reduced_chi_squared(data_linear['tmid'], tmid_linear, data_linear['err'], len(pars_linear)))

print("\nOrbital decay:")
print("Chi-squared", chi_squared(data_decay['tmid'], tmid_decay, data_decay['err']))
print("Reduced chi-squared", reduced_chi_squared(data_decay['tmid'], tmid_decay, data_decay['err'], len(pars_decay)))

print("\nApsidal precession:")
print("Chi-squared", chi_squared(data_precession['tmid'], tmid_precession, data_precession['err']))
print("Reduced chi-squared", reduced_chi_squared(data_precession['tmid'], tmid_precession, data_precession['err'], len(pars_precession)))<jupyter_output>Constant period:
Chi-squared 1563.756832120829
Reduced chi-squared 8.317855490004408

Orbital decay:
Chi-squared 1249.4280224089266
Reduced chi-squared 6.980044817926964

Apsidal precession:
Chi-squared 2577.9444032934603
Reduced chi-squared 13.288373212852887<jupyter_text>As can be seen, orbital decay is the best (or rather, least-worst) fit. More data needs to be included over a longer time-scale. For apsidal precession, further constraints on the parameter initial guesses and uncertainties would be beneficial for Levenberg-Marquadt. A lower DQ on the ETD data may improve the quality at the sacrifice of number of datapoints.  Comparing models: AIC and BIC<jupyter_code>def AIC(chisq, k):
    return chisq + 2*k

def BIC(chisq, k, n):
    return chisq + k*(np.log10(n))<jupyter_output><empty_output><jupyter_text>$\Delta{AIC}$ and $\Delta{BIC}$ for a model is found by the difference between the lowest value of the IC for the three models minus the IC for that model. These can be used in the below functions too:<jupyter_code>def likelihood_ratio(delta_AIC):
    return np.exp(delta_AIC/2)

def bayes_factor(delta_BIC):
    """Approximate Bayes factor"""
    return np.exp(delta_BIC/2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10589835610094929, 'fraction_numerical': 0.02106969205834684, 'mean_word_length': 3.854995082197555, 'pattern_counts': {'":': 6, '<': 94, '<?xml version=': 0, '>': 90, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '262669', 'license': 'MIT', 'n_tokens_mistral': 12084, 'n_tokens_neox': 11196, 'n_words': 3822, 'path': 'TOTEMS_WASP12.ipynb', 'repo_name': 'Ross-Dobson/totems', 'size': '34598'}