Planet_class_LoFo14.py 15.7 KB
Newer Older
1
import os
2
import numpy as np
3
import pandas as pd
4
from scipy.optimize import fsolve
5
import scipy.optimize as optimize
6
7
import astropy.units as u
from astropy import constants as const
8

9
10
11
12
13
14
15
from Lx_evo_and_flux import Lx_evo
from Lx_evo_and_flux import flux_at_planet_earth
from Lx_evo_and_flux import L_xuv_all
from Lx_evo_and_flux import flux_at_planet
# highest level function
from Mass_evolution_function import mass_planet_RK4_forward_LO14

16
17
18

class planet_LoFo14():
    """
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
    Need star and planet dictionary to initialize a LoFo planet object.

    Structure of star_dictionary: {'star_id': "dummySun", 'mass': mass_star,
                                   'radius': radius_star, 'age': age_star,
                                   'L_bol': L_bol, 'Lx_age': Lx_age}
    Structure of planet_dict: {"core_mass": m_c, "fenv": f,
                              "distance": a, "metallicity": metal}

    NOTE: There are actually three options for planet_dict:
         Case 1) An artificial planet with given core mass and envelope mass
                 fraction (M_core & f_env)
                -> in this case we need to calculate the mass and radius of
                the planet (M_pl & R_pl)

         Case 2) An observed planet with a known mass (we have M_pl & R_pl
                 & f_env) -> M_core

         Case 3) An observed planet with radius and mass measurement, plus
                 core mass is specified
                 -> need to calculate/estimate envelope mass fraction

    Additional Input: Lx_sat_info (only needed if you want to scale the
                      1 & 5 Gyr X-ray luminosities for non-solar-mass stars
42
    """
43

44
    def __init__(self, star_dictionary, planet_dict, Lx_sat_info=None):
45
46
47
48
49
50
51
        # initialize stellar parameters
        self.star_id = star_dictionary['star_id']
        self.mass_star = star_dictionary['mass']
        self.radius_star = star_dictionary['radius']
        self.age = star_dictionary['age']
        self.Lbol = star_dictionary['L_bol'] * const.L_sun.cgs.value
        self.Lx_age = star_dictionary['Lx_age']
52
        self.Lx_sat_info = Lx_sat_info
53
54
55

        # initialize planet parameters based on the input dictionary
        # the following 3 parameters are the same for all planets!
56
57
58
        self.distance = planet_dict["distance"]
        self.metallicity = planet_dict["metallicity"]  # solarZ or enhZ
        self.flux = flux_at_planet_earth(self.Lbol, self.distance)
59
60
61

        # flag which tracks if planet has been evolved and result file exists
        self.has_evolved = False
62
        self.planet_id = "None"
63
64

        # the remaining parameters depend on the input dictionary (Case 1, 2, 3)
65
66
        while True:
            try:
67
                # check if envelope mass fraction is specified, then Case 1
68
                planet_dict["fenv"]
69
70
                self.planet_info = "Case 1 - artificial planet"\
                                    + " (with M_core & f_env)"
71
72
                # Case 1: artificial planet with fenv & M_core given, need to
                #         calculate the total mass and radius
73
74
75
76
77
78
                self.fenv = planet_dict["fenv"]
                self.core_mass = planet_dict["core_mass"]
                self.Calculate_core_radius()
                self.Calculate_planet_mass()
                self.Calculate_planet_radius()
                break
79
80
            except KeyError:  # if no f_env is provided, then we are dealing
                              # with Case 2 or 3
81
82
                while True:
                    try:
83
                        # check if planet mass is provided, then Case 3
84
                        planet_dict["mass"]
85
86
87
                        # Case 3: An observed planet with radius and mass
                        #         measurement, plus core mass is specified;
                        #         need to calculate/estimate envelope mass frac.
88
89
90
                        self.planet_info = "Case 3 - obs. planet with radius"\
                        					+ " & mass measurement (and core "\
                        					+ "mass estimate)"
91
92
93
94
                        self.core_mass = planet_dict["core_mass"]
                        self.mass = planet_dict["mass"]
                        self.radius = planet_dict["radius"]
                        self.Calculate_core_radius()
95
96
97
98
99
100
                        self.Solve_for_fenv()  # get for envelope mass fraction

                        # Note to self: add sanity check to make sure the mass
                        # with the calculated fenv matches the input mass!
                        # Note to self: add sanity check to make sure the radius
                        # with the calculated fenv matches the input radius!
101
102
                        break
                    except KeyError:
103
                        # no mass and fenv given -> Case 2
104
105
                        self.planet_info = "Case 2 - obs. planet with"\
                        				    + " radius, but no mass measurement"
106
107
                        # Case 2 - observed planet with a without a mass, but
                        # core mass estimate, need to calculate fenv
108
109
110
                        self.core_mass = planet_dict["core_mass"]
                        self.radius = planet_dict["radius"]
                        self.Calculate_core_radius()
111
                        self.Solve_for_fenv()  # get envelope mass fraction
112
                        self.Calculate_planet_mass()
113
114
                        # Note to self: add sanity check to make sure the radius
                        # with the calculated fenv matches the input radius!
115
116
117
118
119
                        break
                break

    # Class Methods
    def Calculate_planet_mass(self):
120
121
        """ Planet mass determined by core and atmosphere mass
            (specified in terms of envelope mass fraction [in %]). """
122
        self.mass = self.core_mass/(1-(self.fenv/100))
123

124
    def Calculate_planet_core_mass(self):
125
126
        """ Core mass determined by planet mass and envelope mass
            (specified in terms of envelope mass fraction [%]). """
127
        self.core_mass = self.mass*(1-(self.fenv/100))
128

129
130
131
132
133
    def Calculate_core_radius(self):
        """M-R relation for rock/iron Earth-like core. (no envelope)"""
        self.core_radius = (self.core_mass**0.25)

    def Solve_for_fenv(self):
134
135
        """ For known core and planet radius, core mass, age and flux,
            solve for envelope mass fraction."""
136
137
138
139
        if self.radius == self.core_radius:
            self.fenv = 0.0
        else:
            def Calculate_fenv(fenv):
140
141
142
143
144
                age_exponent = {"solarZ": -0.11, "enhZ": -0.18}
                return -self.radius + self.core_radius + (2.06\
                       * (self.core_mass/(1-(fenv/100)))**(-0.21)\
                       * (fenv/5)**0.59 * (self.flux)**0.044\
                       * ((self.age/1e3)/5)**(age_exponent[self.metallicity]))
145
146
147
            f_guess = 0.1
            fenv = optimize.fsolve(Calculate_fenv, x0=f_guess)[0]
            self.fenv = fenv
148
149
150
            # if fenv1 != fenv:
            #    print("Sth went wrong in solving for\
            #          the envelope mass fraction! Check!")
151
152

    def Calculate_R_env(self):
153
154
155
156
        """ Check Planet_models_LoFo14.py for details on input and
            output parameters;
            R_env ~ t**0.18 for *enhanced opacities*;
            R_env ~ t**0.11 for *solar metallicities*
157
        """
158
        age_exponent = {"solarZ": -0.11, "enhZ": -0.18}
159
        R_env = 2.06 * self.mass**(-0.21) * (self.fenv/5)**0.59 * \
160
161
162
            self.flux**0.044 * \
            ((self.age/1e3)/5)**(age_exponent[self.metallicity])
        return R_env  # in units of R_earth
163
164

    def Calculate_planet_radius(self):
165
        """ Check Planet_models_LoFo14.py for details"""
166
        self.radius = self.core_radius + self.Calculate_R_env()
167
168
169
170
171
172

    def set_name(self, t_final, initial_step_size,
                 epsilon, K_on, beta_on, evo_track_dict):
        """ Function to set the right planet name based on the track specified.
        This can then be used to check if a particular planet on a particular
        track has already evolved (i.e. if outpufile exists). """
173
174
175
176
177
178
179
180
        self.planet_id = 'planet_a'\
        				 + str(np.round(self.distance, 3)).replace('.', 'dot')\
                         + '_Mcore'\
                         + str(np.round(self.core_mass, 3)).replace(".", "dot")\
                         + '_fenv'\
                         + str(np.round(self.fenv, 3)) + '_' + self.metallicity\
                         + '_Mstar'\
                         + str(np.round(self.mass_star, 3)).replace(".", "dot")\
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
                         + "_K_" + K_on + "_beta_" + beta_on\
                         + "_track_" + str(evo_track_dict["t_start"]) + "_"\
                         + str(evo_track_dict["t_sat"]) + "_" + str(t_final)\
                         + "_" + str(evo_track_dict["Lx_max"]) + "_"\
                         + str(evo_track_dict["dt_drop"]) + "_"\
                         + str(evo_track_dict["Lx_drop_factor"])

    def evolve_forward(self, t_final, initial_step_size, epsilon, K_on, beta_on,
                       evo_track_dict, path_for_saving, planet_folder_id):
        """ Call this function to make the planet evolve. This also produces the
        output files with the results. This includes a file on the host star
        parameters (host_star.txt), the initial planet parameters (name is the
        planet folder name, e.g. planet_001.txt), a file with the corresponding
        track parameters starting with "track_params_", afile with the
        mass and radius evolution, as well as a file with only the final params.

        See Mass_evolution_function.py for details on the integration."""
198
        # create a planet ID -> used for filname to save the result
199
200
201
202
203
204
205
        self.planet_id = planet_folder_id + "_track_"\
                         + str(evo_track_dict["t_start"]) + "_"\
                         + str(evo_track_dict["t_sat"]) + "_"\
                         + str(t_final) + "_" + str(evo_track_dict["Lx_max"])\
                         + "_" + str(evo_track_dict["dt_drop"]) + "_"\
                         + str(evo_track_dict["Lx_drop_factor"])

206
        if os.path.exists(path_for_saving+self.planet_id+".txt"):
207
            # planet already exists
208
209
            self.has_evolved = True
        else:
210
            #print("Planet: ", self.planet_id+".txt")
211
212
            # Next, create a file (e.g. planet_XXXX.txt) which contains the
            # initial planet params
213
214
215
            if os.path.exists(path_for_saving+planet_folder_id+".txt"):
                pass
            else:
216
                p = open(path_for_saving+planet_folder_id+".txt", "w")
217
218
219
220
221
                planet_params = "a,core_mass,fenv,mass,radius,metallicity,age\n"\
                                + str(self.distance) + "," + str(self.core_mass)\
                                + "," + str(self.fenv) + "," + str(self.mass)\
                                + "," + str(self.radius) + "," + self.metallicity\
                                + "," + str(self.age)
222
                p.write(planet_params)
223
224
225
226
227
228
                p.close()

            # create a file (track_params_planet_....txt) which contains
            # the track parameters
            if os.path.exists(path_for_saving + "track_params_"\
                              + self.planet_id + ".txt"):
229
230
                pass
            else:
231
                t = open(path_for_saving+"track_params_" +
232
233
234
                         self.planet_id+".txt", "w")
                track_params = "t_start,t_sat,t_curr,t_5Gyr,Lx_max,Lx_curr,"\
                                + "Lx_5Gyr,dt_drop,Lx_drop_factor\n"\
235
236
237
238
239
240
241
242
                                + str(evo_track_dict["t_start"]) + ","\
                                + str(evo_track_dict["t_sat"]) + "," \
                                + str(evo_track_dict["t_curr"]) + ","\
                                + str(evo_track_dict["t_5Gyr"]) + ","\
                                + str(evo_track_dict["Lx_max"]) + ","\
                                + str(evo_track_dict["Lx_curr"]) + ","\
                                + str(evo_track_dict["Lx_5Gyr"]) + ","\
                                + str(evo_track_dict["dt_drop"]) + ","\
243
244
245
                                + str(evo_track_dict["Lx_drop_factor"])
                t.write(track_params)
                t.close()
246
247

            # create a file which contains the host star parameters
248
249
250
            if os.path.exists(path_for_saving+"host_star_properties.txt"):
                pass
            else:
251
                s = open(path_for_saving+"host_star_properties.txt", "w")
252
253
254
255
256
                star_params = "star_id,mass_star,radius_star,age,Lbol,Lx_age\n"\
                              + self.star_id + "," + str(self.mass_star) + ","\
                              + str(self.radius_star) + "," + str(self.age)\
                              + "," + str(self.Lbol/const.L_sun.cgs.value)\
                              + "," + str(self.Lx_age)
257
258
                s.write(star_params)
                s.close()
259
260
261
262
263
264
265
266
267
268
269
270

            # call mass_planet_RK4_forward_LO14 to start the integration
            t, M, R, Lx = mass_planet_RK4_forward_LO14(
                                        epsilon=epsilon,
                                        K_on=K_on,
                                        beta_on=beta_on,
                                        planet_object=self,
                                        initial_step_size=initial_step_size,
                                        t_final=t_final,
                                        track_dict=evo_track_dict
                                        )
            # add results to dataframe and save
271
272
            df = pd.DataFrame({"Time": t, "Mass": M, "Radius": R, "Lx": Lx})
            df.to_csv(path_for_saving+self.planet_id+".txt", index=None)
273
274

            # create another file, which contains the final parameters only
275
            if os.path.exists(path_for_saving+\
276
                              self.planet_id+"_final.txt"):
277
278
                pass
            else:
279
                p = open(path_for_saving+self.planet_id+"_final.txt", "w")
280
281
                # get index of last valid entry in the radius array and access
                # its value
282
283
284
285
                index_of_last_entry = df["Radius"][df["Radius"].notna()].index[-1]
                R_final = df["Radius"].loc[index_of_last_entry]
                index_of_last_entry = df["Mass"][df["Mass"].notna()].index[-1]
                M_final = df["Mass"].loc[index_of_last_entry]
286
287
288
289
290
291
                f_env_final = ((M_final-self.core_mass)/M_final)*100  # %
                planet_params = "a,core_mass,fenv,mass,radius,metallicity,track\n"\
                                + str(self.distance) + "," + str(self.core_mass)\
                                + "," + str(f_env_final) + "," + str(M_final)\
                                + "," + str(R_final) + "," + self.metallicity\
                                + "," + self.planet_id
292
                p.write(planet_params)
293
294
295
296
                p.close()

            self.has_evolved = True  # set evolved-flag to True

297
    def read_results(self, file_path):
298
299
300
301
302
303
304
305
306
307
        # planet exists, has been evolved and results can be read in
        if self.has_evolved == True:
            df = pd.read_csv(file_path+self.planet_id+".txt",
                             float_precision='round_trip')
            # 'round_trip': to avoid pandas doing sth. weird to the last digit
            # Info: Pandas uses a dedicated dec 2 bin converter that compromises
            # accuracy in preference to speed. Passing float_precision =
            # 'round_trip' to read_csv fixes this.
            t, M, R, Lx = df["Time"].values, df["Mass"].values,\
                          df["Radius"].values, df["Lx"].values
308
309
310
            return t, M, R, Lx
        else:
            print("Planet has not been evolved & no result file exists.")
311
312
313
314

    def evolve_backwards(self, t_final, initial_step_size, epsilon, K_on,
                         beta_on, evo_track_dict, path_for_saving):
        print("Not implemented yet.")
315