Added cumulative XUV Flux and reproduction of the "Cosmic Shoreline"

Makefile

@@ -30,7 +30,7 @@ legacy:
 	@echo "=========================================================================================================="
 
 debug:
-	-gcc -g -D DEBUG -Wunused-but-set-variable -Wunused-variable -Wfloat-equal -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\"
+	-gcc -g -D DEBUG -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\"
 
 debug_no_AE:
 	-gcc -g -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\"
@@ -52,7 +52,7 @@ cpp:
 	g++ -o bin/vplanet src/*.c -lm -O3 -fopenmp -fpermissive -w -DGITVERSION=\"$(GITVERSION)\"
 
 warnings:
-	clang -Weverything src/*.c -lm -O3 -DGITVERSION=\"$(GITVERSION)\"
+		-gcc -g -D DEBUG -Wunused-but-set-variable -Wunused-variable -Wfloat-equal -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\"
 
 parallel:
 	gcc -o bin/vplanet src/*.c -lm -O3 -fopenmp -DGITVERSION=\"$(GITVERSION)\"

examples/CosmicShoreline/README.rst

@@ -0,0 +1,38 @@
+The "Cosmic Shoreline"
+==========
+
+Overview
+--------
+
+A reproduction of Zahnle & Catling (2018)'s cosmic shoreline,
+a proposed boundary between planet with and without atmosphere.
+
+===================   ============
+**Date**              09/12/2024
+**Author**            Rory Barnes
+**Modules**           Atmesc, STELLAR
+**Approx. runtime**   10 seonds
+===================   ============
+
+
+To run this example
+-------------------
+
+.. code-block:: bash
+
+   python makeplot.py <pdf | png>
+
+
+Expected output
+---------------
+
+.. figure:: CosmicShoreline.png
+   :width: 100%
+   :align: center
+
+The blue line is the cosmic shoreline, which is taken from Fig. 2
+of Zahnle & Catling (2018). The black points are the positions
+of the Solar System planets in this parameter space. The Sun's XUV luminosity
+is assumed to follow the Ribas et al. (2005) model. This figure, while
+not as complete as Zahnle & Catling (2018), is a near exact reproduction
+of their results.

examples/CosmicShoreline/earth.in

@@ -0,0 +1,6 @@
+sName        Earth
+saModules    atmesc
+dMass        3.0018090452e-06
+dRadius      -1.
+dSemi        1.00000321
+saOutputOrder    Time -FXUV -CumulativeXUVFlux

examples/CosmicShoreline/george.in

@@ -0,0 +1,6 @@
+sName        George
+saModules    atmesc
+dMass        -14.54
+dRadius      -4.007
+dSemi        19.19203990
+saOutputOrder    Time -fXUV -CumulativeXUVFlux

examples/CosmicShoreline/jupiter.in

@@ -0,0 +1,6 @@
+sName        Jupiter
+saModules    atmesc
+dMass        -317.83
+dRadius      -11.2
+dSemi        5.20880408
+saOutputOrder    Time -FXUV -CumulativeXUVFlux

examples/CosmicShoreline/makeplot.py

@@ -0,0 +1,73 @@
+import os
+import pathlib
+import subprocess
+import sys
+
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+import vplot as vpl
+from matplotlib.ticker import MaxNLocator
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+import vplanet
+
+path = pathlib.Path(__file__).parents[0].absolute()
+sys.path.insert(1, str(path.parents[0]))
+
+output = vplanet.run(units = False)
+
+# Plot!
+fig = plt.figure(figsize=(8.5, 6))
+
+fxuv_earth = output.log.final.Earth.CumulativeXUVFlux
+
+fxuv = []
+fxuv.append(output.log.final.Mercury.CumulativeXUVFlux/fxuv_earth)
+fxuv.append(output.log.final.Venus.CumulativeXUVFlux/fxuv_earth)
+fxuv.append(output.log.final.Earth.CumulativeXUVFlux/fxuv_earth)
+fxuv.append(output.log.final.Mars.CumulativeXUVFlux/fxuv_earth)
+fxuv.append(output.log.final.Jupiter.CumulativeXUVFlux/fxuv_earth)
+fxuv.append(output.log.final.Saturn.CumulativeXUVFlux/fxuv_earth)
+fxuv.append(output.log.final.George.CumulativeXUVFlux/fxuv_earth)
+fxuv.append(output.log.final.Neptune.CumulativeXUVFlux/fxuv_earth)
+
+escvel = []
+escvel.append(output.log.final.Mercury.EscapeVelocity/1e3)
+escvel.append(output.log.final.Venus.EscapeVelocity/1e3)
+escvel.append(output.log.final.Earth.EscapeVelocity/1e3)
+escvel.append(output.log.final.Mars.EscapeVelocity/1e3)
+escvel.append(output.log.final.Jupiter.EscapeVelocity/1e3)
+escvel.append(output.log.final.Saturn.EscapeVelocity/1e3)
+escvel.append(output.log.final.George.EscapeVelocity/1e3)
+escvel.append(output.log.final.Neptune.EscapeVelocity/1e3)
+
+shorelinex = []
+shorelinex.append(0.2)
+shorelinex.append(60)
+
+shoreliney = []
+shoreliney.append(1e-6)
+shoreliney.append(1e4)
+
+plt.xlabel('Escape Velocity [km/s]')
+plt.ylabel('Normalized Cumulative XUV Flux')
+plt.plot(shorelinex,shoreliney,color=vpl.colors.pale_blue)
+plt.plot(escvel,fxuv,'o',color='k')
+plt.xscale('log')
+plt.yscale('log')
+plt.ylim(1e-6,1e4)
+plt.xlim(0.1,200)
+
+plt.annotate('Mercury',(2.7,9))
+plt.annotate('Venus',(10.2,2.6))
+plt.annotate('Earth',(11.5,0.5))
+plt.annotate('Mars',(5,0.2))
+plt.annotate('Jupiter',(60,0.05))
+plt.annotate('Saturn',(37,0.014))
+plt.annotate('Uranus',(22,0.0034))
+plt.annotate('Neptune',(24,0.0006))
+
+
+# Save figure
+fig.savefig(path / f"CosmicShoreline.png", bbox_inches="tight", dpi=200)

examples/CosmicShoreline/mars.in

@@ -0,0 +1,6 @@
+sName        Mars
+saModules    atmesc
+dMass        -0.107
+dRadius      -0.53202
+dSemi        1.52366290
+saOutputOrder    Time -FXUV -CumulativeXUVFlux

examples/CosmicShoreline/mercury.in

@@ -0,0 +1,6 @@
+sName        Mercury
+saModules    atmesc
+dMass        -0.0553
+dRadius      -0.383
+dSemi        0.38709897
+saOutputOrder    Time -FXUV -CumulativeXUVFlux

examples/CosmicShoreline/neptune.in

@@ -0,0 +1,6 @@
+sName        Neptune
+saModules    atmesc
+dMass        -17.15
+dRadius      -3.883
+dSemi        30.07050641
+saOutputOrder    Time -FXUV -CumulativeXUVFlux

examples/CosmicShoreline/saturn.in

@@ -0,0 +1,6 @@
+sName        Saturn
+saModules    atmesc
+dMass        -95.16
+dRadius      -9.449
+dSemi       9.53999265
+saOutputOrder    Time -FXUV -CumulativeXUVFlux

examples/CosmicShoreline/sun.in

@@ -0,0 +1,18 @@
+# The host star
+sName	        Sun		                    # Body's name
+saModules	    stellar                     # Modules to apply, exact spelling required
+
+# Output
+saOutputOrder Age -Luminosity -LXUVStellar -Radius Temperature
+
+# Physical Parameters
+dMass         1.00                        # Mass, solar masses
+dAge          5e7
+dObliquity    0
+dRotPeriod    -30                         # Rotation period, negative -> days
+dRadGyra      0.5                         # Radius of gyration (moment of inertia constant)
+
+# STELLAR Parameters
+sStellarModel baraffe                     # Stellar evolution model: `baraffe` or `none`
+dSatXUVFrac   1.e-3                       # Saturation level of the XUV luminosity
+dSatXUVTime   -0.1

examples/CosmicShoreline/venus.in

@@ -0,0 +1,7 @@
+# Planet a parameters
+sName         Venus                      # Body's name
+saModules     atmesc
+dMass         -0.815
+dRadius       -0.9499
+dSemi         -0.723
+saOutputOrder Time -FXUV -CumulativeXUVFlux

examples/CosmicShoreline/vpl.in

@@ -0,0 +1,25 @@
+sSystemName   solarsystem                   # System Name
+iVerbose      5                             # Verbosity level
+bOverwrite    1                             # Allow file overwrites?
+
+# List of "body files" that contain body-specific parameters
+saBodyFiles   sun.in mercury.in venus.in earth.in mars.in $
+              jupiter.in saturn.in george.in neptune.in
+
+# Input/Output Units
+sUnitMass      solar                        # Options: gram, kg, Earth, Neptune, Jupiter, solar
+sUnitLength    aU                           # Options: cm, m, km, Earth, Jupiter, solar, AU
+sUnitTime      YEARS                        # Options: sec, day, year, Myr, Gyr
+sUnitAngle     d                            # Options: deg, rad
+
+# Input/Output
+bDoLog         1                            # Write a log file?
+iDigits        6                            # Maximum number of digits to right of decimal
+dMinValue      1e-10                        # Minimum value of eccentricity/obliquity
+
+# Evolution Parameters
+bDoForward    1                             # Perform a forward evolution?
+bVarDt        1                             # Use variable timestepping?
+dEta          0.01                           # Coefficient for variable timestepping
+dStopTime     4.5e9                         # Stop time for evolution
+dOutputTime   1e8                           # Output timesteps (assuming in body files)

src/atmesc.c

@@ -2001,7 +2001,6 @@ void fnPropsAuxAtmEsc(BODY *body, EVOLVE *evolve, IO *io, UPDATE *update,
   }
 }
 
-
 /**
 Assigns functions returning the time-derivatives of each variable
 to the magical matrix of function pointers.

src/atmesc.h

@@ -171,6 +171,7 @@ void FinalizeUpdateMassAtmEsc(BODY *, UPDATE *, int *, int, int, int);
 #define OUT_HDIFFFLUX 1238 /**< The diffusion limited flux (the true flux in the diffusion regime) */
 #define OUT_HREFODRAGMOD 1239 /**< Multiply by H ref flux to get H flux with drag of oxgyen */
 #define OUT_KTIDE 1240 /**< Gravitational enhancement of mass loss */
+#define OUT_CUMULATIVEFXUV 1250 
 
 void InitializeOutputAtmEsc(OUTPUT *, fnWriteOutput[]);
 void InitializeOutputFunctionAtmEsc(OUTPUT *, int, int);
@@ -212,6 +213,7 @@ void LogBodyAtmEsc(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UPDATE *,
 void fnForceBehaviorAtmEsc(BODY *, MODULE *, EVOLVE *, IO *, SYSTEM *, UPDATE *,
                            fnUpdateVariable ***, int, int);
 void fnPropsAuxAtmEsc(BODY *, EVOLVE *, IO *, UPDATE *, int);
+
 double fdDSurfaceWaterMassDt(BODY *, SYSTEM *, int *);
 double fdDEnvelopeMassDt(BODY *, SYSTEM *, int *);
 double fdDEnvelopeMassDtBondiLimited(BODY *, SYSTEM *, int *);

src/body.c

@@ -1528,4 +1528,11 @@ void fdHabitableZoneKopparapu2013(BODY *body, int iNumBodies,
 double fdEffectiveTemperature(BODY *body,int iBody) {
   double dTeff = pow((body[iBody].dLuminosity/(4*PI*SIGMA*body[iBody].dRadius*body[iBody].dRadius)),0.25);
   return dTeff;
+}
+
+double fdEscapeVelocity(BODY *body,int iBody) {
+  double dEscVel;
+
+  dEscVel = sqrt(2*BIGG*body[iBody].dMass/body[iBody].dRadius);
+  return dEscVel;
 }

src/body.h

@@ -110,6 +110,7 @@ double CalcDynEllipEq(BODY *, int);
 void fdHabitableZoneKopparapu2013(BODY *, int, double *);
 
 double fdEffectiveTemperature(BODY*,int);
+double fdEscapeVelocity(BODY *,int);
 
 // RB: Move
 

src/evolve.c

@@ -13,6 +13,11 @@
 
 #include "vplanet.h"
 
+double fdTrapezoidalArea(double dY1,double dY2,double dDeltaX) {
+  double dArea = 0.5*dDeltaX*(dY1 + dY2);
+  return dArea;
+}
+
 void PropsAuxGeneral(BODY *body, CONTROL *control) {
   /* Recompute the mean motion, necessary for most modules */
   int iBody; // Dummy counting variable
@@ -50,6 +55,17 @@ void PropertiesAuxiliary(BODY *body, CONTROL *control, SYSTEM *system,
   }
 }
 
+void fvCumulative(BODY *body,CONTROL *control,SYSTEM *system,double dDt) {
+  int iBody;
+
+  if (body[0].bStellar) {
+    for (iBody=1;iBody<control->Evolve.iNumBodies;iBody++) {
+      fvCumulativeXUVFlux(body,&control->Evolve,system,dDt,iBody);
+    }
+  }
+
+}
+
 void CalculateDerivatives(BODY *body, SYSTEM *system, UPDATE *update,
                           fnUpdateVariable ***fnUpdate, int iNumBodies) {
   int iBody, iVar, iEqn;
@@ -634,6 +650,8 @@ void Evolve(BODY *body, CONTROL *control, FILES *files, MODULE *module,
     dDt = control->Evolve.dTimeStep;
   }
 
+  fvCumulative(body,control,system,dDt);
+
   /* Write out initial conditions */
   WriteOutput(body, control, files, output, system, update, fnWrite);
 
@@ -702,6 +720,8 @@ void Evolve(BODY *body, CONTROL *control, FILES *files, MODULE *module,
        was prior to loop. */
     PropertiesAuxiliary(body, control, system, update);
 
+  fvCumulative(body,control,system,dDt);
+
     // If control->Evolve.bFirstStep hasn't been switched off by now, do so.
     if (control->Evolve.bFirstStep) {
       control->Evolve.bFirstStep = 0;

src/evolve.h

@@ -16,6 +16,7 @@
 
 /* @cond DOXYGEN_OVERRIDE */
 
+double fdTrapezoidalArea(double,double,double);
 void PropertiesAuxiliary(BODY *, CONTROL *, SYSTEM *, UPDATE *);
 void fdGetUpdateInfo(BODY *, CONTROL *, SYSTEM *, UPDATE *,
                      fnUpdateVariable ***);