I made changes to Planet_class_LoFo14, in particular the output files, then I...

I made changes to Planet_class_LoFo14, in particular the output files, then I had to fix some issues with time steps in Lx_evo for the OwWu17 Lx track with const. slope, and I was in mass_evolution_function. These were some important fixes, now I can run the OwWu17 sample on merlot. YAY!

Email_Notifications.ipynb

@@ -64,7 +64,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.7.2"
   }
  },
  "nbformat": 4,

platypos_package/Lx_evo_and_flux.py

@@ -74,12 +74,23 @@ def Lx_evo(t, track_dict):
             print("Make sure t_curr > t_start, t >= t_start")
     
     elif t > t_curr:
-        # this is the regime past 1 Myr
+        # this is the regime past 1 Gyr (the same for all Tu tracks)
         # I use the Tu 2015 slope which is approx the same for all three evolutionary tracks (and based on current solar value);
         # this is mainly for consistency since we approximate our tracks based on the Tu-tracks
-        alpha = (np.log10(Lx_5Gyr/Lx_curr))/(np.log10(t_5Gyr/t_curr))
-        k = 10**(np.log10(Lx_curr) - alpha*np.log10(t_curr))
-        Lx = powerlaw(t, k, alpha)
+        
+        # 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)
+        
         return Lx
     
     elif t_curr > t_start:
@@ -92,6 +103,7 @@ 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)
                     

platypos_package/Planet_class_LoFo14.py

@@ -207,7 +207,7 @@ class planet_LoFo14():
             ###################################
             
             print("Start evolving.")
-            t, M, R, Lx = mass_planet_RK4_forward_LO14(epsilon=epsilon, K_on="yes", beta_on="yes", planet_object=self, 
+            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)
             df = pd.DataFrame({"Time": t, "Mass": M, "Radius": R, "Lx": Lx})
@@ -223,10 +223,11 @@ class planet_LoFo14():
                 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]
-                planet_params = "a,core_mass,mass,radius,metallicity,track\n" + \
+                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(M_final) + "," + str(R_final) + "," + self.metallicity + "," + \
-                                self.planet_id
+                                str(f_env_final) + "," + str(M_final) + "," + str(R_final) + \
+                                "," + self.metallicity + "," + self.planet_id
                 p.write(planet_params)
                 p.close()    
             

platypos_package/Planet_class_LoFo14_PAPER.py

@@ -196,7 +196,7 @@ class planet_LoFo14_PAPER():
             ###################################
         
             print("Start evolving.")
-            t, M, R, Lx = mass_planet_RK4_forward_LO14_PAPER(epsilon=epsilon, K_on="yes", beta_on="yes", planet_object=self, 
+            t, M, R, Lx = mass_planet_RK4_forward_LO14_PAPER(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)
             df = pd.DataFrame({"Time": t, "Mass": M, "Radius": R, "Lx": Lx})

platypos_package/Planet_class_Ot20.py

@@ -156,7 +156,7 @@ class planet_Ot20():
             ###################################
             
             print("Start evolving.")
-            t, M, R, Lx = mass_planet_RK4_forward_Ot14(epsilon=epsilon, K_on="yes", beta_on="yes", planet_object=self, 
+            t, M, R, Lx = mass_planet_RK4_forward_Ot14(epsilon=epsilon, K_on=K_on, beta_on=K_on, planet_object=self, 
                                                        initial_step_size=initial_step_size, t_final=t_final, 
                                                        track_dict=evo_track_dict)
             df = pd.DataFrame({"Time": t, "Mass": M, "Radius": R, "Lx": Lx})

platypos_package/Planet_class_Ot20_PAPER.py

@@ -150,7 +150,7 @@ class planet_Ot20_PAPER():
             ###################################
             
             print("Start evolving.")
-            t, M, R, Lx = mass_planet_RK4_forward_Ot14_PAPER(epsilon=epsilon, K_on="yes", beta_on="yes", planet_object=self, 
+            t, M, R, Lx = mass_planet_RK4_forward_Ot14_PAPER(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)
             df = pd.DataFrame({"Time": t, "Mass": M, "Radius": R, "Lx": Lx})