# highest level 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 import Planet_model_Ot20 as plmoOt20 import Beta_K_functions as bk import astropy.units as u from astropy import constants as const from Mass_loss_rate_function import mass_loss_rate_forward_LO14, mass_loss_rate_forward_Ot20 import numpy as np import math def mass_planet_RK4_forward_LO14(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 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. 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) # 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 = 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_LO14(t, 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(t+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(t+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(t+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 ########################################################################################### # this part is new compared to the one used in the PAPER (there we used a fixed step size!) ########################################################################################### # if radius change is too drastic, decrease step size M_new = M + M_lost M_env_new = M_new - M_core f_env_new = (M_env_new/M_new)*100 # in % R_new = plmo14.calculate_planet_radius(M_core, f_env_new, t, flux_at_planet_earth(planet_object.Lbol, planet_object.distance), planet_object.metallicity) #print("time = ", t) #print(R, R_new, f_env, f_env_new) #print(abs((R-R_new)/R)*100) # 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)): # larger than 1%, make step size smaller dt = dt/10. #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.): # smaller than 0.01%, make step size bigger, but not bigger than 10 Myr dt = dt*10. #print('make step size bigger: ', dt) elif (abs((R-R_new)/R)*100 >= 1.) & (t >= track_dict["t_curr"]) & (dt > 1e-2): # larger than 1%, make step size smaller dt = dt/10. elif (abs((R-R_new)/R)*100 < (0.01)) & (t >= track_dict["t_curr"]) & (dt < 10.): # smaller than 0.01%, make step size bigger, but not bigger than 10 Myr dt = dt*10 else: # if radius change is ok, check if there is still envelope left if ((M + M_lost) - M_core) >= 0: #print("envelope left") # then planet still has some atmosphere left -> continue M = M + M_lost # new planet mass (M_lost is negative) M_arr.append(M) M_env = M - M_core # new envelope mass #print("(t_arr[-1]", t_arr[-1]) t = t_arr[-1] + dt #t_start + i*dt # t_i_plus_1 - update time value # if t > track_dict["t_curr"]: # print(dt, "new time: ", t) t_arr.append(t) # new time t Lx_arr.append(Lx_evo(t=t, track_dict=track_dict)) # Lx at new time t # new envelope mass fraction: f_env = (M_env/M)*100 # in % #print("f_env = ", f_env) # calculate new radius with new planet mass/envelope mass fraction & one time step later R = plmo14.calculate_planet_radius(M_core, f_env, t, flux_at_planet_earth(planet_object.Lbol, planet_object.distance), planet_object.metallicity) 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!") # print("params=", t, M, f_env, R) # print("Lx", Lx_evo(t=t_arr[-1]+dt, track_dict=track_dict)) # 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) print("M_core = ", M_core) R_arr.append(R_core) print("R_core = ", 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_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