diff --git a/src/flare.c b/src/flare.c index 624b47e0f..be9d8d2ea 100644 --- a/src/flare.c +++ b/src/flare.c @@ -26,6 +26,13 @@ void BodyCopyFlare(BODY *dest, BODY *src, int foo, int iNumBodies, int iBody) { dest[iBody].iEnergyBin = src[iBody].iEnergyBin; // dest[iBody].dFlareSlope = src[iBody].dFlareSlope; // dest[iBody].dFlareYInt = src[iBody].dFlareYInt; + + dest[iBody].dFlareSlopeA1 = src[iBody].dFlareSlopeA1; + dest[iBody].dFlareSlopeA2 = src[iBody].dFlareSlopeA2; + dest[iBody].dFlareSlopeA3 = src[iBody].dFlareSlopeA3; + dest[iBody].dFlareYIntB1 = src[iBody].dFlareYIntB1; + dest[iBody].dFlareYIntB2 = src[iBody].dFlareYIntB2; + dest[iBody].dFlareYIntB3 = src[iBody].dFlareYIntB3; } /**************** FLARE options ********************/ @@ -125,7 +132,7 @@ void ReadFlareBandPass(BODY *body, body[iFile - 1].iFlareBandPass = FLARE_UV; } else if (!memcmp(sLower(cTmp), "go", 2)) { body[iFile - 1].iFlareBandPass = FLARE_GOES; - } else if (!memcmp(sLower(cTmp), "sr", 2)) { + } else if (!memcmp(sLower(cTmp), "sx", 2)) { body[iFile - 1].iFlareBandPass = FLARE_SXR; } else if (!memcmp(sLower(cTmp), "te", 2)) { body[iFile - 1].iFlareBandPass = FLARE_TESS_UV; @@ -318,6 +325,109 @@ void ReadFlareSlope(BODY *body, CONTROL *control, FILES *files, } } +void ReadFlareSlopeA1(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dFlareSlopeA1 = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dFlareSlopeA1 = options->dDefault; + } +} + +void ReadFlareSlopeA2(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dFlareSlopeA2 = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dFlareSlopeA2 = options->dDefault; + } +} + +void ReadFlareSlopeA3(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dFlareSlopeA3 = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dFlareSlopeA3 = options->dDefault; + } +} + +void ReadFlareYIntB1(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dFlareYIntB1 = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dFlareYIntB1 = options->dDefault; + } +} + +void ReadFlareYIntB2(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dFlareYIntB2 = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dFlareYIntB2 = options->dDefault; + } +} + +void ReadFlareYIntB3(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dFlareYIntB3 = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dFlareYIntB3 = options->dDefault; + } +} + + // TODO: Include the error in the FFD slopes to calculate the upper and higher // limit of XUV luminosity by flares @@ -618,6 +728,90 @@ void InitializeOptionsFlare(OPTIONS *options, fnReadOption fnRead[]) { "MIN for number of flares per minutes, \n" " HOUR for number of flares per hour, and DAY number of flares per " "for days.\n");*/ + + fvFormattedString(&options[OPT_FLARESLOPEA1].cName, "dFlareSlopeA1"); + fvFormattedString(&options[OPT_FLARESLOPEA1].cDescr, + "a1 coefficient in Davenport (2019) variable-slope FFD: " + "slope = a1*log10(age_Myr) + a2*M_sun + a3"); + fvFormattedString(&options[OPT_FLARESLOPEA1].cDefault, + "-0.07 (Davenport+2019 M-dwarf fit)"); + options[OPT_FLARESLOPEA1].dDefault = -0.07; + options[OPT_FLARESLOPEA1].iType = 2; + options[OPT_FLARESLOPEA1].bMultiFile = 1; + options[OPT_FLARESLOPEA1].iModuleBit = FLARE; + options[OPT_FLARESLOPEA1].bNeg = 0; + fvFormattedString(&options[OPT_FLARESLOPEA1].cDimension, "nd"); + fnRead[OPT_FLARESLOPEA1] = &ReadFlareSlopeA1; + + fvFormattedString(&options[OPT_FLARESLOPEA2].cName, "dFlareSlopeA2"); + fvFormattedString(&options[OPT_FLARESLOPEA2].cDescr, + "a2 coefficient in Davenport (2019) variable-slope FFD; " + "see dFlareSlopeA1"); + fvFormattedString(&options[OPT_FLARESLOPEA2].cDefault, + "0.79 (Davenport+2019 M-dwarf fit)"); + options[OPT_FLARESLOPEA2].dDefault = 0.79; + options[OPT_FLARESLOPEA2].iType = 2; + options[OPT_FLARESLOPEA2].bMultiFile = 1; + options[OPT_FLARESLOPEA2].iModuleBit = FLARE; + options[OPT_FLARESLOPEA2].bNeg = 0; + fvFormattedString(&options[OPT_FLARESLOPEA2].cDimension, "nd"); + fnRead[OPT_FLARESLOPEA2] = &ReadFlareSlopeA2; + + fvFormattedString(&options[OPT_FLARESLOPEA3].cName, "dFlareSlopeA3"); + fvFormattedString(&options[OPT_FLARESLOPEA3].cDescr, + "a3 coefficient in Davenport (2019) variable-slope FFD; " + "see dFlareSlopeA1"); + fvFormattedString(&options[OPT_FLARESLOPEA3].cDefault, + "-1.06 (Davenport+2019 M-dwarf fit)"); + options[OPT_FLARESLOPEA3].dDefault = -1.06; + options[OPT_FLARESLOPEA3].iType = 2; + options[OPT_FLARESLOPEA3].bMultiFile = 1; + options[OPT_FLARESLOPEA3].iModuleBit = FLARE; + options[OPT_FLARESLOPEA3].bNeg = 0; + fvFormattedString(&options[OPT_FLARESLOPEA3].cDimension, "nd"); + fnRead[OPT_FLARESLOPEA3] = &ReadFlareSlopeA3; + + fvFormattedString(&options[OPT_FLAREYINTB1].cName, "dFlareYIntB1"); + fvFormattedString(&options[OPT_FLAREYINTB1].cDescr, + "b1 coefficient in Davenport (2019) variable-slope FFD: " + "yint = b1*log10(age_Myr) + b2*M_sun + b3"); + fvFormattedString(&options[OPT_FLAREYINTB1].cDefault, + "2.01 (Davenport+2019 M-dwarf fit)"); + options[OPT_FLAREYINTB1].dDefault = 2.01; + options[OPT_FLAREYINTB1].iType = 2; + options[OPT_FLAREYINTB1].bMultiFile = 1; + options[OPT_FLAREYINTB1].iModuleBit = FLARE; + options[OPT_FLAREYINTB1].bNeg = 0; + fvFormattedString(&options[OPT_FLAREYINTB1].cDimension, "nd"); + fnRead[OPT_FLAREYINTB1] = &ReadFlareYIntB1; + + fvFormattedString(&options[OPT_FLAREYINTB2].cName, "dFlareYIntB2"); + fvFormattedString(&options[OPT_FLAREYINTB2].cDescr, + "b2 coefficient in Davenport (2019) variable-slope FFD; " + "see dFlareYIntB1"); + fvFormattedString(&options[OPT_FLAREYINTB2].cDefault, + "-25.15 (Davenport+2019 M-dwarf fit)"); + options[OPT_FLAREYINTB2].dDefault = -25.15; + options[OPT_FLAREYINTB2].iType = 2; + options[OPT_FLAREYINTB2].bMultiFile = 1; + options[OPT_FLAREYINTB2].iModuleBit = FLARE; + options[OPT_FLAREYINTB2].bNeg = 0; + fvFormattedString(&options[OPT_FLAREYINTB2].cDimension, "nd"); + fnRead[OPT_FLAREYINTB2] = &ReadFlareYIntB2; + + fvFormattedString(&options[OPT_FLAREYINTB3].cName, "dFlareYIntB3"); + fvFormattedString(&options[OPT_FLAREYINTB3].cDescr, + "b3 coefficient in Davenport (2019) variable-slope FFD; " + "see dFlareYIntB1"); + fvFormattedString(&options[OPT_FLAREYINTB3].cDefault, + "33.99 (Davenport+2019 M-dwarf fit)"); + options[OPT_FLAREYINTB3].dDefault = 33.99; + options[OPT_FLAREYINTB3].iType = 2; + options[OPT_FLAREYINTB3].bMultiFile = 1; + options[OPT_FLAREYINTB3].iModuleBit = FLARE; + options[OPT_FLAREYINTB3].bNeg = 0; + fvFormattedString(&options[OPT_FLAREYINTB3].cDimension, "nd"); + fnRead[OPT_FLAREYINTB3] = &ReadFlareYIntB3; } void ReadOptionsFlare(BODY *body, @@ -1008,6 +1202,20 @@ void WriteLXUVFlare(BODY *body, } } +void WriteFlareSlopeOUT(BODY *body, CONTROL *control, OUTPUT *output, + SYSTEM *system, UNITS *units, UPDATE *update, + int iBody, double *dTmp, char **cUnit) { + *dTmp = fdFlareSlopeOUT(body, iBody); + fvFormattedString(cUnit, ""); +} + +void WriteFlareYIntOUT(BODY *body, CONTROL *control, OUTPUT *output, + SYSTEM *system, UNITS *units, UPDATE *update, + int iBody, double *dTmp, char **cUnit) { + *dTmp = fdFlareYIntOUT(body, iBody); + fvFormattedString(cUnit, ""); +} + // TODO: Include the error in the FFD slopes to calculate the upper and higher // limit of XUV luminosity by flares /* @@ -1455,6 +1663,22 @@ void InitializeOutputFlare(OUTPUT *output, fnWriteOutput fnWrite[]) { output[OUT_LXUVFLARE].iModuleBit = FLARE; fnWrite[OUT_LXUVFLARE] = &WriteLXUVFlare; + fvFormattedString(&output[OUT_FLARESLOPEOUT].cName, "FlareSlopeOUT"); + fvFormattedString(&output[OUT_FLARESLOPEOUT].cDescr, + "Instantaneous Davenport FFD slope at current age and mass"); + output[OUT_FLARESLOPEOUT].bNeg = 0; + output[OUT_FLARESLOPEOUT].iNum = 1; + output[OUT_FLARESLOPEOUT].iModuleBit = FLARE; + fnWrite[OUT_FLARESLOPEOUT] = &WriteFlareSlopeOUT; + + fvFormattedString(&output[OUT_FLAREYINTOUT].cName, "FlareYIntOUT"); + fvFormattedString(&output[OUT_FLAREYINTOUT].cDescr, + "Instantaneous Davenport FFD y-intercept at current age and mass"); + output[OUT_FLAREYINTOUT].bNeg = 0; + output[OUT_FLAREYINTOUT].iNum = 1; + output[OUT_FLAREYINTOUT].iModuleBit = FLARE; + fnWrite[OUT_FLAREYINTOUT] = &WriteFlareYIntOUT; + // TODO: Include the error in the FFD slopes to calculate the upper and higher // limit of XUV luminosity by flares /*fvFormattedString(output[OUT_LXUVFLAREUPPER].cName, "LXUVFlareUpper"); @@ -1642,6 +1866,25 @@ double fdDavenport(double dA1, double dA2, double dA3, double dStarAge, return dA; } + +/* Instantaneous Davenport FFD slope for the current age and mass, evaluated + from the body's variable-slope coefficients. Computed on demand so that + output does not require storing redundant state. */ +double fdFlareSlopeOUT(BODY *body, int iBody) { + return fdDavenport(body[iBody].dFlareSlopeA1, + body[iBody].dFlareSlopeA2, + body[iBody].dFlareSlopeA3, + body[iBody].dAge, body[iBody].dMass); +} + +/* Instantaneous Davenport FFD y-intercept for the current age and mass. */ +double fdFlareYIntOUT(BODY *body, int iBody) { + return fdDavenport(body[iBody].dFlareYIntB1, + body[iBody].dFlareYIntB2, + body[iBody].dFlareYIntB3, + body[iBody].dAge, body[iBody].dMass); +} + // fdFFD calculates FFD of the flares. If LACY mode is choosen, them fdFFD has // to convert from SI to the units that the equation 3 from Davenport et al. // 2019 understand (i.e., flares/day and flares/(day*log10(erg))). @@ -1695,14 +1938,17 @@ double fdLXUVFlare(BODY *body, double dDeltaTime, int iBody) { // slopes(constant)?################################## if (body[iBody].iFlareFFD == FLARE_FFD_DAVENPORT) { - // The coefficient values given here were given by Dr. James Davenport in - // private comunication - dFlareSlope = - fdDavenport(-0.07, 0.79, -1.06, body[iBody].dAge, - body[iBody].dMass); //(-0.07054598,0.81225239,-1.07054511) - dFlareYInt = fdDavenport( - 2.01, -25.15, 33.99, body[iBody].dAge, - body[iBody].dMass); //(2.06012734,-25.79885288,34.44115635) + /* Davenport+2019 variable-slope FFD: slope and y-intercept are each + third-order fits in log10(age_Myr) with a linear mass term. The six + coefficients default to Davenport et al.'s published M-dwarf values. */ + dFlareSlope = fdDavenport(body[iBody].dFlareSlopeA1, + body[iBody].dFlareSlopeA2, + body[iBody].dFlareSlopeA3, + body[iBody].dAge, body[iBody].dMass); + dFlareYInt = fdDavenport(body[iBody].dFlareYIntB1, + body[iBody].dFlareYIntB2, + body[iBody].dFlareYIntB3, + body[iBody].dAge, body[iBody].dMass); } else if (body[iBody].iFlareFFD == FLARE_FFD_LACY) { dFlareSlope = body[iBody].dFlareSlope; dFlareYInt = body[iBody].dFlareYInt; diff --git a/src/flare.h b/src/flare.h index 8d8ca979c..f57c5d34c 100644 --- a/src/flare.h +++ b/src/flare.h @@ -46,6 +46,13 @@ #define OPT_FLAREBANDPASS 2041 #define OPT_LXUVFLARECONST 2042 +#define OPT_FLARESLOPEA1 2043 +#define OPT_FLARESLOPEA2 2044 +#define OPT_FLARESLOPEA3 2045 +#define OPT_FLAREYINTB1 2046 +#define OPT_FLAREYINTB2 2047 +#define OPT_FLAREYINTB3 2048 + /* Output Functinos */ /* FLARE 1900 - 1999 */ @@ -70,6 +77,10 @@ #define OUT_FLAREENERGYMIN 2024 #define OUT_FLAREENERGYMID 2025 #define OUT_FLAREENERGYMAX 2026 + +#define OUT_FLAREYINTOUT 2046 +#define OUT_FLARESLOPEOUT 2047 + /* @cond DOXYGEN_OVERRIDE */ void InitializeControlFlare(CONTROL *); @@ -137,6 +148,10 @@ void WriteFlareEnergyMid(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, int, double *, char**); void WriteFlareEnergyMax(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, int, double *, char**); +void WriteFlareSlopeOUT(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, + UPDATE *, int, double *, char**); +void WriteFlareYIntOUT(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, + UPDATE *, int, double *, char**); /* Logging Functions */ void LogOptionsFlare(CONTROL *, FILE *); @@ -148,6 +163,8 @@ void LogBodyFlare(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UPDATE *, /* FLARE functions */ double fdLXUVFlare(BODY *, double, int); double fdDavenport(double, double, double, double, double); +double fdFlareSlopeOUT(BODY *, int); +double fdFlareYIntOUT(BODY *, int); double fdFFD(BODY *, int, double, double, double); double fdBandPassKepler(BODY *, int, double); double fdBandPassXUV(BODY *, int, double); diff --git a/src/stellar.c b/src/stellar.c index 896746b00..740b661dc 100644 --- a/src/stellar.c +++ b/src/stellar.c @@ -33,6 +33,31 @@ void BodyCopyStellar(BODY *dest, BODY *src, int foo, int iNumBodies, dest[iBody].dLuminosityAmplitude = src[iBody].dLuminosityAmplitude; dest[iBody].dLuminosityFrequency = src[iBody].dLuminosityFrequency; dest[iBody].dLuminosityPhase = src[iBody].dLuminosityPhase; + + dest[iBody].dXUVEngleEarlyA = src[iBody].dXUVEngleEarlyA; + dest[iBody].dXUVEngleEarlyB = src[iBody].dXUVEngleEarlyB; + dest[iBody].dXUVEngleEarlyC = src[iBody].dXUVEngleEarlyC; + dest[iBody].dXUVEngleEarlyD = src[iBody].dXUVEngleEarlyD; + + dest[iBody].dXUVEngleMidLateA = src[iBody].dXUVEngleMidLateA; + dest[iBody].dXUVEngleMidLateB = src[iBody].dXUVEngleMidLateB; + dest[iBody].dXUVEngleMidLateC = src[iBody].dXUVEngleMidLateC; + dest[iBody].dXUVEngleMidLateD = src[iBody].dXUVEngleMidLateD; + + dest[iBody].dRotEngleEarlyA = src[iBody].dRotEngleEarlyA; + dest[iBody].dRotEngleEarlyB = src[iBody].dRotEngleEarlyB; + dest[iBody].dRotEngleEarlyC = src[iBody].dRotEngleEarlyC; + dest[iBody].dRotEngleEarlyD = src[iBody].dRotEngleEarlyD; + + dest[iBody].dRotEngleMidA = src[iBody].dRotEngleMidA; + dest[iBody].dRotEngleMidB = src[iBody].dRotEngleMidB; + dest[iBody].dRotEngleMidC = src[iBody].dRotEngleMidC; + dest[iBody].dRotEngleMidD = src[iBody].dRotEngleMidD; + + dest[iBody].dRotEngleLateA = src[iBody].dRotEngleLateA; + dest[iBody].dRotEngleLateB = src[iBody].dRotEngleLateB; + dest[iBody].dRotEngleLateC = src[iBody].dRotEngleLateC; + dest[iBody].dRotEngleLateD = src[iBody].dRotEngleLateD; } /**************** STELLAR options ********************/ @@ -165,11 +190,18 @@ void ReadMagBrakingModel(BODY *body, CONTROL *control, FILES *files, body[iFile - 1].iMagBrakingModel = STELLAR_DJDT_MA15; } else if (!memcmp(sLower(cTmp), "br", 2)) { body[iFile - 1].iMagBrakingModel = STELLAR_DJDT_BR21; + } else if (!memcmp(sLower(cTmp), "engle23e", 8)) { + body[iFile - 1].iMagBrakingModel = STELLAR_DJDT_ENGLE23EARLY; + } else if (!memcmp(sLower(cTmp), "engle23m", 8)) { + body[iFile - 1].iMagBrakingModel = STELLAR_DJDT_ENGLE23MID; + } else if (!memcmp(sLower(cTmp), "engle23l", 8)) { + body[iFile - 1].iMagBrakingModel = STELLAR_DJDT_ENGLE23LATE; } else { if (control->Io.iVerbose >= VERBERR) { fprintf(stderr, "ERROR: Unknown argument to %s: %s. Options are REINERS, " - "SKUMANICH, MATT, BREIMANN21 or NONE.\n", + "SKUMANICH, MATT, BREIMANN21, ENGLE23EARLY, ENGLE23MID, " + "ENGLE23LATE, or NONE.\n", options->cName, cTmp); } LineExit(files->Infile[iFile].cIn, lTmp); @@ -229,14 +261,19 @@ void ReadXUVModel(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, } else if (!memcmp(sLower(cTmp), "re", 2)) { if (control->Io.iVerbose >= VERBINPUT) { fprintf(stderr, "WARNING: The REINERS XUV model has serious issues. " - "The recommended model is RIBAS.\n"); + "The recommended models are RIBAS, ENGLE24EARLY, or " + "ENGLE24MIDLATE.\n"); } body[iFile - 1].iXUVModel = STELLAR_MODEL_REINERS; + } else if (!memcmp(sLower(cTmp), "engle24e", 8)) { + body[iFile - 1].iXUVModel = STELLAR_MODEL_ENGLE24EARLY; + } else if (!memcmp(sLower(cTmp), "engle24m", 8)) { + body[iFile - 1].iXUVModel = STELLAR_MODEL_ENGLE24MIDLATE; } else { if (control->Io.iVerbose >= VERBERR) { fprintf(stderr, - "ERROR: Unknown argument to %s: %s. Options are RIBAS, REINERS " - "or NONE.\n", + "ERROR: Unknown argument to %s: %s. Options are RIBAS, REINERS, " + "ENGLE24EARLY, ENGLE24MIDLATE, or NONE.\n", options->cName, cTmp); } LineExit(files->Infile[iFile].cIn, lTmp); @@ -396,6 +433,346 @@ void ReadLuminosityPhase(BODY *body, CONTROL *control, FILES *files, body[iFile - 1].dLuminosityPhase = options->dDefault; } +void ReadXUVEngleEarlyA(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dXUVEngleEarlyA = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dXUVEngleEarlyA = options->dDefault; + } +} + +void ReadXUVEngleEarlyB(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dXUVEngleEarlyB = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dXUVEngleEarlyB = options->dDefault; + } +} + +void ReadXUVEngleEarlyC(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dXUVEngleEarlyC = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dXUVEngleEarlyC = options->dDefault; + } +} + +void ReadXUVEngleEarlyD(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dXUVEngleEarlyD = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dXUVEngleEarlyD = options->dDefault; + } +} + +void ReadXUVEngleMidLateA(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dXUVEngleMidLateA = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dXUVEngleMidLateA = options->dDefault; + } +} + +void ReadXUVEngleMidLateB(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dXUVEngleMidLateB = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dXUVEngleMidLateB = options->dDefault; + } +} + +void ReadXUVEngleMidLateC(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dXUVEngleMidLateC = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dXUVEngleMidLateC = options->dDefault; + } +} + +void ReadXUVEngleMidLateD(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dXUVEngleMidLateD = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dXUVEngleMidLateD = options->dDefault; + } +} + +void ReadRotEngleEarlyA(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dRotEngleEarlyA = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dRotEngleEarlyA = options->dDefault; + } +} + +void ReadRotEngleEarlyB(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dRotEngleEarlyB = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dRotEngleEarlyB = options->dDefault; + } +} + +void ReadRotEngleEarlyC(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dRotEngleEarlyC = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dRotEngleEarlyC = options->dDefault; + } +} + +void ReadRotEngleEarlyD(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dRotEngleEarlyD = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dRotEngleEarlyD = options->dDefault; + } +} + +void ReadRotEngleMidA(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dRotEngleMidA = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dRotEngleMidA = options->dDefault; + } +} + +void ReadRotEngleMidB(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dRotEngleMidB = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dRotEngleMidB = options->dDefault; + } +} + +void ReadRotEngleMidC(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dRotEngleMidC = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dRotEngleMidC = options->dDefault; + } +} + +void ReadRotEngleMidD(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dRotEngleMidD = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dRotEngleMidD = options->dDefault; + } +} + +void ReadRotEngleLateA(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dRotEngleLateA = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dRotEngleLateA = options->dDefault; + } +} + +void ReadRotEngleLateB(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dRotEngleLateB = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dRotEngleLateB = options->dDefault; + } +} + +void ReadRotEngleLateC(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dRotEngleLateC = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dRotEngleLateC = options->dDefault; + } +} + +void ReadRotEngleLateD(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dRotEngleLateD = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else if (iFile > 0) { + body[iFile - 1].dRotEngleLateD = options->dDefault; + } +} + /* Halts */ void ReadHaltEndBaraffeGrid(BODY *body, CONTROL *control, FILES *files, @@ -635,6 +1012,286 @@ void InitializeOptionsStellar(OPTIONS *options, fnReadOption fnRead[]) { options[OPT_LUMPHASE].dNeg = DEGRAD; fvFormattedString(&options[OPT_LUMPHASE].cNeg, "Degrees"); fnRead[OPT_LUMPHASE] = &ReadLuminosityPhase; + + iOpt = OPT_XUVENGLEEARLYA; + fvFormattedString(&options[iOpt].cName, "dXUVEngleEarlyA"); + fvFormattedString(&options[iOpt].cDescr, + "Value of a in the Engle (2024) empirical age-XUV " + "luminosity relationship for M0-2 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-0.4896"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -0.4896; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadXUVEngleEarlyA; + + iOpt = OPT_XUVENGLEEARLYB; + fvFormattedString(&options[iOpt].cName, "dXUVEngleEarlyB"); + fvFormattedString(&options[iOpt].cDescr, + "Value of b in the Engle (2024) empirical age-XUV " + "luminosity relationship for M0-2 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-3.2128"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -3.2128; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadXUVEngleEarlyB; + + iOpt = OPT_XUVENGLEEARLYC; + fvFormattedString(&options[iOpt].cName, "dXUVEngleEarlyC"); + fvFormattedString(&options[iOpt].cDescr, + "Value of c in the Engle (2024) empirical age-XUV " + "luminosity relationship for M0-2 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-0.4469"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -0.4469; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadXUVEngleEarlyC; + + iOpt = OPT_XUVENGLEEARLYD; + fvFormattedString(&options[iOpt].cName, "dXUVEngleEarlyD"); + fvFormattedString(&options[iOpt].cDescr, + "Value of d in the Engle (2024) empirical age-XUV " + "luminosity relationship for M0-2 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-0.2985"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -0.2985; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadXUVEngleEarlyD; + + iOpt = OPT_XUVENGLEMIDLATEA; + fvFormattedString(&options[iOpt].cName, "dXUVEngleMidLateA"); + fvFormattedString(&options[iOpt].cDescr, + "Value of a in the Engle (2024) empirical age-XUV " + "luminosity relationship for M2.6-6.5 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-0.1456"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -0.1456; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadXUVEngleMidLateA; + + iOpt = OPT_XUVENGLEMIDLATEB; + fvFormattedString(&options[iOpt].cName, "dXUVEngleMidLateB"); + fvFormattedString(&options[iOpt].cDescr, + "Value of b in the Engle (2024) empirical age-XUV " + "luminosity relationship for M2.6-6.5 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-2.8876"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -2.8876; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadXUVEngleMidLateB; + + iOpt = OPT_XUVENGLEMIDLATEC; + fvFormattedString(&options[iOpt].cName, "dXUVEngleMidLateC"); + fvFormattedString(&options[iOpt].cDescr, + "Value of c in the Engle (2024) empirical age-XUV " + "luminosity relationship for M2.6-6.5 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-1.8187"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -1.8187; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadXUVEngleMidLateC; + + iOpt = OPT_XUVENGLEMIDLATED; + fvFormattedString(&options[iOpt].cName, "dXUVEngleMidLateD"); + fvFormattedString(&options[iOpt].cDescr, + "Value of d in the Engle (2024) empirical age-XUV " + "luminosity relationship for M2.6-6.5 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "0.3545"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = 0.3545; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadXUVEngleMidLateD; + + iOpt = OPT_ROTENGLEEARLYA; + fvFormattedString(&options[iOpt].cName, "dRotEngleEarlyA"); + fvFormattedString(&options[iOpt].cDescr, + "Value of a in the Engle & Guinan (2023) empirical " + "age-rotation relationship for M0-2 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "0.0621"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = 0.0621; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadRotEngleEarlyA; + + iOpt = OPT_ROTENGLEEARLYB; + fvFormattedString(&options[iOpt].cName, "dRotEngleEarlyB"); + fvFormattedString(&options[iOpt].cDescr, + "Value of b in the Engle & Guinan (2023) empirical " + "age-rotation relationship for M0-2 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-1.0437"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -1.0437; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadRotEngleEarlyB; + + iOpt = OPT_ROTENGLEEARLYC; + fvFormattedString(&options[iOpt].cName, "dRotEngleEarlyC"); + fvFormattedString(&options[iOpt].cDescr, + "Value of c in the Engle & Guinan (2023) empirical " + "age-rotation relationship for M0-2 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-0.0528"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -0.0528; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadRotEngleEarlyC; + + iOpt = OPT_ROTENGLEEARLYD; + fvFormattedString(&options[iOpt].cName, "dRotEngleEarlyD"); + fvFormattedString(&options[iOpt].cDescr, + "Value of d in the Engle & Guinan (2023) empirical " + "age-rotation relationship for M0-2 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-23.4933"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -23.4933; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadRotEngleEarlyD; + + iOpt = OPT_ROTENGLEMIDA; + fvFormattedString(&options[iOpt].cName, "dRotEngleMidA"); + fvFormattedString(&options[iOpt].cDescr, + "Value of a in the Engle & Guinan (2023) empirical " + "age-rotation relationship for M2.5-M3.5 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "0.0561"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = 0.0561; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadRotEngleMidA; + + iOpt = OPT_ROTENGLEMIDB; + fvFormattedString(&options[iOpt].cName, "dRotEngleMidB"); + fvFormattedString(&options[iOpt].cDescr, + "Value of b in the Engle & Guinan (2023) empirical " + "age-rotation relationship for M2.5-M3.5 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-0.8900"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -0.8900; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadRotEngleMidB; + + iOpt = OPT_ROTENGLEMIDC; + fvFormattedString(&options[iOpt].cName, "dRotEngleMidC"); + fvFormattedString(&options[iOpt].cDescr, + "Value of c in the Engle & Guinan (2023) empirical " + "age-rotation relationship for M2.5-M3.5 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-0.0521"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -0.0521; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadRotEngleMidC; + + iOpt = OPT_ROTENGLEMIDD; + fvFormattedString(&options[iOpt].cName, "dRotEngleMidD"); + fvFormattedString(&options[iOpt].cDescr, + "Value of d in the Engle & Guinan (2023) empirical " + "age-rotation relationship for M2.5-M3.5 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-24.1888"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -24.1888; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadRotEngleMidD; + + iOpt = OPT_ROTENGLELATEA; + fvFormattedString(&options[iOpt].cName, "dRotEngleLateA"); + fvFormattedString(&options[iOpt].cDescr, + "Value of a in the Engle & Guinan (2023) empirical " + "age-rotation relationship for M4-M6.5 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "0.0251"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = 0.0251; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadRotEngleLateA; + + iOpt = OPT_ROTENGLELATEB; + fvFormattedString(&options[iOpt].cName, "dRotEngleLateB"); + fvFormattedString(&options[iOpt].cDescr, + "Value of b in the Engle & Guinan (2023) empirical " + "age-rotation relationship for M4-M6.5 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-0.1615"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -0.1615; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadRotEngleLateB; + + iOpt = OPT_ROTENGLELATEC; + fvFormattedString(&options[iOpt].cName, "dRotEngleLateC"); + fvFormattedString(&options[iOpt].cDescr, + "Value of c in the Engle & Guinan (2023) empirical " + "age-rotation relationship for M4-M6.5 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-0.0212"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -0.0212; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadRotEngleLateC; + + iOpt = OPT_ROTENGLELATED; + fvFormattedString(&options[iOpt].cName, "dRotEngleLateD"); + fvFormattedString(&options[iOpt].cDescr, + "Value of d in the Engle & Guinan (2023) empirical " + "age-rotation relationship for M4-M6.5 dwarfs"); + fvFormattedString(&options[iOpt].cDefault, "-25.45"); + fvFormattedString(&options[iOpt].cDimension, "nd"); + options[iOpt].dDefault = -25.45; + options[iOpt].iType = 2; + options[iOpt].bMultiFile = 1; + options[iOpt].iModuleBit = STELLAR; + options[iOpt].bNeg = 0; + fnRead[iOpt] = &ReadRotEngleLateD; } void ReadOptionsStellar(BODY *body, CONTROL *control, FILES *files, @@ -1012,6 +1669,14 @@ void fnPropsAuxStellar(BODY *body, EVOLVE *evolve, IO *io, UPDATE *update, body[iBody].dLXUV = body[iBody].dSatXUVFrac * body[iBody].dLuminosity; } + } else if (body[iBody].iXUVModel == STELLAR_MODEL_ENGLE24EARLY) { + body[iBody].dLXUV = fdLXUVEngle( + body, body[iBody].dXUVEngleEarlyA, body[iBody].dXUVEngleEarlyB, + body[iBody].dXUVEngleEarlyC, body[iBody].dXUVEngleEarlyD, iBody); + } else if (body[iBody].iXUVModel == STELLAR_MODEL_ENGLE24MIDLATE) { + body[iBody].dLXUV = fdLXUVEngle( + body, body[iBody].dXUVEngleMidLateA, body[iBody].dXUVEngleMidLateB, + body[iBody].dXUVEngleMidLateC, body[iBody].dXUVEngleMidLateD, iBody); } else { // Constant XUV fraction @@ -1026,6 +1691,13 @@ void fnForceBehaviorStellar(BODY *body, MODULE *module, EVOLVE *evolve, IO *io, if (body[iBody].iStellarModel == STELLAR_MODEL_SINEWAVE) { body[iBody].dLXUV = body[iBody].dSatXUVFrac * body[iBody].dLuminosity; } + + if (body[iBody].iMagBrakingModel == STELLAR_DJDT_ENGLE23EARLY || + body[iBody].iMagBrakingModel == STELLAR_DJDT_ENGLE23MID || + body[iBody].iMagBrakingModel == STELLAR_DJDT_ENGLE23LATE) { + body[iBody].dRotPer = fdRotationPeriodEngle(body, iBody); + body[iBody].dRotRate = fdPerToFreq(body[iBody].dRotPer); + } } void AssignStellarDerivatives(BODY *body, EVOLVE *evolve, UPDATE *update, @@ -1072,11 +1744,184 @@ void NullStellarDerivatives(BODY *body, EVOLVE *evolve, UPDATE *update, } } +void VerifyTwoEngleOptionsNotSet(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, char cModel[OPTLEN], + int iOptionModel, int iOptionConstant, + int iBody) { + if (options[iOptionModel].iLine[iBody + 1] >= 0 && + options[iOptionConstant].iLine[iBody + 1] >= 0) { + if (control->Io.iVerbose >= VERBERR) { + fprintf(stderr, "ERROR: Cannot set %s if %s is set to %s.\n", + options[iOptionConstant].cName, options[iOptionModel].cName, + cModel); + } + DoubleLineExit(files->Infile[iBody + 1].cIn, files->Infile[iBody + 1].cIn, + options[iOptionConstant].iLine[iBody + 1], + options[iOptionModel].iLine[iBody + 1]); + } +} + +void VerifyEngleMassSpectralType(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, char cModel[LINE], + double dMaxMass, double dMinMass, + int iOptionModel, int iBody) { + + if (body[iBody].dMass > dMaxMass * MSUN || + body[iBody].dMass <= dMinMass * MSUN) { + if (control->Io.iVerbose >= VERBINPUT) { + fprintf(stderr, + "ERROR: Stellar mass must be between %.2e and %.2e solar " + "masses to use the %s model.", + dMinMass, dMaxMass, cModel); + } + DoubleLineExit(files->Infile[iBody + 1].cIn, files->Infile[iBody + 1].cIn, + options[iOptionModel].iLine[iBody + 1], + options[OPT_MASS].iLine[iBody + 1]); + } +} + +void VerifyNoRotEngleEarlyOptions(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, char cModel[LINE], + int iBody) { + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_MAGBRAKINGMODEL, OPT_ROTENGLEEARLYA, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_MAGBRAKINGMODEL, OPT_ROTENGLEEARLYB, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_MAGBRAKINGMODEL, OPT_ROTENGLEEARLYC, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_MAGBRAKINGMODEL, OPT_ROTENGLEEARLYD, iBody); +} + +void VerifyNoRotEngleMidOptions(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, char cModel[LINE], + int iBody) { + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_MAGBRAKINGMODEL, OPT_ROTENGLEMIDA, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_MAGBRAKINGMODEL, OPT_ROTENGLEMIDB, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_MAGBRAKINGMODEL, OPT_ROTENGLEMIDC, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_MAGBRAKINGMODEL, OPT_ROTENGLEMIDD, iBody); +} + +void VerifyNoRotEngleLateOptions(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, char cModel[LINE], + int iBody) { + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_MAGBRAKINGMODEL, OPT_ROTENGLELATEA, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_MAGBRAKINGMODEL, OPT_ROTENGLELATEB, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_MAGBRAKINGMODEL, OPT_ROTENGLELATEC, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_MAGBRAKINGMODEL, OPT_ROTENGLELATED, iBody); +} + +void VerifyRotRateEngleEarly(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, int iBody) { + char cModel[LINE] = "ENGLE23EARLY magnetic braking"; + double dMinMass = 0.4; + double dMaxMass = 0.6; + + VerifyEngleMassSpectralType(body, control, files, options, cModel, dMaxMass, + dMinMass, OPT_MAGBRAKINGMODEL, iBody); + VerifyNoRotEngleMidOptions(body, control, files, options, cModel, iBody); + VerifyNoRotEngleLateOptions(body, control, files, options, cModel, iBody); +} + +void VerifyRotRateEngleMid(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, int iBody) { + char cModel[LINE] = "ENGLE23MID magnetic braking"; + double dMinMass = 0.2; + double dMaxMass = 0.4; + + VerifyEngleMassSpectralType(body, control, files, options, cModel, dMaxMass, + dMinMass, OPT_MAGBRAKINGMODEL, iBody); + VerifyNoRotEngleEarlyOptions(body, control, files, options, cModel, iBody); + VerifyNoRotEngleLateOptions(body, control, files, options, cModel, iBody); +} + +void VerifyRotRateEngleLate(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, int iBody) { + char cModel[LINE] = "ENGLE23LATE magnetic braking"; + double dMinMass = 0.1; + double dMaxMass = 0.2; + + VerifyEngleMassSpectralType(body, control, files, options, cModel, dMaxMass, + dMinMass, OPT_MAGBRAKINGMODEL, iBody); + VerifyNoRotEngleEarlyOptions(body, control, files, options, cModel, iBody); + VerifyNoRotEngleMidOptions(body, control, files, options, cModel, iBody); +} + +void VerifyRotRateEngle(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, int iBody) { + if (body[iBody].iMagBrakingModel == STELLAR_DJDT_ENGLE23EARLY) { + VerifyRotRateEngleEarly(body, control, files, options, iBody); + } else if (body[iBody].iMagBrakingModel == STELLAR_DJDT_ENGLE23MID) { + VerifyRotRateEngleMid(body, control, files, options, iBody); + } else if (body[iBody].iMagBrakingModel == STELLAR_DJDT_ENGLE23LATE) { + VerifyRotRateEngleLate(body, control, files, options, iBody); + } +} + +void VerifyXUVEngleEarly(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, int iBody) { + char cModel[LINE] = "ENGLE24EARLY XUV"; + double dMinMass = 0.4; + double dMaxMass = 0.6; + + VerifyEngleMassSpectralType(body, control, files, options, cModel, dMaxMass, + dMinMass, OPT_XUVMODEL, iBody); + + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_XUVMODEL, OPT_XUVENGLEMIDLATEA, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_XUVMODEL, OPT_XUVENGLEMIDLATEB, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_XUVMODEL, OPT_XUVENGLEMIDLATEC, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_XUVMODEL, OPT_XUVENGLEMIDLATED, iBody); +} + +void VerifyXUVEngleMidLate(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, int iBody) { + char cModel[LINE] = "ENGLE24MIDLATE XUV"; + double dMinMass = 0.1; + double dMaxMass = 0.4; + + VerifyEngleMassSpectralType(body, control, files, options, cModel, dMaxMass, + dMinMass, OPT_XUVMODEL, iBody); + + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_XUVMODEL, OPT_XUVENGLEEARLYA, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_XUVMODEL, OPT_XUVENGLEEARLYB, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_XUVMODEL, OPT_XUVENGLEEARLYC, iBody); + VerifyTwoEngleOptionsNotSet(body, control, files, options, cModel, + OPT_XUVMODEL, OPT_XUVENGLEEARLYD, iBody); +} + +void VerifyXUVEngle(BODY *body, CONTROL *control, FILES *files, + OPTIONS *options, int iBody) { + if (body[iBody].iXUVModel == STELLAR_MODEL_ENGLE24EARLY) { + VerifyXUVEngleEarly(body, control, files, options, iBody); + } else if (body[iBody].iXUVModel == STELLAR_MODEL_ENGLE24MIDLATE) { + VerifyXUVEngleMidLate(body, control, files, options, iBody); + } +} + void VerifyStellar(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, OUTPUT *output, SYSTEM *system, UPDATE *update, int iBody, int iModule) { /* Stellar is active for this body if this subroutine is called. */ + VerifyXUVEngle(body, control, files, options, iBody); + VerifyRotRateEngle(body, control, files, options, iBody); + + if (update[iBody].iNumLuminosity > 1) { if (control->Io.iVerbose >= VERBERR) { fprintf(stderr, @@ -2114,3 +2959,61 @@ double fdLuminosityFunctionSineWave(BODY *body, int iBody) { body[iBody].dLuminosityPhase); return dLuminosity; } + +double fdRotationPeriodEngleGeneral(BODY *body, double dA, double dB, double dC, + double dD, int iBody) { + double dRotationPeriod; + + double dMinRotPerInDays = 0.5; + double dRotPerInDays = body[iBody].dRotPer / DAYSEC; + double dLogAgeInGyr = log10(body[iBody].dAge / (1e9 * YEARSEC)); + /* D is negative; -D is the rotation-period threshold (days) that selects + the two-branch fit in Engle & Guinan (2023). */ + if (dRotPerInDays < -dD) { + dRotationPeriod = (dLogAgeInGyr - dB) / dA; + } else { + dRotationPeriod = (dLogAgeInGyr - (dB + dC * dD)) / (dA + dC); + } + + if (dRotationPeriod < dMinRotPerInDays) { + dRotationPeriod = dMinRotPerInDays; + } + + return dRotationPeriod * DAYSEC; +} + +double fdRotationPeriodEngle(BODY *body, int iBody) { + double dRotationPeriod = body[iBody].dRotPer; + + if (body[iBody].iMagBrakingModel == STELLAR_DJDT_ENGLE23EARLY) { + dRotationPeriod = fdRotationPeriodEngleGeneral( + body, body[iBody].dRotEngleEarlyA, body[iBody].dRotEngleEarlyB, + body[iBody].dRotEngleEarlyC, body[iBody].dRotEngleEarlyD, iBody); + } else if (body[iBody].iMagBrakingModel == STELLAR_DJDT_ENGLE23MID) { + dRotationPeriod = fdRotationPeriodEngleGeneral( + body, body[iBody].dRotEngleMidA, body[iBody].dRotEngleMidB, + body[iBody].dRotEngleMidC, body[iBody].dRotEngleMidD, iBody); + } else if (body[iBody].iMagBrakingModel == STELLAR_DJDT_ENGLE23LATE) { + dRotationPeriod = fdRotationPeriodEngleGeneral( + body, body[iBody].dRotEngleLateA, body[iBody].dRotEngleLateB, + body[iBody].dRotEngleLateC, body[iBody].dRotEngleLateD, iBody); + } + + return dRotationPeriod; +} + +double fdLXUVEngle(BODY *body, double dA, double dB, double dC, double dD, + int iBody) { + double dLXUV; + + double dAgeInGyr = body[iBody].dAge / (1.e9 * YEARSEC); + double dLogAge = log10(dAgeInGyr); + if (dLogAge >= dD) { + dLXUV = pow(10., (dA * dLogAge + dB + dC * (dLogAge - dD))); + } else { + dLXUV = pow(10., (dA * dLogAge + dB)); + } + dLXUV *= body[iBody].dLuminosity; + + return dLXUV; +} diff --git a/src/stellar.h b/src/stellar.h index e286e5c90..5e738a407 100644 --- a/src/stellar.h +++ b/src/stellar.h @@ -41,6 +41,10 @@ #define STELLAR_MODEL_RIBAS 4 #define STELLAR_MODEL_PROXIMACEN 5 #define STELLAR_MODEL_SINEWAVE 6 +#define STELLAR_MODEL_ENGLE24EARLY \ + 7 /**< XUV evolution from Engle (2024) for M0-M2 dwarfs */ +#define STELLAR_MODEL_ENGLE24MIDLATE \ + 8 /**< XUV evolution from Engle (2024) for M2.6-M6.5 dwarfs */ #define STELLAR_DJDT_NONE \ 0 /**< No stellar angular momentum loss via magnetic braking */ @@ -49,6 +53,12 @@ 2 /**< dJ/dt according to Skumanich 1972 empirical law */ #define STELLAR_DJDT_MA15 3 /**< dJ/dt according to Matt+2015 */ #define STELLAR_DJDT_BR21 4 /**< dJ/dt according to Breimann+2021 */ +#define STELLAR_DJDT_ENGLE23EARLY \ + 5 /**< Rotation evolution from Engle & Guinan (2023) for M0-M2 dwarfs */ +#define STELLAR_DJDT_ENGLE23MID \ + 6 /**< Rotation evolution from Engle & Guinan (2023) for M2.5-M3.5 dwarfs */ +#define STELLAR_DJDT_ENGLE23LATE \ + 7 /**< Rotation evolution from Engle & Guinan (2023) for M4-M6.5 dwarfs */ #define HZ_MODEL_KOPPARAPU 1 #define DRYRGFLUX 415 /**< W/m^2 from Abe et al. (2011) */ @@ -77,6 +87,31 @@ #define OPT_LUMPERIOD 1555 #define OPT_LUMPHASE 1560 +#define OPT_XUVENGLEEARLYA 1570 +#define OPT_XUVENGLEEARLYB 1571 +#define OPT_XUVENGLEEARLYC 1572 +#define OPT_XUVENGLEEARLYD 1573 + +#define OPT_XUVENGLEMIDLATEA 1574 +#define OPT_XUVENGLEMIDLATEB 1575 +#define OPT_XUVENGLEMIDLATEC 1576 +#define OPT_XUVENGLEMIDLATED 1577 + +#define OPT_ROTENGLEEARLYA 1580 +#define OPT_ROTENGLEEARLYB 1581 +#define OPT_ROTENGLEEARLYC 1582 +#define OPT_ROTENGLEEARLYD 1583 + +#define OPT_ROTENGLEMIDA 1584 +#define OPT_ROTENGLEMIDB 1585 +#define OPT_ROTENGLEMIDC 1586 +#define OPT_ROTENGLEMIDD 1587 + +#define OPT_ROTENGLELATEA 1588 +#define OPT_ROTENGLELATEB 1589 +#define OPT_ROTENGLELATEC 1590 +#define OPT_ROTENGLELATED 1591 + /* Halt Functions */ #define STELLARHALTSYSEND 5 #define STELLARHALTBODYEND 5 @@ -187,6 +222,8 @@ double fdDEDtRotRadGyraStellar(BODY *, SYSTEM *, int *); double fdDEDtRotBrakeStellar(BODY *, SYSTEM *, int *); double fdDEDtStellar(BODY *, SYSTEM *, int *); double fdCranmerSaar2011TauCZ(double); +double fdRotationPeriodEngle(BODY *, int); +double fdLXUVEngle(BODY *, double, double, double, double, int); /* Dummy functions */ double fdSurfEnFluxStellar(BODY *, SYSTEM *, UPDATE *, int, int); diff --git a/src/vplanet.h b/src/vplanet.h index 11c4a80ba..c19688562 100644 --- a/src/vplanet.h +++ b/src/vplanet.h @@ -622,6 +622,30 @@ struct BODY { double dLuminosityFrequency; double dLuminosityPhase; + /* Engle (2024) XUV-evolution coefficients */ + double dXUVEngleEarlyA; + double dXUVEngleEarlyB; + double dXUVEngleEarlyC; + double dXUVEngleEarlyD; + double dXUVEngleMidLateA; + double dXUVEngleMidLateB; + double dXUVEngleMidLateC; + double dXUVEngleMidLateD; + + /* Engle & Guinan (2023) rotation-evolution coefficients */ + double dRotEngleEarlyA; + double dRotEngleEarlyB; + double dRotEngleEarlyC; + double dRotEngleEarlyD; + double dRotEngleMidA; + double dRotEngleMidB; + double dRotEngleMidC; + double dRotEngleMidD; + double dRotEngleLateA; + double dRotEngleLateB; + double dRotEngleLateC; + double dRotEngleLateD; + /* POISE parameters */ int bPoise; /**< Apply POISE module? */ @@ -850,6 +874,14 @@ struct BODY { // coefficient*/ // double dFlareSlopeErrorLower; /**< Lower error of slope /FFD angular // coefficient*/ + double dFlareSlopeA1; /**< a1 in slope = a1*log10(age_Myr) + a2*M_sun + a3 + (Davenport+2019 variable-slope FFD)*/ + double dFlareSlopeA2; /**< a2 coefficient; see dFlareSlopeA1 */ + double dFlareSlopeA3; /**< a3 coefficient; see dFlareSlopeA1 */ + double dFlareYIntB1; /**< b1 in yint = b1*log10(age_Myr) + b2*M_sun + b3 + (Davenport+2019 variable-slope FFD)*/ + double dFlareYIntB2; /**< b2 coefficient; see dFlareYIntB1 */ + double dFlareYIntB3; /**< b3 coefficient; see dFlareYIntB1 */ double dFlareMinEnergy; /**< Flare minimum energy value to calculate the FFD*/ double dFlareMaxEnergy; /**< Flare maximum energy value to calculate the FFD*/ double dFlareFreq1; /**< First value of flare frequency range*/ diff --git a/tests/AtmescFlareStellar/AtmEscFlareDavenportVariableSlope/b.in b/tests/AtmescFlareStellar/AtmEscFlareDavenportVariableSlope/b.in new file mode 100644 index 000000000..c2a2b7a00 --- /dev/null +++ b/tests/AtmescFlareStellar/AtmEscFlareDavenportVariableSlope/b.in @@ -0,0 +1,24 @@ +# Proxima b's parameters +sName b # Body's name +saModules atmesc + +# Physical Properties +dMass -1.27 # Mass, negative -> Earth masses + +# Orbital Properties +dEcc 0.08 # Eccentricity +dSemi 0.0485 # Semi-major axis, negative -> AU + +# ATMESC Properties +dXFrac 1.0 # X-Ray/XUV absorption radius (fraction of planet radius) +dSurfWaterMass -1.5 # Initial surface water (Earth oceans) +dEnvelopeMass -0.001 # 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) +dMinEnvelopeMass -1.e-5 +sPlanetRadiusModel proximacenb +bAtmEscAuto 1 # Is the flow radiation-recombination-limited? + +#Output +saOutputOrder Time -SurfWaterMass -EnvelopeMass -OxygenMass -PlanetRadius -FXUV Instellation diff --git a/tests/AtmescFlareStellar/AtmEscFlareDavenportVariableSlope/star.in b/tests/AtmescFlareStellar/AtmEscFlareDavenportVariableSlope/star.in new file mode 100644 index 000000000..0ac8a41c6 --- /dev/null +++ b/tests/AtmescFlareStellar/AtmEscFlareDavenportVariableSlope/star.in @@ -0,0 +1,27 @@ +sName star # Body's name +saModules stellar flare + +# Physical Parameters +dMass 0.12 +dAge 1e7 + +# STELLAR parameters +sStellarModel baraffe # Stellar evolutin model: `baraffe` or `none` +dSatXUVFrac 1.e-3 # Saturation level of the XUV luminosity +dRadGyra 0.5 + +dFlareMinEnergy -1.0e33 +dFlareMaxEnergy -1.0e36 +sFlareFFD DAVENPORT +sFlareBandPass KEPLER + +# Variable-slope Davenport (2019) coefficients: non-default values chosen +# from the GJ 1132 MCMC posterior mean so this test exercises the new path. +dFlareSlopeA1 -0.1487 +dFlareSlopeA2 0.5167 +dFlareSlopeA3 -0.6160 +dFlareYIntB1 4.7077 +dFlareYIntB2 -16.4409 +dFlareYIntB3 19.3826 + +saOutputOrder Time -LXUVFlare -LXUVStellar FlareSlopeOUT FlareYIntOUT diff --git a/tests/AtmescFlareStellar/AtmEscFlareDavenportVariableSlope/test_AtmEscFlareDavenportVariableSlope.py b/tests/AtmescFlareStellar/AtmEscFlareDavenportVariableSlope/test_AtmEscFlareDavenportVariableSlope.py new file mode 100644 index 000000000..a6f38e482 --- /dev/null +++ b/tests/AtmescFlareStellar/AtmEscFlareDavenportVariableSlope/test_AtmEscFlareDavenportVariableSlope.py @@ -0,0 +1,45 @@ +import astropy.units as u +from benchmark import Benchmark, benchmark + + +@benchmark( + { + "log.initial.system.Age": {"value": 3.1557600000e14, "unit": u.sec}, + "log.initial.system.Time": {"value": 0.0, "unit": u.sec}, + "log.initial.system.TotAngMom": { + "value": 3.7230704732e41, + "unit": (u.kg * u.m**2) / u.sec, + }, + "log.initial.system.TotEnergy": {"value": -7.1876127212e39, "unit": u.Joule}, + "log.initial.star.Mass": {"value": 2.3860992000e29, "unit": u.kg}, + "log.initial.star.Radius": {"value": 3.1661898155e8, "unit": u.m}, + "log.initial.star.RadGyra": {"value": 0.4610332834}, + "log.initial.star.RotPer": {"value": 8.6400000000e4, "unit": u.sec}, + "log.initial.star.Luminosity": { + "value": 6.5779241790e24, + "unit": (u.kg * u.m**2) / u.sec**3, + }, + "log.initial.star.LXUVStellar": {"value": 1.7103286997e-05, "unit": u.LSUN}, + "log.initial.star.LXUVFlare": {"value": 1.0944793096e-05, "unit": u.LSUN}, + "log.initial.star.Temperature": {"value": 3095.3300004448, "unit": u.K}, + "log.initial.star.FlareSlopeOUT": {"value": -0.7026960000}, + "log.initial.star.FlareYIntOUT": {"value": 22.1173920000}, + "log.final.system.Age": {"value": 3.4713360000e15, "unit": u.sec}, + "log.final.system.TotAngMom": { + "value": 3.7295966219e41, + "unit": (u.kg * u.m**2) / u.sec, + }, + "log.final.system.TotEnergy": {"value": -7.1816302434e39, "unit": u.Joule}, + "log.final.star.Luminosity": { + "value": 1.1827853786e24, + "unit": (u.kg * u.m**2) / u.sec**3, + }, + "log.final.star.LXUVStellar": {"value": 2.7351658447e-06, "unit": u.LSUN}, + "log.final.star.LXUVFlare": {"value": 3.4907273326e-06, "unit": u.LSUN}, + "log.final.star.Temperature": {"value": 3070.0997599141, "unit": u.K}, + "log.final.star.FlareSlopeOUT": {"value": -0.8575510923}, + "log.final.star.FlareYIntOUT": {"value": 27.0199563439}, + } +) +class Test_AtmEscFlareDavenportVariableSlope(Benchmark): + pass diff --git a/tests/AtmescFlareStellar/AtmEscFlareDavenportVariableSlope/vpl.in b/tests/AtmescFlareStellar/AtmEscFlareDavenportVariableSlope/vpl.in new file mode 100644 index 000000000..e2b98800f --- /dev/null +++ b/tests/AtmescFlareStellar/AtmEscFlareDavenportVariableSlope/vpl.in @@ -0,0 +1,39 @@ +# Example primary input file for VPLanet +sSystemName flare # System Name +iVerbose 5 # Verbosity level +bOverwrite 1 # Allow file overwrites? + +# All space after a # is ignored, as is white space +# The first lowercase letter(s) denote the cast: b=boolean, i=int, d=double, +# s=string. An "a" indicates an array and multiple arguments are allowed/expected. + +# List of "body files" that contain body-specific parameters +saBodyFiles star.in b.in + +# Array options can continue to the next line with a terminating "$". The $ can be +# at the end of the string or not. Comments are allowed afterwards. + +# 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 +sUnitTemp K + +# Units specified in the primary input file are propagated into the bodies. Otherwise +# specify units on a per body basis in the body files. +# Most string arguments can be in any case and need only be unambiguous. + +# Input/Output +bDoLog 1 # Write a log file? +iDigits 10 # Maximum number of digits to right of decimal +dMinValue 1e-10 # Minimum value of eccentricity/obliquity + +# Option names must be exact in spelling and case. + +# Evolution Parameters +bDoForward 1 # Perform a forward evolution? +bVarDt 1 # Use variable timestepping? +dEta 0.01 # Coefficient for variable timestepping +dStopTime 1e8 +dOutputTime 1e8 diff --git a/tests/Stellar/EngleRotation/star.in b/tests/Stellar/EngleRotation/star.in new file mode 100644 index 000000000..9904f0fde --- /dev/null +++ b/tests/Stellar/EngleRotation/star.in @@ -0,0 +1,20 @@ +# Engle23 Mid magnetic-braking model test — M3-like dwarf + +sName star +saModules stellar + +dMass 0.3 +dRotPeriod -10 +dAge 1e9 + +sStellarModel baraffe +sXUVModel RIBAS +sMagBrakingModel Engle23Mid + +# Engle & Guinan (2023) Mid coefficients — defaults, stated explicitly +dRotEngleMidA 0.0561 +dRotEngleMidB -0.8900 +dRotEngleMidC -0.0521 +dRotEngleMidD -24.1888 + +saOutputOrder Time -Luminosity -RotPer -LXUVStellar diff --git a/tests/Stellar/EngleRotation/test_EngleRotation.py b/tests/Stellar/EngleRotation/test_EngleRotation.py new file mode 100644 index 000000000..06fa8f69b --- /dev/null +++ b/tests/Stellar/EngleRotation/test_EngleRotation.py @@ -0,0 +1,31 @@ +import astropy.units as u +from benchmark import Benchmark, benchmark + + +@benchmark( + { + "log.initial.system.Age": {"value": 3.155760e16, "unit": u.sec}, + "log.initial.system.Time": {"value": 0.0, "unit": u.sec}, + "log.initial.system.TotAngMom": { + "value": 3.729543e40, + "unit": (u.kg * u.m**2) / u.sec, + }, + "log.initial.system.TotEnergy": {"value": -7.038771e40, "unit": u.Joule}, + "log.initial.star.Mass": {"value": 5.965248e29, "unit": u.kg}, + "log.initial.star.Radius": {"value": 2.024487e8, "unit": u.m}, + "log.initial.star.RadGyra": {"value": 0.458000}, + "log.initial.star.RotPer": {"value": 10.0, "unit": u.day}, + "log.initial.star.RotRate": {"value": 7.272205e-6, "unit": 1 / u.sec}, + "log.initial.star.Luminosity": {"value": 0.010328, "unit": u.LSUN}, + "log.initial.star.LXUVStellar": {"value": 6.081350e-7, "unit": u.LSUN}, + "log.initial.star.Temperature": {"value": 3415.002296, "unit": u.K}, + "log.final.system.Age": {"value": 3.155760e17, "unit": u.sec}, + "log.final.system.TotEnergy": {"value": -6.865097e40, "unit": u.Joule}, + "log.final.star.RotPer": {"value": 157.440880, "unit": u.day}, + "log.final.star.Luminosity": {"value": 0.011350, "unit": u.LSUN}, + "log.final.star.LXUVStellar": {"value": 3.935500e-8, "unit": u.LSUN}, + "log.final.star.Temperature": {"value": 3416.001085, "unit": u.K}, + } +) +class Test_EngleRotation(Benchmark): + pass diff --git a/tests/Stellar/EngleRotation/vpl.in b/tests/Stellar/EngleRotation/vpl.in new file mode 100644 index 000000000..e54a57b00 --- /dev/null +++ b/tests/Stellar/EngleRotation/vpl.in @@ -0,0 +1,20 @@ +sSystemName star +iVerbose 0 +bOverwrite 1 + +saBodyFiles star.in + +sUnitMass solar +sUnitLength AU +sUnitTime YEARS +sUnitAngle d + +bDoLog 1 +iDigits 6 +dMinValue 1e-10 + +bDoForward 1 +bVarDt 1 +dEta 0.01 +dStopTime 1e10 +dOutputTime 1e9 diff --git a/tests/Stellar/EngleXUV/star.in b/tests/Stellar/EngleXUV/star.in new file mode 100644 index 000000000..f7cb598fd --- /dev/null +++ b/tests/Stellar/EngleXUV/star.in @@ -0,0 +1,20 @@ +# Engle24 MidLate XUV model test — GJ 1132-like late M dwarf + +sName star +saModules stellar + +dMass 0.181 +dRotPeriod -122 +dAge 1e9 + +sStellarModel baraffe +sXUVModel Engle24MidLate +sMagBrakingModel reiners + +# Engle et al. (2024) MidLate coefficients — defaults, stated explicitly +dXUVEngleMidLateA -0.1456 +dXUVEngleMidLateB -2.8876 +dXUVEngleMidLateC -1.8187 +dXUVEngleMidLateD 0.3545 + +saOutputOrder Time -Luminosity -LXUVStellar -LXUVTot diff --git a/tests/Stellar/EngleXUV/test_EngleXUV.py b/tests/Stellar/EngleXUV/test_EngleXUV.py new file mode 100644 index 000000000..805cc7424 --- /dev/null +++ b/tests/Stellar/EngleXUV/test_EngleXUV.py @@ -0,0 +1,37 @@ +import astropy.units as u +from benchmark import Benchmark, benchmark + + +@benchmark( + { + "log.initial.system.Age": {"value": 3.155760e16, "unit": u.sec}, + "log.initial.system.Time": {"value": 0.0, "unit": u.sec}, + "log.initial.system.TotAngMom": { + "value": 8.694238e38, + "unit": (u.kg * u.m**2) / u.sec, + }, + "log.initial.system.TotEnergy": {"value": -3.754922e40, "unit": u.Joule}, + "log.initial.star.Mass": {"value": 3.599033e29, "unit": u.kg}, + "log.initial.star.Radius": {"value": 1.381422e8, "unit": u.m}, + "log.initial.star.RadGyra": {"value": 0.460833}, + "log.initial.star.RotPer": {"value": 1.054080e7, "unit": u.sec}, + "log.initial.star.RotRate": {"value": 5.960824e-7, "unit": 1 / u.sec}, + "log.initial.star.Luminosity": {"value": 0.003787, "unit": u.LSUN}, + "log.initial.star.LXUVStellar": {"value": 4.906221e-6, "unit": u.LSUN}, + "log.initial.star.LXUVTot": {"value": 4.906221e-6, "unit": u.LSUN}, + "log.initial.star.Temperature": {"value": 3212.579977, "unit": u.K}, + "log.final.system.Age": {"value": 3.108426e17, "unit": u.sec}, + "log.final.system.TotAngMom": { + "value": 8.931191e38, + "unit": (u.kg * u.m**2) / u.sec, + }, + "log.final.system.TotEnergy": {"value": -3.698928e40, "unit": u.Joule}, + "log.final.star.Luminosity": {"value": 0.003956, "unit": u.LSUN}, + "log.final.star.LXUVStellar": {"value": 2.529287e-7, "unit": u.LSUN}, + "log.final.star.LXUVTot": {"value": 2.529287e-7, "unit": u.LSUN}, + "log.final.star.Temperature": {"value": 3219.534123, "unit": u.K}, + "log.final.star.RotPer": {"value": 1.193750e7, "unit": u.sec}, + } +) +class Test_EngleXUV(Benchmark): + pass diff --git a/tests/Stellar/EngleXUV/vpl.in b/tests/Stellar/EngleXUV/vpl.in new file mode 100644 index 000000000..e54a57b00 --- /dev/null +++ b/tests/Stellar/EngleXUV/vpl.in @@ -0,0 +1,20 @@ +sSystemName star +iVerbose 0 +bOverwrite 1 + +saBodyFiles star.in + +sUnitMass solar +sUnitLength AU +sUnitTime YEARS +sUnitAngle d + +bDoLog 1 +iDigits 6 +dMinValue 1e-10 + +bDoForward 1 +bVarDt 1 +dEta 0.01 +dStopTime 1e10 +dOutputTime 1e9