trappist_machine.py

from astropy.io import fits
import pandas as pd
import glob
import os
import numpy as np

def flux_to_lum(flux, distance):
    """Distance should be in parsecs, flux in erg cm^-2 sec^-1"""
    dist_in_cm = 3.086*10**18.0 * distance # cm
    lum = np.multiply(flux, 4 * np.pi * np.power(dist_in_cm, 2.0))
    return lum

def planet_flux(luminosity, distance):
    """Distance should be in AU, luminosity in erg sec^-1"""
    dist_in_cm = (6.685*10**(-14.))**-1.0 * distance # cm /au * cm3.086*10**18.0 * distance # cm
    flux_on_planet = np.true_divide(luminosity, np.multiply(4 * np.pi, np.power(dist_in_cm, 2.0)))
    return flux_on_planet


def tidal_evo(y, t, Q, r, per, dmdtval):
    """Units: earth masses, years, AU
    """
    k = 0.3 # some constant in tidal equations (love number) - doesn't matter if we plot everything in terms of Q' in the end
    G = 0.000118413597733711048158640226628895184135977337110481586 # (*AU^3/(earth mass*year^2),
                                                                # which is G in these weird units*)

    mstar = 27000.0 # mass of trappist-1 star, in earth masses
    rstar = 5.627e-4 # radius of trappist-1 star, in AU
    # also all angles must be in radians
    a,e,m = y
    signval = np.sign(2 * (2*np.pi / 3.3) - 3*(2*np.pi / per))
    qprime_p = Q / ((2/3.) * k)
    qprime_s = 1e5 / ((2/3.) * k) # https://iopscience.iop.org/article/10.3847/1538-4357/aad40e/pdf

    dadt =  (signval * 9./2. * np.sqrt(G/ mstar) * np.power(rstar, 5.) * m * np.power(qprime_s, -1.) - 63./2. * np.multiply(qprime_p, m)**-1*np.sqrt(G * mstar**3.0) *np.power(r,5.0) * np.power(e, 2.0)) * np.power(a, -11/2.) 
    dedt =  (signval * 171./16. * np.sqrt(G/ mstar) * np.power(rstar, 5.) * m * np.power(qprime_s, -1.) - 63./4. * np.multiply(qprime_p, m)**-1*np.sqrt((G * mstar**3.0)) *(r)**5.0) * np.power(a, -13/2.) * e
    #drdt = -1e-18 / (4.0 * np.pi * r**2.) # dm/dt / (4 pi r^2) = dr/dt
    dmdt =  -1 * dmdtval
    dydt = np.asarray([dadt, dedt, dmdt])
    return dydt