worked on evolve_planet, to evolve OwWu17 sample (each planet has a different...

worked on evolve_planet, to evolve OwWu17 sample (each planet has a different host star) along various tracks

platypos_package/Lx_evo_and_flux.py

@@ -87,7 +87,6 @@ def Lx_evo(t, track_dict):
 
         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)
         
@@ -103,7 +102,6 @@ def Lx_evo(t, track_dict):
                 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))
-                    print(alpha)
                     k = 10**(np.log10(Lx_max) - alpha*np.log10(t_sat))
                     Lx = powerlaw(t, k, alpha)
                     
@@ -195,3 +193,18 @@ def undo_what_Lxuv_all_does(L_xuv):
 def calculate_Lx_sat(L_star_bol):
     """ Typical relation (from observations) to estimate the saturation X-ray luminosity given a bolometric luminosity."""
     return 10**(-3)*(L_star_bol*const.L_sun.cgs).value
+
+
+def L_high_energy(t_i, mass_star):
+    """ L_HE (UV to X-ray) in Owen & Wu (2017).
+    @param t_i - in Myr
+    @param mass_star - in solar masses
+    """
+    t_sat = 100. # Myr
+    L_sat = 10**(-3.5)*(mass_star) * (const.L_sun.cgs).value
+    a_0 = 0.5
+    if t_i < t_sat:
+        L_HE = L_sat
+    elif t_i >= t_sat:
+        L_HE = L_sat*(t_i/t_sat)**(-1-a_0)
+    return L_HE # in erg/s
\ No newline at end of file

population_evolution/evolve_planet.py

@@ -8,7 +8,7 @@ sys.path.append('../plaml_package/')
 # Planet Class
 from Planet_class_LoFo14 import planet_LoFo14
 from Planet_class_Ot20 import planet_Ot20
-from Lx_evo_and_flux import undo_what_Lxuv_all_does
+from Lx_evo_and_flux import undo_what_Lxuv_all_does, L_high_energy
 
 
 def evolve_one_planet(pl_folder_pair, t_final, init_step, eps, K_on, beta_on, evo_track_dict_list_or_1, path_for_saving):
@@ -55,19 +55,56 @@ def evolve_one_planet(pl_folder_pair, t_final, init_step, eps, K_on, beta_on, ev
                              path_for_saving=path_for_saving, planet_folder_id=folder)
     
     else: # evolve on planet along more than one track (e.g. to see effects of different tracks on planetary evolution)
-        for track in evo_track_dict_list:
+        # either the dictionary in the list already contains ALL nine parameters (then continue)
+        if len(evo_track_dict_list_or_1[0]) == 9:
             pl = pl_folder_pair[1]
             folder = pl_folder_pair[0]
-            pl.set_name(t_final, init_step, eps, K_on, beta_on, evo_track_dict=track) # set planet name based on specified track
-            # generate a useful planet name for saving the results
-            pl_file_name = folder + "_track" + pl.planet_id.split("_track")[1] + ".txt"
-            if os.path.isdir(path_for_saving+pl_file_name):
-                # skip & don't evolve planet
-                # print("Folder "+pl_file_name+" exists.")
-                pass
-            else:
-                pl.evolve_forward(t_final, init_step, eps, K_on, beta_on, evo_track_dict=track, 
-                                  path_for_saving=path_for_saving, planet_folder_id=folder)
+            for track in evo_track_dict_list:
+                pl.set_name(t_final, init_step, eps, K_on, beta_on, evo_track_dict=track) # set planet name based on specified track
+                # generate a useful planet name for saving the results
+                pl_file_name = folder + "_track" + pl.planet_id.split("_track")[1] + ".txt"
+                if os.path.isdir(path_for_saving+pl_file_name):
+                    # skip & don't evolve planet
+                    # print("Folder "+pl_file_name+" exists.")
+                    pass
+                else:
+                    pl.evolve_forward(t_final, init_step, eps, K_on, beta_on, evo_track_dict=track, 
+                                      path_for_saving=path_for_saving, planet_folder_id=folder)
+
+        else:
+            # calculate missing parameters based on individual planet (host star is different!!)
+            Lx_solar = L_high_energy(1., 1.0) # value with which to shift Lx_1Gyr & Lx_5Gyr using OwWu17 Lx_sat
+            # if I use Tu Lx values, then I should change the solar saturation value to what they use!
+            #Lx_solar = 10**(-3.13)*const.L_sun.cgs.value
+
+            pl = pl_folder_pair[1]
+            folder = pl_folder_pair[0]
+            # get factor with which to scale tracks up or down
+            scaling_factor = pl.Lx_age/Lx_solar
+            #Lx_age_approx = 10**(-3)*(10**mass_luminosity_relation(smass))*const.L_sun.cgs.value
+
+            # create dictionaries with all the values necessary to create the evolutionary paths (tracks converge at 1 Gyr!)
+            Lx_1Gyr = 2.10*10**28 * scaling_factor  # Lx value at 1 Gyr from Tu et al. (2015) model tracks (scaled)
+            Lx_5Gyr = 1.65*10**27 * scaling_factor  # Lx value at 5 Gyr from Tu et al. (2015) model tracks (scaled)
+
+            for t in evo_track_dict_list_or_1:
+                same_track_params = {"t_start": pl.age, "t_curr": 1e3, "t_5Gyr": 5e3, "Lx_curr": Lx_1Gyr, "Lx_5Gyr": Lx_5Gyr, "Lx_max": pl.Lx_age}
+                # combine the same_track params with different track params for a complete track dictionary
+                track = same_track_params.copy()
+                track.update(t)
+                #print(track)
+
+                pl.set_name(t_final, init_step, eps, K_on, beta_on, evo_track_dict=track) # set planet name based on specified track
+                # generate a useful planet name for saving the results
+                pl_file_name = folder + "_track" + pl.planet_id.split("_track")[1] + ".txt"
+
+                if os.path.isdir(path_for_saving+pl_file_name):
+                        # skip & don't evolve planet
+                        # print("Folder "+pl_file_name+" exists.")
+                        pass
+                else:
+                    pl.evolve_forward(t_final, init_step, eps, K_on, beta_on, evo_track_dict=track, 
+                                      path_for_saving=path_for_saving, planet_folder_id=folder)
 
 
 def evolve_ensamble(planets_chunks, t_final, init_step, eps, K_on, beta_on, evo_track_dict_list_or_1, path_save):

population_evolution/get_CKS_sample.py

+
+
+def read_in_CKS_sample(path_to_directory):
+    """Get dataset from CKS (California-Kepler-Survey): https://github.com/California-Planet-Search/cks-website/blob/master/readme.md
+    Column explanations: https://www.astro.caltech.edu/~howard/cks/column-definitions.txt)"""
+    
+    df_cks_orig = pd.read_csv(path_to_directory)
+    print(len(df_cks_orig))
+    df_cks = df_cks_orig.copy()
+    # sample selection based on Fulton 2017
+    mask_confirmed = df_cks["koi_disposition"] == "CONFIRMED" # only select confirmed planets
+    df_cks = df_cks[mask_confirmed]
+    print(len(df_cks))
+    # restrict sample to only the magnitudelimited portion of the larger CKS sample (Kp <14.2):
+    mask_magnitude = df_cks["kic_kmag"] < 14.2
+    df_cks = df_cks[mask_magnitude]
+    print(len(df_cks))
+    # planet-to-star radius ratio (R_pl/R_star) becomes uncertain at high impact parameters (b) due to degeneracies with limbdarkening. 
+    # excluded KOIs with b > 0.7 to minimize the impact of grazing geometries.
+    mask_impactparam = df_cks["koi_impact"] <= 0.7
+    df_cks = df_cks[mask_impactparam]
+    print(len(df_cks))
+    # remove planets with orbital periods longer than 100 days in order to avoid domains of low completeness 
+    # (especially for planets smaller than about 4 R_earth) and low transit probability.
+    mask_period = df_cks["koi_period"] <= 100.
+    df_cks = df_cks[mask_period]
+    print(len(df_cks))
+    # also excised planets orbiting evolved stars since they have somewhat lower detectability and less certain radii. 
+    # implemented using an ad hoc temperature-dependent stellar radius filter:
+    mask_evolved = df_cks["koi_srad"] <= 10**(0.00025*(df_cks["koi_steff"]-5500)+0.20)
+    df_cks = df_cks[mask_evolved]
+    print(len(df_cks))
+    # also restrict sample to planets orbiting stars within the temperature range where we can extract precise stellar parameters 
+    # from our high-resolution optical spectra (6500–4700 K).
+    mask_temperature = (df_cks["koi_steff"] >= 4700) & (df_cks["koi_steff"] <= 6500)
+    df_cks = df_cks[mask_temperature]
+    print(len(df_cks))
+
+    # drop columns which have missing stellar mass, planetary radius, semi-major axis or period
+    df_cks = df_cks.dropna(axis=0, how="any", subset=["koi_sma", "koi_period", "koi_smass", "koi_prad"])
+    df_cks.reset_index(inplace=True)
+
+    return df_cks # final filtered sample
+
+    # # inflate uncertainties on the histogram bin heights by the scaling factors to account for completeness
+    # # corrections -> stellar properties of planet-hosting stars come from a different source and have higher 
+    # # precision than the stellar properties for the full set of Kepler stars
+    # df_R_bin_scaling = pd.read_csv("../supplementary_files/Fulton_Radius_correction.csv", sep="\s+", header=None, names=["bin", "factor"])
+    # df_R_bin_scaling["bin_left"] = np.zeros(len(df_R_bin_scaling))
+    # df_R_bin_scaling["bin_right"] = np.zeros(len(df_R_bin_scaling))
+
+    # for index in df_R_bin_scaling.index:
+    #     # Select rpw by index position using iloc[]
+    #     bin_range = df_R_bin_scaling["bin"].iloc[index].split("–")
+    #     df_R_bin_scaling.loc[index, "bin_left"] = float(bin_range[0])
+    #     df_R_bin_scaling.loc[index, "bin_right"] = float(bin_range[1])
+    #     factor = df_R_bin_scaling["factor"].iloc[index]
+
+    # df_R_bin_scaling.head()    
\ No newline at end of file

population_evolution/mass_luminosity_relation.py

+import numpy as np
+import pandas as pd
+from scipy import interpolate
+
+def M_L_relation():
+    """ main-sequencen mass-luminosity relation""" 
+    # this is from Thomas, I am not sure where exactly it comes from
+    df = pd.read_csv("../supplementary_files/ZAMS_online.dat", header="infer", sep="\s+")
+    logL_from_M = interpolate.interp1d(df["Mass"], np.log10(df["Lbol"]), kind="quadratic")
+    return logL_from_M
\ No newline at end of file

population_evolution/result_files_OwWu17_smass_cks_Mcore5_5Gyr/planet_0001/host_star_properties.txt

+star_id,mass_star,radius_star,age,Lbol,Lx_age
+star_age1.0_mass0.72,0.7184941392412282,None,1.0,1.0,1.392369397981642e+29
\ No newline at end of file

population_evolution/result_files_OwWu17_smass_cks_Mcore5_5Gyr/planet_0001/planet_0001.txt

+a,core_mass,fenv,mass,radius,metallicity,age
+0.1385481957745512,4.266318510202839,21.077576557608428,5.405711487454493,11.691119046005358,solarZ,1.0
\ No newline at end of file

population_evolution/result_files_OwWu17_smass_cks_Mcore5_5Gyr/planet_0001/planet_0001_track_1.0_100.0_5000.0_1.392369397981642e+29_0.0_0.0_final.txt

+a,core_mass,fenv,mass,radius,metallicity,track
+0.1385481957745512,4.266318510202839,15.578695623472305,5.053604112979136,4.846603075263876,solarZ,planet_0001_track_1.0_100.0_5000.0_1.392369397981642e+29_0.0_0.0
\ No newline at end of file

population_evolution/result_files_OwWu17_smass_cks_Mcore5_5Gyr/planet_0001/track_params_planet_0001_track_1.0_100.0_5000.0_1.392369397981642e+29_0.0_0.0.txt

+t_start,t_sat,t_curr,t_5Gyr,Lx_max,Lx_curr,Lx_5Gyr,dt_drop,Lx_drop_factor
+1.0,100.0,5000.0,5000.0,1.392369397981642e+29,3.9382153729178e+26,3.9382153729178e+26,0.0,0.0
\ No newline at end of file

population_evolution/result_files_OwWu17_smass_cks_Mcore5_5Gyr/planet_0002/host_star_properties.txt

+star_id,mass_star,radius_star,age,Lbol,Lx_age
+star_age1.0_mass1.01,1.0097333169215508,None,1.0,1.0,2.047074262504397e+29
\ No newline at end of file

population_evolution/result_files_OwWu17_smass_cks_Mcore5_5Gyr/planet_0002/planet_0002.txt

+a,core_mass,fenv,mass,radius,metallicity,age
+0.2552748377979306,8.327241668651112,1.274857470980825,8.434773002432902,3.389666871505346,solarZ,1.0
\ No newline at end of file