Moved CumulativeXUVFlux to system.c and no longer require AtmEsc to o…

…utput. Small tweaks to examples/ComsicShoreline.
VirtualPlanetaryLaboratory · Sep 13, 2024 · f706ce9 · f706ce9
1 parent 3baad6e
commit f706ce9
Show file tree

Hide file tree

Showing 15 changed files with 66 additions and 82 deletions.
diff --git a/examples/CosmicShoreline/README.rst b/examples/CosmicShoreline/README.rst
@@ -15,8 +15,6 @@ a proposed boundary between planet with and without atmosphere.
 ===================   ============
 
 
-
-
 To run this example
 -------------------
 
@@ -32,6 +30,9 @@ Expected output
    :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.
+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.
diff --git a/examples/CosmicShoreline/earth.in b/examples/CosmicShoreline/earth.in
@@ -3,4 +3,4 @@ saModules    atmesc
 dMass        3.0018090452e-06
 dRadius      -1.
 dSemi        1.00000321
-saOutputOrder    Time -FXUV -CumulativeFXUV
+saOutputOrder    Time -FXUV -CumulativeXUVFlux
diff --git a/examples/CosmicShoreline/george.in b/examples/CosmicShoreline/george.in
@@ -3,4 +3,4 @@ saModules    atmesc
 dMass        -14.54
 dRadius      -4.007
 dSemi        19.19203990
-saOutputOrder    Time -fXUV -CumulativeFXUV
+saOutputOrder    Time -fXUV -CumulativeXUVFlux
diff --git a/examples/CosmicShoreline/jupiter.in b/examples/CosmicShoreline/jupiter.in
@@ -3,4 +3,4 @@ saModules    atmesc
 dMass        -317.83
 dRadius      -11.2
 dSemi        5.20880408
-saOutputOrder    Time -FXUV -CumulativeFXUV
+saOutputOrder    Time -FXUV -CumulativeXUVFlux
diff --git a/examples/CosmicShoreline/mars.in b/examples/CosmicShoreline/mars.in
@@ -3,4 +3,4 @@ saModules    atmesc
 dMass        -0.107
 dRadius      -0.53202
 dSemi        1.52366290
-saOutputOrder    Time -FXUV -CumulativeFXUV
+saOutputOrder    Time -FXUV -CumulativeXUVFlux
diff --git a/examples/CosmicShoreline/mercury.in b/examples/CosmicShoreline/mercury.in
@@ -3,4 +3,4 @@ saModules    atmesc
 dMass        -0.0553
 dRadius      -0.383
 dSemi        0.38709897
-saOutputOrder    Time -FXUV -CumulativeFXUV
+saOutputOrder    Time -FXUV -CumulativeXUVFlux
diff --git a/examples/CosmicShoreline/neptune.in b/examples/CosmicShoreline/neptune.in
@@ -3,4 +3,4 @@ saModules    atmesc
 dMass        -17.15
 dRadius      -3.883
 dSemi        30.07050641
-saOutputOrder    Time -FXUV -CumulativeFXUV
+saOutputOrder    Time -FXUV -CumulativeXUVFlux
diff --git a/examples/CosmicShoreline/saturn.in b/examples/CosmicShoreline/saturn.in
@@ -3,4 +3,4 @@ saModules    atmesc
 dMass        -95.16
 dRadius      -9.449
 dSemi       9.53999265
-saOutputOrder    Time -FXUV -CumulativeFXUV
+saOutputOrder    Time -FXUV -CumulativeXUVFlux
diff --git a/examples/CosmicShoreline/venus.in b/examples/CosmicShoreline/venus.in
@@ -1,29 +1,7 @@
 # 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
+saModules     atmesc
+dMass         -0.815
+dRadius       -0.9499
+dSemi         -0.723
+saOutputOrder Time -FXUV -CumulativeXUVFlux
diff --git a/src/atmesc.c b/src/atmesc.c
@@ -2001,18 +2001,6 @@ 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.
@@ -3465,20 +3453,6 @@ 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.
@@ -3718,16 +3692,6 @@ 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");

diff --git a/src/atmesc.h b/src/atmesc.h
@@ -213,7 +213,6 @@ 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 *);

diff --git a/src/evolve.c b/src/evolve.c
@@ -58,9 +58,9 @@ 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);
+  if (body[0].bStellar) {
+    for (iBody=1;iBody<control->Evolve.iNumBodies;iBody++) {
+      fvCumulativeXUVFlux(body,&control->Evolve,system,dDt,iBody);
     }
   }
 
@@ -650,6 +650,9 @@ 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);
 
@@ -685,8 +688,6 @@ void Evolve(BODY *body, CONTROL *control, FILES *files, MODULE *module,
       }
     }
 
-    fvCumulative(body,control,system,dDt);
-
     fdGetUpdateInfo(body, control, system, update, fnUpdate);
 
     /* Halt? */

diff --git a/src/output.c b/src/output.c
@@ -82,6 +82,20 @@ void WriteCriticalSemi(BODY *body, CONTROL *control, OUTPUT *output,
   }
 }
 
+void WriteCumulativeXUVFlux(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);
+  }
+}
+
 /*
  * D
  */
@@ -1298,6 +1312,15 @@ void InitializeOutputGeneral(OUTPUT *output, fnWriteOutput fnWrite[]) {
         "in unstable orbits. This output parameter prints the instantaneous\n"
         "value of that critical distance.");
 
+  fvFormattedString(&output[OUT_CUMULATIVEFXUV].cName, "CumulativeXUVFlux");
+  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 = 1;
+  fnWrite[OUT_CUMULATIVEFXUV]           = &WriteCumulativeXUVFlux;
+
   /*
    * D
    */

diff --git a/src/system.c b/src/system.c
@@ -1114,4 +1114,19 @@ double fdBondiRadius(BODY *body, int iBody) {
     dBondiRadius = -1;
   }
   return dBondiRadius;
-}
+}
+
+void fvCumulativeXUVFlux(BODY *body,EVOLVE *evolve,SYSTEM *system,double dDt,int iBody) {
+  if (evolve->bFirstStep) {
+    body[iBody].dFXUVCumulative = 0;
+    if (body[iBody].bCalcFXUV) {
+      body[iBody].dFXUV = fdXUVFlux(body, iBody);
+    }
+  } else {
+    if (body[iBody].bCalcFXUV) {
+      body[iBody].dFXUV = fdXUVFlux(body, iBody);
+    }
+    body[iBody].dFXUVCumulative += fdTrapezoidalArea(body[iBody].dFXUV,body[iBody].dFXUVLast,dDt);
+    body[iBody].dFXUVLast = body[iBody].dFXUV;
+  }
+}
diff --git a/src/system.h b/src/system.h
@@ -56,4 +56,7 @@ void fdMergePlanet(BODY *, UPDATE *, fnUpdateVariable ***, int);
 double fdBondiRadius(BODY *, int);
 double fdRocheRadius(BODY *, int, int);
 
+void fvCumulativeXUVFlux(BODY *,EVOLVE *,SYSTEM *,double,int);
+
+
 /* @endcond */