Commit 9a776365 authored by Laura Ketzer's avatar Laura Ketzer
Browse files

fixed multiprocessing issue for now -> newest running code in...

fixed multiprocessing issue for now -> newest running code in Population_evolution_OwWu17__04_06-2020_FAST.ipynb
parent 3cb69307
...@@ -3673,7 +3673,7 @@ ...@@ -3673,7 +3673,7 @@
"name": "python", "name": "python",
"nbconvert_exporter": "python", "nbconvert_exporter": "python",
"pygments_lexer": "ipython3", "pygments_lexer": "ipython3",
"version": "3.7.2" "version": "3.7.3"
} }
}, },
"nbformat": 4, "nbformat": 4,
...@@ -3739,7 +3739,7 @@ ...@@ -3739,7 +3739,7 @@
"name": "python", "name": "python",
"nbconvert_exporter": "python", "nbconvert_exporter": "python",
"pygments_lexer": "ipython3", "pygments_lexer": "ipython3",
"version": "3.7.2" "version": "3.7.3"
} }
}, },
"nbformat": 4, "nbformat": 4,
...@@ -18,8 +18,8 @@ from Mass_loss_rate_function import mass_loss_rate_forward_Ot20 ...@@ -18,8 +18,8 @@ from Mass_loss_rate_function import mass_loss_rate_forward_Ot20
def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
initial_step_size, t_final, track_dict): initial_step_size, t_final, track_dict):
"""USED: 4th order Runge-Kutta as numerical integration method """USED: 4th order Runge-Kutta as numerical integration method
Integrate from the current time (t_start (where planet has R0 and M0) into the future Integrate from the current time (t_start (where planet has R0 and M0) into
taking into account photoevaporative mass loss. the future taking into account photoevaporative mass loss.
Parameters: Parameters:
----------- -----------
...@@ -63,7 +63,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -63,7 +63,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
f_env_0 = f_env = planet_object.fenv f_env_0 = f_env = planet_object.fenv
R0 = R = planet_object.radius R0 = R = planet_object.radius
M0 = M = planet_object.mass M0 = M = planet_object.mass
rho0 = rho = plmo14.density_planet(M0, R0) rho0 = rho = plmoLoFo14.density_planet(M0, R0)
M_env0 = M_env = M0 - planet_object.core_mass M_env0 = M_env = M0 - planet_object.core_mass
M_core = planet_object.core_mass M_core = planet_object.core_mass
R_core = planet_object.core_radius R_core = planet_object.core_radius
...@@ -105,7 +105,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -105,7 +105,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
# complete atmospheric removal # complete atmospheric removal
while t <= t_final: while t <= t_final:
#print("i - t - dt: ", i, t, dt) #print(i, " t= ", t, dt)
# This step (Lx(t) calculation) is just for me to trace Lx and check if # This step (Lx(t) calculation) is just for me to trace Lx and check if
# it is correct. It is NOT required since the Lx(t) calculation is # it is correct. It is NOT required since the Lx(t) calculation is
...@@ -141,7 +141,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -141,7 +141,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
k1 = (dt*Myr_to_sec * Mdot1)/M_earth k1 = (dt*Myr_to_sec * Mdot1)/M_earth
M_05k1 = M + 0.5 * k1 # mass after 1st RK step M_05k1 = M + 0.5 * k1 # mass after 1st RK step
M_env_05k1 = M_05k1 - M_core M_env_05k1 = M_05k1 - M_core
#f_env_05k1 = (M_env_05k1/M_05k1) * 100 # new envelope mass fraction f_env_05k1 = (M_env_05k1/M_05k1) * 100 # new envelope mass fraction
if (i == 1) and (M_05k1 < M_core): if (i == 1) and (M_05k1 < M_core):
# then I am still in the first RK iteration, and the initial # then I am still in the first RK iteration, and the initial
# step size was likely too large -> set step size to minumum # step size was likely too large -> set step size to minumum
...@@ -154,7 +154,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -154,7 +154,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
k2 = (dt*Myr_to_sec * Mdot2)/M_earth k2 = (dt*Myr_to_sec * Mdot2)/M_earth
M_05k2 = M + 0.5 * k2 M_05k2 = M + 0.5 * k2
M_env_05k2 = M_05k2 - M_core M_env_05k2 = M_05k2 - M_core
#f_env_05k2 = (M_env_05k2/M_05k2) * 100 f_env_05k2 = (M_env_05k2/M_05k2) * 100
if (i == 1) and (M_05k2 < M_core): if (i == 1) and (M_05k2 < M_core):
# then I am still in the first RK iteration, and the initial # then I am still in the first RK iteration, and the initial
# step size was likely too large -> set step size to minumum # step size was likely too large -> set step size to minumum
...@@ -167,6 +167,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -167,6 +167,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
k3 = (dt*Myr_to_sec * Mdot3)/M_earth k3 = (dt*Myr_to_sec * Mdot3)/M_earth
M_k3 = M + k3 M_k3 = M + k3
M_env_k3 = M_k3 - M_core M_env_k3 = M_k3 - M_core
f_env_k3 = (M_env_k3/M_k3) * 100
if (i == 1) and (M_k3 < M_core): if (i == 1) and (M_k3 < M_core):
# then I am still in the first RK iteration, and the initial # then I am still in the first RK iteration, and the initial
# step size was likely too large -> set step size to minumum # step size was likely too large -> set step size to minumum
...@@ -186,8 +187,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -186,8 +187,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
# now it is time to check, if atmosphere is gone or # now it is time to check, if atmosphere is gone or
# if planet is close to complete atmosphere removal # if planet is close to complete atmosphere removal
if ((np.isnan(M_lost) == True) or (M_new <= M_core)) if ((np.isnan(M_lost) == True) or (M_new <= M_core)) and (dt == min_step_size):
and (dt == min_step_size):
# if M_lost = nan (i.e. planet evaporates in one of the RK # if M_lost = nan (i.e. planet evaporates in one of the RK
# steps) OR the four RK steps finish and the new planet mass is # steps) OR the four RK steps finish and the new planet mass is
# smaller or equal to the core mass, then the planet counts as # smaller or equal to the core mass, then the planet counts as
...@@ -198,8 +198,8 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -198,8 +198,8 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
# t_i+dt; this terminates the code and returns the final planet # t_i+dt; this terminates the code and returns the final planet
# properties # properties
# print("Atmosphere has evaportated! Only bare rocky core left!\ #print("Atmosphere has evaportated! Only bare rocky core left!"\
# STOP this madness!") # + " STOP this madness!")
# since the the stop criterium is reached, we assume at t_i+1 # since the the stop criterium is reached, we assume at t_i+1
# the planet only consists of the bare rocky core with the # the planet only consists of the bare rocky core with the
...@@ -212,8 +212,8 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -212,8 +212,8 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
envelope_left = False # set flag for complete env. removal envelope_left = False # set flag for complete env. removal
break break
elif ((np.isnan(M_lost) == True) or (M_new <= M_core)) elif ((np.isnan(M_lost) == True) or (M_new <= M_core)) and (dt > min_step_size) and (close_to_evaporation == False):
and (dt > min_step_size) and (close_to_evaporation == False): #print("close to evaporation")
# planet close to evaporation, but since the step size is not # planet close to evaporation, but since the step size is not
# minimum yet, we set it to its minimum value and run the RK # minimum yet, we set it to its minimum value and run the RK
# iteration again (until above stopping condition is fulfilled) # iteration again (until above stopping condition is fulfilled)
...@@ -231,14 +231,16 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -231,14 +231,16 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
# new radius and check how drastic the radius change would be; # new radius and check how drastic the radius change would be;
# adjust the step size if too drastic or too little # adjust the step size if too drastic or too little
f_env_new = (M_env_new/M_new)*100 # in % f_env_new = (M_env_new/M_new)*100 # in %
R_new = plmo14.calculate_planet_radius(M_core, f_env_new, t, R_new = plmoLoFo14.calculate_planet_radius(M_core, f_env_new, t,
flux_at_planet_earth( flux_at_planet_earth(
planet_object.Lbol, planet_object.Lbol,
planet_object.distance), planet_object.distance),
planet_object.metallicity) planet_object.metallicity)
#print("R_new, f_new: ", R_new, f_env_new)
# only adjust step size if planet is not close to complete evaporat. # only adjust step size if planet is not close to complete evaporat.
if (close_to_evaporation == False): if (close_to_evaporation == False):
#print("radius check")
# then do the check on how much the radius changes # then do the check on how much the radius changes
# R(t_i) compared to R(t_i+dt); # R(t_i) compared to R(t_i+dt);
# if radius change is larger than 1%, make step size smaller # if radius change is larger than 1%, make step size smaller
...@@ -250,8 +252,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -250,8 +252,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
R_change = abs((R-R_new)/R)*100 # radius change compared to R_change = abs((R-R_new)/R)*100 # radius change compared to
# previous radius - in percent # previous radius - in percent
if (R_change > 0.1) and (t < track_dict["t_curr"]) if (R_change > 0.5) and (t < track_dict["t_curr"]) and (dt > min_step_size):
and (dt > min_step_size):
dt = dt / 10. dt = dt / 10.
break break
...@@ -259,22 +260,30 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -259,22 +260,30 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
# dt = dt / 10. # dt = dt / 10.
# break # break
elif (R_change < (0.01)) and (t < track_dict["t_curr"]) elif (R_change < (0.01)) and (t < track_dict["t_curr"]) and (dt < max_step_size):
and (dt < max_step_size):
dt = dt * 10. dt = dt * 10.
break break
##################################################################################
# need to make sure that code does not end up in infinite loop
# if I chode my limits too close to each other, in a 0.1 Myr time
# step the radius change might bee too big, but in a 0.01 Myr time
# step it it too small -> then the code is stuck.
# in such situation I could continue with the smaller step size or chose
# sth in between.
##################################################################################
# in principle I can adjust the code such that these hardcoded # in principle I can adjust the code such that these hardcoded
# parameters are different for early planet evolution where # parameters are different for early planet evolution where
# much more is happening typically, and late planet evolution # much more is happening typically, and late planet evolution
# where almost no change is occurring anymore # where almost no change is occurring anymore
elif (R_change > 0.1) and (t >= track_dict["t_curr"]) elif (R_change > 0.5) and (t >= track_dict["t_curr"]) and (dt > min_step_size):
and (dt > min_step_size):
dt = dt / 10. dt = dt / 10.
break break
elif (R_change < (0.01)) and (t >= track_dict["t_curr"]) elif (R_change < (0.01)) and (t >= track_dict["t_curr"]) and (dt < max_step_size):
and (dt < max_step_size):
dt = dt * 10 dt = dt * 10
break break
...@@ -295,7 +304,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -295,7 +304,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
# calculate new radius with new planet mass/envelope # calculate new radius with new planet mass/envelope
# mass fraction & one time step later # mass fraction & one time step later
R = plmo14.calculate_planet_radius(M_core, f_env, t, R = plmoLoFo14.calculate_planet_radius(M_core, f_env, t,
flux_at_planet_earth( flux_at_planet_earth(
planet_object.Lbol, planet_object.Lbol,
planet_object.distance), planet_object.distance),
...@@ -324,7 +333,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -324,7 +333,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
# calculate new radius with new planet mass/envelope mass # calculate new radius with new planet mass/envelope mass
# fraction & one time step later # fraction & one time step later
R = plmo14.calculate_planet_radius(M_core, f_env, t, R = plmoLoFo14.calculate_planet_radius(M_core, f_env, t,
flux_at_planet_earth( flux_at_planet_earth(
planet_object.Lbol, planet_object.Lbol,
planet_object.distance), planet_object.distance),
...@@ -338,7 +347,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, ...@@ -338,7 +347,7 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
if (envelope_left == False): if (envelope_left == False):
# planet has evaporated, so return last planet params for bare core # planet has evaporated, so return last planet params for bare core
return np.array(t_arr), np.array(M_arr), return np.array(t_arr), np.array(M_arr),\
np.array(R_arr), np.array(Lx_arr) np.array(R_arr), np.array(Lx_arr)
#print("Done!") #print("Done!")
...@@ -379,17 +388,17 @@ def mass_planet_RK4_forward_LO14_PAPER(epsilon, K_on, beta_on, planet_object, in ...@@ -379,17 +388,17 @@ def mass_planet_RK4_forward_LO14_PAPER(epsilon, K_on, beta_on, planet_object, in
f_env_0 = f_env = planet_object.fenv f_env_0 = f_env = planet_object.fenv
R0 = R = planet_object.radius # should match observed radius - determined by M_core and f_env R0 = R = planet_object.radius # should match observed radius - determined by M_core and f_env
M0 = M = planet_object.mass M0 = M = planet_object.mass
rho0 = rho = plmo14.density_planet(M0, R0) # initial mean density rho0 = rho = plmoLoFo14.density_planet(M0, R0) # initial mean density
M_env0 = M_env = M0 - planet_object.core_mass # initial envelope mass M_env0 = M_env = M0 - planet_object.core_mass # initial envelope mass
M_core = planet_object.core_mass M_core = planet_object.core_mass
# specify beta0 and K0 # specify beta0 and K0
if beta_on == "yes": # use the beta estimate from Salz et al. 2016 if beta_on == "yes": # use the beta estimate from Salz et al. 2016
beta = beta0 = bk.beta_fct_LO14(M0, Fxuv0, R0) beta = beta0 = bk.beta_fct(M0, Fxuv0, R0)
elif beta_on == "no": elif beta_on == "no":
beta = beta0 = 1. beta = beta0 = 1.
if K_on == "yes": # use K-correction estimation from Erkaev et al. 2007 if K_on == "yes": # use K-correction estimation from Erkaev et al. 2007
K = K0 = bk.K_fct_LO14(planet_object.distance, M0, planet_object.mass_star, R0) K = K0 = bk.K_fct(planet_object.distance, M0, planet_object.mass_star, R0)
elif K_on == "no": elif K_on == "no":
K = K0 = 1. K = K0 = 1.
...@@ -462,7 +471,7 @@ def mass_planet_RK4_forward_LO14_PAPER(epsilon, K_on, beta_on, planet_object, in ...@@ -462,7 +471,7 @@ def mass_planet_RK4_forward_LO14_PAPER(epsilon, K_on, beta_on, planet_object, in
# new envelope mass fraction: # new envelope mass fraction:
f_env = (M_env/M)*100 # in % f_env = (M_env/M)*100 # in %
# calculate new radius with new planet mass/envelope mass fraction & one time step later # calculate new radius with new planet mass/envelope mass fraction & one time step later
R_arr[i+1] = R = plmo14.calculate_planet_radius(M_core, f_env, t, flux_at_planet_earth(planet_object.Lbol, R_arr[i+1] = R = plmoLoFo14.calculate_planet_radius(M_core, f_env, t, flux_at_planet_earth(planet_object.Lbol,
planet_object.distance), planet_object.metallicity) planet_object.distance), planet_object.metallicity)
else: else:
...@@ -534,11 +543,11 @@ def mass_planet_RK4_forward_Ot14(epsilon, K_on, beta_on, planet_object, initial_ ...@@ -534,11 +543,11 @@ def mass_planet_RK4_forward_Ot14(epsilon, K_on, beta_on, planet_object, initial_
# specify beta0 and K0 # specify beta0 and K0
if beta_on == "yes": # use the beta estimate from Salz et al. 2016 if beta_on == "yes": # use the beta estimate from Salz et al. 2016
beta = beta0 = bk.beta_fct_LO14(M0, Fxuv0, R0) beta = beta0 = bk.beta_fct(M0, Fxuv0, R0)
elif beta_on == "no": elif beta_on == "no":
beta = beta0 = 1. beta = beta0 = 1.
if K_on == "yes": # use K-correction estimation from Erkaev et al. 2007 if K_on == "yes": # use K-correction estimation from Erkaev et al. 2007
K = K0 = bk.K_fct_LO14(planet_object.distance, M0, planet_object.mass_star, R0) K = K0 = bk.K_fct(planet_object.distance, M0, planet_object.mass_star, R0)
elif K_on == "no": elif K_on == "no":
K = K0 = 1. K = K0 = 1.
...@@ -682,11 +691,11 @@ def mass_planet_RK4_forward_Ot14_PAPER(epsilon, K_on, beta_on, planet_object, in ...@@ -682,11 +691,11 @@ def mass_planet_RK4_forward_Ot14_PAPER(epsilon, K_on, beta_on, planet_object, in
# specify beta0 and K0 # specify beta0 and K0
if beta_on == "yes": # use the beta estimate from Salz et al. 2016 if beta_on == "yes": # use the beta estimate from Salz et al. 2016
beta = beta0 = bk.beta_fct_LO14(M0, Fxuv0, R0) beta = beta0 = bk.beta_fct(M0, Fxuv0, R0)
elif beta_on == "no": elif beta_on == "no":
beta = beta0 = 1. beta = beta0 = 1.
if K_on == "yes": # use K-correction estimation from Erkaev et al. 2007 if K_on == "yes": # use K-correction estimation from Erkaev et al. 2007
K = K0 = bk.K_fct_LO14(planet_object.distance, M0, planet_object.mass_star, R0) K = K0 = bk.K_fct(planet_object.distance, M0, planet_object.mass_star, R0)
elif K_on == "no": elif K_on == "no":
K = K0 = 1. K = K0 = 1.
......
...@@ -216,7 +216,7 @@ class planet_LoFo14(): ...@@ -216,7 +216,7 @@ class planet_LoFo14():
#print("Saved!\n") #print("Saved!\n")
# make another file, which contains the final parameters # make another file, which contains the final parameters
if os.path.exists(path_for_saving+planet_folder_id+"_"+self.planet_id+"_final.txt"): if os.path.exists(path_for_saving+"_"+self.planet_id+"_final.txt"):
pass pass
else: else:
p = open(path_for_saving+self.planet_id+"_final.txt", "a") p = open(path_for_saving+self.planet_id+"_final.txt", "a")
......
import os
import numpy as np import numpy as np
import pandas as pd
from scipy.optimize import fsolve from scipy.optimize import fsolve
import scipy.optimize as optimize
import astropy.units as u import astropy.units as u
from astropy import constants as const 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 from Lx_evo_and_flux import Lx_evo
from Lx_evo_and_flux import flux_at_planet_earth from Lx_evo_and_flux import flux_at_planet_earth
from Lx_evo_and_flux import L_xuv_all from Lx_evo_and_flux import L_xuv_all
...@@ -65,8 +66,8 @@ class planet_LoFo14(): ...@@ -65,8 +66,8 @@ class planet_LoFo14():
try: try:
# check if envelope mass fraction is specified, then Case 1 # check if envelope mass fraction is specified, then Case 1
planet_dict["fenv"] planet_dict["fenv"]
self.planet_info = "Case 1 - artificial planet \ self.planet_info = "Case 1 - artificial planet"\
(with M_core & f_env)" + " (with M_core & f_env)"
# Case 1: artificial planet with fenv & M_core given, need to # Case 1: artificial planet with fenv & M_core given, need to
# calculate the total mass and radius # calculate the total mass and radius
self.fenv = planet_dict["fenv"] self.fenv = planet_dict["fenv"]
...@@ -84,9 +85,9 @@ class planet_LoFo14(): ...@@ -84,9 +85,9 @@ class planet_LoFo14():
# Case 3: An observed planet with radius and mass # Case 3: An observed planet with radius and mass
# measurement, plus core mass is specified; # measurement, plus core mass is specified;
# need to calculate/estimate envelope mass frac. # need to calculate/estimate envelope mass frac.
self.planet_info = "Case 3 - obs. planet with radius & \ self.planet_info = "Case 3 - obs. planet with radius"\
mass measurement (and core mass \ + " & mass measurement (and core "\
estimate)" + "mass estimate)"
self.core_mass = planet_dict["core_mass"] self.core_mass = planet_dict["core_mass"]
self.mass = planet_dict["mass"] self.mass = planet_dict["mass"]
self.radius = planet_dict["radius"] self.radius = planet_dict["radius"]
...@@ -100,8 +101,8 @@ class planet_LoFo14(): ...@@ -100,8 +101,8 @@ class planet_LoFo14():
break break
except KeyError: except KeyError:
# no mass and fenv given -> Case 2 # no mass and fenv given -> Case 2
self.planet_info = "Case 2 - obs. planet with radius, \ self.planet_info = "Case 2 - obs. planet with"\
but no mass measurement" + " radius, but no mass measurement"
# Case 2 - observed planet with a without a mass, but # Case 2 - observed planet with a without a mass, but
# core mass estimate, need to calculate fenv # core mass estimate, need to calculate fenv
self.core_mass = planet_dict["core_mass"] self.core_mass = planet_dict["core_mass"]
...@@ -169,10 +170,14 @@ class planet_LoFo14(): ...@@ -169,10 +170,14 @@ class planet_LoFo14():
""" Function to set the right planet name based on the track specified. """ 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 This can then be used to check if a particular planet on a particular
track has already evolved (i.e. if outpufile exists). """ track has already evolved (i.e. if outpufile exists). """
self.planet_id = 'planet_a' + str(np.round(self.distance, 3)).replace('.', 'dot')\ self.planet_id = 'planet_a'\
+ '_Mcore' + str(np.round(self.core_mass, 3)).replace(".", "dot")\ + str(np.round(self.distance, 3)).replace('.', 'dot')\
+ '_fenv' + str(np.round(self.fenv, 3)) + '_' + self.metallicity\ + '_Mcore'\
+ '_Mstar' + str(np.round(self.mass_star, 3)).replace(".", "dot")\ + 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\ + "_K_" + K_on + "_beta_" + beta_on\
+ "_track_" + str(evo_track_dict["t_start"]) + "_"\ + "_track_" + str(evo_track_dict["t_start"]) + "_"\
+ str(evo_track_dict["t_sat"]) + "_" + str(t_final)\ + str(evo_track_dict["t_sat"]) + "_" + str(t_final)\
...@@ -202,13 +207,13 @@ class planet_LoFo14(): ...@@ -202,13 +207,13 @@ class planet_LoFo14():
# planet already exists # planet already exists
self.has_evolved = True self.has_evolved = True
else: else:
print("Planet: ", self.planet_id+".txt") #print("Planet: ", self.planet_id+".txt")
# Next, create a file (e.g. planet_XXXX.txt) which contains the # Next, create a file (e.g. planet_XXXX.txt) which contains the
# initial planet params # initial planet params
if os.path.exists(path_for_saving+planet_folder_id+".txt"): if os.path.exists(path_for_saving+planet_folder_id+".txt"):
pass pass
else: else:
p = open(path_for_saving+planet_folder_id+".txt", "a") p = open(path_for_saving+planet_folder_id+".txt", "w")
planet_params = "a,core_mass,fenv,mass,radius,metallicity,age\n"\ planet_params = "a,core_mass,fenv,mass,radius,metallicity,age\n"\
+ str(self.distance) + "," + str(self.core_mass)\ + str(self.distance) + "," + str(self.core_mass)\
+ "," + str(self.fenv) + "," + str(self.mass)\ + "," + str(self.fenv) + "," + str(self.mass)\
...@@ -224,9 +229,9 @@ class planet_LoFo14(): ...@@ -224,9 +229,9 @@ class planet_LoFo14():
pass pass
else: else:
t = open(path_for_saving+"track_params_" + t = open(path_for_saving+"track_params_" +
self.planet_id+".txt", "a") self.planet_id+".txt", "w")
track_params = "t_start,t_sat,t_curr,t_5Gyr,Lx_max,Lx_curr,\ track_params = "t_start,t_sat,t_curr,t_5Gyr,Lx_max,Lx_curr,"\
Lx_5Gyr,dt_drop,Lx_drop_factor\n"\ + "Lx_5Gyr,dt_drop,Lx_drop_factor\n"\
+ str(evo_track_dict["t_start"]) + ","\ + str(evo_track_dict["t_start"]) + ","\
+ str(evo_track_dict["t_sat"]) + "," \ + str(evo_track_dict["t_sat"]) + "," \
+ str(evo_track_dict["t_curr"]) + ","\ + str(evo_track_dict["t_curr"]) + ","\
...@@ -243,7 +248,7 @@ class planet_LoFo14(): ...@@ -243,7 +248,7 @@ class planet_LoFo14():
if os.path.exists(path_for_saving+"host_star_properties.txt"): if os.path.exists(path_for_saving+"host_star_properties.txt"):
pass pass
else: else:
s = open(path_for_saving+"host_star_properties.txt", "a") s = open(path_for_saving+"host_star_properties.txt", "w")
star_params = "star_id,mass_star,radius_star,age,Lbol,Lx_age\n"\ star_params = "star_id,mass_star,radius_star,age,Lbol,Lx_age\n"\
+ self.star_id + "," + str(self.mass_star) + ","\ + self.star_id + "," + str(self.mass_star) + ","\
+ str(self.radius_star) + "," + str(self.age)\ + str(self.radius_star) + "," + str(self.age)\
...@@ -267,11 +272,11 @@ class planet_LoFo14(): ...@@ -267,11 +272,11 @@ class planet_LoFo14():
df.to_csv(path_for_saving+self.planet_id+".txt", index=None) df.to_csv(path_for_saving+self.planet_id+".txt", index=None)
# create another file, which contains the final parameters only # create another file, which contains the final parameters only
if os.path.exists(path_for_saving+planet_folder_id+"_"+\ if os.path.exists(path_for_saving+\
self.planet_id+"_final.txt"): self.planet_id+"_final.txt"):
pass pass
else: else:
p = open(path_for_saving+self.planet_id+"_final.txt", "a") p = open(path_for_saving+self.planet_id+"_final.txt", "w")
# get index of last valid entry in the radius array and access # get index of last valid entry in the radius array and access
# its value # its value
index_of_last_entry = df["Radius"][df["Radius"].notna()].index[-1] index_of_last_entry = df["Radius"][df["Radius"].notna()].index[-1]
...@@ -307,4 +312,4 @@ class planet_LoFo14(): ...@@ -307,4 +312,4 @@ class planet_LoFo14():
def evolve_backwards(self, t_final, initial_step_size, epsilon, K_on, def evolve_backwards(self, t_final, initial_step_size, epsilon, K_on,
beta_on, evo_track_dict, path_for_saving): beta_on, evo_track_dict, path_for_saving):
print("Not implemented yet.") print("Not implemented yet.")
continue
...@@ -20,12 +20,12 @@ def calculate_mass_planet_Ot19(R_pl): ...@@ -20,12 +20,12 @@ def calculate_mass_planet_Ot19(R_pl):
rho_volatile = (M_p_volatile*M_earth) \ rho_volatile = (M_p_volatile*M_earth) \
/ (4/3*np.pi*(R_pl*R_earth)**3) / (4/3*np.pi*(R_pl*R_earth)**3)
if (rho_volatile >= 3.3): if (rho_volatile >= 3.3):
raise Exception("Planet with this radius is too small and likely \ raise Exception("Planet with this radius is too small and"\
rocky; use LoFo14 models instead.") + " likely rocky; use LoFo14 models instead.")
else: else:
if (M_p_volatile >= 120): if (M_p_volatile >= 120):
raise Exception("Planet too massive. M-R relation only valid \ raise Exception("Planet too massive. M-R relation only valid"\
for <120 M_earth.") + " for <120 M_earth.")
else: else:
return M_p_volatile return M_p_volatile
...@@ -61,12 +61,12 @@ def calculate_radius_planet_Ot19(M_pl): ...@@ -61,12 +61,12 @@ def calculate_radius_planet_Ot19(M_pl):
rho_volatile = (M_pl*M_earth) \ rho_volatile = (M_pl*M_earth) \
/ (4/3*np.pi*(R_p_volatile*R_earth)**3) / (4/3*np.pi*(R_p_volatile*R_earth)**3)