Commit 96cee622 authored by Laura Ketzer's avatar Laura Ketzer
Browse files

Post .gitignore changes

parent 08833cfd
.*
!/.gitignore
,planet,Lx_init,t_final_high,t_final_med,t_final_low,M_init,fenv,M_final_high,fenv_high,M_final_med,fenv_med,M_final_low,fenv_low,R_init,R_final_high,R_final_med,R_final_low
0,c,1.44e+30,51.88,85.3,927.9,5.4,7.42,5.0,0.0,5.0,0.0,5.0,0.0,5.59,1.5,1.5,1.5
1,c,1.3e+30,54.9,95.9,1035.6,5.4,7.42,5.0,0.0,5.0,0.0,5.0,0.0,5.59,1.5,1.5,1.5
2,c,1.16e+30,58.64,110.1,1165.9,5.4,7.42,5.0,0.0,5.0,0.0,5.0,0.0,5.59,1.5,1.5,1.5
3,c,1.44e+30,176.0,5006.4,5007.1,10.92,8.45,10.0,0.0,10.1,1.03,10.34,3.26,5.59,1.78,2.4,3.0
4,c,1.3e+30,192.8,5005.8,5000.3,10.92,8.45,10.0,0.0,10.12,1.16,10.35,3.39,5.59,1.78,2.44,3.02
5,c,1.16e+30,213.8,5004.7,5002.4,10.92,8.45,10.0,0.0,10.13,1.3,10.36,3.52,5.59,1.78,2.49,3.05
6,c,1.44e+30,5006.0,5007.0,5008.0,26.39,,23.35,,25.49,,25.88,,5.59,5.17,5.47,5.52
7,c,1.3e+30,5007.0,5003.0,5007.0,26.39,,23.56,,25.52,,25.89,,5.59,5.2,5.47,5.52
8,c,1.16e+30,5007.0,5010.0,5006.0,26.39,,23.79,,25.55,,25.9,,5.59,5.23,5.48,5.52
9,d,1.44e+30,84.21,460.8,5001.9,5.6,10.67,5.0,0.0,5.0,0.0,5.03,0.5,6.41,1.5,1.5,1.95
10,d,1.3e+30,90.8,577.3,5008.8,5.6,10.67,5.0,0.0,5.0,0.0,5.03,0.59,6.41,1.5,1.5,2.0
11,d,1.16e+30,99.0,741.8,5006.5,5.6,10.67,5.0,0.0,5.0,0.0,5.03,0.68,6.41,1.5,1.5,2.05
12,d,1.44e+30,1471.5,5008.5,5003.9,11.42,12.43,10.0,0.0,10.49,4.69,10.85,7.84,6.41,1.78,3.25,3.75
13,d,1.3e+30,5008.5,5003.8,5008.2,11.42,12.43,10.01,0.08,10.52,4.92,10.87,7.99,6.41,1.91,3.29,3.77
14,d,1.16e+30,5000.6,5007.7,5001.4,11.42,12.43,10.03,0.27,10.54,5.16,10.89,8.14,6.41,2.05,3.33,3.79
15,d,1.44e+30,5008.0,5005.0,5003.0,32.76,,30.74,,32.17,,32.43,,6.41,6.16,6.34,6.37
16,d,1.3e+30,5003.0,5003.0,5003.0,32.76,,30.89,,32.2,,32.44,,6.41,6.18,6.34,6.37
17,d,1.16e+30,5003.0,5003.0,5003.0,32.76,,31.04,,32.22,,32.45,,6.41,6.19,6.34,6.37
18,b,1.44e+30,280.7,5001.5,5006.9,7.56,33.88,5.0,0.0,5.27,5.16,5.92,15.5,10.27,1.5,3.22,4.72
19,b,1.3e+30,336.2,5001.5,5007.9,7.56,33.88,5.0,0.0,5.3,5.69,5.96,16.13,10.27,1.5,3.32,4.79
20,b,1.16e+30,467.0,5001.1,5002.8,7.56,33.88,5.0,0.0,5.34,6.31,6.01,16.8,10.27,1.5,3.43,4.86
21,b,1.44e+30,5004.0,5006.0,5006.0,17.65,43.33,13.06,23.46,16.22,38.34,16.97,41.08,10.27,5.26,6.23,6.37
22,b,1.3e+30,5006.0,5001.0,5005.0,17.65,43.33,13.36,25.15,16.29,38.6,17.0,41.18,10.27,5.39,6.24,6.37
23,b,1.16e+30,5006.0,5005.0,5004.0,17.65,43.33,13.68,26.88,16.36,38.87,17.03,41.28,10.27,5.52,6.25,6.38
24,b,1.44e+30,5003.0,5003.0,5003.0,69.0,,67.6,,68.59,,68.77,,10.27,10.14,10.23,10.25
25,b,1.3e+30,5003.0,5003.0,5003.0,69.0,,67.7,,68.61,,68.77,,10.27,10.15,10.23,10.25
26,b,1.16e+30,5003.0,5003.0,5003.0,69.0,,67.81,,68.62,,68.78,,10.27,10.16,10.23,10.25
27,e,1.44e+30,5002.2,5005.4,5009.3,6.73,25.69,5.26,5.02,6.03,17.04,6.36,21.41,8.74,3.11,4.71,5.14
28,e,1.3e+30,5000.2,5006.2,5009.6,6.73,25.69,5.31,5.76,6.05,17.41,6.38,21.58,8.74,3.24,4.75,5.15
29,e,1.16e+30,5006.8,5006.8,5008.9,6.73,25.69,5.35,6.62,6.08,17.81,6.39,21.75,8.74,3.38,4.79,5.17
30,e,1.44e+30,5005.0,5003.0,5010.0,14.63,31.66,13.64,26.71,14.34,30.28,14.5,31.01,8.74,5.31,5.55,5.59
31,e,1.3e+30,5005.0,5010.0,5009.0,14.63,31.66,13.72,27.11,14.36,30.35,14.5,31.04,8.74,5.34,5.55,5.59
32,e,1.16e+30,5005.0,5006.0,5008.0,14.63,31.66,13.8,27.51,14.37,30.42,14.51,31.07,8.74,5.37,5.55,5.59
33,e,1.44e+30,5003.0,5003.0,5003.0,53.47,,53.16,,53.38,,53.42,,8.74,8.71,8.73,8.73
34,e,1.3e+30,5003.0,5003.0,5003.0,53.47,,53.18,,53.39,,53.42,,8.74,8.71,8.73,8.73
35,e,1.16e+30,5003.0,5003.0,5003.0,53.47,,53.21,,53.39,,53.42,,8.74,8.71,8.73,8.73
# beta and K functions
import astropy.units as u
from astropy import constants as const
import numpy as np
# rename to beta_fct & K_fct
def beta_fct_LO14(M_p, F_xuv, R_p):
"""Function estimates the beta parameter, which is a correction to the planetary absorption radius in XUV
-> planet appers larger than in optical; this approximation is used in Lalitha et al. 2018 and comes
originally from Salz M., Schneider P., Czesla S., Schmitt J., 2016a, Astronomy & Astrophysics, 585, L2
-> NOTE from paper: 'The atmospheric expansion can be neglected for massive hot Jupiters, but in the range
of superEarth-sized planets the expansion causes mass-loss rates that are higher by a factor of four.'
Sidenote: fct works for arrays of M_p/Fxuv or single values.
Parameters:
----------
@param M_p: (float) An float for input; current mass of the planet
@param F_xuv: (float)
@param R_p: (float)
"""
M_earth = const.M_earth.cgs.value
R_earth = const.R_earth.cgs.value
if (type(F_xuv)==float) or (type(F_xuv)==np.float64): # if F_xuv is single value
grav_pot = -const.G.cgs.value * (M_p*M_earth) / (R_p*R_earth)
log_beta = max(0.0, -0.185 * np.log10(-grav_pot) + 0.021 * np.log10(F_xuv) + 2.42)
beta = 10**log_beta
#print("beta = ", beta)
return beta
elif len(F_xuv) > 1: # if F_xuv is array
betas = []
for i in range(len(F_xuv)):
grav_pot_i = -const.G.cgs.value * (M_p[i]*M_earth) / (R_p[i]*R_earth)
log_beta_i = max(0.0, -0.185 * np.log10(-grav_pot_i) + 0.021 * np.log10(F_xuv[i]) + 2.42)
beta_i = 10**log_beta_i
betas.append(beta_i)
betas = np.array(betas)
return betas
def K_fct_LO14(a_pl, M_pl, M_star, R_pl):
"""K correction factor to take into account tidal influences of the host star on the planetary atmosphere (from Erkaev et al. 2007)
Sidenote: fct works for arrays of M_p/Fxuv or single values.
Parameters:
----------
a_pl:
M_pl:
M_star:
R_pl:
"""
AU = const.au.cgs.value
M_sun = const.M_sun.cgs.value
M_earth = const.M_earth.cgs.value
R_earth = const.R_earth.cgs.value
#print("M_pl = ", M_pl)
if (type(M_pl)==float) or (type(M_pl)==np.float64): # if M_pl is single value
R_roche = (a_pl*AU) * ((M_pl*M_earth)/(3*(M_star*M_sun)))**(1./3)
K = 1. - 3.*(R_pl*R_earth)/(2.*R_roche) + (R_pl*R_earth)**3./(2.*R_roche**3.)
#print("K = ", K)
return K
elif len(M_pl) > 1: # if F_xuv is array
Ks = []
for i in range(len(M_pl)):
R_roche_i = (a_pl*AU)*((M_pl[i]*M_earth)/(3*(M_star*M_sun)))**(1./3)
K_i = 1. - 3.*(R_pl[i]*R_earth)/(2.*R_roche_i) + (R_pl[i]*R_earth)**3./(2.*R_roche_i**3.)
Ks.append(K_i)
Ks = np.array(Ks)
return Ks
# this file contains functions related to Lx:
#
# test
#
####################################################################################################################
import numpy as np
import astropy.units as u
from astropy import constants as const
# re-write this function to pass a dictionary
def Lx_evo(t, track_dict):
"""
Function to calculate the X-ray luminosity, Lx, of the host star at a given time t.
NOTE: I tried to set up this function such that it works going
forward (for young system) and backwards (for older system) in time;
-> so t_curr(ent) is a little confusing! -> 2 we have two cases!!
Parameters Case 1 - a young system which you want to evolve forward in time
(i.e. estimate the mass loss the planet will undergo in the future):
---------------------------------------------------------------------------
@param t_start: (float)
Starting time for the forward simulation in Myr (e.g. 23 Myr for V1298 Tau)
@param Lx_max: (float)
Measured Lx value at the system age (i.e. t_start)
@param t_curr: (1-d numpy array or integer)
Set to 1 Gyr (with Lx_curr = Lx at 1 Gyr), which is where, by design, the three
activity tracks converge (based on Tu et al.(2015)/ Johnstone et al.(2015) models);
@param Lx_curr
@param t_5Gyr: (float)
Defines, together with Lx_5Gyr, the slope common for all three
activity evolution tracks past 1 Gyr (also based on Tu et al.(2015));
@param Lx_5Gyr: (float)
Defines, together with t_5Gyr, the slope common for all three
activity evolution tracks past 1 Gyr (also based on Tu et al.(2015));
Parameters Case 2 - an older system, which you want to evolve backwards in time
(i.e. estimate the mass loss the planet has undergone until now):
---------------------------------------------------------------------------
@param t_start is the youngest age you want to calculate backwards
(i.e. this shouldbe sth. like ~10 Myr, or close to the disk clearing timescale);
@param Lx_max: (float) Some saturation X-ray luminosity (e.g. Lx,sat/Lbol ~ 10^(-3)),
@param t_curr: (float) Needs to be set to the system age (e.g. t_curr=700 Myr for K2-136b);
@param Lx_curr: (float) If a measured Lx value is available for the system, otherwise need to estimate;
@param t_5Gyr: (float)
Defines, together with Lx_5Gyr, the slope common for all three
activity evolution tracks past 1 Gyr (also based on Tu et al.(2015));
@param Lx_5Gyr: (float)
Defines, together with t_5Gyr, the slope common for all three
activity evolution tracks past 1 Gyr (also based on Tu et al.(2015));
NOTE: t_5Gyr (and Lx_5Gyr) could be changed to create whatever slope desired (for a future evolution).
Sidenote: this function is not very pretty, needs improvements...
@return Lx_at_time_t: (float)
Pass one t-value, and get the corresponding Lx
"""
# read in the parameters from the provided dictionary -> these are all the parameters required to define the tracks
t_start, t_sat, t_curr, t_5Gyr = track_dict["t_start"], track_dict["t_sat"], track_dict["t_curr"], track_dict["t_5Gyr"]
Lx_max, Lx_curr, Lx_5Gyr = track_dict["Lx_max"], track_dict["Lx_curr"], track_dict["Lx_5Gyr"]
dt_drop, Lx_drop_factor = track_dict["dt_drop"], track_dict["Lx_drop_factor"]
########################################################################################################################
# Define function for calculating a power law
powerlaw = lambda x, k, index: k * (x**index)
# t_curr has to be greater than t_start
if (t_curr < t_start) or (t < t_start):
if t <= t_start: # this allows me (for the V1298 Tau system) to produce a saturation regime before 23 Myr with Lx_max
Lx = Lx_max
else:
print("Make sure t_curr > t_start, t >= t_start")
elif t > t_curr:
# this is the regime past 1 Myr
# I use the Tu 2015 slope which is approx the same for all three evolutionary tracks (and based on current solar value);
# this is mainly for consistency since we approximate our tracks based on the Tu tracks
alpha = (np.log10(Lx_5Gyr/Lx_curr))/(np.log10(t_5Gyr/t_curr))
k = 10**(np.log10(Lx_curr) - alpha*np.log10(t_curr))
Lx = powerlaw(t, k, alpha)
return Lx
elif t_curr > t_start:
###################################################################################
# dt_drop==0 means we create a track with only the saturation regime and
# a drop to the converging Lx at 1 Gyr
###################################################################################
if dt_drop==0: # then t_sat == t_drop
if t_start >= t_sat: # then t_sat <= t_start < t_curr
#print("t_sat <= t_start < t_curr")
if t >= t_sat: #and t <= t_curr:
# only data in power law regime
alpha = (np.log10(Lx_curr/Lx_max))/(np.log10(t_curr/t_sat))
k = 10**(np.log10(Lx_max) - alpha*np.log10(t_sat))
Lx = powerlaw(t, k, alpha)
#print(Lx)
elif t_start < t_sat: # then t_start < t_sat < t_curr
if t > t_sat: #and t <= t_curr:
alpha = (np.log10(Lx_curr/Lx_max))/(np.log10(t_curr/t_sat))
k = 10**(np.log10(Lx_max) - alpha*np.log10(t_sat))
Lx = powerlaw(t, k, alpha)
elif t <= t_sat:
Lx = Lx_max
elif dt_drop > 0: # then t_sat != t_drop\
###################################################################################
# dt_drop>0 means we create a more fancy track with the saturation regime and then
# a drop with slope change to the converging Lx at 1 Gyr
###################################################################################
t_drop = t_sat + dt_drop # t_sat < t_drop
if t <= t_sat:
Lx = Lx_max
elif (t > t_sat) and (t <= t_drop): # first of the two slopes
alpha_drop1 = (np.log10((Lx_max/Lx_drop_factor)/Lx_max))/(np.log10(t_drop/t_sat))
k_drop1 = 10**(np.log10(Lx_max) - alpha_drop1*np.log10(t_sat))
Lx = powerlaw(t, k_drop1, alpha_drop1)
elif t > t_drop:
alpha_drop2 = (np.log10(Lx_curr/(Lx_max/Lx_drop_factor)))/(np.log10(t_curr/t_drop))
k_drop2 = 10**(np.log10((Lx_max/Lx_drop_factor)) - alpha_drop2*np.log10(t_drop))
Lx = powerlaw(t, k_drop2, alpha_drop2)
return Lx
#####################################################################################################################
# define flux
def flux_at_planet_earth(L, a_p):
"""Function calculates the flux that the planet recieves at a given distance normalized with Earth's flux
-> which I approximate with the semi-major axis
(i.e. I ignore the eccentricity of the planetary orbits)"""
A = (4.*np.pi*(a_p*const.au.cgs)**2)
F = (L*u.erg/u.s)/A
flux_earth = 1373*1e7*(u.erg/u.s)/(100*u.cm*100*u.cm)
F_earth = (F/flux_earth).value
return F_earth
#####################################################################################################################
# define flux
def flux_at_planet(L, a_p):
"""Function calculates the flux that the planet recieves at a given distance -> which I approximate with the semi-major axis
(i.e. I ignore the eccentricity of the planetary orbits)"""
A = (4.*np.pi*(a_p*const.au.cgs)**2)
F = ((L*u.erg/u.s)/A).value
return F
#####################################################################################################################
def L_xuv_all(Lx):
"""Function to estimate the EUV luminosity using the scaling relations given by Sanz-Forcada et al. (2011);
the measuredradius_planet_volatilety can be extrapolated to the total high-energy radiation as given below.
To get the combined high-energy flux, we add L_x and L_euv to get L_xuv.
THEY USE: Lx(0.1-2.4 keV) ROSAT band
(LX and LEUV are the X-ray (λλ 5−100 Å) and EUV (λλ 100−920 Å) luminosities, in erg s−1.)
"""
log_L_euv = 0.860*np.log10(Lx) + 4.80
Leuv = 10.**log_L_euv
log_L_xuv = np.log10(Lx+Leuv)
return 10.**log_L_xuv#*u.erg/u.s
#####################################################################################################################
#####################################################################################################################
def Lx_relation_Booth(t, R_star):
"""R_star in terms of solar units"""
log_Lx = 54.65 - 2.8*np.log10(t) + 2*np.log10(R_star)
return 10.**log_Lx
#####################################################################################################################
def calculate_Lx_sat(L_star_bol):
"""*NOTE*: it is not clear yet, if L_x is directly related to L_bol
-> or rather that the size of star (radius) is correlated with the
number of active regions, i.e. bigger surface, bigger L_bol, brighter L_x"""
return 10**(-3)*(L_star_bol*const.L_sun.cgs).value
\ No newline at end of file
This diff is collapsed.
# mass loss rate function
from Lx_evo_and_flux import Lx_evo, flux_at_planet_earth, L_xuv_all, flux_at_planet
import Planet_models_LoFo14 as plmo14 #rename plmoLoFo14
import Planet_model_Ot20 as plmoOt20
import Beta_K_functions as bk
import astropy.units as u
from astropy import constants as const
def mass_loss_rate_forward_LO14(t_, epsilon, K_on, beta_on, planet_object, f_env, track_dict):
"""I use this to calculate the updated mass-loss rates at any given time step of my integration.
Input is a given time t_; this fct. calculates the corresponding Lx(t_) (then Lxuv & Fxuv given a_p),
then uses the radius at t_ to estimate a mass [based on a given mass-radius relation, which can be specified ->
e.g. mass_planet=mass_planet_CK and radius_planet=radius_planet_CK];
then we have a density, can calculate the beta and K parameter (if beta_on and K_on are set to "yes",
otherwise at each time step they beta=1 and/or K=1);
finally put it all into the mass-loss rate equation to get Mdot at t_
Parameters:
----------
t_:
Lx_evo:
calculate_planet_radius:
epsilon:
K_on:
beta_on:
planet_object:
f_env: f_env at a given time
"""
Lx = Lx_evo(t=t_, track_dict=track_dict)
Fxuv = flux_at_planet(L_xuv_all(Lx), planet_object.distance) # get flux at orbital separation a_p
# get planet density
R_p = plmo14.calculate_planet_radius(planet_object.core_mass, f_env, t_, flux_at_planet_earth(planet_object.Lbol,
planet_object.distance),
planet_object.metallicity)
M_p = plmo14.calculate_planet_mass(planet_object.core_mass, f_env) # planetary mass (mass core + mass envelope)
rho_p = rho = plmo14.density_planet(M_p, R_p) # initial approx. density
# specify beta and K
if beta_on == "yes": # use the beta estimate from Salz et al. 2016
beta = bk.beta_fct_LO14(M_p, Fxuv, R_p)
elif beta_on == "no":
beta = 1.
if K_on == "yes": # use K-correction estimation from Erkaev et al. 2007
K = bk.K_fct_LO14(planet_object.distance, M_p, planet_object.mass_star , R_p)
elif K_on == "no":
K = 1.
num = 3. * beta**2. * epsilon * Fxuv
denom = 4. * const.G.cgs * K * rho_p
M_dot = num / denom
################################################################################################################
### I define my mass loss rate to be negative, this way the planet loses mass as I integrate forward in time.###
### my RK4 function returns the mass evolution of the planet, not the mass lost over time. ###
################################################################################################################
return -(M_dot.decompose(bases=u.cgs.bases)).value # return the mass loss rate in cgs units [g/s]
def mass_loss_rate_forward_Ot20(t_, epsilon, K_on, beta_on, planet_object, radius_at_t, track_dict):
"""I use this to calculate the updated mass-loss rates at any given time step of my integration.
Input is a given time t_; this fct. calculates the corresponding Lx(t_) (then Lxuv & Fxuv given a_p),
then uses the radius at t_ to estimate a mass [based on a given mass-radius relation, which can be specified ->
e.g. mass_planet=mass_planet_CK and radius_planet=radius_planet_CK];
then we have a density, can calculate the beta and K parameter (if beta_on and K_on are set to "yes",
otherwise at each time step they beta=1 and/or K=1);
finally put it all into the mass-loss rate equation to get Mdot at t_
Parameters:
----------
t_:
Lx_evo:
calculate_planet_radius:
epsilon:
K_on:
beta_on:
planet_object:
f_env: f_env at a given time
"""
Lx = Lx_evo(t=t_, track_dict=track_dict)
Fxuv = flux_at_planet(L_xuv_all(Lx), planet_object.distance) # get flux at orbital separation a_p
# get planet density
M_p = plmoOt20.calculate_mass_planet_Ot19(radius_at_t) # estimate planetary mass using a mass-radius relation
rho_p = rho = plmoOt20.density_planet(M_p, radius_at_t) # initial approx. density
# specify beta and K
if beta_on == "yes": # use the beta estimate from Salz et al. 2016
beta = bk.beta_fct_LO14(M_p, Fxuv, radius_at_t)
elif beta_on == "no":
beta = 1.
if K_on == "yes": # use K-correction estimation from Erkaev et al. 2007
K = bk.K_fct_LO14(planet_object.distance, M_p, planet_object.mass_star , radius_at_t)
elif K_on == "no":
K = 1.
num = 3. * beta**2. * epsilon * Fxuv
denom = 4. * const.G.cgs * K * rho_p
M_dot = num / denom
################################################################################################################
### I define my mass loss rate to be negative, this way the planet loses mass as I integrate forward in time.###
### my RK4 function returns the mass evolution of the planet, not the mass lost over time. ###
################################################################################################################
return -(M_dot.decompose(bases=u.cgs.bases)).value # return the mass loss rate in cgs units [g/s]
# file with my planet class
import numpy as np
from scipy.optimize import fsolve
import astropy.units as u
from astropy import constants as const
import scipy.optimize as optimize
import os
import pandas as pd
from Lx_evo_and_flux import Lx_evo, flux_at_planet_earth, L_xuv_all, flux_at_planet
from Mass_evolution_function import mass_planet_RK4_forward_LO14 # highest level function
class planet_LoFo14():
"""
Structure of star_dictionary: {'star_id': "dummySun", 'mass': mass_star, 'radius': radius_star, 'age': age_star,
'L_bol': L_bol, 'Lx_age': Lx_age}
Structure of planet_dict: {"core_mass": m_c, "fenv": f, "distance": a, "metallicity": metal}
"""
def __init__(self, star_dictionary, planet_dict):
# initialize stellar parameters
self.star_id = star_dictionary['star_id']
self.mass_star = star_dictionary['mass']
self.radius_star = star_dictionary['radius']
self.age = star_dictionary['age']
self.Lbol = star_dictionary['L_bol'] * const.L_sun.cgs.value
self.Lx_age = star_dictionary['Lx_age']
# initialize planet parameters
# I have three options for planet_dict:
# 1) artificial planet with M_core & f_env -> need to calculate M_pl & R_pl
# 2) observed planet with a mass: M_pl & R_pl & f_env -> M_core
# 3) observed planet w/o mass: M_core, R_pl -> calulate f_env & M_pl
# following 3 params are the same for all planets!
self.distance = planet_dict["distance"]
self.metallicity = planet_dict["metallicity"] # solarZ or enhZ
self.flux = flux_at_planet_earth(self.Lbol, self.distance)
self.has_evolved = False # flag which tracks if planet has been evolved and result file exists
self.planet_id = "None"
# the remaining parameters depend on if you choose a case 1, 2, or 3 planet
while True:
try:
planet_dict["fenv"]
#print("Case 1 - artificial planet (with M_core & f_env)") # Case 1
self.planet_info = "Case 1 - artificial planet (with M_core & f_env)"
# Case 1: artificial planet with fenv & M_core given -> need to calculate the planet mass
self.fenv = planet_dict["fenv"]
self.core_mass = planet_dict["core_mass"]
self.Calculate_core_radius()
self.Calculate_planet_mass()
self.Calculate_planet_radius()
break
except KeyError: # Case 2 or 3
#print("No fenv given.")
while True:
try:
planet_dict["mass"]
# Case 3 - observed planet with radius but no mass & M_core specified -> need to calculate/estimate fenv & M_pl
#print("Case 3 - obs. planet with radius & mass measurement (with R_pl, M_pl, M_core)")
self.planet_info = "Case 3 - obs. planet with radius & mass measurement (with R_pl, M_pl, M_core)"
self.core_mass = planet_dict["core_mass"]
self.mass = planet_dict["mass"]
self.radius = planet_dict["radius"]
self.Calculate_core_radius()
self.Solve_for_fenv() # function to solve for the envelope mass fraction
break
except KeyError:
#print("Case 2 - obs. planet with radius, but no mass (with R_pl, M_core)")
self.planet_info = "Case 2 - obs. planet with radius, but no mass (with R_pl, M_core)"
# Case 2 - observed planet with a mass, radius & core mass specified -> need to calculate fenv
self.core_mass = planet_dict["core_mass"]
self.radius = planet_dict["radius"]
self.Calculate_core_radius()
self.Solve_for_fenv() # function to solve for the envelope mass fraction
self.Calculate_planet_mass()
break
break
# Class Methods
def Calculate_planet_mass(self):
""" Planet mass determined by core mass and atmosphere mass (specified in terms of atm. mass fraction [%]). """
#M_core = (self.core_mass/const.M_earth.cgs).value
self.mass = self.core_mass/(1-(self.fenv/100))
def Calculate_planet_core_mass(self):
""" Planet core mass determined by planet mass and atmosphere mass (specified in terms of atm. mass fraction [%]). """
self.core_mass = self.mass*(1-(self.fenv/100))
def Calculate_core_radius(self):
"""M-R relation for rock/iron Earth-like core. (no envelope)"""
self.core_radius = (self.core_mass**0.25)
def Solve_for_fenv(self):
if self.radius == self.core_radius:
self.fenv = 0.0
else:
def Calculate_fenv(fenv):
age_exponent = {"solarZ": -0.11, "enhZ": -0.18}
return -self.radius + self.core_radius + (2.06 * (self.core_mass/(1-(fenv/100)))**(-0.21) \
* (fenv/5)**0.59 * (self.flux)**0.044 * \
((self.age/1e3)/5)**(age_exponent[self.metallicity]))
f_guess = 0.1
fenv1 = fsolve(Calculate_fenv, f_guess)[0]
#print("fenv1 = ", fenv1)
fenv = optimize.fsolve(Calculate_fenv, x0=f_guess)[0]
#print("fenv = ", fenv)
self.fenv = fenv
if fenv1 != fenv:
print("Sth went wrong in solving for the envelope mass fraction! Check!")
def Calculate_R_env(self):
""" M_p in units of Earth masses, f_env in percent, Flux in earth units, age in Gyr
R_env ~ t**0.18 for *enhanced opacities*;
R_env ~ t**0.11 for *solar metallicities*
"""
age_exponent = {"solarZ": -0.11, "enhZ": -0.18}
R_env = 2.06 * self.mass**(-0.21) * (self.fenv/5)**0.59 * \
self.flux**0.044 * ((self.age/1e3)/5)**(age_exponent[self.metallicity])
return R_env # in units of R_earth
def Calculate_planet_radius(self):
""" description: """
self.radius = self.core_radius + self.Calculate_R_env()
def set_name(self, t_final, initial_step_size, epsilon, K_on, beta_on, evo_track_dict):
""" function to set the right planet name based on the track specified. This can then be used to check if
a particular planet on a particular track has already evolved. """
self.planet_id = 'planet_a'+str(np.round(self.distance, 3)).replace('.', 'dot') + \
'_Mcore'+str(np.round(self.core_mass,3)).replace(".", "dot") + '_fenv' + \
str(np.round(self.fenv,3)) + '_' + self.metallicity + \
'_Mstar' + str(np.round(self.mass_star,3)).replace(".", "dot") + "_K_" + K_on + "_beta_" + beta_on + \
"_track_" + str(evo_track_dict["t_start"]) + "_" + str(evo_track_dict["t_sat"]) + \
"_" + str(t_final) + "_" + str(evo_track_dict["Lx_max"]) + \
"_" + str(evo_track_dict["dt_drop"]) + "_" + str(evo_track_dict["Lx_drop_factor"])
def evolve_forward(self, t_final, initial_step_size, epsilon, K_on, beta_on, evo_track_dict, path_for_saving, planet_folder_id):
""" Call this function to make the planet evolve. """
# create a planet ID -> used for filname to save the result
# self.planet_id = 'planet_a'+str(np.round(self.distance, 3)).replace('.', 'dot') + \
# '_Mcore'+str(np.round(self.core_mass,3)).replace(".", "dot") + '_fenv' + \
# str(np.round(self.fenv,3)) + '_' + self.metallicity + \
# '_Mstar' + str(np.round(self.mass_star,3)).replace(".", "dot") + "_K_" + K_on + "_beta_" + beta_on + \
# "_track_" + str(evo_track_dict["t_start"]) + "_" + str(evo_track_dict["t_sat"]) + \
# "_" + str(t_final) + "_" + str(evo_track_dict["Lx_max"]) + \
# "_" + str(evo_track_dict["dt_drop"]) + "_" + str(evo_track_dict["Lx_drop_factor"])
# write info with planet parameters to file!!
# also track params!
self.planet_id = planet_folder_id + "_track_" + str(evo_track_dict["t_start"]) + "_" + str(evo_track_dict["t_sat"]) + \
"_" + str(t_final) + "_" + str(evo_track_dict["Lx_max"]) + \
"_" + str(evo_track_dict["dt_drop"]) + "_" + str(evo_track_dict["Lx_drop_factor"])
if os.path.exists(path_for_saving+self.planet_id+".txt"):
#print("Planet already existis!")
self.has_evolved = True
else:
print("Planet: ", self.planet_id+".txt")
###################################
# create a file: planet_XXXX.txt which contains the initial planet params
if os.path.exists(path_for_saving+planet_folder_id+".txt"):
pass
else:
p = open(path_for_saving+planet_folder_id+".txt", "a")
planet_params = "a,core_mass,fenv,metallicity\n"+ str(self.distance) + "," + str(self.core_mass) + "," \
+ str(self.fenv) + "," + self.metallicity
p.write(planet_params)
p.close()
###################################
# create a file: planet_XXXX.txt which contains the track params
if os.path.exists(path_for_saving+"track_params_"+self.planet_id+".txt"):
pass
else:
t = open(path_for_saving+"track_params_"+self.planet_id+".txt", "a")
track_params = "t_start,t_sat,t_curr,t_5Gyr,Lx_max,Lx_curr,Lx_5Gyr,dt_drop,Lx_drop_factor\n" \
+ str(evo_track_dict["t_start"]) + "," + str(evo_track_dict["t_sat"]) + "," \
+ str(evo_track_dict["t_curr"]) + "," + str(evo_track_dict["t_5Gyr"]) + "," \
+ str(evo_track_dict["Lx_max"]) + "," + str(evo_track_dict["Lx_curr"]) + "," \
+ str(evo_track_dict["Lx_5Gyr"]) + "," + str(evo_track_dict["dt_drop"]) + "," \
+ str(evo_track_dict["Lx_drop_factor"])
t.write(track_params)
t.close()
###################################
# create a file: host_star_properties.txt which contains the host star params
if os.path.exists(path_for_saving+"host_star_properties.txt"):
pass
else:
s = open(path_for_saving+"host_star_properties.txt", "a")
star_params = "star_id,mass_star,radius_star,age,Lbol,Lx_age\n" + self.star_id + "," + \
str(self.mass_star) + "," + str(self.radius_star) + "," + str(self.age) + "," + \
str(self.Lbol/const.L_sun.cgs.value) + "," + str(self.Lx_age)
s.write(star_params)
s.close()
###################################