Commit f032e478 authored by Laura Ketzer's avatar Laura Ketzer
Browse files

are hidden files still there?

parent 96cee622
%% Cell type:markdown id: tags:
 
With this notebook you can reproduce the results from the paper <br> **X-ray irradiation and evaporation of the four young planets around V1298 Tau - Poppenhaeger, Ketzer, Mallonn (2020)**
.
 
%% Cell type:markdown id: tags:
 
# Import
 
%% Cell type:code id: tags:
 
``` python
import sys
sys.path.append('../platypos_package/')
 
# Planet Class
from Planet_class_LoFo14_PAPER import planet_LoFo14_PAPER # this is the code with fixed step size
from Planet_class_LoFo14 import planet_LoFo14 # this is the code with variable step size
from Planet_class_Ot20_PAPER import planet_Ot20_PAPER # this is the code with fixed step size
from Planet_class_Ot20 import planet_Ot20 # this is the code with variable step size
import Planet_models_LoFo14 as plmo14
import Planet_model_Ot20 as plmoOt20
from Lx_evo_and_flux import Lx_evo, flux_at_planet_earth
 
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
import matplotlib.ticker as ticker
import os
from astropy import constants as const
from astroquery.nasa_exoplanet_archive import NasaExoplanetArchive
from PyAstronomy import pyasl
 
p = "../supplementary_files/"
# Tu et al. (2015) - model tracks
blueTu15 = pd.read_csv(p+'Lx_blue_track.csv')
redTu15 = pd.read_csv(p+'Lx_red_track.csv')
greenTu15 = pd.read_csv(p+'Lx_green_track.csv')
 
# Jackson et al. (2012) - Lx sample
jack12 = 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
```
 
%% Cell type:markdown id: tags:
 
## Present V1298 Tau parameters, $L_x$ evolutionary tracks, and planet models
First we need to define all the necessary system parameters. <br>
This includes the host star parameters, parameters to set the shape of the assumed future $L_x$ evolutionary tracks, and the planets themselves. <br>
To model the radius evolution of the planets we use the results from *Lopez & Fortney (2014)* and *Otegi et al. (2020)*.
 
%% Cell type:code id: tags:
 
``` python
# Stellar Parameters:
L_sun = const.L_sun
L_bol = 0.934 # David et al. 2019
mass_star, radius_star = 1.101, 1.345 # solar units
age_star = 23. # Myr
Lx_age = Lx_chandra = 1.3e30 # erg/s in energy band: (0.1-2.4 keV), error +- 1.4e29
# use dictionary to store star-parameters
star_V1298Tau = {'star_id': 'V1298Tau', 'mass': mass_star, 'radius': radius_star, 'age': age_star, 'L_bol': L_bol, 'Lx_age': Lx_age}
age = star_V1298Tau["age"]
 
# Lx evolutionary tracks:
# 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": age, "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_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": 16.}
 
# Create planet objects using either LoFO14 or Ot20 results:
# 'fluffy' LoFo14 planets with 5 M_earth core
planet_c = {"core_mass": 5.0, "radius": 5.59, "distance": 0.0825, "metallicity": "solarZ"}
planet_d = {"core_mass": 5.0, "radius": 6.41, "distance": 0.1083, "metallicity": "solarZ"}
planet_b = {"core_mass": 5.0, "radius": 10.27, "distance": 0.1688, "metallicity": "solarZ"}
planet_e = {"core_mass": 5.0, "radius": 8.74, "distance": 0.308, "metallicity": "solarZ"}
 
# fixed step size
pl_c_5_PAPER = planet_LoFo14_PAPER(star_dictionary=star_V1298Tau, planet_dict=planet_c)
pl_d_5_PAPER = planet_LoFo14_PAPER(star_dictionary=star_V1298Tau, planet_dict=planet_d)
pl_b_5_PAPER = planet_LoFo14_PAPER(star_dictionary=star_V1298Tau, planet_dict=planet_b)
pl_e_5_PAPER = planet_LoFo14_PAPER(star_dictionary=star_V1298Tau, planet_dict=planet_e)
 
# variable step size
pl_c_5 = planet_LoFo14(star_dictionary=star_V1298Tau, planet_dict=planet_c)
pl_d_5 = planet_LoFo14(star_dictionary=star_V1298Tau, planet_dict=planet_d)
pl_b_5 = planet_LoFo14(star_dictionary=star_V1298Tau, planet_dict=planet_b)
pl_e_5 = planet_LoFo14(star_dictionary=star_V1298Tau, planet_dict=planet_e)
 
# 'fluffy' LoFo14 planets with 10 M_earth core
planet_c = {"core_mass": 10.0, "radius": 5.59, "distance": 0.0825, "metallicity": "solarZ"}
planet_d = {"core_mass": 10.0, "radius": 6.41, "distance": 0.1083, "metallicity": "solarZ"}
planet_b = {"core_mass": 10.0, "radius": 10.27, "distance": 0.1688, "metallicity": "solarZ"}
planet_e = {"core_mass": 10.0, "radius": 8.74, "distance": 0.308, "metallicity": "solarZ"}
 
# fixed step size
pl_c_10_PAPER = planet_LoFo14_PAPER(star_dictionary=star_V1298Tau, planet_dict=planet_c)
pl_d_10_PAPER = planet_LoFo14_PAPER(star_dictionary=star_V1298Tau, planet_dict=planet_d)
pl_b_10_PAPER = planet_LoFo14_PAPER(star_dictionary=star_V1298Tau, planet_dict=planet_b)
pl_e_10_PAPER = planet_LoFo14_PAPER(star_dictionary=star_V1298Tau, planet_dict=planet_e)
 
# variable step size
pl_c_10 = planet_LoFo14(star_dictionary=star_V1298Tau, planet_dict=planet_c)
pl_d_10 = planet_LoFo14(star_dictionary=star_V1298Tau, planet_dict=planet_d)
pl_b_10 = planet_LoFo14(star_dictionary=star_V1298Tau, planet_dict=planet_b)
pl_e_10 = planet_LoFo14(star_dictionary=star_V1298Tau, planet_dict=planet_e)
 
# 'high-density' Ot20 planets
planet_c = {"radius": 5.59, "distance": 0.0825}
planet_d = {"radius": 6.41, "distance": 0.1083}
planet_b = {"radius": 10.27, "distance": 0.1688}
planet_e = {"radius": 8.74, "distance": 0.308}
 
# fixed step size
pl_c_Ot_PAPER = planet_Ot20_PAPER(star_dictionary=star_V1298Tau, planet_dict=planet_c)
pl_d_Ot_PAPER = planet_Ot20_PAPER(star_dictionary=star_V1298Tau, planet_dict=planet_d)
pl_b_Ot_PAPER = planet_Ot20_PAPER(star_dictionary=star_V1298Tau, planet_dict=planet_b)
pl_e_Ot_PAPER = planet_Ot20_PAPER(star_dictionary=star_V1298Tau, planet_dict=planet_e)
 
# variable step size
pl_c_Ot = planet_Ot20(star_dictionary=star_V1298Tau, planet_dict=planet_c)
pl_d_Ot = planet_Ot20(star_dictionary=star_V1298Tau, planet_dict=planet_d)
pl_b_Ot = planet_Ot20(star_dictionary=star_V1298Tau, planet_dict=planet_b)
pl_e_Ot = planet_Ot20(star_dictionary=star_V1298Tau, planet_dict=planet_e)
```
 
%% Cell type:markdown id: tags:
 
## Plot current V1298 Tau $L_x$ & evolutionary tracks
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(figsize=(10,6))
#ax.set_title('$L_x$ evolution for different "evolutionary tracks"')
 
# plot Tu15 tracks (for a Sun-like star!)
ax.plot(blueTu15["time"], blueTu15["Lx"], marker="None", linestyle=":", color="blue", linewidth=2.5, alpha=0.6, label="__nolabel__")#, label="fast rot. (solar model)")
ax.plot(redTu15["time"], redTu15["Lx"], marker="None", linestyle=":", color="red", linewidth=2.5, alpha=0.6, label="__nolabel__")#, label="slow rot. (solar model)")
#ax.plot(greenTu15["time"], greenTu15["Lx"], marker="None", linestyle=":", color="lime", linewidth=2.5, alpha=0.5, label="__nolabel__")#, label="interm. rot. (solar model)")
 
# plot approximated tracks
step_size, t_track_start, t_track_end = 1., age, 5000. # Myr
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="high activity track", lw=2.2)
#####
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:grey", zorder=3, lw=2.2, alpha=1., label="medium activity track")
#####
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=1, ls="-", lw=2.2)
 
# plot current X-ray luminosity of V1298 Tau as measured with Chandra & the assumed X-ray luminosity at 5 Gyr
ax.scatter(age, Lx_chandra, marker='*', c='xkcd:pale yellow', edgecolors='black', linewidths=1.1, s=500, alpha=1, zorder=4, label="__nolabel__")#, label="today"
ax.scatter(5000., Lx_5Gyr, marker='*', c='white', edgecolors='black', linewidths=1.2, s=350, zorder=4, label="__nolabel__")#, label="at 5 Gyr"
 
ax.loglog()
ax.set_xlabel("Time [Myr]", fontsize=22)
ax.set_ylabel("L$_\mathrm{x}$ [erg/s]", fontsize=22)
ax.set_xticks([10, 100, 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=6500)
ylim = ax.get_ylim()
ax.set_ylim(abs(ylim[0]), ylim[1])
ax.set_ylim(10.**27, ylim[1])
ax.tick_params(direction="in", which="both", labelsize=18)
ax.legend(loc="best", fontsize=15)
#plt.savefig("./Plots_PAPER/Activity_tracks_v1298Tau_largelabels.jpg", dpi=300)
#plt.savefig("./Plots_PAPER/Fig8_largelabels.jpg", dpi=300)
 
plt.show()
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(figsize=(12,8))
ax.set_title('$L_x$ evolution evolutionary tracks with Jackson12 sample')
 
# plot Tu15 tracks (for a Sun-like star!)
ax.plot(blueTu15["time"], blueTu15["Lx"], marker="None", linestyle=":", color="blue", linewidth=2.5, alpha=0.6, label="__nolabel__")#, label="fast rot. (solar model)")
ax.plot(redTu15["time"], redTu15["Lx"], marker="None", linestyle=":", color="red", linewidth=2.5, alpha=0.6, label="__nolabel__")#, label="slow rot. (solar model)")
ax.plot(greenTu15["time"], greenTu15["Lx"], marker="None", linestyle=":", color="lime", linewidth=2.5, alpha=0.5, label="__nolabel__")#, 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, label="cluster stars from \nJackson et al. (2012)")
 
# plot approximated tracks
step_size, t_track_start, t_track_end = 1., age, 5000. # Myr
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)
#####
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")
#####
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:
 
# Evolve LoFo14_PAPER planets (fixed step size)
 
%% Cell type:markdown id: tags:
 
### Evolve LoFo14_PAPER planets w. 10 Mcore
 
%% Cell type:code id: tags:
 
``` python
%%time
try:
p = os.getcwd()+"/LoFo14_Results_PAPER/10_Mcore/"
os.mkdir(p)
curr_path = p
print(curr_path)
except:
try:
os.mkdir(os.getcwd()+"/LoFo14_Results_PAPER/")
os.mkdir(os.getcwd()+"/LoFo14_Results_PAPER/10_Mcore/")
curr_path = os.getcwd()+"/LoFo14_Results_PAPER/10_Mcore/"
print(curr_path)
except:
print("Path exists")
curr_path = os.getcwd()+"/LoFo14_Results_PAPER/10_Mcore/"
print(curr_path)
 
folders = ["planet_c", "planet_d", "planet_b", "planet_e"]
for f in folders:
if os.path.isdir(curr_path+f) == True:
pass
else:
os.mkdir(curr_path+f)
 
# set epsilon, the step size, and how far into the future you want to evolve the planet
eps, step_size, t_final = 0.1, 0.1, 5000.
# planet c
folder = folders[0]
pl_c_10_PAPER.evolve_forward(t_final=t_final, initial_step_size=step_size, epsilon=eps, K_on="yes", beta_on="yes",
evo_track_dict=track1, path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_c_10P, M1_c_10P, R1_c_10P, Lx1_c_10 = pl_c_10_PAPER.read_results(curr_path+folder+"/")
pl_c_10_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_c_10P, M2_c_10P, R2_c_10P, Lx2_c_10 = pl_c_10_PAPER.read_results(curr_path+folder+"/")
pl_c_10_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_c_10P, M3_c_10P, R3_c_10P, Lx3_c_10 = pl_c_10_PAPER.read_results(curr_path+folder+"/")
 
# planet d
folder = folders[1]
pl_d_10_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_d_10P, M1_d_10P, R1_d_10P, Lx1_d_10P = pl_d_10_PAPER.read_results(curr_path+folder+"/")
pl_d_10_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_d_10P, M2_d_10P, R2_d_10P, Lx2_d_10P = pl_d_10_PAPER.read_results(curr_path+folder+"/")
pl_d_10_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_d_10P, M3_d_10P, R3_d_10P, Lx3_d_10P = pl_d_10_PAPER.read_results(curr_path+folder+"/")
 
# planet b
folder = folders[2]
pl_b_10_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_b_10P, M1_b_10P, R1_b_10P, Lx1_b_10P_2 = pl_b_10_PAPER.read_results(curr_path+folder+"/")
pl_b_10_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_b_10P, M2_b_10P, R2_b_10P, Lx2_b_10P_2 = pl_b_10_PAPER.read_results(curr_path+folder+"/")
pl_b_10_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_b_10P, M3_b_10P, R3_b_10P, Lx3_b_10P_2 = pl_b_10_PAPER.read_results(curr_path+folder+"/")
 
# planet e
folder = folders[3]
pl_e_10_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_e_10P, M1_e_10P, R1_e_10P, Lx1_e_10P_2 = pl_e_10_PAPER.read_results(curr_path+folder+"/")
pl_e_10_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_e_10P, M2_e_10P, R2_e_10P, Lx2_e_10P_2 = pl_e_10_PAPER.read_results(curr_path+folder+"/")
pl_e_10_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_e_10P, M3_e_10P, R3_e_10P, Lx3_e_10P_2 = pl_e_10_PAPER.read_results(curr_path+folder+"/")
```
 
%% Output
 
/media/laura/SSD2lin/PhD/work/gitlab/platypos/example_V1298Tau/LoFo14_Results_PAPER/10_Mcore/
Planet: planet_c_track_23.0_240.0_5000.0_1.3e+30_0.0_0.0.txt
Start evolving.
stepsize= 0.10000200928288694
 
../platypos_package/Planet_models_LoFo14.py:20: RuntimeWarning: invalid value encountered in double_scalars
R_env = 2.06 * (M_p)**(-0.21) * (fenv/5)**0.59 * (F_p)**0.044 * ((age/1e3)/5)**(age_exponent[metallicity]) # R_earth
 
Atmosphere has evaportated! Only bare rocky core left! STOP this madness!
192.80341176234526
192.80341176234526
Saved!
Planet: planet_c_track_23.0_23.0_5000.0_1.3e+30_0.0_0.0.txt
Start evolving.
stepsize= 0.10000200928288694
Done!
Saved!
Planet: planet_c_track_23.0_23.0_5000.0_1.3e+30_20.0_16.0.txt
Start evolving.
stepsize= 0.10000200928288694
 
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
<timed exec> in <module>()
/media/laura/SSD2lin/PhD/work/gitlab/platypos/platypos_package/Planet_class_LoFo14_PAPER.py in evolve_forward(self, t_final, initial_step_size, epsilon, K_on, beta_on, evo_track_dict, path_for_saving, planet_folder_id)
207 t, M, R, Lx = mass_planet_RK4_forward_LO14_PAPER(epsilon=epsilon, K_on="yes", beta_on="yes", planet_object=self,
208 initial_step_size=initial_step_size, t_final=t_final,
--> 209 track_dict=evo_track_dict)
210 df = pd.DataFrame({"Time": t, "Mass": M, "Radius": R, "Lx": Lx})
211 df.to_csv(path_for_saving+self.planet_id+".txt", index=None)
/media/laura/SSD2lin/PhD/work/gitlab/platypos/platypos_package/Mass_evolution_function.py in mass_planet_RK4_forward_LO14_PAPER(epsilon, K_on, beta_on, planet_object, initial_step_size, t_final, track_dict)
345 f_env_05k2 = (M_env_05k2/M_05k2)*100 # new mass fraction
346
--> 347 Mdot3 = mass_loss_rate_forward_LO14(times[i]+0.5*dt, epsilon, K_on, beta_on, planet_object, f_env_05k2, track_dict)
348 # R_p=radius_planet(M + 0.5*k2))
349 k3 = (dt*Myr_to_sec * Mdot3)/M_earth
/media/laura/SSD2lin/PhD/work/gitlab/platypos/platypos_package/Mass_loss_rate_function.py in mass_loss_rate_forward_LO14(t_, epsilon, K_on, beta_on, planet_object, f_env, track_dict)
46 beta = 1.
47 if K_on == "yes": # use K-correction estimation from Erkaev et al. 2007
---> 48 K = bk.K_fct_LO14(planet_object.distance, M_p, planet_object.mass_star , R_p)
49 elif K_on == "no":
50 K = 1.
/media/laura/SSD2lin/PhD/work/gitlab/platypos/platypos_package/Beta_K_functions.py in K_fct_LO14(a_pl, M_pl, M_star, R_pl)
52 AU = const.au.cgs.value
53 M_sun = const.M_sun.cgs.value
---> 54 M_earth = const.M_earth.cgs.value
55 R_earth = const.R_earth.cgs.value
56
~/anaconda3/lib/python3.7/site-packages/astropy/constants/constant.py in cgs(self)
208 """
209
--> 210 return self._instance_or_super('cgs')
211
212 def __array_finalize__(self, obj):
~/anaconda3/lib/python3.7/site-packages/astropy/constants/constant.py in _instance_or_super(self, key)
192 return inst
193 else:
--> 194 return getattr(super(), key)
195
196 @property
~/anaconda3/lib/python3.7/site-packages/astropy/units/quantity.py in cgs(self)
768 cgs_unit = self.unit.cgs
769 return self._new_view(self.value * cgs_unit.scale,
--> 770 cgs_unit / cgs_unit.scale)
771
772 @property
~/anaconda3/lib/python3.7/site-packages/astropy/units/quantity.py in _new_view(self, obj, unit)
570 quantity_subclass = Quantity
571 else:
--> 572 unit = Unit(unit)
573 quantity_subclass = getattr(unit, '_quantity_class', Quantity)
574 if isinstance(self, quantity_subclass):
~/anaconda3/lib/python3.7/site-packages/astropy/units/core.py in __call__(self, s, represents, format, namespace, doc, parse_strict)
1791 bases=s.unit.bases,
1792 powers=s.unit.powers,
-> 1793 _error_check=False)
1794
1795 # now decide what we really need to do; define derived Unit?
~/anaconda3/lib/python3.7/site-packages/astropy/units/core.py in __init__(self, scale, bases, powers, decompose, decompose_bases, _error_check)
2031 if power == 1:
2032 scale *= unit.scale
-> 2033 self._bases = unit.bases
2034 self._powers = unit.powers
2035 elif power == 0:
KeyboardInterrupt:
 
%% Cell type:markdown id: tags:
 
### Evolve LoFo14_PAPER planets w. 5 Mcore
 
%% Cell type:code id: tags:
 
``` python
%%time
try:
p = os.getcwd()+"/LoFo14_Results_PAPER/5_Mcore/"
os.mkdir(p)
curr_path = p
print(curr_path)
except:
try:
os.mkdir(os.getcwd()+"/LoFo14_Results_PAPER/")
os.mkdir(os.getcwd()+"/LoFo14_Results_PAPER/5_Mcore/")
curr_path = os.getcwd()+"/LoFo14_Results_PAPER/5_Mcore/"
print(curr_path)
except:
print("Path exists")
curr_path = os.getcwd()+"/LoFo14_Results_PAPER/5_Mcore/"
print(curr_path)
 
folders = ["planet_c", "planet_d", "planet_b", "planet_e"]
for f in folders:
if os.path.isdir(curr_path+f) == True:
pass
else:
os.mkdir(curr_path+f)
 
# set epsilon, the step size, and how far into the future you want to evolve the planet
eps, step_size, t_final = 0.1, 0.1, 5000.
# planet c
folder = folders[0]
pl_c_5_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_c_5P, M1_c_5P, R1_c_5P, Lx1_c_5P = pl_c_5_PAPER.read_results(curr_path+folder+"/")
pl_c_5_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_c_5P, M2_c_5P, R2_c_5P, Lx2_c_5P = pl_c_5_PAPER.read_results(curr_path+folder+"/")
pl_c_5_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_c_5P, M3_c_5P, R3_c_5P, Lx3_c_5P = pl_c_5_PAPER.read_results(curr_path+folder+"/")
 
# planet d
folder = folders[1]
pl_d_5_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_d_5P, M1_d_5P, R1_d_5P, Lx1_d_5P = pl_d_5_PAPER.read_results(curr_path+folder+"/")
pl_d_5_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_d_5P, M2_d_5P, R2_d_5P, Lx2_d_5P = pl_d_5_PAPER.read_results(curr_path+folder+"/")
pl_d_5_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_d_5P, M3_d_5P, R3_d_5P, Lx3_d_5P = pl_d_5_PAPER.read_results(curr_path+folder+"/")
 
# planet b
folder = folders[2]
pl_b_5_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_b_5P, M1_b_5P, R1_b_5P, Lx1_b_5P_2 = pl_b_5_PAPER.read_results(curr_path+folder+"/")
pl_b_5_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_b_5P, M2_b_5P, R2_b_5P, Lx2_b_5P_2 = pl_b_5_PAPER.read_results(curr_path+folder+"/")
pl_b_5_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_b_5P, M3_b_5P, R3_b_5P, Lx3_b_5P_2 = pl_b_5_PAPER.read_results(curr_path+folder+"/")
 
# planet e
folder = folders[3]
pl_e_5_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_e_5P, M1_e_5P, R1_e_5P, Lx1_e_5P_2 = pl_e_5_PAPER.read_results(curr_path+folder+"/")
pl_e_5_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_e_5P, M2_e_5P, R2_e_5P, Lx2_e_5P_2 = pl_e_5_PAPER.read_results(curr_path+folder+"/")
pl_e_5_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_e_5P, M3_e_5P, R3_e_5P, Lx3_e_5P_2 = pl_e_5_PAPER.read_results(curr_path+folder+"/")
```
 
%% Output
 
Path exists
/media/laura/SSD2lin/PhD/work/gitlab/plaml/example_V1298Tau/LoFo14_Results_PAPER/5_Mcore/
CPU times: user 293 ms, sys: 36 ms, total: 329 ms
Wall time: 330 ms
 
%% Cell type:markdown id: tags:
 
### Evolve Ot20_PAPER Planets
 
%% Cell type:code id: tags:
 
``` python
%%time
try:
p = os.getcwd()+"/Otegi_results_PAPER/"
os.mkdir(p)
curr_path = p
print(curr_path)
except:
print("Path exists")
curr_path = os.getcwd()+"/Otegi_results_PAPER/"
print(curr_path)
 
folders = ["planet_c", "planet_d", "planet_b", "planet_e"]
for f in folders:
if os.path.isdir(curr_path+f) == True:
pass
else:
os.mkdir(curr_path+f)
 
# set epsilon, the step size, and how far into the future you want to evolve the planet
eps, step_size, t_final = 0.1, 1., 5000.
# planet c
folder = folders[0]
pl_c_Ot_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_c_OtP, M1_c_OtP, R1_c_OtP, Lx1_c_OtP = pl_c_Ot_PAPER.read_results(curr_path+folder+"/")
pl_c_Ot_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_c_OtP, M2_c_OtP, R2_c_OtP, Lx2_c_OtP = pl_c_Ot_PAPER.read_results(curr_path+folder+"/")
pl_c_Ot_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_c_OtP, M3_c_OtP, R3_c_OtP, Lx3_c_OtP = pl_c_Ot_PAPER.read_results(curr_path+folder+"/")
 
# planet d
folder = folders[1]
pl_d_Ot_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_d_OtP, M1_d_OtP, R1_d_OtP, Lx1_d_OtP = pl_d_Ot_PAPER.read_results(curr_path+folder+"/")
pl_d_Ot_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_d_OtP, M2_d_OtP, R2_d_OtP, Lx2_d_OtP = pl_d_Ot_PAPER.read_results(curr_path+folder+"/")
pl_d_Ot_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_d_OtP, M3_d_OtP, R3_d_OtP, Lx3_d_OtP = pl_d_Ot_PAPER.read_results(curr_path+folder+"/")
 
# planet b
folder = folders[2]
pl_b_Ot_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_b_OtP, M1_b_OtP, R1_b_OtP, Lx1_b_OtP = pl_b_Ot_PAPER.read_results(curr_path+folder+"/")
pl_b_Ot_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_b_OtP, M2_b_OtP, R2_b_OtP, Lx2_b_OtP = pl_b_Ot_PAPER.read_results(curr_path+folder+"/")
pl_b_Ot_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_b_OtP, M3_b_OtP, R3_b_OtP, Lx3_b_OtP = pl_b_Ot_PAPER.read_results(curr_path+folder+"/")
 
# planet e
folder = folders[3]
pl_e_Ot_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_e_OtP, M1_e_OtP, R1_e_OtP, Lx1_e_OtP = pl_e_Ot_PAPER.read_results(curr_path+folder+"/")
pl_e_Ot_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_e_OtP, M2_e_OtP, R2_e_OtP, Lx2_e_OtP = pl_e_Ot_PAPER.read_results(curr_path+folder+"/")
pl_e_Ot_PAPER.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_e_OtP, M3_e_OtP, R3_e_OtP, Lx3_e_OtP = pl_e_Ot_PAPER.read_results(curr_path+folder+"/")
```
 
%% Output
 
Path exists
/media/laura/SSD2lin/PhD/work/gitlab/platypos/example_V1298Tau/Otegi_results_PAPER/
Planet: planet_c_track_23.0_240.0_5000.0_1.3e+30_0.0_0.0.txt
Start evolving.
Done!
Saved!
CPU times: user 20.6 s, sys: 7.63 ms, total: 20.6 s
Wall time: 20.6 s
 
%% Cell type:markdown id: tags:
 
# Evolve LoFo14 planets (code with variable step size)
 
%% Cell type:markdown id: tags:
 
### Evolve LoFo14 planets w. 10 Mcore
 
%% Cell type:code id: tags:
 
``` python
%%time
try:
p = os.getcwd()+"/LoFo14_Results_varstep/10_Mcore/"
os.mkdir(p)
curr_path = p
print(curr_path)
except:
try:
os.mkdir(os.getcwd()+"/LoFo14_Results_varstep/")
os.mkdir(os.getcwd()+"/LoFo14_Results_varstep/10_Mcore/")
curr_path = os.getcwd()+"/LoFo14_Results_varstep/10_Mcore/"
print(curr_path)
except:
print("Path exists")
curr_path = os.getcwd()+"/LoFo14_Results_varstep/10_Mcore/"
print(curr_path)
 
folders = ["planet_c", "planet_d", "planet_b", "planet_e"]
for f in folders:
if os.path.isdir(curr_path+f) == True:
pass
else:
os.mkdir(curr_path+f)
 
# set epsilon, the step size, and how far into the future you want to evolve the planet
eps, step_size, t_final = 0.1, 0.1, 5000.
# planet c
folder = folders[0]
pl_c_10.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_c_10, M1_c_10, R1_c_10, Lx1_c_10 = pl_c_10.read_results(curr_path+folder+"/")
pl_c_10.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_c_10, M2_c_10, R2_c_10, Lx2_c_10 = pl_c_10.read_results(curr_path+folder+"/")
pl_c_10.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_c_10, M3_c_10, R3_c_10, Lx3_c_10 = pl_c_10.read_results(curr_path+folder+"/")
 
# planet d
folder = folders[1]
pl_d_10.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_d_10, M1_d_10, R1_d_10, Lx1_d_10 = pl_d_10.read_results(curr_path+folder+"/")
pl_d_10.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_d_10, M2_d_10, R2_d_10, Lx2_d_10 = pl_d_10.read_results(curr_path+folder+"/")
pl_d_10.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_d_10, M3_d_10, R3_d_10, Lx3_d_10 = pl_d_10.read_results(curr_path+folder+"/")
 
# planet b
folder = folders[2]
pl_b_10.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_b_10, M1_b_10, R1_b_10, Lx1_b_10_2 = pl_b_10.read_results(curr_path+folder+"/")
pl_b_10.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_b_10, M2_b_10, R2_b_10, Lx2_b_10_2 = pl_b_10.read_results(curr_path+folder+"/")
pl_b_10.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_b_10, M3_b_10, R3_b_10, Lx3_b_10_2 = pl_b_10.read_results(curr_path+folder+"/")
 
# planet e
folder = folders[3]
pl_e_10.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_e_10, M1_e_10, R1_e_10, Lx1_e_10_2 = pl_e_10.read_results(curr_path+folder+"/")
pl_e_10.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_e_10, M2_e_10, R2_e_10, Lx2_e_10_2 = pl_e_10.read_results(curr_path+folder+"/")
pl_e_10.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_e_10, M3_e_10, R3_e_10, Lx3_e_10_2 = pl_e_10.read_results(curr_path+folder+"/")
```
 
%% Output
 
Path exists
/media/laura/SSD2lin/PhD/work/gitlab/plaml/example_V1298Tau/LoFo14_Results_new/10_Mcore/
CPU times: user 66.8 ms, sys: 0 ns, total: 66.8 ms
Wall time: 66.6 ms
 
%% Cell type:markdown id: tags:
 
### Evolve LoFo14 planets w. 5 Mcore
 
%% Cell type:code id: tags:
 
``` python
%%time
try:
p = os.getcwd()+"/LoFo14_Results_varstep/5_Mcore/"
os.mkdir(p)
curr_path = p
print(curr_path)
except:
try:
os.mkdir(os.getcwd()+"/LoFo14_Results_varstep/")
os.mkdir(os.getcwd()+"/LoFo14_Results_varstep/5_Mcore/")
curr_path = os.getcwd()+"/LoFo14_Results_varstep/5_Mcore/"
print(curr_path)
except:
print("Path exists")
curr_path = os.getcwd()+"/LoFo14_Results_varstep/5_Mcore/"
print(curr_path)
 
folders = ["planet_c", "planet_d", "planet_b", "planet_e"]
for f in folders:
if os.path.isdir(curr_path+f) == True:
pass
else:
os.mkdir(curr_path+f)
 
# set epsilon, the step size, and how far into the future you want to evolve the planet
eps, step_size, t_final = 0.1, 0.1, 5000.
# planet c
folder = folders[0]
pl_c_5.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_c_5, M1_c_5, R1_c_5, Lx1_c_5 = pl_c_5.read_results(curr_path+folder+"/")
pl_c_5.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_c_5, M2_c_5, R2_c_5, Lx2_c_5 = pl_c_5.read_results(curr_path+folder+"/")
pl_c_5.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_c_5, M3_c_5, R3_c_5, Lx3_c_5 = pl_c_5.read_results(curr_path+folder+"/")
 
# planet d
folder = folders[1]
pl_d_5.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_d_5, M1_d_5, R1_d_5, Lx1_d_5 = pl_d_5.read_results(curr_path+folder+"/")
pl_d_5.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_d_5, M2_d_5, R2_d_5, Lx2_d_5 = pl_d_5.read_results(curr_path+folder+"/")
pl_d_5.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_d_5, M3_d_5, R3_d_5, Lx3_d_5 = pl_d_5.read_results(curr_path+folder+"/")
 
# planet b
folder = folders[2]
pl_b_5.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_b_5, M1_b_5, R1_b_5, Lx1_b_5_2 = pl_b_5.read_results(curr_path+folder+"/")
pl_b_5.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_b_5, M2_b_5, R2_b_5, Lx2_b_5_2 = pl_b_5.read_results(curr_path+folder+"/")
pl_b_5.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_b_5, M3_b_5, R3_b_5, Lx3_b_5_2 = pl_b_5.read_results(curr_path+folder+"/")
 
# planet e
folder = folders[3]
pl_e_5.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_e_5, M1_e_5, R1_e_5, Lx1_e_5_2 = pl_e_5.read_results(curr_path+folder+"/")
pl_e_5.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_e_5, M2_e_5, R2_e_5, Lx2_e_5_2 = pl_e_5.read_results(curr_path+folder+"/")
pl_e_5.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_e_5, M3_e_5, R3_e_5, Lx3_e_5_2 = pl_e_5.read_results(curr_path+folder+"/")
```
 
%% Output
 
Path exists
/media/laura/SSD2lin/PhD/work/gitlab/plaml/example_V1298Tau/LoFo14_Results_new/5_Mcore/
CPU times: user 50.8 ms, sys: 7.93 ms, total: 58.7 ms
Wall time: 58.4 ms
 
%% Cell type:markdown id: tags:
 
### Evolve Ot20 Planets
 
%% Cell type:code id: tags:
 
``` python
%%time
try:
p = os.getcwd()+"/Otegi_results_varstep/"
os.mkdir(p)
curr_path = p
print(curr_path)
except:
print("Path exists")
curr_path = os.getcwd()+"/Otegi_results_varstep/"
print(curr_path)
 
folders = ["planet_c", "planet_d", "planet_b", "planet_e"]
for f in folders:
if os.path.isdir(curr_path+f) == True:
pass
else:
os.mkdir(curr_path+f)
 
# set epsilon, the step size, and how far into the future you want to evolve the planet
eps, step_size, t_final = 0.1, 1., 5000.
# planet c
folder = folders[0]
pl_c_Ot.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_c_Ot, M1_c_Ot, R1_c_Ot, Lx1_c_Ot = pl_c_Ot.read_results(curr_path+folder+"/")
pl_c_Ot.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_c_Ot, M2_c_Ot, R2_c_Ot, Lx2_c_Ot = pl_c_Ot.read_results(curr_path+folder+"/")
pl_c_Ot.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_c_Ot, M3_c_Ot, R3_c_Ot, Lx3_c_Ot = pl_c_Ot.read_results(curr_path+folder+"/")
 
# planet d
folder = folders[1]
pl_d_Ot.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_d_Ot, M1_d_Ot, R1_d_Ot, Lx1_d_Ot = pl_d_Ot.read_results(curr_path+folder+"/")
pl_d_Ot.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_d_Ot, M2_d_Ot, R2_d_Ot, Lx2_d_Ot = pl_d_Ot.read_results(curr_path+folder+"/")
pl_d_Ot.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_d_Ot, M3_d_Ot, R3_d_Ot, Lx3_d_Ot = pl_d_Ot.read_results(curr_path+folder+"/")
 
# planet b
folder = folders[2]
pl_b_Ot.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_b_Ot, M1_b_Ot, R1_b_Ot, Lx1_b_Ot = pl_b_Ot.read_results(curr_path+folder+"/")
pl_b_Ot.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_b_Ot, M2_b_Ot, R2_b_Ot, Lx2_b_Ot = pl_b_Ot.read_results(curr_path+folder+"/")
pl_b_Ot.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_b_Ot, M3_b_Ot, R3_b_Ot, Lx3_b_Ot = pl_b_Ot.read_results(curr_path+folder+"/")
 
# planet e
folder = folders[3]
pl_e_Ot.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track1,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t1_e_Ot, M1_e_Ot, R1_e_Ot, Lx1_e_Ot = pl_e_Ot.read_results(curr_path+folder+"/")
pl_e_Ot.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track2,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t2_e_Ot, M2_e_Ot, R2_e_Ot, Lx2_e_Ot = pl_e_Ot.read_results(curr_path+folder+"/")
pl_e_Ot.evolve_forward(t_final, step_size, eps, K_on="yes", beta_on="yes", evo_track_dict=track3,
path_for_saving=curr_path+folder+"/", planet_folder_id=folder)
t3_e_Ot, M3_e_Ot, R3_e_Ot, Lx3_e_Ot = pl_e_Ot.read_results(curr_path+folder+"/")
```
 
%% Output
 
Path exists
/media/laura/SSD2lin/PhD/work/gitlab/plaml/example_V1298Tau/Otegi_results_new/
CPU times: user 36 ms, sys: 0 ns, total: 36 ms
Wall time: 35.6 ms
 
%% Cell type:markdown id: tags:
 
# Plots from the Paper
 
%% Cell type:markdown id: tags:
 
## Plot mass & radius evolution planet c
 
%% Cell type:code id: tags:
 
``` python
fig, axs = plt.subplots(3,2, figsize=(13,10), sharex=True, sharey=False)
 
fig.suptitle("Planet c")
#ax1.set_title("Planet c")
fig.subplots_adjust(top=0.93)
 
# Otegi/heavy
axs[0, 0].plot(t1_c_OtP, (M1_c_OtP), label="high activity track", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[0, 0].plot(t2_c_OtP, (M2_c_OtP), label="medium activity track", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[0, 0].plot(t3_c_OtP, (M3_c_OtP), label="low activity track", ls="-", color="xkcd:red", lw=1.5, zorder=1)
 
# 10Mcore
axs[1, 0].plot(t1_c_10P, (M1_c_10P), label=r"fast, step = 0.1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[1, 0].plot(t2_c_10P, (M2_c_10P), label=r"med, step = 0.1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[1, 0].plot(t3_c_10P, (M3_c_10P), label=r"slow, step = 0.1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
axs[1, 0].hlines(pl_c_10_PAPER.core_mass, t1_c_10P[0], 5000., linestyle="--", color="k", lw=0.9)
 
# 5Mcore
axs[2, 0].plot(t1_c_5P, (M1_c_5P), label=r"fast, step = 0.1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[2, 0].plot(t2_c_5P, (M2_c_5P), label=r"med, step = 0.1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[2, 0].plot(t3_c_5P, (M3_c_5P), label=r"slow, step = 0.1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
axs[2, 0].hlines(pl_c_5_PAPER.core_mass, t1_c_5[0], 5000., linestyle="--", color="k", lw=0.9)
axs[2, 0].text(1600, 5.006, "core mass", fontsize=11.5)
 
#radii
# Otegi/heavy
axs[0, 1].plot(t1_c_OtP, (R1_c_OtP), label=r"fast, step = 1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[0, 1].plot(t2_c_OtP, (R2_c_OtP), label=r"med, step = 1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[0, 1].plot(t3_c_OtP, (R3_c_OtP), label=r"slow, step = 1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
 
# 10Mcore
axs[1, 1].plot(t1_c_10P, (R1_c_10P), label=r"fast, step = 0.1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[1, 1].plot(t2_c_10P, (R2_c_10P), label=r"med, step = 0.1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[1, 1].plot(t3_c_10P, (R3_c_10P), label=r"slow, step = 0.1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
age_arr = np.linspace(23., 5000., 1000)
axs[1, 1].plot(age_arr, plmo14.calculate_planet_radius(pl_c_10_PAPER.core_mass, pl_c_10_PAPER.fenv, age_arr, pl_c_10_PAPER.flux, pl_c_10_PAPER.metallicity), color="k", ls=":", lw=1)
axs[1, 1].hlines(plmo14.calculate_core_radius(pl_c_10_PAPER.core_mass), pl_c_10_PAPER.age, 5000., linestyle="--", color="k", lw=0.9)
 
dy = 1.25
axs[2, 1].text(575, 3.19+dy, "thermal contraction", fontsize=11.5, rotation=-12.5)
axs[2, 1].text(570, 2.82+dy, "without mass loss", fontsize=11.5, rotation=-12.5)
#xkcd:goldenrod
 
# 5 mcore
axs[2, 1].plot(t1_c_5P, (R1_c_5P), label=r"fast track", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[2, 1].plot(t2_c_5P, (R2_c_5P), label=r"medim track", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[2, 1].plot(t3_c_5P, (R3_c_5P), label=r"slow track", ls="-", color="xkcd:red", lw=1.5, zorder=1)
axs[2, 1].plot(age_arr, plmo14.calculate_planet_radius(pl_c_5_PAPER.core_mass, pl_c_5_PAPER.fenv, age_arr, pl_c_5_PAPER.flux, pl_c_5_PAPER.metallicity), color="k", ls=":", lw=1)
axs[2, 1].hlines(plmo14.calculate_core_radius(pl_c_5_PAPER.core_mass), pl_c_5_PAPER.age, 5000., linestyle="--", color="k", lw=0.9)
axs[2, 1].text(1475, 1.55, "core radius", fontsize=11.5)
 
axs[0, 0].legend(fontsize=10, loc=3)
for ax in [axs[0, 0], axs[1, 0], axs[2, 0]]:
ax.set_xscale("log")
ax.set_xticks([23, 100, 300, 1000, 5000])
ax.set_xlim(right= 6000)
ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.0f}'))
ax.tick_params(direction="in", which="both", labelsize=16)
ax.set_xlabel("Time [Myr]", fontsize=16, labelpad=8)
ax.set_ylabel('M [M$_\oplus$]', fontsize=16, labelpad=8)
ax.xaxis.set_major_formatter(FormatStrFormatter('%.0f')) # No decimal places
ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f')) # No decimal places
 
ylim = axs[2, 0].get_ylim()
axs[2, 0].set_ylim(top=ylim[1]+0.02)
ylim = axs[1, 0].get_ylim()
axs[1, 0].set_ylim(top=11.02)
 
for ax in [axs[0, 1], axs[1, 1], axs[2, 1]]:
ax.set_xscale("log")
ax.set_xticks([23, 100, 300, 1000, 5000])
ax.set_xlim(right= 6000)
ax.tick_params(direction="in", which="both", labelsize=16)
ax.set_xlabel("Time [Myr]", fontsize=16, labelpad=8)
ax.set_ylabel('R [R$_\oplus$]', fontsize=16, labelpad=8)
ax.xaxis.set_major_formatter(FormatStrFormatter('%.0f')) # No decimal places
ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f')) # No decimal places
 
plt.subplots_adjust(hspace=0, wspace=0.25)
fig.align_ylabels(axs[:, 0])
fig.align_ylabels(axs[:, 1])
#plt.tight_layout()
#plt.savefig("./planet_c_EVO_eps01_Zsun.jpg", dpi=300)
plt.show()
```
 
%% Output
 
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-32-c74756b68ab6> in <module>()
6
7 # Otegi/heavy
----> 8 axs[0, 0].plot(t1_c_Ot, (M1_c_Ot), label="high activity track", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
9 axs[0, 0].plot(t2_c_Ot, (M2_c_Ot), label="medium activity track", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
10 axs[0, 0].plot(t3_c_Ot, (M3_c_Ot), label="low activity track", ls="-", color="xkcd:red", lw=1.5, zorder=1)
NameError: name 't1_c_Ot' is not defined
 
 
%% Cell type:markdown id: tags:
 
### Planet c - larger labels for in-text plot
 
%% Cell type:code id: tags:
 
``` python
fig, axs = plt.subplots(3,2, figsize=(13,10), sharex=True, sharey=False)
fig.suptitle("Planet c")
fig.subplots_adjust(top=0.93)
 
# Otegi/heavy
axs[0, 0].plot(t1_c_OtP, (M1_c_OtP), label="high activity track", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[0, 0].plot(t2_c_OtP, (M2_c_OtP), label="medium activity track", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[0, 0].plot(t3_c_OtP, (M3_c_OtP), label="low activity track", ls="-", color="xkcd:red", lw=1.5, zorder=1)
 
# 10Mcore
axs[1, 0].plot(t1_c_10P, (M1_c_10P), label=r"fast, step = 0.1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[1, 0].plot(t2_c_10P, (M2_c_10P), label=r"med, step = 0.1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[1, 0].plot(t3_c_10P, (M3_c_10P), label=r"slow, step = 0.1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
axs[1, 0].hlines(pl_c_10_PAPER.core_mass, t1_c_10P[0], 5000., linestyle="--", color="k", lw=0.9)
 
# 5Mcore
axs[2, 0].plot(t1_c_5P, (M1_c_5P), label=r"high activity track", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[2, 0].plot(t2_c_5P, (M2_c_5P), label=r"medium activity track", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[2, 0].plot(t3_c_5P, (M3_c_5P), label=r"low activity track", ls="-", color="xkcd:red", lw=1.5, zorder=1)
axs[2, 0].hlines(pl_c_5_PAPER.core_mass, t1_c_5[0], 5000., linestyle="--", color="k", lw=0.9)
axs[2, 0].text(1200, 5.006, "core mass", fontsize=15)
 
#radii
# Otegi/heavy
axs[0, 1].plot(t1_c_OtP, (R1_c_OtP), label=r"fast, step = 1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[0, 1].plot(t2_c_OtP, (R2_c_OtP), label=r"med, step = 1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[0, 1].plot(t3_c_OtP, (R3_c_OtP), label=r"slow, step = 1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
 
# 10Mcore
axs[1, 1].plot(t1_c_10P, (R1_c_10P), label=r"fast, step = 0.1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[1, 1].plot(t2_c_10P, (R2_c_10P), label=r"med, step = 0.1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[1, 1].plot(t3_c_10P, (R3_c_10P), label=r"slow, step = 0.1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
age_arr = np.linspace(23., 5000., 1000)
axs[1, 1].plot(age_arr, plmo14.calculate_planet_radius(pl_c_10_PAPER.core_mass, pl_c_10_PAPER.fenv, age_arr, pl_c_10_PAPER.flux, pl_c_10_PAPER.metallicity), color="k", ls=":", lw=1)
axs[1, 1].hlines(plmo14.calculate_core_radius(pl_c_10_PAPER.core_mass), pl_c_10_PAPER.age, 5000., linestyle="--", color="k", lw=0.9)
 
dy = 1.3
axs[2, 1].text(345, 3.31+dy, "thermal contraction", fontsize=15, rotation=-12.5)
axs[2, 1].text(340, 2.9+dy, "without mass loss", fontsize=15, rotation=-12.5)
#xkcd:goldenrod
 
# 5 mcore
axs[2, 1].plot(t1_c_5P, (R1_c_5P), label=r"fast track", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[2, 1].plot(t2_c_5P, (R2_c_5P), label=r"medim track", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[2, 1].plot(t3_c_5P, (R3_c_5P), label=r"slow track", ls="-", color="xkcd:red", lw=1.5, zorder=1)
axs[2, 1].plot(age_arr, plmo14.calculate_planet_radius(pl_c_5_PAPER.core_mass, pl_c_5_PAPER.fenv, age_arr, pl_c_5_PAPER.flux, pl_c_5_PAPER.metallicity), color="k", ls=":", lw=1)
axs[2, 1].hlines(plmo14.calculate_core_radius(pl_c_5_PAPER.core_mass), pl_c_5_PAPER.age, 5000., linestyle="--", color="k", lw=0.9)
axs[2, 1].text(1075, 1.55, "core radius", fontsize=15)
 
axs[2, 0].legend(fontsize=15)#, loc=3)
for ax in [axs[0, 0], axs[1, 0], axs[2, 0]]:
ax.set_xscale("log")
ax.set_xticks([23, 100, 300, 1000, 5000])
ax.set_xlim(right= 6000)
ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.0f}'))
ax.tick_params(direction="in", which="both", labelsize=18)
ax.set_xlabel("Time [Myr]", fontsize=22, labelpad=8)
ax.set_ylabel('M [M$_\oplus$]', fontsize=22, labelpad=8)
ax.xaxis.set_major_formatter(FormatStrFormatter('%.0f')) # No decimal places
ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f')) # No decimal places
 
ylim = axs[2, 0].get_ylim()
axs[2, 0].set_ylim(top=ylim[1]+0.02)
ylim = axs[1, 0].get_ylim()
axs[1, 0].set_ylim(top=11.02)
 
for ax in [axs[0, 1], axs[1, 1], axs[2, 1]]:
ax.set_xscale("log")
ax.set_xticks([23, 100, 300, 1000, 5000])
ax.set_xlim(right= 6000)
ax.tick_params(direction="in", which="both", labelsize=18)
ax.set_xlabel("Time [Myr]", fontsize=22, labelpad=8)
ax.set_ylabel('R [R$_\oplus$]', fontsize=22, labelpad=8)
ax.xaxis.set_major_formatter(FormatStrFormatter('%.0f')) # No decimal places
ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f')) # No decimal places
 
plt.subplots_adjust(hspace=0, wspace=0.25)
fig.align_ylabels(axs[:, 0])
fig.align_ylabels(axs[:, 1])
#plt.tight_layout()
plt.savefig("./Plots_PAPER/planet_c_EVO_eps01_Zsun_largelabels.jpg", dpi=300)
plt.savefig("./Plots_PAPER/Fig9_largelabels.jpg", dpi=300)
plt.show()
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
## Plot mass & radius evolution planet d
 
%% Cell type:code id: tags:
 
``` python
fig, axs = plt.subplots(3,2, figsize=(13,10), sharex=True, sharey=False)
 
fig.suptitle("Planet d")
#ax1.set_title("Planet c")
fig.subplots_adjust(top=0.93)
 
# Otegi/heavy
axs[0, 0].plot(t1_d_OtP, (M1_d_OtP), label="high activity track", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[0, 0].plot(t2_d_OtP, (M2_d_OtP), label="medium activity track", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[0, 0].plot(t3_d_OtP, (M3_d_OtP), label="low activity track", ls="-", color="xkcd:red", lw=1.5, zorder=1)
 
# 10Mcore
axs[1, 0].plot(t1_d_10P, (M1_d_10P), label=r"fast, step = 1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[1, 0].plot(t2_d_10P, (M2_d_10P), label=r"med, step = 1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[1, 0].plot(t3_d_10P, (M3_d_10P), label=r"slow, step = 1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
axs[1, 0].hlines(pl_d_10_PAPER.core_mass, t1_d_10P[0], 5000., linestyle="--", color="k", lw=0.9)
axs[1, 0].text(25, 10.025, "core mass", fontsize=11.5)
 
# 5Mcore
axs[2, 0].plot(t1_d_5P, (M1_d_5P), label=r"fast, step = 0.1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[2, 0].plot(t2_d_5P, (M2_d_5P), label=r"med, step = 0.1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[2, 0].plot(t3_d_5P, (M3_d_5P), label=r"slow, step = 0.1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
axs[2, 0].hlines(pl_d_5_PAPER.core_mass, t1_d_5P[0], 5000., linestyle="--", color="k", lw=0.9)
 
#radii
# Otegi/heavy
axs[0, 1].plot(t1_d_OtP, (R1_d_OtP), label=r"fast, step = 1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[0, 1].plot(t2_d_OtP, (R2_d_OtP), label=r"med, step = 1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[0, 1].plot(t3_d_OtP, (R3_d_OtP), label=r"slow, step = 1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
 
# 10Mcore
axs[1, 1].plot(t1_d_10P, (R1_d_10P), label=r"fast, step = 1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[1, 1].plot(t2_d_10P, (R2_d_10P), label=r"med, step = 1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[1, 1].plot(t3_d_10P, (R3_d_10P), label=r"slow, step = 1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
age_arr = np.linspace(23., 5000., 1000)
axs[1, 1].plot(age_arr, plmo14.calculate_planet_radius(pl_d_10_PAPER.core_mass, pl_d_10_PAPER.fenv, age_arr, pl_d_10_PAPER.flux, pl_d_10_PAPER.metallicity), color="k", ls=":", lw=1)
axs[1, 1].hlines(plmo14.calculate_core_radius(pl_d_10_PAPER.core_mass), pl_d_10_PAPER.age, 5000., linestyle="--", color="k", lw=0.9)
axs[1, 1].text(25, 1.9, "core radius", fontsize=11.5)
 
# 5 mcore
axs[2, 1].plot(t1_d_5P, (R1_d_5P), label=r"fast track", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[2, 1].plot(t2_d_5P, (R2_d_5P), label=r"medim track", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[2, 1].plot(t3_d_5P, (R3_d_5P), label=r"slow track", ls="-", color="xkcd:red", lw=1.5, zorder=1)
axs[2, 1].plot(age_arr, plmo14.calculate_planet_radius(pl_d_5_PAPER.core_mass, pl_d_5_PAPER.fenv, age_arr, pl_d_5_PAPER.flux, pl_d_5_PAPER.metallicity), color="k", ls=":", lw=1)
axs[2, 1].hlines(plmo14.calculate_core_radius(pl_d_5_PAPER.core_mass), pl_d_5_PAPER.age, 5000., linestyle="--", color="k", lw=0.9)
dy = 1.82
axs[2, 1].text(575, 3.25+dy, "thermal contraction", fontsize=11.5, rotation=-13)
axs[2, 1].text(570, 2.8+dy, "without mass loss", fontsize=11.5, rotation=-13)
 
axs[0, 0].legend(fontsize=10, loc=3)
for ax in [axs[0, 0], axs[1, 0], axs[2, 0]]:
ax.set_xscale("log")
ax.set_xticks([23, 100, 300, 1000, 5000])
ax.set_xlim(right= 6000)
ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.0f}'))
ax.tick_params(direction="in", which="both", labelsize=16)
ax.set_xlabel("Time [Myr]", fontsize=16, labelpad=8)
ax.set_ylabel('M [M$_\oplus$]', fontsize=16, labelpad=8)
ax.xaxis.set_major_formatter(FormatStrFormatter('%.0f')) # No decimal places
ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f')) # No decimal places
 
ylim = axs[2, 0].get_ylim()
axs[2, 0].set_ylim(top=ylim[1]+0.02)
ylim = axs[1, 0].get_ylim()
axs[1, 0].set_ylim(top=11.55)
 
for ax in [axs[0, 1], axs[1, 1], axs[2, 1]]:
ax.set_xscale("log")
ax.set_xticks([23, 100, 300, 1000, 5000])
ax.set_xlim(right= 6000)
ax.tick_params(direction="in", which="both", labelsize=16)
ax.set_xlabel("Time [Myr]", fontsize=16, labelpad=8)
ax.set_ylabel('R [R$_\oplus$]', fontsize=16, labelpad=8)
ax.xaxis.set_major_formatter(FormatStrFormatter('%.0f')) # No decimal places
ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f')) # No decimal places
 
axs[0,1].yaxis.set_major_formatter(FormatStrFormatter('%.2f')) # No decimal places
plt.subplots_adjust(hspace=0, wspace=0.25)
fig.align_ylabels(axs[:, 0])
fig.align_ylabels(axs[:, 1])
#plt.tight_layout()
#plt.savefig("./planet_d_EVO_eps01_Zsun.jpg", dpi=300)
plt.show()
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
### Planet d - larger labels for in-text plot
 
%% Cell type:code id: tags:
 
``` python
fig, axs = plt.subplots(3,2, figsize=(13,10), sharex=True, sharey=False)
fig.suptitle("Planet d")
fig.subplots_adjust(top=0.93)
 
# Otegi/heavy
axs[0, 0].plot(t1_d_OtP, (M1_d_OtP), label="high activity track", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[0, 0].plot(t2_d_OtP, (M2_d_OtP), label="medium activity track", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[0, 0].plot(t3_d_OtP, (M3_d_OtP), label="low activity track", ls="-", color="xkcd:red", lw=1.5, zorder=1)
 
# 10Mcore
axs[1, 0].plot(t1_d_10P, (M1_d_10P), label=r"fast, step = 1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[1, 0].plot(t2_d_10P, (M2_d_10P), label=r"med, step = 1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[1, 0].plot(t3_d_10P, (M3_d_10P), label=r"slow, step = 1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
axs[1, 0].hlines(pl_d_10_PAPER.core_mass, t1_d_10P[0], 5000., linestyle="--", color="k", lw=0.9)
axs[1, 0].text(25, 10.025, "core mass", fontsize=15)
 
# 5Mcore
axs[2, 0].plot(t1_d_5P, (M1_d_5P), label=r"high activity track", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[2, 0].plot(t2_d_5P, (M2_d_5P), label=r"medium activity track", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[2, 0].plot(t3_d_5P, (M3_d_5P), label=r"low activity track", ls="-", color="xkcd:red", lw=1.5, zorder=1)
axs[2, 0].hlines(pl_d_5_PAPER.core_mass, t1_d_5P[0], 5000., linestyle="--", color="k", lw=0.9)
 
#radii
# Otegi/heavy
axs[0, 1].plot(t1_d_OtP, (R1_d_OtP), label=r"fast, step = 1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[0, 1].plot(t2_d_OtP, (R2_d_OtP), label=r"med, step = 1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[0, 1].plot(t3_d_OtP, (R3_d_OtP), label=r"slow, step = 1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
 
# 10Mcore
axs[1, 1].plot(t1_d_10P, (R1_d_10P), label=r"fast, step = 1 Myr", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[1, 1].plot(t2_d_10P, (R2_d_10P), label=r"med, step = 1 Myr", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[1, 1].plot(t3_d_10P, (R3_d_10P), label=r"slow, step = 1 Myr", ls="-", color="xkcd:red", lw=1.5, zorder=1)
age_arr = np.linspace(23., 5000., 1000)
axs[1, 1].plot(age_arr, plmo14.calculate_planet_radius(pl_d_10_PAPER.core_mass, pl_d_10_PAPER.fenv, age_arr, pl_d_10_PAPER.flux, pl_d_10_PAPER.metallicity), color="k", ls=":", lw=1)
axs[1, 1].hlines(plmo14.calculate_core_radius(pl_d_10_PAPER.core_mass), pl_d_10_PAPER.age, 5000., linestyle="--", color="k", lw=0.9)
axs[1, 1].text(25, 1.9, "core radius", fontsize=15)
 
# 5 mcore
axs[2, 1].plot(t1_d_5P, (R1_d_5P), label=r"high activity track", ls="-", color="xkcd:royal blue", lw=1, zorder=3)
axs[2, 1].plot(t2_d_5P, (R2_d_5P), label=r"medium activity track", ls="-", color="xkcd:grey", lw=2.5, zorder=2)
axs[2, 1].plot(t3_d_5P, (R3_d_5P), label=r"low activity track", ls="-", color="xkcd:red", lw=1.5, zorder=1)
axs[2, 1].plot(age_arr, plmo14.calculate_planet_radius(pl_d_5_PAPER.core_mass, pl_d_5_PAPER.fenv, age_arr, pl_d_5_PAPER.flux, pl_d_5_PAPER.metallicity), color="k", ls=":", lw=1)
axs[2, 1].hlines(plmo14.calculate_core_radius(pl_d_5_PAPER.core_mass), pl_d_5_PAPER.age, 5000., linestyle="--", color="k", lw=0.9)
 
dy = 1.88
axs[2, 1].text(345, 3.33+dy, "thermal contraction", fontsize=15, rotation=-12.5)
axs[2, 1].text(340, 2.87+dy, "without mass loss", fontsize=15, rotation=-12.5)
 
axs[2, 0].legend(fontsize=15)
for ax in [axs[0, 0], axs[1, 0], axs[2, 0]]:
ax.set_xscale("log")
ax.set_xticks([23, 100, 300, 1000, 5000])
ax.set_xlim(right= 6000)
ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.0f}'))
ax.tick_params(direction="in", which="both", labelsize=18)
ax.set_xlabel("Time [Myr]", fontsize=22, labelpad=8)
ax.set_ylabel('M [M$_\oplus$]', fontsize=22, labelpad=8)
ax.xaxis.set_major_formatter(FormatStrFormatter('%.0f')) # No decimal places
ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f')) # No decimal places
 
ylim = axs[2, 0].get_ylim()
axs[2, 0].set_ylim(top=ylim[1]+0.02)
ylim = axs[1, 0].get_ylim()
axs[1, 0].set_ylim(top=11.55)
 
for ax in [axs[0, 1], axs[1, 1], axs[2, 1]]:
ax.set_xscale("log")
ax.set_xticks([23, 100, 300, 1000, 5000])
ax.set_xlim(right= 6000)
ax.tick_params(direction="in", which="both", labelsize=18)
ax.set_xlabel("Time [Myr]", fontsize=22, labelpad=8)
ax.set_ylabel('R [R$_\oplus$]', fontsize=22, labelpad=8)
ax.xaxis.set_major_formatter(FormatStrFormatter('%.0f')) # No decimal places
ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f')) # No decimal places
 
axs[0,1].yaxis.set_major_formatter(FormatStrFormatter('%.1f')) # No decimal places
plt.subplots_adjust(hspace=0, wspace=0.25)
fig.align_ylabels(axs[:, 0])
fig.align_ylabels(axs[:, 1])
#plt.tight_layout()
plt.savefig("./Plots_PAPER/planet_d_EVO_eps01_Zsun_largelabels.jpg", dpi=300)
plt.savefig("./Plots_PAPER/Fig10_largelabels.jpg", dpi=300)
plt.show()
```
 
%% Output