Skip to content

Adding the Engle 2023/2024 XUV and gyrochronology models #332

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
sengle01 wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Adding the Engle 2023/2024 XUV and gyrochronology models #332

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6254afb
Add files via upload
sengle01 May 29, 2025
3b47068
Add files via upload
sengle01 Jul 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
205 changes: 203 additions & 2 deletions src/stellar.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
\rst
This module implements the Baraffe stellar evolution tracks
:cite:`Baraffe15`, the Ribas XUV evolution model :cite:`Ribas05`,
the Engle XUV evolution models,
and a simple wind/magnetic braking model.
\endrst

Expand Down Expand Up @@ -129,6 +130,10 @@ void ReadStellarModel(BODY *body, CONTROL *control, FILES *files,
body[iFile - 1].iStellarModel = STELLAR_MODEL_PROXIMACEN;
} else if (!memcmp(sLower(cTmp), "si", 2)) {
body[iFile - 1].iStellarModel = STELLAR_MODEL_SINEWAVE;
} else if (!memcmp(sLower(cTmp), "engle23e", 8)) {
body[iFile - 1].iStellarModel = STELLAR_MODEL_ENGLE24EARLY;
} else if (!memcmp(sLower(cTmp), "engle23m", 8)) {
body[iFile - 1].iStellarModel = STELLAR_MODEL_ENGLE24MIDLATE;
} else {
if (control->Io.iVerbose >= VERBERR) {
fprintf(stderr,
Expand Down Expand Up @@ -165,11 +170,20 @@ 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;
printf("Magnetic Braking Model set to ENGLE23EARLY for body %d\n", iFile);
} else if (!memcmp(sLower(cTmp), "engle23m", 8)) {
body[iFile - 1].iMagBrakingModel = STELLAR_DJDT_ENGLE23MIDLATE;
printf("Magnetic Braking Model set to ENGLE23MIDLATE for body %d\n", iFile);
} else if (!memcmp(sLower(cTmp), "engle23l", 8)) {
body[iFile - 1].iMagBrakingModel = STELLAR_DJDT_ENGLE23LATE;
printf("Magnetic Braking Model set to ENGLE23LATE for body %d\n", iFile);
} 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/MIDLATE/LATE, or NONE.\n",
options->cName, cTmp);
}
LineExit(files->Infile[iFile].cIn, lTmp);
Expand Down Expand Up @@ -224,6 +238,10 @@ void ReadXUVModel(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
control->Io.iVerbose);
if (!memcmp(sLower(cTmp), "ri", 2)) {
body[iFile - 1].iXUVModel = STELLAR_MODEL_RIBAS;
} 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 (!memcmp(sLower(cTmp), "no", 2)) {
body[iFile - 1].iXUVModel = STELLAR_MODEL_NONE;
} else if (!memcmp(sLower(cTmp), "re", 2)) {
Expand Down Expand Up @@ -521,7 +539,8 @@ void InitializeOptionsStellar(OPTIONS *options, fnReadOption fnRead[]) {
fvFormattedString(&options[OPT_XUVMODEL].cName, "sXUVModel");
fvFormattedString(&options[OPT_XUVMODEL].cDescr, "XUV Evolution Model");
fvFormattedString(&options[OPT_XUVMODEL].cDefault, "RIBAS");
fvFormattedString(&options[OPT_XUVMODEL].cValues, "RIBAS REINERS NONE");
fvFormattedString(&options[OPT_XUVMODEL].cValues,
"RIBAS ENGLE24EARLY ENGLE24MIDLATE REINERS NONE");
options[OPT_XUVMODEL].iType = 3;
options[OPT_XUVMODEL].bMultiFile = 1;
options[OPT_XUVMODEL].iModuleBit = STELLAR;
Expand Down Expand Up @@ -1012,6 +1031,35 @@ 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) {

// Engle segmented power-law decay model
double dAge = body[iBody].dAge / (1.e9 * YEARSEC);
double logAge = log10(dAge);
double shift = -0.2985;
//double dTMin = body[iBody].dSatXUVTime / (1.e9 * YEARSEC);
if (logAge >= shift) {
body[iBody].dLXUV = pow(10., (-0.4896 * logAge - 3.2128 - 0.4469 *
(logAge - shift))) * body[iBody].dLuminosity;
} else {
body[iBody].dLXUV = pow(10., (-0.4896 * logAge - 3.2128)) *
body[iBody].dLuminosity;
}

} else if (body[iBody].iXUVModel == STELLAR_MODEL_ENGLE24MIDLATE) {

// Engle segmented power-law decay model
double dAge = body[iBody].dAge / (1.e9 * YEARSEC);
double logAge = log10(dAge);
double shift = 0.3545;
if (logAge >= shift) {
body[iBody].dLXUV = pow(10., (-0.1456 * logAge - 2.8876 - 1.8187 *
(logAge - shift))) * body[iBody].dLuminosity;
} else {
body[iBody].dLXUV = pow(10., (-0.1456 * logAge - 2.8876)) *
body[iBody].dLuminosity;
}

} else {

// Constant XUV fraction
Expand Down Expand Up @@ -1356,6 +1404,78 @@ void WriteDRotPerDtStellar(BODY *body, CONTROL *control, OUTPUT *output,
}
}

void WriteEmpirRotPerStellar(BODY *body, CONTROL *control, OUTPUT *output,
SYSTEM *system, UNITS *units, UPDATE *update,
int iBody, double *dTmp, char **cUnit) {

if (body[iBody].iMagBrakingModel == STELLAR_DJDT_ENGLE23EARLY) {
// Engle 2023 Early M-dwarf braking law
double logAge = log10(body[iBody].dAge / (1e9 * YEARSEC)); // Age in Gyr
double dEmpirRotPer;

if (logAge < (0.0621 * 23.4933 - 1.0437)) {
dEmpirRotPer = (logAge + 1.0437) / 0.0621;
} else {
dEmpirRotPer = 23.4933; // Initial guess
double tolerance = 1e-6;
double dRotPerOld;

do {
dRotPerOld = dEmpirRotPer;
dEmpirRotPer = (logAge + 1.0437 + 0.0528 * (dRotPerOld - 23.4933)) / 0.0621;
} while (fabs(dEmpirRotPer - dRotPerOld) > tolerance);
}
*dTmp = (dEmpirRotPer <= 1.0 ? 1.0 : dEmpirRotPer);
fsUnitsTime(1, cUnit);

}

else if (body[iBody].iMagBrakingModel == STELLAR_DJDT_ENGLE23MIDLATE) {
// Engle 2023 Mid-Late M-dwarf braking law
double logAge = log10(body[iBody].dAge / (1e9 * YEARSEC)); // Age in Gyr
double dEmpirRotPer;

if (logAge < (0.0561 * 24.1888 - 0.8900)) {
dEmpirRotPer = (logAge + 0.8900) / 0.0561;
} else {
dEmpirRotPer = 24.1888; // Initial guess
double tolerance = 1e-6;
double dRotPerOld;

do {
dRotPerOld = dEmpirRotPer;
dEmpirRotPer = (logAge + 0.8900 + 0.0521 * (dRotPerOld - 24.1888)) / 0.0561;
} while (fabs(dEmpirRotPer - dRotPerOld) > tolerance);
}
*dTmp = (dEmpirRotPer <= 1.0 ? 1.0 : dEmpirRotPer);
fsUnitsTime(1, cUnit);

}

else if (body[iBody].iMagBrakingModel == STELLAR_DJDT_ENGLE23LATE) {
// Engle 2023 Late M-dwarf braking law
double logAge = log10(body[iBody].dAge / (1e9 * YEARSEC)); // Age in Gyr
double dEmpirRotPer;

if (logAge < (0.0251 * 25.45 - 0.1615)) {
dEmpirRotPer = (logAge + 0.1615) / 0.0251;
} else {
dEmpirRotPer = 25.45; // Initial guess
double tolerance = 1e-6;
double dRotPerOld;

do {
dRotPerOld = dEmpirRotPer;
dEmpirRotPer = (logAge + 0.1615 + 0.0212 * (dRotPerOld - 25.45)) / 0.0251;
} while (fabs(dEmpirRotPer - dRotPerOld) > tolerance);
}
*dTmp = (dEmpirRotPer <= 1.0 ? 1.0 : dEmpirRotPer);
fsUnitsTime(1, cUnit);

}

}

void InitializeOutputStellar(OUTPUT *output, fnWriteOutput fnWrite[]) {

fvFormattedString(&output[OUT_LUMINOSITY].cName, "Luminosity");
Expand Down Expand Up @@ -1415,6 +1535,15 @@ void InitializeOutputStellar(OUTPUT *output, fnWriteOutput fnWrite[]) {
output[OUT_DROTPERDTSTELLAR].iNum = 1;
output[OUT_DROTPERDTSTELLAR].iModuleBit = STELLAR;
fnWrite[OUT_DROTPERDTSTELLAR] = &WriteDRotPerDtStellar;

fvFormattedString(&output[OUT_EMPIRROTPERSTELLAR].cName, "EmpirRotPer");
fvFormattedString(&output[OUT_EMPIRROTPERSTELLAR].cDescr,
"Rotation period according to empirical Engle relationships");
fvFormattedString(&output[OUT_EMPIRROTPERSTELLAR].cNeg, "days");
output[OUT_EMPIRROTPERSTELLAR].bNeg = 0;
output[OUT_EMPIRROTPERSTELLAR].iNum = 1;
output[OUT_EMPIRROTPERSTELLAR].iModuleBit = STELLAR;
fnWrite[OUT_EMPIRROTPERSTELLAR] = &WriteEmpirRotPerStellar;
}

/************ STELLAR Logging Functions **************/
Expand Down Expand Up @@ -2114,3 +2243,75 @@ double fdLuminosityFunctionSineWave(BODY *body, int iBody) {
body[iBody].dLuminosityPhase);
return dLuminosity;
}
iError == STELLAR_ERR_OUTOFBOUNDS_LO) {
fprintf(stderr,
"ERROR: Temperature out of bounds (low) in fdBaraffe().\n");
} else if (iError == STELLAR_ERR_FILE) {
fprintf(stderr,
"ERROR: File access error in temperature routine fdBaraffe().\n");
} else if (iError == STELLAR_ERR_BADORDER) {
fprintf(stderr,
"ERROR: Bad temperature interpolation order in routine "
"fdBaraffe().\n");
} else {
fprintf(stderr, "ERROR: Undefined temperature error in fdBaraffe().\n");
}
exit(EXIT_INT);
}
}

double fdLuminosityFunctionProximaCen(double dAge, double dMass) {
int iError;
double L = fdProximaCenStellar(PROXIMACEN_L, dAge, dMass, &iError);
if (iError == PROXIMACEN_ERROR) {
return NAN;
} else {
return L;
}
}

double fdTemperatureFunctionProximaCen(double dAge, double dMass) {
int iError;
double L = fdProximaCenStellar(PROXIMACEN_T, dAge, dMass, &iError);
if (iError == PROXIMACEN_ERROR) {
return NAN;
} else {
return L;
}
}

double fdRadiusFunctionProximaCen(double dAge, double dMass) {
int iError;
double L = fdProximaCenStellar(PROXIMACEN_R, dAge, dMass, &iError);
if (iError == PROXIMACEN_ERROR) {
return NAN;
} else {
return L;
}
}

/* Compute the convective turnover timescale (in seconds) based on Equation 36
from Cranmer and Saar 2011. This equation is valid for 3300 <= Teff <= 7000
K, which pretty much brackets all stellar evolution models considered in
stellar.
*/
double fdCranmerSaar2011TauCZ(double Teff) {
double tau =
314.24 * exp(-(Teff / 1952.5) - pow((Teff / 6250.0), 18)) + 0.002;
return tau * 86400.0;
}

double fdSurfEnFluxStellar(BODY *body, SYSTEM *system, UPDATE *update,
int iBody, int iFoo) {
// This is silly, but necessary!
// RORY: I don't think so! -- This function should be set to ReturnOutputZero
return 0;
}

double fdLuminosityFunctionSineWave(BODY *body, int iBody) {
double dLuminosity =
body[iBody].dLuminosityInitial + body[iBody].dLuminosityAmplitude *
sin(body[iBody].dAge * body[iBody].dLuminosityFrequency +
body[iBody].dLuminosityPhase);
return dLuminosity;
}
13 changes: 13 additions & 0 deletions src/stellar.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,12 @@
#define STELLAR_MODEL_PROXIMACEN 5
#define STELLAR_MODEL_SINEWAVE 6

#define STELLAR_DJDT_ENGLE23EARLY 7
#define STELLAR_DJDT_ENGLE23MIDLATE 8
#define STELLAR_DJDT_ENGLE23LATE 9
#define STELLAR_MODEL_ENGLE24EARLY 10
#define STELLAR_MODEL_ENGLE24MIDLATE 11

#define STELLAR_DJDT_NONE \
0 /**< No stellar angular momentum loss via magnetic braking */
#define STELLAR_DJDT_RM12 1 /**< dJ/dt according to Reiners & Mohanty 2012 */
Expand Down Expand Up @@ -95,6 +101,7 @@
#define OUT_ROSSBYNUMBER 1514
#define OUT_DROTPERDTSTELLAR 1515
#define OUT_WINDTORQUE 1516
#define OUT_EMPIRROTPERSTELLAR 1517

/* @cond DOXYGEN_OVERRIDE */

Expand Down Expand Up @@ -191,4 +198,10 @@ double fdCranmerSaar2011TauCZ(double);
/* Dummy functions */
double fdSurfEnFluxStellar(BODY *, SYSTEM *, UPDATE *, int, int);

/* @endcond */
fdCranmerSaar2011TauCZ(double);

/* Dummy functions */
double fdSurfEnFluxStellar(BODY *, SYSTEM *, UPDATE *, int, int);

/* @endcond */