Skip to content
GitLab
Commit 811557e4 authored by Laura Ketzer's avatar Laura Ketzer
Browse files

minor edits in results notebook

parent d8ba4ca3
Expand all Hide whitespace changes
Inline Side-by-side
population_evolution/Population_evolution_OwWu17__04_06-2020.ipynb
View file @ 811557e4
%% Cell type:code id: tags:
 
``` python
# last worked on 4.6.2020
# last worked on 5.6.2020
```
 
%% Cell type:markdown id: tags:
 
# Import
 
%% Cell type:code id: tags:
 
``` python
import sys
import importlib
sys.path.append('../platypos_package/')
# Planet Class
from Planet_class_LoFo14 import planet_LoFo14
importlib.reload(sys.modules['Planet_class_LoFo14'])
from Planet_class_Ot20 import planet_Ot20
from Lx_evo_and_flux import Lx_evo, flux_at_planet_earth, L_xuv_all
import keplers_3rd_law as kepler3
import notify
 
#importlib.reload(sys.modules['Mass_evolution_function'])
from Lx_evo_and_flux import undo_what_Lxuv_all_does
importlib.reload(sys.modules['Lx_evo_and_flux'])
 
 
sys.path.append('../platypos_package/population_evolution/')
 
 
 
from evolve_planet import evolve_one_planet, evolve_ensamble
importlib.reload(sys.modules['evolve_planet'])
from create_planet_chunks import create_planet_chunks
from create_summary_files import create_summary_files_with_final_planet_parameters
 
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib
matplotlib.rcParams.update({'font.size': 18, 'legend.fontsize': 14})
mpl.rcParams['axes.linewidth'] = 1.1 #set the value globally
from matplotlib.ticker import ScalarFormatter, FormatStrFormatter, FuncFormatter
import matplotlib.ticker as ticker
import os
import math
from astroquery.nasa_exoplanet_archive import NasaExoplanetArchive
from PyAstronomy import pyasl
from astropy import constants as const
# Use multiple cores for parallel processing
import multiprocessing
import multiprocessing as mp
print("Number of processors: ", mp.cpu_count())
import time
import scipy
from scipy import stats
from sklearn.neighbors import KernelDensity
import seaborn as sns
import scipy.optimize as optimize
from scipy.optimize import fsolve
 
import random
random.seed(42)
 
#################################################################################################
p = "../supplementary_files/"
# Tu et al. (2015) model tracks
blue = pd.read_csv(p+'Lx_blue_track.csv', header=None).sort_values(0)
red = pd.read_csv(p+'Lx_red_track.csv', header=None).sort_values(0)
green = pd.read_csv(p+'Lx_green_track.csv', header=None).sort_values(0).reset_index(drop=True)
# Lx sample Jackson et al. (2012)
jack = pd.read_csv(p+"Jackson2012_Lx_clean.csv")
 
def read_results_file(path, filename):
# read in results files
df = pd.read_csv(path+filename)
t, M, R, Lx = df["Time"].values, df["Mass"].values, df["Radius"].values, df["Lx"].values
return t, M, R, Lx
```
 
%% Output
 
Number of processors: 4
 
%% Cell type:markdown id: tags:
 
# Tu evolutionary tracks
 
%% Cell type:code id: tags:
 
``` python
# stellar parameters
mass_star = 1.0 # M_sun
radius_star = 1.0 # R_sun
age_star = 10.
L_bol = 1.0 # L_sun
Lx_age = 2.*10**30 # calculate_Lx_sat(L_bol) # erg/s
 
# use dictionary to store star-parameters
star = {'star_id': "dummySun", 'mass': mass_star, 'radius': radius_star, 'age': age_star,
'L_bol': L_bol, 'Lx_age': Lx_age}
age = star['age']
Lx_age = star["Lx_age"]
 
# create dictionaries with all the values necessary to create the evolutionary paths
Lx_1Gyr = 2.10*10**28 # Lx value at 1 Gyr from Tu et al. (2015) model tracks
Lx_5Gyr = 1.65*10**27 # Lx value at 5 Gyr from Tu et al. (2015) model tracks
 
track1 = {"t_start": age, "t_sat": 240., "t_curr": 1000., "t_5Gyr": 5000., "Lx_max": Lx_age, "Lx_curr": Lx_1Gyr, "Lx_5Gyr": Lx_5Gyr, "dt_drop": 0., "Lx_drop_factor": 0.}
track2 = {"t_start": age, "t_sat": 40., "t_curr": 1000., "t_5Gyr": 5000., "Lx_max": Lx_age, "Lx_curr": Lx_1Gyr, "Lx_5Gyr": Lx_5Gyr, "dt_drop": 20., "Lx_drop_factor": 7.}
track2_2 = {"t_start": age, "t_sat": 70., "t_curr": 1000., "t_5Gyr": 5000., "Lx_max": Lx_age, "Lx_curr": Lx_1Gyr, "Lx_5Gyr": Lx_5Gyr, "dt_drop": 20., "Lx_drop_factor": 5.}
track2_3 = {"t_start": age, "t_sat": 100., "t_curr": 1000., "t_5Gyr": 5000., "Lx_max": Lx_age, "Lx_curr": Lx_1Gyr, "Lx_5Gyr": Lx_5Gyr, "dt_drop": 0., "Lx_drop_factor": 0.}
track3 = {"t_start": age, "t_sat": age, "t_curr": 1000., "t_5Gyr": 5000., "Lx_max": Lx_age, "Lx_curr": Lx_1Gyr, "Lx_5Gyr": Lx_5Gyr, "dt_drop": 20., "Lx_drop_factor": 18.}
track3_1 = {"t_start": age, "t_sat": 20, "t_curr": 1000., "t_5Gyr": 5000., "Lx_max": Lx_age, "Lx_curr": Lx_1Gyr, "Lx_5Gyr": Lx_5Gyr, "dt_drop": 20., "Lx_drop_factor": 10.}
track_list = [track1, track2, track2_2, track2_3, track3, track3_1]
```
 
%% Cell type:code id: tags:
 
``` python
# PLOT
########################################################################################################
 
fig, ax = plt.subplots(figsize=(12,8))
ax.set_title('$L_x$ evolution evolutionary tracks')
 
ax.plot(blue[0], blue[1], marker="None", linestyle=":", color="blue", linewidth=2.5, alpha=0.6, label='_nolegend_')#, label="fast rot. (solar model)")
ax.plot(red[0], red[1], marker="None", linestyle=":", color="red", linewidth=2.5, alpha=0.6, label='_nolegend_')#, label="slow rot. (solar model)")
ax.plot(green[0], green[1], marker="None", linestyle=":", color="green", linewidth=2.5, alpha=0.5, label='_nolegend_')#, label="interm. rot. (solar model)")
ax.plot(jack["age"]/1e6, 10**jack["logLx_cgs"], ls="None", marker="o", color="grey", mec="k", alpha=0.3, zorder=1)
#####################################################################
step_size = 1.
t_track_start = 10.
t_track_end = 10000.
t_arr = np.arange(t_track_start, t_track_end+step_size, step_size)
Lx_arr = np.array([Lx_evo(t=t_i, track_dict=track1) for t_i in t_arr])
ax.plot(t_arr, Lx_arr, color="xkcd:royal blue", ls="-", zorder=2, label="fast activity track", lw=2.2)
#####
t_arr = np.arange(t_track_start, t_track_end+step_size, step_size)
Lx_arr = np.array([Lx_evo(t=t_i, track_dict=track2) for t_i in t_arr])
ax.plot(t_arr, Lx_arr, color="xkcd:green", zorder=3, lw=2.2, alpha=1., label="medium activity track 1")
#####
t_arr = np.arange(t_track_start, t_track_end+step_size, step_size)
Lx_arr = np.array([Lx_evo(t=t_i, track_dict=track2_2) for t_i in t_arr])
ax.plot(t_arr, Lx_arr, color="green", zorder=3, lw=2.2, alpha=1., label="medium activity track 2")
#####
t_arr = np.arange(t_track_start, t_track_end+step_size, step_size)
Lx_arr = np.array([Lx_evo(t=t_i, track_dict=track2_3) for t_i in t_arr])
ax.plot(t_arr, Lx_arr, color="xkcd:orange", zorder=1, lw=2.2, alpha=1., label="100 Myr saturation")
#####
t_arr = np.arange(t_track_start, t_track_end+step_size, step_size)
Lx_arr = np.array([Lx_evo(t=t_i, track_dict=track3_1) for t_i in t_arr])
ax.plot(t_arr, Lx_arr, color="xkcd:dark red", zorder=2, label="lower activity track", alpha=0.9, ls="-", lw=2.2)
#####
t_arr = np.arange(t_track_start, t_track_end+step_size, step_size)
Lx_arr = np.array([Lx_evo(t=t_i, track_dict=track3) for t_i in t_arr])
ax.plot(t_arr, Lx_arr, color="xkcd:red", zorder=2, label="low activity track", alpha=0.9, ls="-", lw=2.2)
#####################################################################
 
ax.loglog()
ax.set_xlabel("Time [Myr]")
ax.set_ylabel("L$_\mathrm{x}$ [erg/s]")
ax.set_xticks([5, 10, 20, 50, 100, 300, 1000, 5000])
ax.set_yticks([10**27., 10**28., 10**29., 10**30., 10**31.])
ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.0f}'))
ax.set_xlim(left=4.9, right=11000)
ylim = ax.get_ylim()
ax.set_ylim(abs(ylim[0]), ylim[1])
ax.legend(loc="best", fontsize=12)
#plt.savefig("./tracks_v1298Tau.png", dpi=300)
plt.show()
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
# Create a Planet Population
Valley detected in Kepler (also K2) sample -> choose stellar & planetary population such that it matches the observed
- orbital period distro similar to Kepler targets
 
<font color=red> NEED: </font>
- distribution of semi-major axes/orbital periods/equilibrium temperatures
- core mass distribution
- initial envelope mass fraction distribution <br>
 
 
**Lopez & Fortney 2014**
- f_env = envelope mass in % of total planet mass -> the planet models I use! Take different definition into account when creating the underlying distributions
 
**Owen & Wu 2017**
- orbital period distro similar to Kepler targets (from fitting Kepler planet sample & correcting for transit probabilities): <br>
<font color=blue> dN/dlogP = const. (P>7.6 days) & P^1.9 (P<=7.6 days) </font>
- logarithmically flat initial envelope mass fraction [f=M_env/M_c] distribution (<font color=blue>f=10^(-5) to 1 </font>) <br>
(NOTE: f<1, such that core mass dominates & self-gravity of envelope can be neglected!)
- Rayleigh distro for core masses (1 to 11 M_earth cores): <br>
<font color=blue> dN/dM_c ~ M_c*exp(-M_c^2/(2*sig_M^2), where sig_M = 3 M_earth </font>
 
-> models which fall into the gap (R=1.8 +/-0.2 R_earth) after a few Gyr (due to cooling contraction & photoevaporation) are disvavored by observations and thus flagged! <br>
**Note**: Their final models DO include many planets within the gap! They just exclude them and conclude the initial population must exist of either planets with M_pl >2 M_earth and envelopes of several %, or lower mass planets with essentially no atmospheres
 
**Ginzburg et al. 2018**
- orbital period distibution: <br>
<font color=blue> dN/dlog(T_eq) = const. (500<T_eq [K]<1000) & T_eq^(-6) (1000<T_eq [K]<2000) </font>
- correlate f_env [f=M_atm/M_c] with core mass according to theoretical planet formation models: <br>
<font color=blue> f ~ M_c^(-1/2) </font>
- broken power-law core mass distro (reason: matches high mass tail of observed M_c distro (RV mass measurements) better <br>
<font color=blue> dN/dM_c = const. (for M_c <= 5 M_earth) & M_c^(-2) (for M_c > 5 M_earth) </font>
 
%% Cell type:markdown id: tags:
 
**--------------------------------------------------------------------------------------------------------------------------------------------------------------**
 
%% Cell type:markdown id: tags:
 
# Population based on **Owen & Wu (2017)**
 
%% Cell type:markdown id: tags:
 
## Stellar Parameters & Evo Tracks
- Gaussian distribution in stellar mass with center at 1.3 M$_\oplus$ and variance $\sigma$ of 0.3 M$_\oplus$ <br>
- L$_{XUV}$ = L$_{sat}$ for t<100 Myr & L$_{sat}$*(t/t$_{sat}$)$^{-1-a_0}$ for t>=100 Myr
- $a_0=0.5$, t$_{sat}=100$ Myr, L$_{sat}=10^{-3.5}L_{\oplus}(M_{*}/M_\oplus$)
 
%% Cell type:code id: tags:
 
``` python
# get dataset from CKS (California-Kepler-Survey): https://california-planet-search.github.io/cks-website/
# column explanations: https://www.astro.caltech.edu/~howard/cks/column-definitions.txt)
df_cks_orig = pd.read_csv("../supplementary_files/cks_physical_merged.csv")
print(len(df_cks_orig))
df_cks = df_cks_orig.copy()
# sample selection based on Fulton 2017
mask_confirmed = df_cks["koi_disposition"] == "CONFIRMED" # only select confirmed planets
df_cks = df_cks[mask_confirmed]
print(len(df_cks))
# restrict sample to only the magnitudelimited portion of the larger CKS sample (Kp <14.2):
mask_magnitude = df_cks["kic_kmag"] < 14.2
df_cks = df_cks[mask_magnitude]
print(len(df_cks))
# planet-to-star radius ratio (R_pl/R_star) becomes uncertain at high impact parameters (b) due to degeneracies with limbdarkening.
# excluded KOIs with b > 0.7 to minimize the impact of grazing geometries.
mask_impactparam = df_cks["koi_impact"] <= 0.7
df_cks = df_cks[mask_impactparam]
print(len(df_cks))
# remove planets with orbital periods longer than 100 days in order to avoid domains of low completeness
# (especially for planets smaller than about 4 R_earth) and low transit probability.
mask_period = df_cks["koi_period"] <= 100.
df_cks = df_cks[mask_period]
print(len(df_cks))
# also excised planets orbiting evolved stars since they have somewhat lower detectability and less certain radii.
# implemented using an ad hoc temperature-dependent stellar radius filter:
mask_evolved = df_cks["koi_srad"] <= 10**(0.00025*(df_cks["koi_steff"]-5500)+0.20)
df_cks = df_cks[mask_evolved]
print(len(df_cks))
# also restrict sample to planets orbiting stars within the temperature range where we can extract precise stellar parameters
# from our high-resolution optical spectra (6500–4700 K).
mask_temperature = (df_cks["koi_steff"] >= 4700) & (df_cks["koi_steff"] <= 6500)
df_cks = df_cks[mask_temperature]
print(len(df_cks))
 
# drop columns which have missing stellar mass, planetary radius, semi-major axis or period
df_cks = df_cks.dropna(axis=0, how="any", subset=["koi_sma", "koi_period", "koi_smass", "koi_prad"])
df_cks.reset_index(inplace=True)
 
 
# Gaussian kernel density estimation
# actual data -> transform to log space
log3P = np.log(df_cks["koi_period"])/np.log(3)
log2R = np.log(df_cks["koi_prad"])/np.log(2)
 
# make a grid (in log space - to match the data)
n = np.arange(log3P.min(),log3P.max()+1,1)
m = np.arange(log2R.min(),log2R.max()+1,1)
nn, mm = np.mgrid[n.min():n.max():200j, m.min():m.max():200j]
 
nm_sample = np.vstack([nn.ravel(), mm.ravel()]).T # now I have a grid of points covering my radii and periods
nm_train = np.vstack([log3P, log2R]).T # this is my data in grid-form
 
# construct a Gaussian kernel density estimate of the distribution
kde = KernelDensity(bandwidth=0.25, kernel='gaussian')
kde.fit(nm_train) # fit/load real dataset
 
# score_samples() returns the log-likelihood of the samples
prob_density = np.exp(kde.score_samples(nm_sample))
# nm_sample is the range/grid over which I compute the density based on the data (2-D)
prob_density = np.reshape(prob_density, nn.shape)
```
 
%% Output
 
2025
1298
1290
1074
1017
983
921
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(figsize=(10,7))
ax.set_facecolor('xkcd:off white')
ax.plot(df_cks_orig["koi_period"], df_cks_orig["koi_prad"], ls="None", marker=".", color="k", alpha=0.1)
ax.plot(df_cks["koi_period"], df_cks["koi_prad"], ls="None", marker=".", ms=12, color="k", mec="white")
 
# plot gaussian kernel density of known exoplanet population
ax.contourf(3**nn, 2**mm, prob_density, cmap="YlOrBr", levels=15, alpha=0.5)
 
ax.set(xlabel="Period [d]", ylabel="Radius [R$_\oplus$]", xlim=(1,100), ylim=(0.3,20), title="Planets in CKS sample")
ax.set_ylim(1,10)
ax.set_xlim(0.7,100)
ax.loglog(basex=3, basey=2)
ax.set_yticks([1,2,4,6,10])
ax.set_xticks([1,3,10,30,100])
ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.0f}'))
ax.yaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.0f}'))
plt.show()
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7,6))
mask_radius = df_cks["koi_prad"] <= 20.
hist, bins, _ = ax1.hist(df_cks["koi_prad"][mask_radius], bins=18, range=(0.5, 6.5), histtype="step", density=True)
 
# histogram on log scale.
# Use non-equal bin sizes, such that they look equal on log scale.
logbins = np.logspace(np.log10(bins[0]), np.log10(bins[-1]), len(bins))
ax2.hist(df_cks["koi_prad"][mask_radius], bins=logbins, histtype="step", density=True)
 
ax2.vlines(1.8, ax2.get_ylim()[0], ax2.get_ylim()[1], color="k", linestyle="--")
ax2.set_xscale('log')
ax2.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.0f}'))
ax2.set_xticks([1,2,4,6])
ax1.set_title("Radius Histogram")
ax2.set_xlabel("Radius [R$_\oplus$]")
plt.show()
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(figsize=(10,6))
#hist, bins, _ = ax.hist(df_cks["koi_smass"][mask_radius], bins=10, histtype="step", density=True)
 
# 1D KDE
sns.distplot(df_cks["koi_smass"], hist=True, kde=True,
norm_hist = True,
color = 'darkblue',
hist_kws={'edgecolor':'black'},
kde_kws={'linewidth': 4})
 
ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.1f}'))
ax.set_title("Stellar Mass Histogram & KDE", fontsize=15)
ax.set_xlabel("Stellar Mass [M$_\odot$]")
plt.show()
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(figsize=(10,6))
#hist, bins, _ = ax.hist(df_cks["koi_smass"][mask_radius], bins=10, histtype="step", density=True)
 
# 1D KDE
sns.distplot(df_cks["koi_srad"], hist=True, kde=True,
norm_hist = True,
color = 'darkblue',
hist_kws={'edgecolor':'black'},
kde_kws={'linewidth': 4})
 
ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.1f}'))
ax.set_title("Stellar Radius Histogram & KDE", fontsize=15)
ax.set_xlabel("Stellar Radius [R$_\odot$]")
plt.show()
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
# HR-diagram
fig, ax = plt.subplots()
ax.set_yscale("log", basey=1.2)
ax.plot(df_cks_orig["koi_steff"], df_cks_orig["koi_srad"], ls="None", marker=".", color="k", alpha=0.2)
ax.plot(df_cks["koi_steff"], df_cks["koi_srad"], ls="None", marker=".")
ax.set(ylim=(0.5, 6.0), xlim=(6600, 4500))
ax.yaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.1f}'))
ax.set_yticks([0.5, 0.7, 1.0, 1.5, 2.1, 2.9, 4.2, 6.0])
plt.show()
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
# creating a distribution for the stellar mass
# info: df_cks is the
cks_mean = df_cks["koi_smass"].mean()
cks_std = df_cks["koi_smass"].std()
from scipy.stats import norm
smass_dist_cks = norm(cks_mean, cks_std) # Gaussian
 
# sample probabilities for a range of outcomes
starmass_arr = np.linspace(df_cks["koi_smass"].min(), df_cks["koi_smass"].max(), 100)
probabilities = [smass_dist_cks.pdf(value) for value in starmass_arr]
 
# Owen & Wu 2017 paper:
# Section 3.2. "We adopt host star masses similar to those in the CKS sample (Fulton et al. 2017),
# a Gaussian distribution in mass, centered at 1.3 M_sun with a variance of 0.3 M_sun."
mean_OW17, variance_OW17 = 1.3, 0.3
std_OW17 = np.sqrt(variance_OW17)
smass_dist_OW17 = norm(mean_OW17, std_OW17)
smass_dist_OW17_2 = norm(mean_OW17, variance_OW17)
 
starmass_arr_OW17 = np.linspace(0,5,100)
probabilities_OW17 = [smass_dist_OW17.pdf(value) for value in starmass_arr_OW17]
probabilities_OW17_2 = [smass_dist_OW17_2.pdf(value) for value in starmass_arr_OW17]
 
#plot
fig, ax = plt.subplots(figsize=(10,5))
ax.hist(df_cks["koi_smass"], bins=20, histtype="step", density=True, lw=2)
 
ax.plot(starmass_arr, probabilities, ls="-", lw=2, label=r"$\mu$={:.2f}, $\sigma$={:.2f}".format(cks_mean, cks_std))
ax.plot(starmass_arr_OW17, probabilities_OW17, ls="-", label=r"$\mu$={}, $\sigma$={:.2f}".format(mean_OW17, std_OW17))
ax.plot(starmass_arr_OW17, probabilities_OW17_2, ls="-", label=r"$\mu$={}, $\sigma$={}".format(mean_OW17, variance_OW17))
ylim = ax.get_ylim()
ax.vlines(cks_mean, ylim[0], ylim[1], linestyle=":", color="k", lw=1.3)
ax.vlines(mean_OW17, ylim[0], ylim[1], linestyle=":", color="k", lw=1.3)
ax.set_xlim(0,2.5)
ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.1f}'))
ax.set_title("Stellar Mass Histogram", fontsize=15)
ax.set_xlabel("Stellar Mass [M$_\odot$]")
ax.legend()
plt.show()
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
def L_high_energy(t_i, mass_star):
""" L_HE (UV to X-ray) in Owen & Wu (2017).
@param t_i - in Myr
@param mass_star - in solar masses
"""
t_sat = 100. # Myr
L_sat = 10**(-3.5)*(mass_star) * (const.L_sun.cgs).value
a_0 = 0.5
if t_i < t_sat:
L_HE = L_sat
elif t_i >= t_sat:
L_HE = L_sat*(t_i/t_sat)**(-1-a_0)
return L_HE # in erg/s
 
# make test plot
step_size, t_track_start, t_track_end = 1., 1., 3000.
fig, ax = plt.subplots()
for m in [0.7, 1.0, 1.3]:
t_arr = np.arange(t_track_start, t_track_end+step_size, step_size)
Lx_arr = np.array([L_high_energy(t_i, m) for t_i in t_arr])
ax.plot(t_arr, Lx_arr, label="{} M$_\odot$".format(m))
ax.loglog()
ax.legend()
ax.set(xlabel="Time [Myr]", ylabel="L$_{\, \mathrm{high\ energy}}$")
plt.show()
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
# draw N times out of pre-defined distribution of stellar masses which resembles the observed data
N = 1000
masses_cks = dist_cks.rvs(N)
# calculate corresonding saturation luminosity
L_HE_cks = L_high_energy(1., masses_cks)
 
# plot L_HE distro
fig, ax = plt.subplots(figsize=(10,6))
# 1D KDE
sns.distplot(L_HE_cks, hist=True, kde=True,
norm_hist = True,
color = 'darkblue',
hist_kws={'edgecolor':'black'},
kde_kws={'linewidth': 4})
 
ax.set_xticks([6*10**29, 10**30, 2*10**30])
ax.set_xscale("log")
ax.set_xlabel("L$_{HE}$ [erg/s]")
 
plt.show()
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
## Planetary Parameters
 
%% Cell type:markdown id: tags:
 
- orbital period distro similar to Kepler targets (from fitting Kepler planet sample & correcting for transit probabilities): <br>
<font color=blue> dN/dlogP = const. (P>7.6 days) & P^1.9 (P<=7.6 days) </font>
- logarithmically flat initial envelope mass fraction [f=M_env/M_c] distribution (<font color=blue>f=10^(-5) to 1 </font>) <br>
(NOTE: f<1, such that core mass dominates & self-gravity of envelope can be neglected!)
- Rayleigh distro for core masses (1 to 11 M_earth cores): <br>
<font color=blue> dN/dM_c ~ M_c*exp(-M_c^2/(2*sig_M^2), where sig_M = 3 M_earth </font>
 
-> models which fall into the gap (R=1.8 +/-0.2 R_earth) after a few Gyr (due to cooling contraction & photoevaporation) are disvavored by observations and thus flagged! <br>
**Note**: Their final models DO include many planets within the gap! They just exclude them and conclude the initial population must exist of either planets with M_pl >2 M_earth and envelopes of several %, or lower mass planets with essentially no atmospheres
 
 
%% Cell type:markdown id: tags:
 
### Semi-major axis distribution
 
%% Cell type:code id: tags:
 
``` python
N = 1000
deltaP = 0.01
period_values = np.arange(0.01, 100, deltaP)
```
 
%% Cell type:code id: tags:
 
``` python
# this is dN/dlogP
 
class logP_pdf_unnormalized(stats.rv_continuous):
def _pdf(self, P):
if P < 7.6:
return P**1.9
elif P >= 7.6:
return float(7.6**1.9)
 
rv_logP = logP_pdf_unnormalized(a=0, b=100, name='logP_pdf')
pdfs = np.array([rv_logP.pdf(p) for p in period_values])
area_under_curve_logP = scipy.integrate.trapz(pdfs, period_values)
 
class logP_pdf_normalized(stats.rv_continuous):
def _pdf(self, P):
if P < 7.6:
return (P**1.9)/area_under_curve_logP
elif P >= 7.6:
return float(7.6**1.9)/area_under_curve_logP
 
rv_logP = logP_pdf_normalized(a=0, b=100, name='logP_pdf')
pdfs = np.array([rv_logP.pdf(p) for p in period_values])
print(scipy.integrate.trapz(pdfs, dx=deltaP))
 
rvs = rv_logP.rvs(size=N)
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(period_values, pdfs)
ax.hist(rvs, bins=10**np.linspace(np.log10(np.min(period_values)), np.log10(np.max(period_values)), 50, endpoint=False), density=True, rwidth=0.9)
ax.set(xlabel="logP", ylabel="Probability density")
ax.set_xscale("log")
plt.show()
 
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(period_values, pdfs)
ax.hist(rvs, bins=10**np.linspace(np.log10(np.min(period_values)), np.log10(np.max(period_values)), 50, endpoint=False), density=True, rwidth=0.9)
ax.set(xlabel="logP", ylabel="Probability density")
plt.show()
 
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(period_values, pdfs/period_values)
#ax.hist(rvs, bins=10**np.linspace(np.log10(np.min(period_values)), np.log10(np.max(period_values)), 50, endpoint=False), density=True, rwidth=0.9)
ax.set(xlabel="Period", ylabel="Probability density")
plt.show()
```
 
%% Output
 
1.0000000000000002
 
 
 
 
%% Cell type:code id: tags:
 
``` python
# this is dN/dP (dN/dlogP*dP/dP = dN/dP*dlogP/dP = dN/dP*P)
 
class P_pdf_unnormalized(stats.rv_continuous):
def _pdf(self, P):
if P < 7.6:
return P**1.9/P
elif P >= 7.6:
return float(7.6**1.9)/P
 
rv_P = P_pdf_unnormalized(a=0, b=100, name='P_pdf')
pdfs = np.array([rv_P.pdf(p) for p in period_values])
area_under_curve_P = scipy.integrate.trapz(pdfs, period_values)
 
class P_pdf_normalized(stats.rv_continuous):
def _pdf(self, P):
if P < 7.6:
return (P**1.9)/P/area_under_curve_P
elif P >= 7.6:
return float(7.6**1.9)/P/area_under_curve_P
 
rv_P = P_pdf_normalized(a=0, b=100, name='P_pdf')
pdfs = np.array([rv_P.pdf(p) for p in period_values])
print(scipy.integrate.trapz(pdfs, dx=deltaP))
 
rvs = rv_P.rvs(size=N)
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(period_values, pdfs)
ax.hist(rvs, bins=10**np.linspace(np.log10(np.min(period_values)), np.log10(np.max(period_values)), 50, endpoint=False), density=True, rwidth=0.9)
ax.set(xlabel="Period [d]", ylabel="Probability density")
plt.show()
 
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(period_values, pdfs)
ax.hist(rvs, bins=10**np.linspace(np.log10(np.min(period_values)), np.log10(np.max(period_values)), 50, endpoint=False), density=True, rwidth=0.9)
ax.set(xlabel="Period [d]", ylabel="Probability density")
ax.set_xscale("log")
plt.show()
```
 
%% Output
 
1.0
 
 
 
%% Cell type:code id: tags:
 
``` python
# if P <= 7.6, probability density is distributed like P^1.9
 
class my_pdf_unnormalized(stats.rv_continuous):
def _pdf(self,P):
return P**1.9
 
my_rv_unnormalized = my_pdf_unnormalized(a=0, b=7.6, name='my_pdf')
area_under_unnormalized_pdf = scipy.integrate.trapz(my_rv_unnormalized.pdf(period_values), period_values)
 
class my_pdf_normalized(stats.rv_continuous):
def _pdf(self,P):
return P**1.9/area_under_unnormalized_pdf
 
my_rv_normalized = my_pdf_normalized(a=0, b=7.6, name='my_pdf')
 
# PLOT
fig, axs = plt.subplots(1,2, figsize=(14,4))
axs[0].set_title("PDF")
axs[0].plot(period_values, my_rv_normalized.pdf(period_values), label="PDF")
axs[1].set_title("CDF")
axs[1].plot(period_values, my_rv_normalized.cdf(period_values), label="CDF")
 
for ax in axs:
ax.set_xlim(0,25)
ylim = ax.get_ylim()
ax.vlines(7.6, ylim[0], ylim[1], linestyle="--", color="k")
 
plt.tight_layout()
#plt.show()
plt.close()
 
print(scipy.integrate.trapz(my_rv_normalized.pdf(period_values), period_values))
```
 
%% Output
 
1.0
 
%% Cell type:code id: tags:
 
``` python
#if P > 7.6, probability density is const. (uniformly distributed)
 
rv_uniform = stats.uniform(loc=7.6, scale=100)
 
# PLOT
fig, axs = plt.subplots(1,2, figsize=(10,4))
axs[0].set_title("PDF")
axs[0].plot(np.log10(period_values), rv_uniform.pdf(np.log10(period_values)), label="PDF")
axs[1].set_title("CDF")
axs[1].plot(np.log10(period_values), rv_uniform.cdf(np.log10(period_values)), label="CDF")
# for ax in axs:
# ax.set_xlim(0,120)
# ax.set_xlabel("Period [d]")
# ylim = ax.get_ylim()
# ax.vlines(7.6, ylim[0], ylim[1], linestyle="--", color="k")
plt.tight_layout()
#plt.show()
plt.close()
 
print(scipy.integrate.trapz(rv_uniform.pdf(period_values), dx=0.01))
```
 
%% Output
 
1.0000000000000002
 
%% Cell type:markdown id: tags:
 
### Core-mass distribution
 
%% Cell type:code id: tags:
 
``` python
from scipy.stats import rayleigh
 
rv_Mcore = rayleigh(scale=3)
# The probability density above is defined in the “standardized” form.
# To shift and/or scale the distribution use the loc and scale parameters.
# Specifically, rayleigh.pdf(x, loc, scale) is identically equivalent to rayleigh.pdf(y) / scale with y = (x - loc) / scale.
 
core_masses = np.arange(0,20,0.01)
# generate N random variables
rvs_ray = rv_Mcore.rvs(size=N)
 
# PLOT
fig, axs = plt.subplots(1,2, figsize=(10,4))
axs[0].set_title("PDF")
axs[0].plot(core_masses, rv_Mcore.pdf(core_masses), label="PDF")
axs[0].hist(rvs_ray, bins=20, density=True, rwidth=0.9)
axs[0].set_xlabel("Core mass [M$_\oplus$]")
 
axs[1].set_title("CDF")
axs[1].plot(core_masses, rv_Mcore.cdf(core_masses), label="CDF")
axs[1].set_xlabel("Core mass [M$_\oplus$]")
 
plt.tight_layout()
plt.show()
 
print(scipy.integrate.trapz(rv_Mcore.pdf(core_masses), dx=0.01))
print(rv_Mcore.stats(moments="mvsk"))
```
 
%% Output
 
 
0.9999990738451694
(array(3.75994241), array(3.86283306), array(0.63111066), array(0.2450893))
 
%% Cell type:markdown id: tags:
 
### Initial envelope mass fraction distribution
logarithmically flat initial envelope mass fraction [X=M_env/M_c] distribution; from X=0.01-0.3
 
%% Cell type:code id: tags:
 
``` python
deltaX = 0.001
X_values = np.arange(0.01, 0.3, deltaX)
```
 
%% Cell type:code id: tags:
 
``` python
class logX_pdf_unnormalized(stats.rv_continuous):
def _pdf(self, X):
return 0.1
 
rv_logX = logX_pdf_unnormalized(a=0.01, b=0.3, name='logX_pdf')
pdfs = np.array([rv_logX.pdf(x) for x in X_values])
area_under_curve_logX = scipy.integrate.trapz(pdfs, X_values)
print(area_under_curve_logX)
 
class logX_pdf_normalized(stats.rv_continuous):
def _pdf(self, X):
return 0.1/area_under_curve_logX
 
rv_logX = logX_pdf_normalized(a=0.01, b=0.3, name='logX_pdf')
pdfs = np.array([rv_logX.pdf(x) for x in X_values])
print(scipy.integrate.trapz(pdfs, dx=deltaX))
 
rvs = rv_logX.rvs(size=N)
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(X_values, pdfs)
ax.hist(rvs, bins=10**np.linspace(np.log10(np.min(X_values)), np.log10(np.max(X_values)), 50, endpoint=False), density=True, rwidth=0.9)
ax.set(xlabel="Period [d]", ylabel="Probability density")
ax.set_xscale("log")
plt.show()
 
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(X_values, pdfs)
ax.hist(rvs, bins=10**np.linspace(np.log10(np.min(X_values)), np.log10(np.max(X_values)), 50, endpoint=False), density=True, rwidth=0.9)
ax.set(xlabel="Envelope mass fraction X", ylabel="Probability density")
plt.show()
 
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(X_values, pdfs/X_values)
#ax.hist(rvs, bins=10**np.linspace(np.log10(np.min(period_values)), np.log10(np.max(period_values)), 50, endpoint=False), density=True, rwidth=0.9)
ax.set(xlabel="Envelope mass fraction X", ylabel="Probability density")
plt.show()
```
 
%% Output
 
0.028899999999999974
1.000000000000001
 
 
 
 
%% Cell type:code id: tags:
 
``` python
class X_pdf_unnormalized(stats.rv_continuous):
def _pdf(self, X):
return 0.1/X
 
rv_X = X_pdf_unnormalized(a=0.01, b=0.3, name='X_pdf')
pdfs = np.array([rv_X.pdf(x) for x in X_values])
area_under_curve_X = scipy.integrate.trapz(pdfs, X_values)
print(area_under_curve_X)
 
class X_pdf_normalized(stats.rv_continuous):
def _pdf(self, X):
return 0.1/X/area_under_curve_X
 
rv_X = X_pdf_normalized(a=0.01, b=0.3, name='X_pdf')
pdfs = np.array([rv_X.pdf(x) for x in X_values])
print(scipy.integrate.trapz(pdfs, dx=deltaX))
 
rvs = rv_X.rvs(size=N)
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(X_values, pdfs)
ax.hist(rvs, bins=10**np.linspace(np.log10(np.min(X_values)), np.log10(np.max(X_values)), 50, endpoint=False), density=True, rwidth=0.9)
ax.set(xlabel="Envelope mass fraction X", ylabel="Probability density")
plt.show()
 
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(X_values, pdfs)
ax.hist(rvs, bins=10**np.linspace(np.log10(np.min(X_values)), np.log10(np.max(X_values)), 50, endpoint=False), density=True, rwidth=0.9)
ax.set(xlabel="Envelope mass fraction X", ylabel="Probability density")
ax.set_xscale("log")
plt.show()
```
 
%% Output
 
0.33986900521952396
1.0000000000000009