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

updated Beta_K_functions, Planet_models_LoFo14, Planet_models_Ot20 and...

updated Beta_K_functions, Planet_models_LoFo14, Planet_models_Ot20 and Mass_evolution_function to be compliant with PEP8 formatting standards.
parent e8be8feb
......@@ -3,67 +3,85 @@ import astropy.units as u
from astropy import constants as const
import numpy as np
# TO DO: 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 comes from a study by Salz et al. (2016)
-> 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.'
def beta_fct(M_p, F_xuv, R_p):
"""Function estimates the beta parameter, which is a correction to the
planetary absorption radius in XUV, as planet appers larger than in
optical; this approximation comes from a study by Salz et al. (2016)
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."
Parameters:
----------
@param M_p: (float) mass of the planet in Earth units
@param F_xuv: (float) XUV flux recieved by planet
@param R_p: (float) radius mass of the planet in Earth units
Output:
beta parameter
M_p (float): mass of the planet in Earth units
F_xuv (float or list): XUV flux recieved by planet
R_p (float): radius mass of the planet in Earth units
Returns:
--------
beta (float or array): beta parameter
"""
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
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)
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
elif len(F_xuv) > 1:
# if F_xuv is a list
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)
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)
def K_fct(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).
Parameters:
----------
a_pl: star-planet separation in A.U.
M_pl: mass of the planet in Earth units
M_star: mass of the star in solar units
R_pl: radius of the planet in Earth units
-----------
a_pl (float): star-planet separation in A.U.
M_pl (float or list): mass of the planet in Earth units
M_star (float): mass of the star in solar units
R_pl (float): radius of the planet in Earth units
Returns:
--------
K (float) or K (array)
"""
# define some constants
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
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)
K = 1. - 3.*(R_pl*R_earth)/(2.*R_roche) \
+ (R_pl*R_earth)**3./(2.*R_roche**3.)
return K
elif len(M_pl) > 1: # if F_xuv is array
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.)
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 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 plmoLoFo14
import Planet_model_Ot20 as plmoOt20
......
......@@ -2,30 +2,39 @@ import numpy as np
import astropy.units as u
from astropy import constants as const
def calculate_mass_planet_Ot19(R_pl):
""" Estimate a planet mass (in Earth masses) based on the mass-radius
relation derived for the volatile rich regime in Otegi et al. 2020
(radius planet larger than 2.115 Earth radii, or density < 3.3 g/cc).
Mass has to be smaller than 120 Earth masses. Radius should be given
in Earth radii.
"""
I only use the volatile rich regime! This means the radius needs
to be greater than ~2.115 R_earth (radius with density > 3.3)
"""
M_earth = const.M_earth.cgs.value
R_earth = const.R_earth.cgs.value
if (type(R_pl)==int) or (type(R_pl)==float) or (type(R_pl)==np.float64): # if R is single value
M_p_volatile = 1.74*(R_pl)**1.58 # if rho < 3.3 g/cm^3
rho_volatile = M_p_volatile*M_earth/(4/3*np.pi*(R_pl*R_earth)**3)
if (type(R_pl) == int) or (type(R_pl) == float) \
or (type(R_pl) == np.float64): # if R is single value
M_p_volatile = 1.74*(R_pl)**1.58 # if rho < 3.3 g/cm^3
rho_volatile = (M_p_volatile*M_earth) \
/ (4/3*np.pi*(R_pl*R_earth)**3)
if (rho_volatile >= 3.3):
raise Exception("Planet with this radius is too small and likely rocky; use LoFo14 models instead.")
raise Exception("Planet with this radius is too small and likely \
rocky; use LoFo14 models instead.")
else:
if (M_p_volatile >= 120):
raise Exception("Planet too massive. M-R relation only valid for <120 M_earth.")
raise Exception("Planet too massive. M-R relation only valid \
for <120 M_earth.")
else:
return M_p_volatile
elif len(R_pl) > 1: # if R is array
elif len(R_pl) > 1: # if R is array
Ms = []
for i in range(len(R_pl)):
M_p_volatile_i = 1.74*(R_pl[i])**1.58 # if rho < 3.3 g/cm^3
rho_volatile_i = M_p_volatile_i*M_earth/(4/3*np.pi*(R_pl[i]*R_earth)**3)
M_p_volatile_i = 1.74*(R_pl[i])**1.58 # if rho < 3.3 g/cm^3
rho_volatile_i = (M_p_volatile_i * M_earth) \
/ (4/3*np.pi*(R_pl[i]*R_earth)** 3)
if (rho_volatile_i >= 3.3) or (M_p_volatile_i >= 120):
M_i = np.nan
else:
......@@ -34,33 +43,40 @@ def calculate_mass_planet_Ot19(R_pl):
Ms = np.array(Ms)
return Ms
def calculate_radius_planet_Ot19(M_pl):
""" Estimate a planet radius (in Earth radii) based on the mass-radius
relation derived for the volatile rich regime in Otegi et al. 2020.
Mass needs to be between ~5.7 and 120 Earth masses and should be given
in Earth masses.
"""
I only use the volatile rich regime! This means the mass needs
to be bewteen ~5.7 and 120 M_earth
"""
M_earth = const.M_earth.cgs.value
R_earth = const.R_earth.cgs.value
if (type(M_pl)==int) or (type(M_pl)==float) or (type(M_pl)==np.float64): # if M is single value
#R_p_volatile = 0.7*M_pl**0.63 # if rho < 3.3 g/cm^3
if (type(M_pl) == int) or (type(M_pl) == float) \
or (type(M_pl) == np.float64): # if M is single value
# R_p_volatile = 0.7*M_pl**0.63 # if rho < 3.3 g/cm^3
R_p_volatile = (M_pl/1.74)**(1./1.58)
rho_volatile = M_pl*M_earth/(4/3*np.pi*(R_p_volatile*R_earth)**3)
rho_volatile = (M_pl*M_earth) \
/ (4/3*np.pi*(R_p_volatile*R_earth)**3)
if (rho_volatile >= 3.3):
raise Exception("Planet with this mass/radius is too small and \
likely rocky; use LoFo14 models instead.")
else:
if (M_pl >= 120):
raise Exception("Planet too massive. M-R relation only valid for <120 M_earth.")
raise Exception("Planet too massive. M-R relation only \
valid for <120 M_earth.")
else:
return R_p_volatile
elif len(M_pl) > 1: # if M is array
elif len(M_pl) > 1: # if M is array
Rs = []
for i in range(len(M_pl)):
#R_p_volatile_i = 0.7*M_pl[i]**0.63 # if rho < 3.3 g/cm^3
# R_p_volatile_i = 0.7*M_pl[i]**0.63 # if rho < 3.3 g/cm^3
R_p_volatile_i = (M_pl[i]/1.74)**(1./1.58)
rho_volatile_i = M_pl[i]*M_earth/(4/3*np.pi*(R_p_volatile_i*R_earth)**3)
rho_volatile_i = (M_pl[i]*M_earth)\
/ (4/3*np.pi*(R_p_volatile_i*R_earth)**3)
if (rho_volatile_i >= 3.3) or (M_pl[i] >= 120):
R_i = np.nan
else:
......@@ -68,8 +84,22 @@ def calculate_radius_planet_Ot19(M_pl):
Rs.append(R_i)
Rs = np.array(Rs)
return Rs
def density_planet(M_p, R_p):
"""once I have the radius and a mass estimate for the planet and can estimate a density"""
rho = (M_p*const.M_earth.cgs/(4./3*np.pi*(R_p*const.R_earth.cgs)**3)).cgs
return rho.value
\ No newline at end of file
""" Calculate the mean density in cgs units using the radius and
mass of the planet.
Parameters:
-----------
M_p (float): mass of the planet (in Earth masses)
R_p (float): radius of the planet (in Earth radii)
Returns:
--------
rho (float): density in cgs units (g/ccm)
"""
rho = (M_p*const.M_earth.cgs) \
/ (4./3*np.pi*(R_p*const.R_earth.cgs)**3)
return rho.value
# Lopez & Fortney 2014 models
import astropy.units as u
from astropy import constants as const
import numpy as np
import warnings
import astropy.units as u
from astropy import constants as const
# import warnings
# warnings.filterwarnings("ignore", category=RuntimeWarning)
# ignoring this warining is certainly not the best way, but it is currently
# thrown when one of the "in-between" time steps inside the Runge-Kutta
# integration method results in a planet which has no atmosphere left.
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)
# ignoring this warining is certainly not the best way, but it is currently thrown when one of the "in-between"
# time steps inside the Runge-Kutta integration method results in a planet which has no atmosphere left.
def calculate_core_radius(M_core):
""" M-R relation for rock/iron Earth-like core. (no envelope) """
""" M-R relation for rock/iron Earth-like core. (no envelope)
(see Lopez & Fortney 2014 for details.)
"""
R_core = (M_core**0.25)
return R_core
def calculate_planet_mass(M_core, fenv):
""" Planet mass determined by core mass and atmosphere mass (specified in terms of atm. mass fraction [%]). """
""" Planet mass determined by core mass and atmosphere mass
(specified in terms of envelope mass fraction [in % !]). """
M_pl = M_core/(1-(fenv/100))
return M_pl
def calculate_R_env(M_p, fenv, F_p, age, metallicity):
""" M_p in Earth masses, f_env in percent, Flux in earth units, age in Gyr
R_env ~ t**0.18 for *enhanced opacities* """
""" Calculate planetary envelope radius using the parametrized
results from Lopez & Fortney 2014 for planets with rocky core
and H/He gaseous envelope.
Parameters:
-----------
M_p (float): planet mass (in Earth masses)
f_env (float): envelope mass fraction (given in percent!)
F_p (float): bolometric flux revcieved by the planet (in units
of Earth's insolation)
age (float): age of the planet (in Myr)
metallicity (str): chose models with solar or enhanced metallicity
(set to "solarZ" or "enhZ")
(check Lopez & Fortney 2014 for further details and range of model
validity).
Returns:
--------
R_env (float): radius of the envelope (in Earth radii)
"""
age_exponent = {"solarZ": -0.11, "enhZ": -0.18}
R_env = 2.06 * (M_p)**(-0.21) * (fenv/5)**0.59 * (F_p)**0.044 * ((age/1e3)/5)**(age_exponent[metallicity]) # R_earth
return R_env
R_env = 2.06 * (M_p)**(-0.21) * (fenv/5)**0.59 * (F_p)**0.044 \
* ((age/1e3)/5)**(age_exponent[metallicity])
return R_env # R_earth
def calculate_planet_radius(M_core, fenv, age,
flux_solar, metallicity):
""" Calculate planetary radius (core + envelope) using the
parametrized results from Lopez & Fortney 2014 for planets with
rocky core and H/He gaseous envelope.
Parameters:
-----------
M_core (float): core mass (in Earth masses)
f_env (float): envelope mass fraction (given in percent!)
flux_solar (float): bolometric flux revcieved by the planet
(in units of Earth's flux)
age (float): age of the planet (in Myr)
metallicity (str): chose models with solar or enhanced metallicity
(set to "solarZ" or "enhZ")
(check Lopez & Fortney 2014 for further details and range of model
validity).
Returns:
--------
R_pl (float): radius of the planet (in Earth radii)
"""
def calculate_planet_radius(M_core, fenv, age, flux_solar, metallicity):
""" description: enhanced opacities; age in Gyr, flux in solar units """
R_core = calculate_core_radius(M_core)
M_pl = calculate_planet_mass(M_core, fenv)
R_env = calculate_R_env(M_pl, fenv, flux_solar, age, metallicity)
R_pl = R_core + R_env
return R_pl
def density_planet(M_p, R_p):
"""once I have the radius and a mass estimate for the planet and can estimate a density"""
rho = (M_p*const.M_earth.cgs/(4./3*np.pi*(R_p*const.R_earth.cgs)**3))
return rho.value
\ No newline at end of file
""" Calculate the mean density in cgs units using the radius and
mass of the planet.
Parameters:
-----------
M_p (float): mass of the planet (in Earth masses)
R_p (float): radius of the planet (in Earth radii)
Returns:
--------
rho (float): density in cgs units (g/ccm)
"""
rho = (M_p*const.M_earth.cgs) \
/ (4./3*np.pi*(R_p*const.R_earth.cgs)**3)
return rho.value
{
"cells": [
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"a = 5\n",
"b = True"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"step 1\n",
"step 2\n",
"final result\n",
"update planet params & check radius change\n"
]
}
],
"source": [
"for a in [0,5,4,2]:\n",
" while (b==True):\n",
" print(\"step 1\")\n",
" if a == 1:\n",
" print(\"break out of 1\")\n",
" break\n",
" print(\"step 2\")\n",
" if a == 2:\n",
" print(\"break out of 2\")\n",
" break\n",
" print(\"final result\")\n",
" print(\"update planet params & check radius change\")\n",
" b = False\n",
" break"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"False"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"nan nan\n",
"a == b: False\n",
"True True\n"
]
}
],
"source": [
"import numpy as np\n",
"\n",
"nan1 = np.nan\n",
"nan2 = float(\"nan\")\n",
"print(nan1, nan2)\n",
"print(\"nan1 == nan2: \", nan1 == nan2)\n",
"print(np.isnan(nan1), np.isnan(nan2))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
star_id,mass_star,radius_star,age,Lbol,Lx_age
star_age5.0_mass0.73,0.7269788466496349,None,5.0,0.22435540668462725,4.980171083198093e+29
\ No newline at end of file
a,core_mass,fenv,mass,radius,metallicity,age
0.17107652861345243,4.771653826371514,8.491315192780585,5.214427282420153,6.416076972774628,solarZ,5.0
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment