Added Cumulative XUV flux, but somehow tests are not passing

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,37 @@
+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 and the blue points
+represent the Solar System planets. The Sun's XUV luminosity
+is assumed to follow the Ribas et al. (2005) model.

examples/CosmicShoreline/earth.in

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

examples/CosmicShoreline/george.in

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

examples/CosmicShoreline/jupiter.in

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

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.CumulativeFXUV
+
+fxuv = []
+fxuv.append(output.log.final.Mercury.CumulativeFXUV/fxuv_earth)
+fxuv.append(output.log.final.Venus.CumulativeFXUV/fxuv_earth)
+fxuv.append(output.log.final.Earth.CumulativeFXUV/fxuv_earth)
+fxuv.append(output.log.final.Mars.CumulativeFXUV/fxuv_earth)
+fxuv.append(output.log.final.Jupiter.CumulativeFXUV/fxuv_earth)
+fxuv.append(output.log.final.Saturn.CumulativeFXUV/fxuv_earth)
+fxuv.append(output.log.final.George.CumulativeFXUV/fxuv_earth)
+fxuv.append(output.log.final.Neptune.CumulativeFXUV/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 -CumulativeFXUV

examples/CosmicShoreline/mercury.in

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

examples/CosmicShoreline/neptune.in

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

examples/CosmicShoreline/saturn.in

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

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,29 @@
+# Planet a parameters
+sName         Venus                      # Body's name
+saModules     atmesc                      # Modules to apply, exact spelling required
+
+# Physical Properties
+dMass         -0.815                      # Mass, negative -> Earth masses
+dRadius       -0.9499                     # Radius, negative -> Earth radii
+dRotPeriod    -243.                       # Rotation period, negative -> days
+dObliquity    180.                        # Retrograde rotation
+dRadGyra      0.5                         # Radius of gyration (moment of inertia constant)
+
+# ATMESC Properties
+dXFrac                    1.0             # X-Ray/XUV absorption radius (fraction of planet radius)
+dSurfWaterMass            -3.0            # Initial surface water (Earth oceans)
+dEnvelopeMass             0               # Initial envelope mass (Earth masses)
+bHaltSurfaceDesiccated    0               # Halt when dry?
+bHaltEnvelopeGone         0               # Halt when evaporated?
+dMinSurfWaterMass         -1.e-5          # Planet is desiccated when water content drops below this (Earth oceans)
+sWaterLossModel           lbexact
+sPlanetRadiusModel        none
+bInstantO2Sink            1
+sAtmXAbsEffH2OModel       bolmont16
+
+# Orbital Properties
+dSemi         -0.723                      # Semi-major axis, negative -> AU
+dEcc          0.006772                    # Eccentricity
+
+# Output
+saOutputOrder Time -FXUV -CumulativeFXUV

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

@@ -2002,6 +2002,17 @@ void fnPropsAuxAtmEsc(BODY *body, EVOLVE *evolve, IO *io, UPDATE *update,
 }
 
 
+void fvCumulativeAtmEsc(BODY *body,EVOLVE *evolve,SYSTEM *system,double dDt,int iBody) {
+  if (evolve->bFirstStep) {
+    body[iBody].dFXUVCumulative = 0;
+    body[iBody].dFXUVLast = fdXUVFlux(body,iBody);
+  } else {
+    body[iBody].dFXUV = fdXUVFlux(body,iBody);
+    body[iBody].dFXUVCumulative += fdTrapezoidalArea(body[iBody].dFXUV,body[iBody].dFXUVLast,dDt);
+    body[iBody].dFXUVLast = body[iBody].dFXUV;
+  }
+}
+
 /**
 Assigns functions returning the time-derivatives of each variable
 to the magical matrix of function pointers.
@@ -3454,6 +3465,20 @@ void WriteHRefODragMod(BODY *body, CONTROL *control, OUTPUT *output,
   fvFormattedString(cUnit, "");
 }
 
+void WriteCumulativeFXUV(BODY *body, CONTROL *control, OUTPUT *output,
+                       SYSTEM *system, UNITS *units, UPDATE *update, int iBody,
+                       double *dTmp, char **cUnit) {
+  *dTmp = body[iBody].dFXUVCumulative;
+
+    if (output->bDoNeg[iBody]) {
+    *dTmp *= output->dNeg;
+    fvFormattedString(cUnit, output->cNeg);
+  } else {
+    *dTmp *= fdUnitsEnergyFlux(units->iTime,units->iMass,units->iLength);
+    fsUnitsEnergyFlux(units,cUnit);
+  }
+}
+
 
 /**
 Logs the molecular oxygen mixing ratio.
@@ -3693,6 +3718,16 @@ void InitializeOutputAtmEsc(OUTPUT *output, fnWriteOutput fnWrite[]) {
   output[OUT_FXUV].iModuleBit = ATMESC;
   fnWrite[OUT_FXUV]           = &WriteFXUV;
 
+  fvFormattedString(&output[OUT_CUMULATIVEFXUV].cName, "CumulativeFXUV");
+  fvFormattedString(&output[OUT_CUMULATIVEFXUV].cDescr, "Cumulative XUV flux");
+  fvFormattedString(&output[OUT_CUMULATIVEFXUV].cNeg, "W/m^2");
+  output[OUT_CUMULATIVEFXUV].bNeg       = 1;
+  output[OUT_CUMULATIVEFXUV].dNeg       = 1;
+  output[OUT_CUMULATIVEFXUV].iNum       = 1;
+  output[OUT_CUMULATIVEFXUV].iModuleBit = ATMESC;
+  fnWrite[OUT_CUMULATIVEFXUV]           = &WriteCumulativeFXUV;
+
+
   fvFormattedString(&output[OUT_HESCAPEREGIME].cName, "HEscapeRegime");
   fvFormattedString(&output[OUT_HESCAPEREGIME].cDescr,
           "Integer flag for H envelope escape regime");

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,8 @@ 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);
+void fvCumulativeAtmEsc(BODY *,EVOLVE *,SYSTEM *,double,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;
+
+  for (iBody=0;iBody<control->Evolve.iNumBodies;iBody++) {
+    if (body[iBody].bAtmEsc) {
+      fvCumulativeAtmEsc(body,&control->Evolve,system,dDt,iBody);
+    }
+  }
+
+}
+
 void CalculateDerivatives(BODY *body, SYSTEM *system, UPDATE *update,
                           fnUpdateVariable ***fnUpdate, int iNumBodies) {
   int iBody, iVar, iEqn;
@@ -669,6 +685,8 @@ void Evolve(BODY *body, CONTROL *control, FILES *files, MODULE *module,
       }
     }
 
+    fvCumulative(body,control,system,dDt);
+
     fdGetUpdateInfo(body, control, system, update, fnUpdate);
 
     /* Halt? */

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 ***);