Merge pull request #7 from FormingWorlds/new_fxuv_w_mors

Add a Fxuv function using mors tracks

examples/comparison_work/escape_formulation.py

@@ -26,7 +26,6 @@
 Sun_Lxuv    = Sun.Tracks['Lx'] + Sun.Tracks['Leuv']                         # XUV luminosity            [erg s-1]
 Sun_Fxuv    = (Sun_Lxuv/(4 * np.pi * a_earth*au2cm **2)) * ergcm2stoWm2     # XUV flux                  [W m-2]
 
-
 ########################### Escape computations ####################################
 
 for mass in planet_masses : 

examples/comparison_work/test_fxuv_mors.py

@@ -0,0 +1,25 @@
+import matplotlib.pyplot as plt
+
+from zephyrus.constants import *
+from zephyrus.planets_parameters import *
+from zephyrus.XUV_flux import Fxuv_mors
+
+########################### Fxuv computation ####################################
+
+star_age_sun, star_fxuv_sun = Fxuv_mors(1.0, 1.0, 1.0)
+star_age_2, star_fxuv_2 = Fxuv_mors(1.0, 0.5, 1.0)
+
+########################### Plot ####################################
+
+plt.figure(figsize=(8, 6))
+
+plt.loglog(star_age_sun,star_fxuv_sun, color='gold', linestyle='-', label='1.0 AU, 1.0 Msun, 1.0 OmegaSun')
+plt.loglog(star_age_2,star_fxuv_2, color='red', linestyle='-', label='1.0 AU, 0.5 Msun, 1.0 OmegaSun')
+
+plt.xlabel('Age [Myr]', fontsize=15)
+plt.ylabel(r'XUV Flux received by the planet [W $m^{-2}$]', fontsize=15)
+plt.legend()
+plt.grid(alpha=0.5)
+plt.title('Test of the Fxuv_mors() function, using Star()', fontsize=15)
+
+plt.savefig('output/test_Fxuv_mors_function.pdf',dpi=180)

src/zephyrus/XUV_flux.py

@@ -7,6 +7,7 @@
 import pathlib
 from zephyrus.constants import *
 from zephyrus.planets_parameters import *
+import mors
 
 ########################### IsoFATE Fxuv functions ###########################
 #### The following functions are taken from the IsoFATE code written by Colin Cherubim 
@@ -181,3 +182,25 @@ def Fxuv_Johnstone_Sun(t,d):
     Fxuv_SI     = Fxuv * ergcm2stoWm2                                       # [W m-2]
 
     return np.interp(t,age,Fxuv_SI)
+
+
+def Fxuv_mors(M_star,Omega_star,d):
+    '''
+    Computes the incident XUV flux received by a planet at a distance d for a given star mass and rotation rate
+    Using stellar evolution files from MORS (Star() function)
+
+    Inputs:
+        - M_star : Mass of the star               [Msun]   
+        - Omega_star : Rotation rate of the star  [Omega sun]
+        - t : Time/age                            [s]
+        - d : Orbital distance of the planet      [AU]
+
+    Output: 
+        - Age of the star                         [Myr]
+        - Incident XUV flux                       [W m-2]
+    '''
+    star_data   = mors.Star(Mstar=M_star, Omega=Omega_star)                                                     # Extract luminosities using the mors.Star() function
+    Star_age    = star_data.Tracks['Age']                                                                       # Age of the Star            [Myr]
+    Star_Fxuv   = ((star_data.Tracks['Lx']+star_data.Tracks['Leuv'])/(4*np.pi*(d*au2cm)**2)) * ergcm2stoWm2     # XUV flux                   [W m-2]
+
+    return Star_age,Star_Fxuv