Mass_evolution_function.py 30.8 KB
Newer Older
1
2
3
4
5
6
7
8
# 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
9
from Mass_loss_rate_function import mass_loss_rate_forward_LO14, mass_loss_rate_forward_Ot20
10
11
12
13
import numpy as np
import math


14

15
def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, initial_step_size, t_final, track_dict):
16
17
18
19
20
    """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:
21
    ----------
22
23
24
25
26
27
28
29
30
31
    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.]
32
    t_final: final time of simulation
33
34
35
36
37
    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
38
    """ 
39
    
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
    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

111
        ###########################################################################################
112
113
        M_lost = (k1 + 2*k2 + 2*k3 + k4)/6. # after time-step dt
        ###########################################################################################
114
        
115
116
117
118
119
120
121
122
123
        # 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)
124
125
126
        #print("time = ", t)
        #print(R, R_new, f_env, f_env_new)
        #print(abs((R-R_new)/R)*100)
127
128
        
        # adjust step size, if radius change is too drastic or too little
129
        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
130
131
132
            dt = dt/10.
            #print('make step size smaller: ', dt)
            # don't write anything to file, do iteration again with new step size
133
        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
134
135
136
            dt = dt*10.
            #print('make step size bigger: ', dt)
        
137
        elif (abs((R-R_new)/R)*100 >= 1.) & (t >= track_dict["t_curr"]) & (dt > 1e-2): # larger than 1%, make step size smaller
138
139
            dt = dt/10.
        
140
        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
141
142
143
144
            dt = dt*10
        
        else: # if radius change is ok, check if there is still envelope left
            if ((M + M_lost) - M_core) >= 0:
145
                #print("envelope left")
146
147
148
149
                # 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
150
                #print("(t_arr[-1]", t_arr[-1])
151
                t = t_arr[-1] + dt #t_start + i*dt # t_i_plus_1 - update time value
152
153
#                 if t > track_dict["t_curr"]:
#                     print(dt, "new time: ", t)
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
                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!")
170
171
172
#                 print("params=", t, M, f_env, R)
#                 print("Lx", Lx_evo(t=t_arr[-1]+dt, track_dict=track_dict))
                
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
                # 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):
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
    """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
    """ 
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
    
    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)
247
    #print('stepsize=', step_size2)
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
    
    # 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
        
294
        ###############################################################################################################
295
296
        M_lost = (k1 + 2*k2 + 2*k3 + k4)/6. # after time-step dt 
        ###############################################################################################################
297
298
299
300
301
302
303
304
305
306
307
308
309
        
        # 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,
310
                                                            planet_object.distance), planet_object.metallicity)
311
312
313
314
315
316
317
318
            
        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)
319
            #print(t_arr[i]+dt)
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
            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):
345
346
347
348
349
    """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:
350
    ----------
351
352
353
354
355
356
357
358
359
360
    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.]
361
    t_final: final time of simulation
362
363
364
365
366
    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
367
    """ 
368
    
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
    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

438
        ###############################################################################
439
440
        M_lost = (k1 + 2*k2 + 2*k3 + k4)/6. # after time-step dt
        ###############################################################################
441
        
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
        # 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)


##################################################################################################################################
##################################################################################################################################
##################################################################################################################################


501
def mass_planet_RK4_forward_Ot14_PAPER(epsilon, K_on, beta_on, planet_object, initial_step_size, t_final, track_dict):
502
503
504
505
506
    """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:
507
    ----------
508
509
510
511
512
513
    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
514
    t_final: final time of simulation
515
    track_dict: dictionary with Lx evolutionary track parameters
516
    
517
518
519
520
    Output:
    ----------
    t_arr, M_arr, R_arr, Lx_arr: array of time, mass, radius and Lx values from t_start to t_final
    """ 
521
    
522
523
524
525
526
527
528
529
530
531
532
533
    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
534
    rho0 = rho = plmoOt20.density_planet(M0, R0)  # initial approx. density
535
536
537
538
539
540
541
542
543
544
545
    
    # 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.

546
547
548
549
550
551
    # 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)
552
    
553
554
555
556
557
558
559
560
561
    # 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
562
563
    
    # CRITERION when to stop the mass-loss
564
565
566
    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)
567
    
568
569
570
571
572
573
    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)
574
575
        
        ###############################################################################################################
576
        Mdot1 = mass_loss_rate_forward_Ot20(times[i], epsilon, K_on, beta_on, planet_object, R, track_dict) # mass M, radius R
577
        k1 = (dt*Myr_to_sec * Mdot1)/M_earth # mass lost in one timestep in earth masses
578
579
        M_05k1 = M + 0.5*k1     
        R_05k1 = plmoOt20.calculate_radius_planet_Ot19(M_05k1)
580
        
581
        Mdot2 = mass_loss_rate_forward_Ot20(times[i]+0.5*dt, epsilon, K_on, beta_on, planet_object, R_05k1, track_dict)
582
583
        k2 = (dt*Myr_to_sec * Mdot2)/M_earth
        M_05k2 = M + 0.5*k2
584
        R_05k2 = plmoOt20.calculate_radius_planet_Ot19(M_05k2)
585
        
586
        Mdot3 = mass_loss_rate_forward_Ot20(times[i]+0.5*dt, epsilon, K_on, beta_on, planet_object, R_05k2, track_dict) 
587
588
        k3 = (dt*Myr_to_sec * Mdot3)/M_earth
        M_05k3 = M + 0.5*k3
589
        R_05k3 = plmoOt20.calculate_radius_planet_Ot19(M_05k3)
590
        
591
        Mdot4 = mass_loss_rate_forward_Ot20(times[i]+dt, epsilon, K_on, beta_on, planet_object, R_05k3, track_dict) 
592
593
        k4 = (dt*Myr_to_sec * Mdot4)/M_earth

594
        ###############################################################################################################
595
596
597
598
599
600
        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
601
            
602
603
604
605
606
607
608
609
            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!")
610

611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
            # 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)
629
    print("Done!")
630
631
    return t_arr, M_arr, R_arr, Lx_arr