Skip to content

Commit d41457b

Browse files
brcekadambrcekadam
committed
Updated Shikauchi coefficients, fixed surface and core abundances for HeMS and HeHG
1 parent ecb9af3 commit d41457b

File tree

4 files changed

+30
-12
lines changed

4 files changed

+30
-12
lines changed

‎src/HeHG.h‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -76,6 +76,12 @@ class HeHG: virtual public BaseStar, public HeMS {
7676

7777
double CalculateHeCoreMassAtPhaseEnd() const { return CalculateHeCoreMassOnPhase(); } // Same as on phase
7878
double CalculateHeCoreMassOnPhase() const { return m_Mass; } // NO-OP
79+
80+
double CalculateHeliumAbundanceCoreAtPhaseEnd() const { return 0.0; }
81+
double CalculateHeliumAbundanceCoreOnPhase() const { return 0.0; }
82+
83+
double CalculateHydrogenAbundanceCoreAtPhaseEnd() const { return 0.0; }
84+
double CalculateHydrogenAbundanceCoreOnPhase() const { return 0.0; }
7985

8086
double CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) const;
8187
double CalculateLambdaNanjingEnhanced(const int p_MassIndex, const STELLAR_POPULATION p_StellarPop) const { return CalculateLambdaNanjingStarTrack(0.0, 0.0); } // 0.0 are dummy values that are not used

‎src/HeMS.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,10 @@ class HeMS: virtual public BaseStar, public TPAGB {
6161
protected:
6262

6363
void Initialise() {
64+
// initialise surface abundances
65+
m_HydrogenAbundanceSurface = 0.0;
66+
m_HeliumAbundanceSurface = 1.0 - m_Metallicity;
67+
6468
CalculateTimescales();
6569
// JR: Age for HeMS is partially calculated before switching -
6670
// can get here from various places in ResolveEnvelopeLoss(),

‎src/MainSequence.cpp‎

Lines changed: 15 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -291,8 +291,7 @@ double MainSequence::CalculateLuminosityOnPhase(const double p_Time, const doubl
291291
// If BRCEK core prescription is used, return luminosity from Shikauchi et al. (2024) during core hydrogen burning (valid for MZAMS >= 15 Msol) or
292292
// luminosity that smoothly connects MS and HG during MS hook (valid for MZAMS >= BRCEK_LOWER_MASS_LIMIT)
293293
if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0)) {
294-
double tMS = timescales(tMS);
295-
if (utils::Compare(p_Time, 0.99 * tMS) > 0) // star in MS hook?
294+
if (utils::Compare(p_Time, 0.99 * timescales(tMS)) > 0) // star in MS hook?
296295
return CalculateLuminosityTransitionToHG(p_Mass, p_Time, p_LZAMS);
297296
else {
298297
if (utils::Compare(m_MZAMS, 15.0) >= 0) // use Shikauchi luminosity if MZAMS >= 15 Msun
@@ -340,9 +339,16 @@ double MainSequence::CalculateLuminosityOnPhase(const double p_Time, const doubl
340339
*/
341340
double MainSequence::CalculateLuminosityShikauchi(const double p_CoreMass, const double p_HeliumAbundanceCore) const {
342341
DBL_VECTOR L_COEFFICIENTS = std::get<2>(SHIKAUCHI_COEFFICIENTS);
343-
double logMixingCoreMass = std::log10(p_CoreMass);
344342

345-
double logL = L_COEFFICIENTS[0] * logMixingCoreMass + L_COEFFICIENTS[1] * p_HeliumAbundanceCore + L_COEFFICIENTS[2] * logMixingCoreMass * p_HeliumAbundanceCore + L_COEFFICIENTS[3] * logMixingCoreMass * logMixingCoreMass + L_COEFFICIENTS[4] * p_HeliumAbundanceCore * p_HeliumAbundanceCore + L_COEFFICIENTS[5] * logMixingCoreMass * logMixingCoreMass * logMixingCoreMass + L_COEFFICIENTS[6] * p_HeliumAbundanceCore * p_HeliumAbundanceCore * p_HeliumAbundanceCore + L_COEFFICIENTS[7] * logMixingCoreMass * logMixingCoreMass * p_HeliumAbundanceCore + L_COEFFICIENTS[8] * logMixingCoreMass * p_HeliumAbundanceCore * p_HeliumAbundanceCore + L_COEFFICIENTS[9];
343+
// common factors
344+
double logMixingCoreMass = std::log10(p_CoreMass);
345+
double logMixingCoreMass_2 = logMixingCoreMass * logMixingCoreMass;
346+
double logMixingCoreMass_3 = logMixingCoreMass_2 * logMixingCoreMass;
347+
348+
double heliumAbundanceCore_2 = p_HeliumAbundanceCore * p_HeliumAbundanceCore;
349+
double heliumAbundanceCore_3 = heliumAbundanceCore_2 * p_HeliumAbundanceCore;
350+
351+
double logL = L_COEFFICIENTS[0] * logMixingCoreMass + L_COEFFICIENTS[1] * p_HeliumAbundanceCore + L_COEFFICIENTS[2] * logMixingCoreMass * p_HeliumAbundanceCore + L_COEFFICIENTS[3] * logMixingCoreMass_2 + L_COEFFICIENTS[4] * heliumAbundanceCore_2 + L_COEFFICIENTS[5] * logMixingCoreMass_3 + L_COEFFICIENTS[6] * heliumAbundanceCore_3 + L_COEFFICIENTS[7] * logMixingCoreMass_2 * p_HeliumAbundanceCore + L_COEFFICIENTS[8] * logMixingCoreMass * heliumAbundanceCore_2 + L_COEFFICIENTS[9] * logMixingCoreMass_3 * logMixingCoreMass + L_COEFFICIENTS[10] * heliumAbundanceCore_3 * p_HeliumAbundanceCore + L_COEFFICIENTS[11] * logMixingCoreMass * heliumAbundanceCore_3 + L_COEFFICIENTS[12] * logMixingCoreMass_2 * heliumAbundanceCore_2 + L_COEFFICIENTS[13] * logMixingCoreMass_3 * p_HeliumAbundanceCore + L_COEFFICIENTS[14];
346352

347353
return PPOW(10.0, logL);
348354
}
@@ -764,11 +770,11 @@ DBL_DBL MainSequence::CalculateMainSequenceCoreMassBrcek(const double p_Dt, cons
764770

765771
auto fmix = [&](double mass) { return FMIX_COEFFICIENTS[0] + FMIX_COEFFICIENTS[1] * std::exp(-mass / FMIX_COEFFICIENTS[2]); }; // Shikauchi et al. (2024), eq (A3)
766772
double alpha = PPOW(10.0, std::max(-2.0, ALPHA_COEFFICIENTS[1] * m_MainSequenceCoreMass + ALPHA_COEFFICIENTS[2])) + ALPHA_COEFFICIENTS[0]; // ibid, eq (A2)
767-
double g = -0.0044 * m_MZAMS + 0.27; // ibid, eq (A7)
773+
double g = SHIKAUCHI_DELTA_COEFFICIENTS[1] * m_MainSequenceCoreMass + SHIKAUCHI_DELTA_COEFFICIENTS[2]; // ibid, eq (A7)
768774

769775
double delta;
770776
if (p_MassLossRate <= 0.0)
771-
delta = std::min(PPOW(10.0, -(m_HeliumAbundanceCore - m_InitialHeliumAbundance) / (1.0 - m_InitialHeliumAbundance - m_Metallicity) + g), 1.0); // ibid, eq (A6)
777+
delta = std::min(PPOW(10.0, -SHIKAUCHI_DELTA_COEFFICIENTS[0] * (m_HeliumAbundanceCore - m_InitialHeliumAbundance) / (1.0 - m_InitialHeliumAbundance - m_Metallicity) + g), 1.0); // ibid, eq (A6)
772778
else
773779
delta = PPOW(2.0, -(m_HeliumAbundanceCore - m_InitialHeliumAbundance) / (1.0 - m_InitialHeliumAbundance - m_Metallicity)); // updated prescription for mass gain
774780

@@ -1324,7 +1330,7 @@ std::tuple <DBL_VECTOR, DBL_VECTOR, DBL_VECTOR> MainSequence::InterpolateShikauc
13241330

13251331
DBL_VECTOR alphaCoeff(3, 0.0);
13261332
DBL_VECTOR fmixCoeff(3, 0.0);
1327-
DBL_VECTOR lCoeff(10, 0.0);
1333+
DBL_VECTOR lCoeff(15, 0.0);
13281334

13291335
// Skip calculation if BRCEK core prescription is not used
13301336
if (OPTIONS->MainSequenceCoreMassPrescription() != CORE_MASS_PRESCRIPTION::BRCEK)
@@ -1357,7 +1363,7 @@ std::tuple <DBL_VECTOR, DBL_VECTOR, DBL_VECTOR> MainSequence::InterpolateShikauc
13571363
alphaCoeff[i] = (SHIKAUCHI_ALPHA_COEFFICIENTS[0][i] * middle_logZ + SHIKAUCHI_ALPHA_COEFFICIENTS[1][i] * logZ_low) / middle_low;
13581364
fmixCoeff[i] = (SHIKAUCHI_FMIX_COEFFICIENTS[0][i] * middle_logZ + SHIKAUCHI_FMIX_COEFFICIENTS[1][i] * logZ_low) / middle_low;
13591365
}
1360-
for (size_t i = 0; i < 10; i++)
1366+
for (size_t i = 0; i < 15; i++)
13611367
lCoeff[i] = (SHIKAUCHI_L_COEFFICIENTS[0][i] * middle_logZ + SHIKAUCHI_L_COEFFICIENTS[1][i] * logZ_low) / middle_low;
13621368
}
13631369
// Linear interpolation between metallicity middle and high
@@ -1366,7 +1372,7 @@ std::tuple <DBL_VECTOR, DBL_VECTOR, DBL_VECTOR> MainSequence::InterpolateShikauc
13661372
alphaCoeff[i] = (SHIKAUCHI_ALPHA_COEFFICIENTS[1][i] * high_logZ + SHIKAUCHI_ALPHA_COEFFICIENTS[2][i] * logZ_middle) / high_middle;
13671373
fmixCoeff[i] = (SHIKAUCHI_FMIX_COEFFICIENTS[1][i] * high_logZ + SHIKAUCHI_FMIX_COEFFICIENTS[2][i] * logZ_middle) / high_middle;
13681374
}
1369-
for (size_t i = 0; i < 10; i++)
1375+
for (size_t i = 0; i < 15; i++)
13701376
lCoeff[i] = (SHIKAUCHI_L_COEFFICIENTS[1][i] * high_logZ + SHIKAUCHI_L_COEFFICIENTS[2][i] * logZ_middle) / high_middle;
13711377
}
13721378
// Linear extrapolation (constant) for metallicity equal to solar or higher

‎src/constants.h‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -3755,6 +3755,8 @@ const std::vector<std::vector<std::vector<LoveridgeCoefficients>>> LOVERIDGE_COE
37553755

37563756
// Coefficients for determining Main Sequence core mass
37573757
// from Shikauchi et al. (2024), https://arxiv.org/abs/2409.00460
3758+
// Section A.4
3759+
const DBL_VECTOR SHIKAUCHI_DELTA_COEFFICIENTS = {0.54491412, -0.00900365, 0.08936248};
37583760
// Table 2
37593761
const std::vector<DBL_VECTOR> SHIKAUCHI_ALPHA_COEFFICIENTS = {
37603762
{0.45, -0.0557105, -0.86589929}, // 0.1*Z_Sun
@@ -3769,9 +3771,9 @@ const std::vector<DBL_VECTOR> SHIKAUCHI_FMIX_COEFFICIENTS = {
37693771
};
37703772
// Table 4
37713773
const std::vector<DBL_VECTOR> SHIKAUCHI_L_COEFFICIENTS = {
3772-
{3.2555795, 1.84666823, -0.79986388, -0.75728099, -0.38831172, 0.08223542, 0.49543834, 0.31314176, -0.36705796, 1.72200581}, // 0.1*Z_Sun
3773-
{3.35622529, 1.96904931, -0.88894808, -0.81112488, -0.47925922, 0.09056925, 0.53094768, 0.33971972, -0.35581284, 1.65390003}, // 1/3*Z_Sun
3774-
{3.27883249, 1.79370338, -0.71413866, -0.77019351, -0.3898752, 0.07499563, 0.5920458, 0.33846556, -0.49649838, 1.71263853} // Solar metallicity Z_Sun
3774+
{3.38627891, 1.13599187, -0.97389238, -0.87675442, 1.65386007, 0.07661174, -1.78737297, 0.622451, -0.47511355, 0.02483567, 0.94243277, -0.06798225, 0.11086108, -0.14859538, 1.78029915}, // 0.1*Z_Sun
3775+
{3.45464814, 0.94880846, -1.11409154, -0.86672079, 2.38986855, 0.04448855, -2.74913945, 0.60905625, -0.27648361, 0.03514139, 1.37569819, -0.19184532, 0.12816567, -0.14392935, 1.76390159}, // 1/3*Z_Sun
3776+
{3.80166901, 0.37407948, -1.29904749, -1.34541622, 3.70934166, 0.28320469, -3.92327169, 0.92444477, -0.40146717, -0.00821364, 1.80297947, -0.15776603, 0.09205681, -0.21913557, 1.78496679} // Solar metallicity Z_Sun
37753777
};
37763778

37773779
#endif // __constants_h__

0 commit comments

Comments
 (0)