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

updated exampleV1298Tau notebook

parent 22170acc
......@@ -3,21 +3,21 @@ import astropy.units as u
from astropy import constants as const
import numpy as np
# rename to beta_fct & K_fct
# 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 is used in Lalitha et al. 2018 and comes
originally from Salz M., Schneider P., Czesla S., Schmitt J., 2016a, Astronomy & Astrophysics, 585, L2
-> 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.'
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)
@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_earth = const.M_earth.cgs.value
R_earth = const.R_earth.cgs.value
......@@ -40,14 +40,13 @@ def beta_fct_LO14(M_p, F_xuv, R_p):
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:
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
"""
AU = const.au.cgs.value
M_sun = const.M_sun.cgs.value
......
# 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!!
-> so t_curr(ent) is a little confusing! -> 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):
......@@ -26,10 +17,11 @@ def Lx_evo(t, track_dict):
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)
@param t_curr: (float)
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
activity tracks converge (based on Tu et al. (2015)/ Johnstone et al.(2015) models);
@param Lx_curr: (float)
Lx value at 1 Gyr taken from Tu et al. (2015)
@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));
......@@ -37,6 +29,7 @@ def Lx_evo(t, track_dict):
Defines, together with t_5Gyr, the slope common for all three
activity evolution tracks past 1 Gyr (also based on Tu et al.(2015));
# IGNORE FOR NOW CASE 2 - not implemented yet
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):
---------------------------------------------------------------------------
......@@ -73,7 +66,6 @@ def Lx_evo(t, track_dict):
# 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:
......@@ -82,7 +74,7 @@ def Lx_evo(t, track_dict):
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
# 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)
......@@ -95,24 +87,21 @@ def Lx_evo(t, track_dict):
###################################################################################
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\
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
......@@ -136,8 +125,7 @@ def Lx_evo(t, track_dict):
# 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)"""
-> the semi-major axis is used for the distance (the eccentricity of the planetary orbits is ignored)"""
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)
......@@ -147,16 +135,17 @@ def flux_at_planet_earth(L, a_p):
#####################################################################################################################
# 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)"""
"""Function calculates the flux that the planet recieves at a given distance (in erg/s/cm^2)
-> the semi-major axis is used for the distance (the eccentricity of the planetary orbits is ignored)"""
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.
the measured X-ray luminosity can be extrapolated to the total high-energy radiation as given below.
To get the combined high-energy luminosity (L_xuv), we add L_x and L_euv.
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.)
"""
......@@ -165,17 +154,13 @@ def L_xuv_all(Lx):
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 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"""
""" Typical relation (from observations) to estimate the saturation X-ray luminosity given a bolometric luminosity."""
return 10**(-3)*(L_star_bol*const.L_sun.cgs).value
\ No newline at end of file
# 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_models_LoFo14 as 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.
"""Calculate the updated mass-loss rates at any given time step (of the integration) for a
Lopez & Fortney (2014) planet with rocky core and H/He envelope.
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 uses the radius at t_ to estimate a mass (given the current envelope mass fraction fenv);
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);
otherwise at each time step they are beta=1 and/or K=1);
finally put it all into the mass-loss rate equation to get Mdot at t_
Parameters:
......@@ -32,12 +32,11 @@ def mass_loss_rate_forward_LO14(t_, epsilon, K_on, beta_on, planet_object, f_env
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
R_p = plmoLoFo14.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 = plmoLoFo14.calculate_planet_mass(planet_object.core_mass, f_env) # planetary mass (mass core + mass envelope)
rho_p = rho = plmoLoFo14.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
......@@ -56,16 +55,16 @@ def mass_loss_rate_forward_LO14(t_, epsilon, K_on, beta_on, planet_object, f_env
### 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]
return -(M_dot.decompose(bases=u.cgs.bases)).value # 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.
"""Calculate the updated mass-loss rates at any given time step (of the integration) for a
planet which follows the "mature" mass-radius relation as given in Otegi et al. (2020).
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 uses the radius at t_ to estimate a mass based on the volatile-regime mass-radius relation;
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);
otherwise at each time step they are set to beta=1 and/or K=1);
finally put it all into the mass-loss rate equation to get Mdot at t_
Parameters:
......
......@@ -3,15 +3,23 @@ import math
import numpy as np
def create_planet_chunks(curr_path, folder_name, list_planets, chunk_size):
""" """
# check if directroy for saving the results existst, if not create
""" Function creates the directory & subdirectory structure for saving the results
for all the planets in "list_planets". In addition, it divides list_planets into smaller chunks.
This is needed if "list_planets" is very long, and avoids the multiprocessing-action in the evolve_ensamble-function from crashing.
Return:
------
path_save: path specifying where to save the results
planets_chunks: chunked-up list of planets, which can be passed to "evolve_ensemble" to evolve
"""
path_save = curr_path+folder_name
if os.path.isdir(path_save):
# check if directroy for saving the results existst, if not create
print("Folder -> "+folder_name+" <- exists.")
pass
else:
os.mkdir(path_save)
# create folders, one for each planet in the population
# (what this means: two planets with the exact same parametes will get two different folders)
......
......@@ -2,7 +2,7 @@ import os
import pandas as pd
def create_summary_files_with_final_planet_parameters(path_save):
""" Function NOT finished. In principle it should create a summary file of initial and final parameters."""
files = os.listdir(path_save)
non_empty_folders = []
for f in files:
......
......@@ -3,7 +3,7 @@ import multiprocessing as mp
import multiprocessing
import time
import sys
sys.path.append('../plaml_package/')
sys.path.append('../platypos_package/')
# Planet Class
from Planet_class_LoFo14 import planet_LoFo14
from Planet_class_Ot20 import planet_Ot20
......@@ -12,37 +12,39 @@ from Planet_class_Ot20_PAPER import planet_Ot20_PAPER
def evolve_one_planet(pl_folder_pair, t_final, init_step, eps, K_on, beta_on, evo_track_dict_list, path_for_saving):
"""Evolves one planet, pl, at a time, but through all evo tracks in evo_track_dict_list.
So all calculations for one planet belong to each other, but the planets are independent of each other.
Each time this function is called represents an independent process. """
"""Function to evolve one planet at a time, but through all evo tracks specified in evo_track_dict_list
(dictionary with info about evolutionary tracks).
So all calculations for this one planet belong to each other, but are independent from another planet.
Each function call to "evolve_one_planet" represents an independent process."""
for track in evo_track_dict_list:
pl = pl_folder_pair[1]
folder = pl_folder_pair[0]
pl = pl_folder_pair[1] # get the planet object
folder = pl_folder_pair[0] # and the correcponding folder name where results are saved
pl.set_name(t_final, init_step, eps, K_on, beta_on, evo_track_dict=track) # set planet name based on specified track
# generate a useful planet name for saving the results
pl_file_name = folder + "_track" + pl.planet_id.split("_track")[1] + ".txt"
if os.path.isdir(path_for_saving+pl_file_name):
# skip & don't evolve planet
# print("Folder "+pl_file_name+" exists.")
# if file exists, skip & don't evolve planet
pass
else:
# evolve the planet
pl.evolve_forward(t_final, init_step, eps, K_on, beta_on, evo_track_dict=track,
path_for_saving=path_for_saving, planet_folder_id=folder)
def evolve_ensamble(planets_chunks, t_final, initial_step_size, epsilon, K_on, beta_on, evo_track_dict_list, path_save):
"""Function which parallelizes the multiprocessing_func.
@param planets_chunks: (array of lists) Array of lists of [folder, planet] pairs.
Each list (i.e. list of planets) in the array will be run seperately.
@param t_final: (float)
@param initial_step_size: (float)
@param epsilon: (float)
@param K_on: (str)
@param beta_on: (str)
@param evo_track_dict_list: (list) List of track_dictionaries along which to evolve each planet
"""Function which parallelizes the multiprocessing_func "evolve_one_planet".
@param planets_chunks: (array of lists) Array of lists of [folder, planet] pairs. Each chunk of planets
in the array will be run seperately, one after the other. This avoids computer crashes.
@param t_final: (float) Time at which to end the simulation (e.g. 5 or 10 Gyrs)
@param initial_step_size: (float) Initial step size for integration (old code uses a fixed step size, newer version uses an adptive one)
@param epsilon: (float) Evaporation efficieny parameter
@param K_on: (str) Tidal-effects correction
@param beta_on: (str) X-ray cross section correction
@param evo_track_dict_list: (list) List of track dictionaries along which to evolve each planet
@param path_save: (str) Filepath to the folder in which to save the planet-folders
@result: Each planet will evolve along all the specified tracks.
Each track-planet evolution result will be saved as a csv file in the planet-folder.
Each track-planet evolution result will be saved as a .txt file in the planet-folder.
"""
print("start")
#if __name__ == '__main__':
......
......@@ -2,7 +2,7 @@ import os
import pandas as pd
def read_results_file(path, filename):
"""function to read in the results file"""
"""Function to read in the results file for an individual track. """
df = pd.read_csv(path+filename, float_precision='round_trip') # to avoid pandas doing sth. weird to the last digit
# Pandas uses a dedicated dec 2 bin converter that compromises accuracy in preference to speed.
# Passing float_precision='round_trip' to read_csv fixes this.
......@@ -10,7 +10,7 @@ def read_results_file(path, filename):
return df#t, M, R, Lx
def read_in_PLATYPOS_results(path_to_results, N_tracks):
""" call this function to read in the results.
""" Call this function to read in ALL the results.
Input:
------
path_to_results: path to the folder which containts all the results (i.e. all the folders for the individual planets)
......@@ -44,10 +44,7 @@ def read_in_PLATYPOS_results(path_to_results, N_tracks):
# Next, read in the results for each planet
planet_df_dict = {} # dictionary of results, with planet names as keys, and corresponding results-dataframes as values
# track_name_dict = {} # dictionary with planet names as keys, and list of track names within each planet folder as values
# track_evolved_off_dict = {} # same structure as track_name_dict, but with True or False for each track indicating
# whether planet has evolved off grid (True) or not (False)
tracks_dict = {} # combined track_name_dict & track_evolved_off_dict
tracks_dict = {} # dictionary with planet names as keys, and list of track infos for each planet folder as values
planet_init_dict = {} # dictionary of initial planet parameters,with planet names are keys, parameters the values
for i, f in enumerate(non_empty_folders):
......@@ -96,27 +93,8 @@ def read_in_PLATYPOS_results(path_to_results, N_tracks):
# NEXT, read in track names in case I need this info later, for example for knowing the properties of each track
# track number corresponds to t1,2,3, etc..
track_dict = {}
# flag planets which have moved off (important for MESA planets for now)
# flag planets which have moved off (only important for MESA planets for now)
list_planets_evolved_off = ["track"+file.rstrip(".txt").split("track")[1] for file in all_files_in_f if "evolved_off" in file]
##########################################
# old code (contains unnecessary steps)
# track_evooff_dict = {}
# for i, file in enumerate(result_files_sorted):
# track_dict[i+1] = file[len(f+"_"):].rstrip(".txt") # track name
# if "track"+file.rstrip(".txt").split("track")[1] in list_planets_evolved_off:
# track_evooff_dict[i+1] = True # if track name in list of planets which has evolved off, set flag to True
# else:
# track_evooff_dict[i+1] = False
# # now I have a dataframe for each planet, which contains all the evolutions for each track
# planet_df_dict[f] = df_all_tracks # add to dictionary with planet name as key
# track_name_dict[f] = track_dict # add tracks list to dictionary
# track_evolved_off_dict[f] = track_evooff_dict # add flags to dictionary
# # combine the two trac dicts into one
# tracks_dict[f] = list(zip(track_name_dict[f].keys(), track_name_dict[f].values(), track_evolved_off_dict[f].values()))
##########################################
track_info_list = []
for i, file in enumerate(result_files_sorted):
if "track"+file.rstrip(".txt").split("track")[1] in list_planets_evolved_off:
......@@ -133,18 +111,13 @@ def read_in_PLATYPOS_results(path_to_results, N_tracks):
planet_init_df = pd.DataFrame.from_dict(planet_init_dict, orient='index', columns=df_pl.columns)
# now I can easily access the planet name, semi-major axies, core mass, initial mass and radius
print("\nTotal number of planets to analyze: ", len(planet_init_df))
# convert semi-major axis to period (host-star properties known)
#P_p_arr = get_period_from_a(1.0, a_p_arr)
# **In principle planet_df_dict & planet_init_df contain all info for analyzing the results**
# **(initial and final mass and radius, final age)**
return planet_df_dict, planet_init_df, tracks_dict
def read_in_PLATYPOS_results_dataframe(path_to_results, N_tracks):
"""
does some more re-aranging to the data to make it easier to handle
Calls read_in_PLATYPOS_results & then does some more re-aranging to the data to make it easier to handle.
"""
# call read_in_PLATYPOS_results
planet_df_dict, planet_init_df, tracks_dict = read_in_PLATYPOS_results(path_to_results, N_tracks)
......@@ -183,5 +156,4 @@ def read_in_PLATYPOS_results_dataframe(path_to_results, N_tracks):
for number, evooff in zip(track_number, track_evooff):
planet_all_df.at[key_pl, "track"+str(number)] = evooff
#planet_all_df.head()
return planet_all_df, tracks_dict
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