Skip to content

Commit ba854fd

Browse files
RoryBarnesclaude
RoryBarnes
and
claude
committed
Add Engle stellar rotation/XUV models and Davenport variable-slope flare option
Extends the STELLAR module with Engle et al. rotation and XUV evolution tracks, and adds a variable-slope option to the Davenport flare model in the FLARE module. New tests cover EngleRotation, EngleXUV, and AtmEscFlareDavenportVariableSlope. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
1 parent 9086d64 commit ba854fd

15 files changed

Lines changed: 1531 additions & 13 deletions

File tree

‎src/flare.c‎

Lines changed: 255 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,13 @@ void BodyCopyFlare(BODY *dest, BODY *src, int foo, int iNumBodies, int iBody) {
2626
dest[iBody].iEnergyBin = src[iBody].iEnergyBin;
2727
// dest[iBody].dFlareSlope = src[iBody].dFlareSlope;
2828
// dest[iBody].dFlareYInt = src[iBody].dFlareYInt;
29+
30+
dest[iBody].dFlareSlopeA1 = src[iBody].dFlareSlopeA1;
31+
dest[iBody].dFlareSlopeA2 = src[iBody].dFlareSlopeA2;
32+
dest[iBody].dFlareSlopeA3 = src[iBody].dFlareSlopeA3;
33+
dest[iBody].dFlareYIntB1 = src[iBody].dFlareYIntB1;
34+
dest[iBody].dFlareYIntB2 = src[iBody].dFlareYIntB2;
35+
dest[iBody].dFlareYIntB3 = src[iBody].dFlareYIntB3;
2936
}
3037

3138
/**************** FLARE options ********************/
@@ -125,7 +132,7 @@ void ReadFlareBandPass(BODY *body,
125132
body[iFile - 1].iFlareBandPass = FLARE_UV;
126133
} else if (!memcmp(sLower(cTmp), "go", 2)) {
127134
body[iFile - 1].iFlareBandPass = FLARE_GOES;
128-
} else if (!memcmp(sLower(cTmp), "sr", 2)) {
135+
} else if (!memcmp(sLower(cTmp), "sx", 2)) {
129136
body[iFile - 1].iFlareBandPass = FLARE_SXR;
130137
} else if (!memcmp(sLower(cTmp), "te", 2)) {
131138
body[iFile - 1].iFlareBandPass = FLARE_TESS_UV;
@@ -318,6 +325,109 @@ void ReadFlareSlope(BODY *body, CONTROL *control, FILES *files,
318325
}
319326
}
320327

328+
void ReadFlareSlopeA1(BODY *body, CONTROL *control, FILES *files,
329+
OPTIONS *options, SYSTEM *system, int iFile) {
330+
int lTmp = -1;
331+
double dTmp;
332+
333+
AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp,
334+
control->Io.iVerbose);
335+
if (lTmp >= 0) {
336+
NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp,
337+
control->Io.iVerbose);
338+
body[iFile - 1].dFlareSlopeA1 = dTmp;
339+
UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile);
340+
} else if (iFile > 0) {
341+
body[iFile - 1].dFlareSlopeA1 = options->dDefault;
342+
}
343+
}
344+
345+
void ReadFlareSlopeA2(BODY *body, CONTROL *control, FILES *files,
346+
OPTIONS *options, SYSTEM *system, int iFile) {
347+
int lTmp = -1;
348+
double dTmp;
349+
350+
AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp,
351+
control->Io.iVerbose);
352+
if (lTmp >= 0) {
353+
NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp,
354+
control->Io.iVerbose);
355+
body[iFile - 1].dFlareSlopeA2 = dTmp;
356+
UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile);
357+
} else if (iFile > 0) {
358+
body[iFile - 1].dFlareSlopeA2 = options->dDefault;
359+
}
360+
}
361+
362+
void ReadFlareSlopeA3(BODY *body, CONTROL *control, FILES *files,
363+
OPTIONS *options, SYSTEM *system, int iFile) {
364+
int lTmp = -1;
365+
double dTmp;
366+
367+
AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp,
368+
control->Io.iVerbose);
369+
if (lTmp >= 0) {
370+
NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp,
371+
control->Io.iVerbose);
372+
body[iFile - 1].dFlareSlopeA3 = dTmp;
373+
UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile);
374+
} else if (iFile > 0) {
375+
body[iFile - 1].dFlareSlopeA3 = options->dDefault;
376+
}
377+
}
378+
379+
void ReadFlareYIntB1(BODY *body, CONTROL *control, FILES *files,
380+
OPTIONS *options, SYSTEM *system, int iFile) {
381+
int lTmp = -1;
382+
double dTmp;
383+
384+
AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp,
385+
control->Io.iVerbose);
386+
if (lTmp >= 0) {
387+
NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp,
388+
control->Io.iVerbose);
389+
body[iFile - 1].dFlareYIntB1 = dTmp;
390+
UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile);
391+
} else if (iFile > 0) {
392+
body[iFile - 1].dFlareYIntB1 = options->dDefault;
393+
}
394+
}
395+
396+
void ReadFlareYIntB2(BODY *body, CONTROL *control, FILES *files,
397+
OPTIONS *options, SYSTEM *system, int iFile) {
398+
int lTmp = -1;
399+
double dTmp;
400+
401+
AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp,
402+
control->Io.iVerbose);
403+
if (lTmp >= 0) {
404+
NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp,
405+
control->Io.iVerbose);
406+
body[iFile - 1].dFlareYIntB2 = dTmp;
407+
UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile);
408+
} else if (iFile > 0) {
409+
body[iFile - 1].dFlareYIntB2 = options->dDefault;
410+
}
411+
}
412+
413+
void ReadFlareYIntB3(BODY *body, CONTROL *control, FILES *files,
414+
OPTIONS *options, SYSTEM *system, int iFile) {
415+
int lTmp = -1;
416+
double dTmp;
417+
418+
AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp,
419+
control->Io.iVerbose);
420+
if (lTmp >= 0) {
421+
NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp,
422+
control->Io.iVerbose);
423+
body[iFile - 1].dFlareYIntB3 = dTmp;
424+
UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile);
425+
} else if (iFile > 0) {
426+
body[iFile - 1].dFlareYIntB3 = options->dDefault;
427+
}
428+
}
429+
430+
321431

322432
// TODO: Include the error in the FFD slopes to calculate the upper and higher
323433
// limit of XUV luminosity by flares
@@ -618,6 +728,90 @@ void InitializeOptionsFlare(OPTIONS *options, fnReadOption fnRead[]) {
618728
"MIN for number of flares per minutes, \n"
619729
" HOUR for number of flares per hour, and DAY number of flares per "
620730
"for days.\n");*/
731+
732+
fvFormattedString(&options[OPT_FLARESLOPEA1].cName, "dFlareSlopeA1");
733+
fvFormattedString(&options[OPT_FLARESLOPEA1].cDescr,
734+
"a1 coefficient in Davenport (2019) variable-slope FFD: "
735+
"slope = a1*log10(age_Myr) + a2*M_sun + a3");
736+
fvFormattedString(&options[OPT_FLARESLOPEA1].cDefault,
737+
"-0.07 (Davenport+2019 M-dwarf fit)");
738+
options[OPT_FLARESLOPEA1].dDefault = -0.07;
739+
options[OPT_FLARESLOPEA1].iType = 2;
740+
options[OPT_FLARESLOPEA1].bMultiFile = 1;
741+
options[OPT_FLARESLOPEA1].iModuleBit = FLARE;
742+
options[OPT_FLARESLOPEA1].bNeg = 0;
743+
fvFormattedString(&options[OPT_FLARESLOPEA1].cDimension, "nd");
744+
fnRead[OPT_FLARESLOPEA1] = &ReadFlareSlopeA1;
745+
746+
fvFormattedString(&options[OPT_FLARESLOPEA2].cName, "dFlareSlopeA2");
747+
fvFormattedString(&options[OPT_FLARESLOPEA2].cDescr,
748+
"a2 coefficient in Davenport (2019) variable-slope FFD; "
749+
"see dFlareSlopeA1");
750+
fvFormattedString(&options[OPT_FLARESLOPEA2].cDefault,
751+
"0.79 (Davenport+2019 M-dwarf fit)");
752+
options[OPT_FLARESLOPEA2].dDefault = 0.79;
753+
options[OPT_FLARESLOPEA2].iType = 2;
754+
options[OPT_FLARESLOPEA2].bMultiFile = 1;
755+
options[OPT_FLARESLOPEA2].iModuleBit = FLARE;
756+
options[OPT_FLARESLOPEA2].bNeg = 0;
757+
fvFormattedString(&options[OPT_FLARESLOPEA2].cDimension, "nd");
758+
fnRead[OPT_FLARESLOPEA2] = &ReadFlareSlopeA2;
759+
760+
fvFormattedString(&options[OPT_FLARESLOPEA3].cName, "dFlareSlopeA3");
761+
fvFormattedString(&options[OPT_FLARESLOPEA3].cDescr,
762+
"a3 coefficient in Davenport (2019) variable-slope FFD; "
763+
"see dFlareSlopeA1");
764+
fvFormattedString(&options[OPT_FLARESLOPEA3].cDefault,
765+
"-1.06 (Davenport+2019 M-dwarf fit)");
766+
options[OPT_FLARESLOPEA3].dDefault = -1.06;
767+
options[OPT_FLARESLOPEA3].iType = 2;
768+
options[OPT_FLARESLOPEA3].bMultiFile = 1;
769+
options[OPT_FLARESLOPEA3].iModuleBit = FLARE;
770+
options[OPT_FLARESLOPEA3].bNeg = 0;
771+
fvFormattedString(&options[OPT_FLARESLOPEA3].cDimension, "nd");
772+
fnRead[OPT_FLARESLOPEA3] = &ReadFlareSlopeA3;
773+
774+
fvFormattedString(&options[OPT_FLAREYINTB1].cName, "dFlareYIntB1");
775+
fvFormattedString(&options[OPT_FLAREYINTB1].cDescr,
776+
"b1 coefficient in Davenport (2019) variable-slope FFD: "
777+
"yint = b1*log10(age_Myr) + b2*M_sun + b3");
778+
fvFormattedString(&options[OPT_FLAREYINTB1].cDefault,
779+
"2.01 (Davenport+2019 M-dwarf fit)");
780+
options[OPT_FLAREYINTB1].dDefault = 2.01;
781+
options[OPT_FLAREYINTB1].iType = 2;
782+
options[OPT_FLAREYINTB1].bMultiFile = 1;
783+
options[OPT_FLAREYINTB1].iModuleBit = FLARE;
784+
options[OPT_FLAREYINTB1].bNeg = 0;
785+
fvFormattedString(&options[OPT_FLAREYINTB1].cDimension, "nd");
786+
fnRead[OPT_FLAREYINTB1] = &ReadFlareYIntB1;
787+
788+
fvFormattedString(&options[OPT_FLAREYINTB2].cName, "dFlareYIntB2");
789+
fvFormattedString(&options[OPT_FLAREYINTB2].cDescr,
790+
"b2 coefficient in Davenport (2019) variable-slope FFD; "
791+
"see dFlareYIntB1");
792+
fvFormattedString(&options[OPT_FLAREYINTB2].cDefault,
793+
"-25.15 (Davenport+2019 M-dwarf fit)");
794+
options[OPT_FLAREYINTB2].dDefault = -25.15;
795+
options[OPT_FLAREYINTB2].iType = 2;
796+
options[OPT_FLAREYINTB2].bMultiFile = 1;
797+
options[OPT_FLAREYINTB2].iModuleBit = FLARE;
798+
options[OPT_FLAREYINTB2].bNeg = 0;
799+
fvFormattedString(&options[OPT_FLAREYINTB2].cDimension, "nd");
800+
fnRead[OPT_FLAREYINTB2] = &ReadFlareYIntB2;
801+
802+
fvFormattedString(&options[OPT_FLAREYINTB3].cName, "dFlareYIntB3");
803+
fvFormattedString(&options[OPT_FLAREYINTB3].cDescr,
804+
"b3 coefficient in Davenport (2019) variable-slope FFD; "
805+
"see dFlareYIntB1");
806+
fvFormattedString(&options[OPT_FLAREYINTB3].cDefault,
807+
"33.99 (Davenport+2019 M-dwarf fit)");
808+
options[OPT_FLAREYINTB3].dDefault = 33.99;
809+
options[OPT_FLAREYINTB3].iType = 2;
810+
options[OPT_FLAREYINTB3].bMultiFile = 1;
811+
options[OPT_FLAREYINTB3].iModuleBit = FLARE;
812+
options[OPT_FLAREYINTB3].bNeg = 0;
813+
fvFormattedString(&options[OPT_FLAREYINTB3].cDimension, "nd");
814+
fnRead[OPT_FLAREYINTB3] = &ReadFlareYIntB3;
621815
}
622816

623817
void ReadOptionsFlare(BODY *body,
@@ -1008,6 +1202,20 @@ void WriteLXUVFlare(BODY *body,
10081202
}
10091203
}
10101204

1205+
void WriteFlareSlopeOUT(BODY *body, CONTROL *control, OUTPUT *output,
1206+
SYSTEM *system, UNITS *units, UPDATE *update,
1207+
int iBody, double *dTmp, char **cUnit) {
1208+
*dTmp = fdFlareSlopeOUT(body, iBody);
1209+
fvFormattedString(cUnit, "");
1210+
}
1211+
1212+
void WriteFlareYIntOUT(BODY *body, CONTROL *control, OUTPUT *output,
1213+
SYSTEM *system, UNITS *units, UPDATE *update,
1214+
int iBody, double *dTmp, char **cUnit) {
1215+
*dTmp = fdFlareYIntOUT(body, iBody);
1216+
fvFormattedString(cUnit, "");
1217+
}
1218+
10111219
// TODO: Include the error in the FFD slopes to calculate the upper and higher
10121220
// limit of XUV luminosity by flares
10131221
/*
@@ -1455,6 +1663,22 @@ void InitializeOutputFlare(OUTPUT *output, fnWriteOutput fnWrite[]) {
14551663
output[OUT_LXUVFLARE].iModuleBit = FLARE;
14561664
fnWrite[OUT_LXUVFLARE] = &WriteLXUVFlare;
14571665

1666+
fvFormattedString(&output[OUT_FLARESLOPEOUT].cName, "FlareSlopeOUT");
1667+
fvFormattedString(&output[OUT_FLARESLOPEOUT].cDescr,
1668+
"Instantaneous Davenport FFD slope at current age and mass");
1669+
output[OUT_FLARESLOPEOUT].bNeg = 0;
1670+
output[OUT_FLARESLOPEOUT].iNum = 1;
1671+
output[OUT_FLARESLOPEOUT].iModuleBit = FLARE;
1672+
fnWrite[OUT_FLARESLOPEOUT] = &WriteFlareSlopeOUT;
1673+
1674+
fvFormattedString(&output[OUT_FLAREYINTOUT].cName, "FlareYIntOUT");
1675+
fvFormattedString(&output[OUT_FLAREYINTOUT].cDescr,
1676+
"Instantaneous Davenport FFD y-intercept at current age and mass");
1677+
output[OUT_FLAREYINTOUT].bNeg = 0;
1678+
output[OUT_FLAREYINTOUT].iNum = 1;
1679+
output[OUT_FLAREYINTOUT].iModuleBit = FLARE;
1680+
fnWrite[OUT_FLAREYINTOUT] = &WriteFlareYIntOUT;
1681+
14581682
// TODO: Include the error in the FFD slopes to calculate the upper and higher
14591683
// limit of XUV luminosity by flares
14601684
/*fvFormattedString(output[OUT_LXUVFLAREUPPER].cName, "LXUVFlareUpper");
@@ -1642,6 +1866,25 @@ double fdDavenport(double dA1, double dA2, double dA3, double dStarAge,
16421866

16431867
return dA;
16441868
}
1869+
1870+
/* Instantaneous Davenport FFD slope for the current age and mass, evaluated
1871+
from the body's variable-slope coefficients. Computed on demand so that
1872+
output does not require storing redundant state. */
1873+
double fdFlareSlopeOUT(BODY *body, int iBody) {
1874+
return fdDavenport(body[iBody].dFlareSlopeA1,
1875+
body[iBody].dFlareSlopeA2,
1876+
body[iBody].dFlareSlopeA3,
1877+
body[iBody].dAge, body[iBody].dMass);
1878+
}
1879+
1880+
/* Instantaneous Davenport FFD y-intercept for the current age and mass. */
1881+
double fdFlareYIntOUT(BODY *body, int iBody) {
1882+
return fdDavenport(body[iBody].dFlareYIntB1,
1883+
body[iBody].dFlareYIntB2,
1884+
body[iBody].dFlareYIntB3,
1885+
body[iBody].dAge, body[iBody].dMass);
1886+
}
1887+
16451888
// fdFFD calculates FFD of the flares. If LACY mode is choosen, them fdFFD has
16461889
// to convert from SI to the units that the equation 3 from Davenport et al.
16471890
// 2019 understand (i.e., flares/day and flares/(day*log10(erg))).
@@ -1695,14 +1938,17 @@ double fdLXUVFlare(BODY *body, double dDeltaTime, int iBody) {
16951938
// slopes(constant)?##################################
16961939

16971940
if (body[iBody].iFlareFFD == FLARE_FFD_DAVENPORT) {
1698-
// The coefficient values given here were given by Dr. James Davenport in
1699-
// private comunication
1700-
dFlareSlope =
1701-
fdDavenport(-0.07, 0.79, -1.06, body[iBody].dAge,
1702-
body[iBody].dMass); //(-0.07054598,0.81225239,-1.07054511)
1703-
dFlareYInt = fdDavenport(
1704-
2.01, -25.15, 33.99, body[iBody].dAge,
1705-
body[iBody].dMass); //(2.06012734,-25.79885288,34.44115635)
1941+
/* Davenport+2019 variable-slope FFD: slope and y-intercept are each
1942+
third-order fits in log10(age_Myr) with a linear mass term. The six
1943+
coefficients default to Davenport et al.'s published M-dwarf values. */
1944+
dFlareSlope = fdDavenport(body[iBody].dFlareSlopeA1,
1945+
body[iBody].dFlareSlopeA2,
1946+
body[iBody].dFlareSlopeA3,
1947+
body[iBody].dAge, body[iBody].dMass);
1948+
dFlareYInt = fdDavenport(body[iBody].dFlareYIntB1,
1949+
body[iBody].dFlareYIntB2,
1950+
body[iBody].dFlareYIntB3,
1951+
body[iBody].dAge, body[iBody].dMass);
17061952
} else if (body[iBody].iFlareFFD == FLARE_FFD_LACY) {
17071953
dFlareSlope = body[iBody].dFlareSlope;
17081954
dFlareYInt = body[iBody].dFlareYInt;

‎src/flare.h‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,13 @@
4646
#define OPT_FLAREBANDPASS 2041
4747
#define OPT_LXUVFLARECONST 2042
4848

49+
#define OPT_FLARESLOPEA1 2043
50+
#define OPT_FLARESLOPEA2 2044
51+
#define OPT_FLARESLOPEA3 2045
52+
#define OPT_FLAREYINTB1 2046
53+
#define OPT_FLAREYINTB2 2047
54+
#define OPT_FLAREYINTB3 2048
55+
4956
/* Output Functinos */
5057

5158
/* FLARE 1900 - 1999 */
@@ -70,6 +77,10 @@
7077
#define OUT_FLAREENERGYMIN 2024
7178
#define OUT_FLAREENERGYMID 2025
7279
#define OUT_FLAREENERGYMAX 2026
80+
81+
#define OUT_FLAREYINTOUT 2046
82+
#define OUT_FLARESLOPEOUT 2047
83+
7384
/* @cond DOXYGEN_OVERRIDE */
7485

7586
void InitializeControlFlare(CONTROL *);
@@ -137,6 +148,10 @@ void WriteFlareEnergyMid(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *,
137148
UPDATE *, int, double *, char**);
138149
void WriteFlareEnergyMax(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *,
139150
UPDATE *, int, double *, char**);
151+
void WriteFlareSlopeOUT(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *,
152+
UPDATE *, int, double *, char**);
153+
void WriteFlareYIntOUT(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *,
154+
UPDATE *, int, double *, char**);
140155

141156
/* Logging Functions */
142157
void LogOptionsFlare(CONTROL *, FILE *);
@@ -148,6 +163,8 @@ void LogBodyFlare(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UPDATE *,
148163
/* FLARE functions */
149164
double fdLXUVFlare(BODY *, double, int);
150165
double fdDavenport(double, double, double, double, double);
166+
double fdFlareSlopeOUT(BODY *, int);
167+
double fdFlareYIntOUT(BODY *, int);
151168
double fdFFD(BODY *, int, double, double, double);
152169
double fdBandPassKepler(BODY *, int, double);
153170
double fdBandPassXUV(BODY *, int, double);

0 commit comments

Comments
 (0)