Commit 8f1cd3df authored by Laura Ketzer's avatar Laura Ketzer
Browse files

started updating Mass_loss_evolution_function to be compliant with PEP8

parent d27803b2
......@@ -310,3 +310,435 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object,
#print("Done!")
return np.array(t_arr), np.array(M_arr), np.array(R_arr), np.array(Lx_arr)
def mass_planet_RK4_forward_LO14_PAPER(epsilon, K_on, beta_on, planet_object, initial_step_size, t_final, track_dict):
"""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
taking into account photoevaporative mass loss.
Input:
----------
epsilon: evaporation efficiency
K_on: set use of K parameter on or off
beta_on: set use of beta parameter on or off
planet_object: object of planet class which contains also stellar parameters
and info about stellar evo track
step_size: initial step_size, fixed
t_final: final time of simulation
track_dict: dictionary with Lx evolutionary track parameters
Output:
----------
t_arr, M_arr, R_arr, Lx_arr: array of time, mass, radius and Lx values from t_start to t_final
"""
M_earth = const.M_earth.cgs.value
R_earth = const.R_earth.cgs.value
Myr_to_sec = 1e6*365*86400
# initialize the starting values for Lxuv(t_start), mass, density, beta, K
Lx0 = Lx_evo(t=track_dict["t_start"], track_dict=track_dict)
Lxuv0 = L_xuv_all(Lx0) # use Sanz-Forcada2010 scaling law to get total XUV luminosity
Fxuv0 = flux_at_planet(Lxuv0, planet_object.distance) # get flux at orbital separation a_p
# initial planet parameters at t_start
f_env_0 = f_env = planet_object.fenv
R0 = R = planet_object.radius # should match observed radius - determined by M_core and f_env
M0 = M = planet_object.mass
rho0 = rho = plmo14.density_planet(M0, R0) # initial mean density
M_env0 = M_env = M0 - planet_object.core_mass # initial envelope mass
M_core = planet_object.core_mass
# specify beta0 and K0
if beta_on == "yes": # use the beta estimate from Salz et al. 2016
beta = beta0 = bk.beta_fct_LO14(M0, Fxuv0, R0)
elif beta_on == "no":
beta = beta0 = 1.
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)
elif K_on == "no":
K = K0 = 1.
t_start = track_dict["t_start"]
t_max = t_final
step_size = initial_step_size
# create time array for integration (with user-specified step size)
number = math.ceil((t_max-t_start)/step_size)
times, step_size2 = np.linspace(t_start, t_max, number, endpoint=True, retstep=True)
#print('stepsize=', step_size2)
# here I make lists of all the values I would like to track & output in the end:
M_arr = [0]*len(times)
M_arr[0] = M = M0 # inital value for Mdot at t_curr
R_arr = [0]*len(times)
R_arr[0] = R = R0
t_arr = [0]*len(times)
t_arr[0] = t_start # inital value for t
Lx_arr = [0]*len(times)
#Lx_arr[0] = Lx0
# CRITERION when to stop the mass-loss
# for the Lopez planets I have a specified core mass and thus a fixed core radius (bare rocky core)
R_core = planet_object.core_radius # stop when this radius is reached!
dt = step_size2
#for i in tqdm(range(1, len(times))):
for i in range(0, len(times)-1):
#print(t_arr[i])
################
# this is just for me to return the Lx(t) evolution to check if it is correct (not required since the Lx(t)
# calculation is embedded in the mass_loss_rate_fancy function)
Lx_i = Lx_evo(t=t_arr[i], track_dict=track_dict)
###############################################################################################################
Mdot1 = mass_loss_rate_forward_LO14(times[i], epsilon, K_on, beta_on, planet_object, f_env, track_dict)
k1 = (dt*Myr_to_sec * Mdot1)/M_earth # mass lost in one timestep in earth masses
M_05k1 = M + 0.5*k1
M_env_05k1 = M_05k1 - M_core
f_env_05k1 = (M_env_05k1/M_05k1)*100 # new mass fraction
Mdot2 = mass_loss_rate_forward_LO14(times[i]+0.5*dt, epsilon, K_on, beta_on, planet_object, f_env_05k1, track_dict)
k2 = (dt*Myr_to_sec * Mdot2)/M_earth
M_05k2 = M + 0.5*k2
M_env_05k2 = M_05k2 - M_core
f_env_05k2 = (M_env_05k2/M_05k2)*100 # new mass fraction
Mdot3 = mass_loss_rate_forward_LO14(times[i]+0.5*dt, epsilon, K_on, beta_on, planet_object, f_env_05k2, track_dict)
k3 = (dt*Myr_to_sec * Mdot3)/M_earth
M_k3 = M + k3
M_env_k3 = M_k3 - M_core
f_env_k3 = (M_env_k3/M_k3)*100 # new mass fraction
Mdot4 = mass_loss_rate_forward_LO14(times[i]+dt, epsilon, K_on, beta_on, planet_object, f_env_k3, track_dict)
k4 = (dt*Myr_to_sec * Mdot4)/M_earth
###############################################################################################################
M_lost = (k1 + 2*k2 + 2*k3 + k4)/6. # after time-step dt
###############################################################################################################
# check if planet with new mass does still have some atmosphere:
if ((M + M_lost) - M_core) >= 0:
# then planet still has some atmosphere left -> continue
M_arr[i+1] = M = M + M_lost # new planet mass (assume envelope mass is lost)
M_env = M - M_core # new envelope mass
#M_arr[i] = M = M + (k1 + 2*k2 + 2*k3 + k4)/6. # M_i_plus_1 - update Mass value
t_arr[i+1] = t = t_arr[i] + dt #t_start + i*dt # t_i_plus_1 - update time value
# new envelope mass fraction:
f_env = (M_env/M)*100 # in %
# 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,
planet_object.distance), planet_object.metallicity)
else:
# all atmosphere is gone -> terminate
print("Atmosphere has evaportated! Only bare rocky core left! STOP this madness!")
# if the stop criterium is reached, I add the core mass & core radius to the array at t_i+1
t_arr = np.trim_zeros(t_arr)
print(t_arr[-1]+dt)
#print(t_arr[i]+dt)
t_arr = np.append(np.array(t_arr), t_arr[-1]+dt)
M_arr = np.trim_zeros(M_arr)
#M_arr = np.append(np.array([mass.value for mass in M_arr]), M_core.value)*u.g
M_arr = np.append(np.array(M_arr), M_core)
R_arr = np.trim_zeros(R_arr)
#R_arr = np.append(np.array([radius.value for radius in R_arr]), R_core.value)*u.cm
R_arr = np.append(np.array(R_arr), R_core)
Lx_arr = np.trim_zeros(Lx_arr)
Lx_arr = np.array(Lx_arr.append(Lx_evo(t=t_arr[-1]+dt, track_dict=track_dict)))
return t_arr, M_arr, R_arr, Lx_arr
Lx_arr[i+1] = Lx_evo(t=t_arr[-1]+dt, track_dict=track_dict)
Lx_arr = np.array(Lx_arr)
t_arr = np.array(t_arr)
M_arr = np.array(M_arr)
R_arr = np.array(R_arr)
print("Done!")
return t_arr, M_arr, R_arr, Lx_arr
def mass_planet_RK4_forward_Ot14(epsilon, K_on, beta_on, planet_object, initial_step_size, t_final, track_dict):
"""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
taking into account photoevaporative mass loss.
Input:
----------
epsilon: evaporation efficiency
K_on: set use of K parameter on or off
beta_on: set use of beta parameter on or off
planet_object: object of planet class which contains also stellar parameters
and info about stellar evo track
step_size: initial step_size, variable
[NOTE: the implementation of a variable step size is somewhat preliminary. The step size is adjusted
(made smaller or bigger depending how fast or slow the mass/radius changes) until the final time step
greater than t_final. This means that if the step size in the end is e.g. 10 Myr, and the integration
is at 4999 Myr, then last time entry will be 4999+10 -> 5009 Myr.]
t_final: final time of simulation
track_dict: dictionary with Lx evolutionary track parameters
Output:
----------
t_arr, M_arr, R_arr, Lx_arr: array of time, mass, radius and Lx values from t_start to t_final
"""
M_earth = const.M_earth.cgs.value
R_earth = const.R_earth.cgs.value
Myr_to_sec = 1e6*365*86400
# initialize the starting values for Lxuv(t_start), mass, density, beta, K
Lx0 = Lx_evo(t=track_dict["t_start"], track_dict=track_dict)
Lxuv0 = L_xuv_all(Lx0) # use scaling law to get total XUV luminosity
Fxuv0 = flux_at_planet(Lxuv0, planet_object.distance) # get flux at orbital separation a_p
# "make" initial planet at t_start
R0 = R = planet_object.radius # should match observed radius
M0 = M = planet_object.mass # planetary mass estimate based on mass-radius relation
rho0 = rho = plmoOt20.density_planet(M0, R0) # initial approx. density
# specify beta0 and K0
if beta_on == "yes": # use the beta estimate from Salz et al. 2016
beta = beta0 = bk.beta_fct_LO14(M0, Fxuv0, R0)
elif beta_on == "no":
beta = beta0 = 1.
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)
elif K_on == "no":
K = K0 = 1.
M_arr = []
M_arr.append(M0)
R_arr = []
R_arr.append(R0)
t_arr = []
t0 = t = track_dict["t_start"]
t_arr.append(t0)
Lx_arr = []
Lx_arr.append(Lx0)
#print(M0, R0, t0)
# CRITERION when to stop the mass-loss
R_core = 2.15 # stop when this radius is reached! (this is the minimum radius for which the volatile regime is valid)
M_core = plmoOt20.calculate_mass_planet_Ot19(R_core)
dt = initial_step_size
i = 1 # counter
while t <= t_final:
#print(i, ' - ', t, "- dt = ", dt)
################
# this is just for me to return the Lx(t) evolution to check if it is correct (not required since the Lx(t)
# calculation is embedded in the mass_loss_rate_fancy function)
Lx_i = Lx_evo(t=t, track_dict=track_dict)
################
# for first step use the parameters initialized above (beta=beta0, K=K0, M=M0)
###############################################################################################################
Mdot1 = mass_loss_rate_forward_Ot20(t, epsilon, K_on, beta_on, planet_object, R, track_dict) # mass M, radius R
k1 = (dt*Myr_to_sec * Mdot1)/M_earth # mass lost in one timestep in earth masses
M_05k1 = M + 0.5*k1
R_05k1 = plmoOt20.calculate_radius_planet_Ot19(M_05k1)
Mdot2 = mass_loss_rate_forward_Ot20(t+0.5*dt, epsilon, K_on, beta_on, planet_object, R_05k1, track_dict)
k2 = (dt*Myr_to_sec * Mdot2)/M_earth
M_05k2 = M + 0.5*k2
R_05k2 = plmoOt20.calculate_radius_planet_Ot19(M_05k2)
Mdot3 = mass_loss_rate_forward_Ot20(t+0.5*dt, epsilon, K_on, beta_on, planet_object, R_05k2, track_dict)
k3 = (dt*Myr_to_sec * Mdot3)/M_earth
M_05k3 = M + 0.5*k3
R_05k3 = plmoOt20.calculate_radius_planet_Ot19(M_05k3)
Mdot4 = mass_loss_rate_forward_Ot20(t+dt, epsilon, K_on, beta_on, planet_object, R_05k3, track_dict)
k4 = (dt*Myr_to_sec * Mdot4)/M_earth
###############################################################################
M_lost = (k1 + 2*k2 + 2*k3 + k4)/6. # after time-step dt
###############################################################################
# this part is new compared to the one used in the PAPER
###############################################################################
# if radius change is too drastic, decrease step size
M_new = M+M_lost
R_new = plmoOt20.calculate_radius_planet_Ot19(M_new)
# adjust step size, if radius change is too drastic or too little
if (abs((R-R_new)/R)*100 >= 1.) & (t < track_dict["t_curr"]) & (dt > 1e-2):#10**(-1)): # smaller than 1%
dt = dt/10.
#print("radius change: ", abs((R-R_new)/R)*100)
#print('make step size smaller: ', dt)
# don't write anything to file, do iteration again with new step size
elif (abs((R-R_new)/R)*100 < (0.01)) & (t < track_dict["t_curr"]) & (dt < 10.):
dt = dt*10.
#print('make step size bigger: ', dt)
elif (abs((R-R_new)/R)*100 >= 1) & (t >= track_dict["t_curr"]) & (dt > 1e-1): # smaller than 1%
dt = dt/10.
#print('make step size smaller: ', dt)
elif (abs((R-R_new)/R)*100 < (0.01)) & (t >= track_dict["t_curr"]) & (dt < 10.):
dt = dt*10
#print('make step size bigger: ', dt)
else:#
if ((M + M_lost) - M_core) >= 0:
# then planet still has some atmosphere left -> continue
M = M + M_lost # new planet mass (M_lost is negative)
M_arr.append(M)
t = t_arr[-1] + dt #t_start + i*dt # t_i_plus_1 - update time value
t_arr.append(t) # new time t
Lx_arr.append(Lx_evo(t=t, track_dict=track_dict)) # Lx at new time t
# calculate new radius with new planet mass & one time step later
R = plmoOt20.calculate_radius_planet_Ot19(M)
R_arr.append(R)
i = i+1 # update step to i+1
else:
# all atmosphere is gone -> terminate
print("Atmosphere has evaportated! Only bare rocky core left! STOP this madness!")
# if the stop criterium is reached, I add the core mass & core radius to the array at t_i+1
print("t = ", t_arr[-1]+dt)
t_arr.append(t_arr[-1]+dt)
M_arr.append(M_core)
R_arr.append(R_core)
Lx_arr.append(Lx_evo(t=t_arr[-1]+dt, track_dict=track_dict))
return np.array(t_arr), np.array(M_arr), np.array(R_arr), np.array(Lx_arr)
print("Done!")
return np.array(t_arr), np.array(M_arr), np.array(R_arr), np.array(Lx_arr)
def mass_planet_RK4_forward_Ot14_PAPER(epsilon, K_on, beta_on, planet_object, initial_step_size, t_final, track_dict):
"""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
taking into account photoevaporative mass loss.
Input:
----------
epsilon: evaporation efficiency
K_on: set use of K parameter on or off
beta_on: set use of beta parameter on or off
planet_object: object of planet class which contains also stellar parameters
and info about stellar evo track
step_size: initial step_size, fixed
t_final: final time of simulation
track_dict: dictionary with Lx evolutionary track parameters
Output:
----------
t_arr, M_arr, R_arr, Lx_arr: array of time, mass, radius and Lx values from t_start to t_final
"""
M_earth = const.M_earth.cgs.value
R_earth = const.R_earth.cgs.value
Myr_to_sec = 1e6*365*86400
# initialize the starting values for Lxuv(t_start), mass, density, beta, K
Lx0 = Lx_evo(t=track_dict["t_start"], track_dict=track_dict)
Lxuv0 = L_xuv_all(Lx0) # use scaling law to get total XUV luminosity
Fxuv0 = flux_at_planet(Lxuv0, planet_object.distance) # get flux at orbital separation a_p
# "make" initial planet at t_start
R0 = R = planet_object.radius # should match observed radius
M0 = M = planet_object.mass # planetary mass estimate based on mass-radius relation
rho0 = rho = plmoOt20.density_planet(M0, R0) # initial approx. density
# specify beta0 and K0
if beta_on == "yes": # use the beta estimate from Salz et al. 2016
beta = beta0 = bk.beta_fct_LO14(M0, Fxuv0, R0)
elif beta_on == "no":
beta = beta0 = 1.
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)
elif K_on == "no":
K = K0 = 1.
# create time array for integration (with user-specified step size)
t_start, t_max = track_dict["t_start"], t_final
step_size = initial_step_size
number = math.ceil((t_max-t_start)/step_size)
times, step_size2 = np.linspace(t_start, t_max, number, endpoint=True, retstep=True)
#print('stepsize=', step_size2)
# make lists of all the values wwe want to track & output in the end:
M_arr = [0]*len(times)
M_arr[0] = M = M0 # inital value for Mdot at t_curr
R_arr = [0]*len(times)
R_arr[0] = R = R0
t_arr = [0]*len(times)
t_arr[0] = t_start # inital value for t
Lx_arr = [0]*len(times)
#Lx_arr[0] = Lx0
# CRITERION when to stop the mass-loss
R_core = 2.15 # stop when this radius is reached! (this is the minimum radius
# for which the volatile regime is valid)
M_core = plmoOt20.calculate_mass_planet_Ot19(R_core)
dt = step_size2
for i in range(0, len(times)-1):
# this is just for me to return the Lx(t) evolution to check if it is correct
# (not required since the Lx(t) calculation is embedded in
# the mass_loss_rate_fancy function)
Lx_i = Lx_evo(t=t_arr[i], track_dict=track_dict)
###############################################################################################################
Mdot1 = mass_loss_rate_forward_Ot20(times[i], epsilon, K_on, beta_on, planet_object, R, track_dict) # mass M, radius R
k1 = (dt*Myr_to_sec * Mdot1)/M_earth # mass lost in one timestep in earth masses
M_05k1 = M + 0.5*k1
R_05k1 = plmoOt20.calculate_radius_planet_Ot19(M_05k1)
Mdot2 = mass_loss_rate_forward_Ot20(times[i]+0.5*dt, epsilon, K_on, beta_on, planet_object, R_05k1, track_dict)
k2 = (dt*Myr_to_sec * Mdot2)/M_earth
M_05k2 = M + 0.5*k2
R_05k2 = plmoOt20.calculate_radius_planet_Ot19(M_05k2)
Mdot3 = mass_loss_rate_forward_Ot20(times[i]+0.5*dt, epsilon, K_on, beta_on, planet_object, R_05k2, track_dict)
k3 = (dt*Myr_to_sec * Mdot3)/M_earth
M_05k3 = M + 0.5*k3
R_05k3 = plmoOt20.calculate_radius_planet_Ot19(M_05k3)
Mdot4 = mass_loss_rate_forward_Ot20(times[i]+dt, epsilon, K_on, beta_on, planet_object, R_05k3, track_dict)
k4 = (dt*Myr_to_sec * Mdot4)/M_earth
###############################################################################################################
M_lost = (k1 + 2*k2 + 2*k3 + k4)/6. # mass lost after time-step dt
###############################################################################################################
# check if planet with new mass does still have some atmosphere
if ((M + M_lost) - M_core) >= 0:
# then planet still has some atmosphere left -> continue
M_arr[i+1] = M = M + M_lost # new planet mass (M_lost is negative)
t_arr[i+1] = t = t_arr[i] + dt #t_start + i*dt # t_i_plus_1 - update time value
# calculate new radius with new planet mass
R_arr[i+1] = R = plmoOt20.calculate_radius_planet_Ot19(M)
else:
# all atmosphere is gone (based on criterion set at the top)-> terminate
print("Atmosphere has evaportated! Only bare rocky core left! STOP this madness!")
# if the stop criterium is reached, I add the core mass & core radius to the array at t_i+1
t_arr = np.trim_zeros(t_arr)
print(t_arr[-1]+dt)
t_arr = np.append(np.array(t_arr), t_arr[-1]+dt)
M_arr = np.trim_zeros(M_arr)
M_arr = np.append(np.array(M_arr), M_core)
R_arr = np.trim_zeros(R_arr)
R_arr = np.append(np.array(R_arr), R_core)
Lx_arr = np.trim_zeros(Lx_arr)
Lx_arr = np.array(Lx_arr.append(Lx_evo(t=t_arr[-1]+dt, track_dict=track_dict)))
return t_arr, M_arr, R_arr, Lx_arr
# if planet survives, output the final arrays
Lx_arr[i+1] = Lx_evo(t=t_arr[-1]+dt, track_dict=track_dict)
Lx_arr = np.array(Lx_arr)
t_arr = np.array(t_arr)
M_arr = np.array(M_arr)
R_arr = np.array(R_arr)
print("Done!")
return t_arr, M_arr, R_arr, Lx_arr
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