diff --git a/examples/CosmicShoreline/README.rst b/examples/CosmicShoreline/README.rst index ff41549f..628f8ebb 100644 --- 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 index d0374ceb..31db62cd 100644 --- 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 index 61850837..28e6c178 100644 --- 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 index d3c34bc2..24e07eec 100644 --- 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 index 67b9ead8..7a9c9e0d 100644 --- 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 index d7a13d95..4debc018 100644 --- 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 index e7a47758..8da0ec95 100644 --- 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 index 979218df..de7e74ae 100644 --- 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 index 7fef281d..5254830d 100644 --- 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 index 69dfa656..deee3500 100644 --- 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 index 89a0e740..2eaabd1d 100644 --- 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 index e4cf87ef..435b0b2e 100644 --- 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;iBodyEvolve.iNumBodies;iBody++) { - if (body[iBody].bAtmEsc) { - fvCumulativeAtmEsc(body,&control->Evolve,system,dDt,iBody); + if (body[0].bStellar) { + for (iBody=1;iBodyEvolve.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 index 354386f2..34528cd3 100644 --- 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 index 21a59810..321a0333 100644 --- a/src/system.c +++ b/src/system.c @@ -1114,4 +1114,19 @@ double fdBondiRadius(BODY *body, int iBody) { dBondiRadius = -1; } return dBondiRadius; -} \ No newline at end of file +} + +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 index c38dabb4..686bca3d 100644 --- 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 */