Lx_evo_and_flux.py 10.1 KB
 Laura Ketzer committed May 07, 2020 1 2 3 ``````import numpy as np import astropy.units as u from astropy import constants as const `````` 4 5 ``````import scipy.optimize as optimize from scipy.optimize import fsolve `````` Laura Ketzer committed May 07, 2020 6 7 8 9 10 11 12 `````` def Lx_evo(t, track_dict): """ Function to calculate the X-ray luminosity, Lx, of the host star at a given time t. NOTE: I tried to set up this function such that it works going forward (for young system) and backwards (for older system) in time; `````` Laura Ketzer committed May 08, 2020 13 `````` -> so t_curr(ent) is a little confusing! -> we have two cases!! `````` Laura Ketzer committed May 07, 2020 14 15 16 17 18 19 20 21 `````` Parameters Case 1 - a young system which you want to evolve forward in time (i.e. estimate the mass loss the planet will undergo in the future): --------------------------------------------------------------------------- @param t_start: (float) Starting time for the forward simulation in Myr (e.g. 23 Myr for V1298 Tau) @param Lx_max: (float) Measured Lx value at the system age (i.e. t_start) `````` Laura Ketzer committed May 08, 2020 22 `````` @param t_curr: (float) `````` Laura Ketzer committed May 07, 2020 23 `````` Set to 1 Gyr (with Lx_curr = Lx at 1 Gyr), which is where, by design, the three `````` Laura Ketzer committed May 08, 2020 24 25 26 `````` activity tracks converge (based on Tu et al. (2015)/ Johnstone et al.(2015) models); @param Lx_curr: (float) Lx value at 1 Gyr taken from Tu et al. (2015) `````` Laura Ketzer committed May 07, 2020 27 28 29 30 31 32 33 `````` @param t_5Gyr: (float) Defines, together with Lx_5Gyr, the slope common for all three activity evolution tracks past 1 Gyr (also based on Tu et al.(2015)); @param Lx_5Gyr: (float) Defines, together with t_5Gyr, the slope common for all three activity evolution tracks past 1 Gyr (also based on Tu et al.(2015)); `````` Laura Ketzer committed May 08, 2020 34 `````` # IGNORE FOR NOW CASE 2 - not implemented yet `````` Laura Ketzer committed May 07, 2020 35 36 37 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 `````` Parameters Case 2 - an older system, which you want to evolve backwards in time (i.e. estimate the mass loss the planet has undergone until now): --------------------------------------------------------------------------- @param t_start is the youngest age you want to calculate backwards (i.e. this shouldbe sth. like ~10 Myr, or close to the disk clearing timescale); @param Lx_max: (float) Some saturation X-ray luminosity (e.g. Lx,sat/Lbol ~ 10^(-3)), @param t_curr: (float) Needs to be set to the system age (e.g. t_curr=700 Myr for K2-136b); @param Lx_curr: (float) If a measured Lx value is available for the system, otherwise need to estimate; @param t_5Gyr: (float) Defines, together with Lx_5Gyr, the slope common for all three activity evolution tracks past 1 Gyr (also based on Tu et al.(2015)); @param Lx_5Gyr: (float) Defines, together with t_5Gyr, the slope common for all three activity evolution tracks past 1 Gyr (also based on Tu et al.(2015)); NOTE: t_5Gyr (and Lx_5Gyr) could be changed to create whatever slope desired (for a future evolution). Sidenote: this function is not very pretty, needs improvements... @return Lx_at_time_t: (float) Pass one t-value, and get the corresponding Lx """ # read in the parameters from the provided dictionary -> these are all the parameters required to define the tracks t_start, t_sat, t_curr, t_5Gyr = track_dict["t_start"], track_dict["t_sat"], track_dict["t_curr"], track_dict["t_5Gyr"] Lx_max, Lx_curr, Lx_5Gyr = track_dict["Lx_max"], track_dict["Lx_curr"], track_dict["Lx_5Gyr"] dt_drop, Lx_drop_factor = track_dict["dt_drop"], track_dict["Lx_drop_factor"] ######################################################################################################################## # Define function for calculating a power law powerlaw = lambda x, k, index: k * (x**index) # t_curr has to be greater than t_start if (t_curr < t_start) or (t < t_start): if t <= t_start: # this allows me (for the V1298 Tau system) to produce a saturation regime before 23 Myr with Lx_max Lx = Lx_max else: print("Make sure t_curr > t_start, t >= t_start") elif t > t_curr: `````` Laura Ketzer committed Jun 05, 2020 77 `````` # this is the regime past 1 Gyr (the same for all Tu tracks) `````` Laura Ketzer committed May 07, 2020 78 `````` # I use the Tu 2015 slope which is approx the same for all three evolutionary tracks (and based on current solar value); `````` Laura Ketzer committed May 08, 2020 79 `````` # this is mainly for consistency since we approximate our tracks based on the Tu-tracks `````` Laura Ketzer committed Jun 05, 2020 80 81 82 83 84 85 86 87 88 89 90 91 92 93 `````` # exception for the OwWu17 case, where the slope past t_sat is constant for all t > t_sat if t_5Gyr == t_curr: # then we are dealing with the OwWu17 case #print("in if") alpha = (np.log10(Lx_curr/Lx_max))/(np.log10(t_curr/t_sat)) k = 10**(np.log10(Lx_max) - alpha*np.log10(t_sat)) Lx = powerlaw(t, k, alpha) else: # normal case, i.e. slope past t_curr given by input parameters alpha = (np.log10(Lx_5Gyr/Lx_curr))/(np.log10(t_5Gyr/t_curr)) print(alpha) k = 10**(np.log10(Lx_curr) - alpha*np.log10(t_curr)) Lx = powerlaw(t, k, alpha) `````` Laura Ketzer committed May 07, 2020 94 95 96 97 98 99 100 101 102 103 104 105 `````` return Lx elif t_curr > t_start: ################################################################################### # dt_drop==0 means we create a track with only the saturation regime and # a drop to the converging Lx at 1 Gyr ################################################################################### if dt_drop==0: # then t_sat == t_drop if t_start >= t_sat: # then t_sat <= t_start < t_curr if t >= t_sat: #and t <= t_curr: # only data in power law regime alpha = (np.log10(Lx_curr/Lx_max))/(np.log10(t_curr/t_sat)) `````` Laura Ketzer committed Jun 05, 2020 106 `````` print(alpha) `````` Laura Ketzer committed May 07, 2020 107 108 109 110 111 112 113 114 115 116 117 `````` k = 10**(np.log10(Lx_max) - alpha*np.log10(t_sat)) Lx = powerlaw(t, k, alpha) elif t_start < t_sat: # then t_start < t_sat < t_curr if t > t_sat: #and t <= t_curr: alpha = (np.log10(Lx_curr/Lx_max))/(np.log10(t_curr/t_sat)) k = 10**(np.log10(Lx_max) - alpha*np.log10(t_sat)) Lx = powerlaw(t, k, alpha) elif t <= t_sat: Lx = Lx_max `````` Laura Ketzer committed May 08, 2020 118 `````` elif dt_drop > 0: # then t_sat != t_drop `````` Laura Ketzer committed May 07, 2020 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 `````` ################################################################################### # dt_drop>0 means we create a more fancy track with the saturation regime and then # a drop with slope change to the converging Lx at 1 Gyr ################################################################################### t_drop = t_sat + dt_drop # t_sat < t_drop if t <= t_sat: Lx = Lx_max elif (t > t_sat) and (t <= t_drop): # first of the two slopes alpha_drop1 = (np.log10((Lx_max/Lx_drop_factor)/Lx_max))/(np.log10(t_drop/t_sat)) k_drop1 = 10**(np.log10(Lx_max) - alpha_drop1*np.log10(t_sat)) Lx = powerlaw(t, k_drop1, alpha_drop1) elif t > t_drop: alpha_drop2 = (np.log10(Lx_curr/(Lx_max/Lx_drop_factor)))/(np.log10(t_curr/t_drop)) k_drop2 = 10**(np.log10((Lx_max/Lx_drop_factor)) - alpha_drop2*np.log10(t_drop)) Lx = powerlaw(t, k_drop2, alpha_drop2) return Lx ##################################################################################################################### # define flux def flux_at_planet_earth(L, a_p): """Function calculates the flux that the planet recieves at a given distance normalized with Earth's flux `````` Laura Ketzer committed May 08, 2020 142 `````` -> the semi-major axis is used for the distance (the eccentricity of the planetary orbits is ignored)""" `````` Laura Ketzer committed May 07, 2020 143 144 145 146 147 148 149 150 151 `````` A = (4.*np.pi*(a_p*const.au.cgs)**2) F = (L*u.erg/u.s)/A flux_earth = 1373*1e7*(u.erg/u.s)/(100*u.cm*100*u.cm) F_earth = (F/flux_earth).value return F_earth ##################################################################################################################### # define flux def flux_at_planet(L, a_p): `````` Laura Ketzer committed May 08, 2020 152 153 `````` """Function calculates the flux that the planet recieves at a given distance (in erg/s/cm^2) -> the semi-major axis is used for the distance (the eccentricity of the planetary orbits is ignored)""" `````` Laura Ketzer committed May 07, 2020 154 155 156 `````` A = (4.*np.pi*(a_p*const.au.cgs)**2) F = ((L*u.erg/u.s)/A).value return F `````` Laura Ketzer committed May 08, 2020 157 `````` `````` Laura Ketzer committed May 07, 2020 158 159 160 ``````##################################################################################################################### def L_xuv_all(Lx): """Function to estimate the EUV luminosity using the scaling relations given by Sanz-Forcada et al. (2011); `````` Laura Ketzer committed May 08, 2020 161 162 `````` the measured X-ray luminosity can be extrapolated to the total high-energy radiation as given below. To get the combined high-energy luminosity (L_xuv), we add L_x and L_euv. `````` Laura Ketzer committed May 07, 2020 163 164 165 166 167 168 169 170 `````` THEY USE: Lx(0.1-2.4 keV) ROSAT band (LX and LEUV are the X-ray (λλ 5−100 Å) and EUV (λλ 100−920 Å) luminosities, in erg s−1.) """ log_L_euv = 0.860*np.log10(Lx) + 4.80 Leuv = 10.**log_L_euv log_L_xuv = np.log10(Lx+Leuv) return 10.**log_L_xuv#*u.erg/u.s `````` 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 ``````##################################################################################################################### def undo_what_Lxuv_all_does(L_xuv): """ If you have LXUV given, this takes the Sanz-Forcada et al. (2011) scaling relation and reverses it to estimate Lx.""" def Calculate_Lx(Lx): return Lx + 10**(0.86*np.log10(Lx)+4.8) - L_xuv if (L_xuv > 1.*10**29): f_guess = L_xuv elif (L_xuv <= 1.*10**29) and (L_xuv > 1.*10**26): f_guess = 1.*10**25 elif (L_xuv <= 1.*10**26): f_guess = 1.*10**22 #print(f_guess) Lx = optimize.fsolve(Calculate_Lx, x0=f_guess)[0] return Lx `````` Laura Ketzer committed May 07, 2020 188 ``````##################################################################################################################### `````` Laura Ketzer committed May 08, 2020 189 190 191 192 ``````# def Lx_relation_Booth(t, R_star): # """R_star in terms of solar units""" # log_Lx = 54.65 - 2.8*np.log10(t) + 2*np.log10(R_star) # return 10.**log_Lx `````` Laura Ketzer committed May 07, 2020 193 194 195 `````` ##################################################################################################################### def calculate_Lx_sat(L_star_bol): `````` Laura Ketzer committed May 08, 2020 196 `````` """ Typical relation (from observations) to estimate the saturation X-ray luminosity given a bolometric luminosity.""" `````` 197 `` return 10**(-3)*(L_star_bol*const.L_sun.cgs).value``