Calculate and store orbital period

src/proteus/escape/wrapper.py

@@ -96,8 +96,6 @@ def RunZEPHYRUS(config, hf_row, stellar_track):
                     hf_row["R_xuv"], #XUV optically thick planetary radius [m]
                     Fxuv_star_SI)   # [kg s-1]
 
-    log.info('Zephyrus escape computation done :)')
-
     return mlr
 
 def SpeciesEscapeFromTotalEscape(hf_row:dict, dt:float):

src/proteus/orbit/__init__.py

@@ -0,0 +1,3 @@
+from __future__ import annotations
+
+__all__ = []

src/proteus/orbit/wrapper.py

@@ -0,0 +1,53 @@
+# Generic orbital dynamics stuff
+from __future__ import annotations
+
+import logging
+
+import numpy as np
+
+from proteus.utils.constants import AU, const_G
+
+log = logging.getLogger("fwl."+__name__)
+
+def update_separation(hf_row:dict, sma:float, ecc:float):
+    '''
+    Calculate time-averaged orbital separation on an elliptical path.
+
+    Converts from AU to metres.
+    https://physics.stackexchange.com/a/715749
+
+    Parameters
+    -------------
+        hf_row: dict
+            Current helpfile row
+        sma: float
+            Semimajor axis [AU]
+        ecc: float
+            Eccentricity
+    '''
+
+    hf_row["separation"] = sma * AU *  (1 + 0.5*ecc*ecc)
+
+def update_period(hf_row:dict, sma:float):
+    '''
+    Calculate orbital period on an elliptical path.
+
+    Assuming that M_volatiles << M_star + M_mantle + M_core.
+    https://en.wikipedia.org/wiki/Elliptic_orbit#Orbital_period
+
+    Parameters
+    -------------
+        hf_row: dict
+            Current helpfile row
+        sma: float
+            Semimajor axis [AU]
+    '''
+
+    # Standard gravitational parameter neglecting volatile mass.
+    mu = const_G * (hf_row["M_star"] + hf_row["M_mantle"] + hf_row["M_core"])
+
+    # Semimajor axis in SI units
+    a = sma * AU
+
+    # Orbital period [seconds]
+    hf_row["period"] = 2 * np.pi * (a*a*a/mu)**0.5

src/proteus/proteus.py

@@ -12,6 +12,7 @@
 from proteus.config import read_config_object
 from proteus.escape.wrapper import RunEscape
 from proteus.interior import run_interior
+from proteus.orbit.wrapper import update_period, update_separation
 from proteus.outgas.wrapper import calc_target_elemental_inventories, run_outgassing
 from proteus.star.wrapper import (
     get_new_spectrum,
@@ -21,7 +22,6 @@
     write_spectrum,
 )
 from proteus.utils.constants import (
-    AU,
     M_earth,
     R_earth,
     const_G,
@@ -263,16 +263,6 @@ def start(self, *, resume: bool = False):
                 )
             )
 
-            ############### ORBIT AND TIDES
-
-            # Calculate time-averaged orbital separation (and convert from AU to metres)
-            # https://physics.stackexchange.com/a/715749
-            self.hf_row["separation"] = self.config.orbit.semimajoraxis * AU * \
-                                        (1 + 0.5 * self.config.orbit.eccentricity**2.0)
-
-
-            ############### / ORBIT AND TIDES
-
             ############### INTERIOR
 
             # Run interior model
@@ -285,6 +275,18 @@ def start(self, *, resume: bool = False):
 
             ############### / INTERIOR
 
+            ############### ORBIT AND TIDES
+
+            # Update orbital separation
+            update_separation(self.hf_row,
+                                self.config.orbit.semimajoraxis,
+                                self.config.orbit.eccentricity)
+
+            # Update orbital period
+            update_period(self.hf_row, self.config.orbit.semimajoraxis)
+
+            ############### / ORBIT AND TIDES
+
             ############### STELLAR FLUX MANAGEMENT
             PrintHalfSeparator()
             log.info("Stellar flux management...")

src/proteus/star/wrapper.py

@@ -8,7 +8,7 @@
 
 import numpy as np
 
-from proteus.utils.constants import AU, R_sun, const_sigma
+from proteus.utils.constants import AU, M_sun, R_sun, const_sigma
 
 log = logging.getLogger("fwl."+__name__)
 
@@ -184,9 +184,10 @@ def write_spectrum(wl_arr, fl_arr, hf_row:dict, output_dir:str):
 
 def update_stellar_quantities(hf_row:dict, config:Config, stellar_track=None):
 
-    # Update value for star's radius
-    log.info("Update stellar radius")
+    # Update value for star's radius and mass
+    log.info("Update stellar radius and mass")
     update_stellar_radius(hf_row, config, stellar_track)
+    update_stellar_mass(hf_row, config)
 
     # Update value for instellation flux
     log.info("Update instellation")
@@ -200,6 +201,12 @@ def update_stellar_quantities(hf_row:dict, config:Config, stellar_track=None):
     # Assuming a grey stratosphere in radiative eqm (https://doi.org/10.5194/esd-7-697-2016)
     hf_row["T_skin"] = hf_row["T_eqm"] * (0.5**0.25)
 
+def update_stellar_mass(hf_row:dict, config:Config):
+    '''
+    Update stellar mass in hf_row, stored in SI units.
+    '''
+    hf_row["M_star"] = config.star.mass * M_sun
+
 def update_stellar_radius(hf_row:dict, config:Config, stellar_track=None):
     '''
     Update stellar radius in hf_row, stored in SI units.

src/proteus/utils/coupler.py

@@ -243,8 +243,12 @@ def GetHelpfileKeys():
 
     # Basic keys
     keys = [
-            # Model tracking and basic parameters
-            "Time", "separation", # [yr], [m]
+            # Model tracking
+            "Time", # [yr]
+
+            # Orbital dynamics
+            "separation", # [m]
+            "period", # [s]
 
             # Input parameters (converted to SI)
             "R_int", "M_int", # [m], [kg]
@@ -261,7 +265,7 @@ def GetHelpfileKeys():
             "M_mantle_solid", "M_mantle_liquid", # all [kg]
 
             # Stellar
-            "R_star", "age_star", # [m], [yr]
+            "M_star", "R_star", "age_star", # [kg], [m], [yr]
 
             # Observational (from infinity)
             "z_obs", # observed height relative to R_int [m]
@@ -270,7 +274,7 @@ def GetHelpfileKeys():
             "bond_albedo", # bond albedo [1]
 
             # Escape
-            "esc_rate_total", "p_xuv", "z_xuv", "R_xuv", # [kg s-1], [m], [m]
+            "esc_rate_total", "p_xuv", "z_xuv", "R_xuv", # [kg s-1], [bar], [m], [m]
 
             # Atmospheric composition
             "M_atm", "P_surf", "atm_kg_per_mol", # [kg], [bar], [kg mol-1]