Lx_evo_and_flux.py 10.1 KB
Newer Older
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
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;
13
    -> so t_curr(ent) is a little confusing! -> we have two cases!!
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)
22
    @param t_curr: (float) 
23
        Set to 1 Gyr (with Lx_curr = Lx at 1 Gyr), which is where, by design, the three 
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)
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));
    
34
    # IGNORE FOR NOW CASE 2 - not implemented yet
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:
77
        # this is the regime past 1 Gyr (the same for all Tu tracks)
78
        # I use the Tu 2015 slope which is approx the same for all three evolutionary tracks (and based on current solar value);
79
        # this is mainly for consistency since we approximate our tracks based on the Tu-tracks
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)
        
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))
106
                    print(alpha)
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
                    
118
        elif dt_drop > 0: # then t_sat != t_drop
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
142
    -> the semi-major axis is used for the distance (the eccentricity of the planetary orbits is ignored)"""
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):
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)"""
154
155
156
    A = (4.*np.pi*(a_p*const.au.cgs)**2)
    F = ((L*u.erg/u.s)/A).value
    return F
157

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);
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.
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

188
#####################################################################################################################
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
193
194
195

#####################################################################################################################
def calculate_Lx_sat(L_star_bol):
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