Skip to content

Commit 30ff94c

Browse files
ilyamandelilyamandel
authored
Merge pull request #1406 from brcekadam/fix_zsol_reference
Fix for issue #400
2 parents 348395b + 09bef2d commit 30ff94c

File tree

7 files changed

+52
-43
lines changed

7 files changed

+52
-43
lines changed

‎src/BaseStar.cpp‎

Lines changed: 24 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -91,10 +91,10 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed,
9191

9292
// calculate coefficients, constants etc.
9393

94-
CalculateRCoefficients(LogMetallicityXi(), m_RCoefficients);
95-
CalculateLCoefficients(LogMetallicityXi(), m_LCoefficients);
94+
CalculateRCoefficients(LogMetallicityXiHurley(), m_RCoefficients);
95+
CalculateLCoefficients(LogMetallicityXiHurley(), m_LCoefficients);
9696

97-
CalculateMassCutoffs(m_Metallicity, LogMetallicityXi(), m_MassCutoffs);
97+
CalculateMassCutoffs(m_Metallicity, LogMetallicityXiHurley(), m_MassCutoffs);
9898

9999
CalculateAnCoefficients(m_AnCoefficients, m_LConstants, m_RConstants, m_GammaConstants);
100100
CalculateBnCoefficients(m_BnCoefficients);
@@ -506,7 +506,7 @@ void BaseStar::CalculateAnCoefficients(DBL_VECTOR &p_AnCoefficients,
506506
#define GammaConstants(x) p_GammaConstants[static_cast<int>(GAMMA_CONSTANTS::x)] // for convenience and readability - undefined at end of function
507507

508508
double Z = m_Metallicity;
509-
double xi = LogMetallicityXi();
509+
double xi = LogMetallicityXiHurley();
510510
double sigma = LogMetallicitySigma();
511511

512512
// pow() is slow - use multiplication
@@ -606,7 +606,7 @@ void BaseStar::CalculateBnCoefficients(DBL_VECTOR &p_BnCoefficients) {
606606

607607

608608
double Z = m_Metallicity;
609-
double xi = LogMetallicityXi();
609+
double xi = LogMetallicityXiHurley();
610610
double sigma = LogMetallicitySigma();
611611
double rho = LogMetallicityRho();
612612

@@ -853,8 +853,8 @@ void BaseStar::CalculateMassCutoffs(const double p_Metallicity, const double p_L
853853
massCutoffs(MHook) = 1.0185 + (0.16015 * p_LogMetallicityXi) + (0.0892 * xi_2); // MHook - Hurley et al. 2000, eq 1
854854
massCutoffs(MHeF) = 1.995 + (0.25 * p_LogMetallicityXi) + (0.087 * xi_2); // MHeF - Hurley et al. 2000, eq 2
855855

856-
double top = 13.048 * PPOW((p_Metallicity / ZSOL), 0.06);
857-
double bottom = 1.0 + (0.0012 * PPOW((ZSOL / p_Metallicity), 1.27));
856+
double top = 13.048 * PPOW((p_Metallicity / ZSOL_HURLEY), 0.06);
857+
double bottom = 1.0 + (0.0012 * PPOW((ZSOL_HURLEY / p_Metallicity), 1.27));
858858
massCutoffs(MFGB) = top / bottom; // MFGB - Hurley et al. 2000, eq 3
859859

860860
massCutoffs(MCHE) = 100.0; // MCHE - Mandel/Butler - CHE calculation
@@ -879,7 +879,7 @@ void BaseStar::CalculateMassCutoffs(const double p_Metallicity, const double p_L
879879
double BaseStar::CalculateGBRadiusXExponent() const {
880880

881881
// pow()is slow - use multiplication
882-
double xi = LogMetallicityXi();
882+
double xi = LogMetallicityXiHurley();
883883
double xi_2 = xi * xi;
884884
double xi_3 = xi_2 * xi;
885885
double xi_4 = xi_2 * xi_2;
@@ -1557,7 +1557,7 @@ double BaseStar::CalculateMassLossRateNieuwenhuijzenDeJager() const {
15571557

15581558
if (utils::Compare(m_Luminosity, NJ_MINIMUM_LUMINOSITY) > 0) { // check for minimum luminosity
15591559
double smoothTaper = min(1.0, (m_Luminosity - 4000.0) / 500.0); // smooth taper between no mass loss and mass loss
1560-
rate = std::sqrt((m_Metallicity / ZSOL)) * smoothTaper * 9.6E-15 * PPOW(m_Radius, 0.81) * PPOW(m_Luminosity, 1.24) * PPOW(m_Mass, 0.16);
1560+
rate = std::sqrt((m_Metallicity / ZSOL_HURLEY)) * smoothTaper * 9.6E-15 * PPOW(m_Radius, 0.81) * PPOW(m_Luminosity, 1.24) * PPOW(m_Mass, 0.16);
15611561
} else {
15621562
rate = 0.0;
15631563
}
@@ -1803,7 +1803,7 @@ double BaseStar::CalculateMassLossRateWolfRayetZDependent(const double p_Mu) con
18031803
// TW - Haven't seen StarTrack but I think H&K gives the original equation and V&dK gives the Z dependence
18041804
double rate = 0.0;
18051805
if (utils::Compare(p_Mu, 1.0) < 0) {
1806-
rate = 1.0E-13 * PPOW(m_Luminosity, 1.5) * PPOW(m_Metallicity / ZSOL, 0.86) * (1.0 - p_Mu);
1806+
rate = 1.0E-13 * PPOW(m_Luminosity, 1.5) * PPOW(m_Metallicity / ZSOL_ANDERS, 0.86) * (1.0 - p_Mu);
18071807
}
18081808
return rate;
18091809
}
@@ -1827,14 +1827,14 @@ double BaseStar::CalculateMassLossRateOBVink2001() const {
18271827
double teff = m_Temperature * TSOL;
18281828

18291829
if (utils::Compare(teff, VINK_MASS_LOSS_MINIMUM_TEMP) >= 0 && utils::Compare(teff, VINK_MASS_LOSS_BISTABILITY_TEMP) <= 0) {
1830-
double v = 1.3; // v_inf/v_esc
1831-
v = v * PPOW(m_Metallicity / ZSOL, OPTIONS->ScaleTerminalWindVelocityWithMetallicityPower()); // Scale Vinf with metallicity
1830+
double v = 1.3; // v_inf/v_esc
1831+
v = v * PPOW(m_Metallicity / ZSOL_ANDERS, OPTIONS->ScaleTerminalWindVelocityWithMetallicityPower()); // Scale Vinf with metallicity
18321832

18331833
double logMdotOB = -6.688 +
18341834
(2.210 * log10(m_Luminosity / 1.0E5)) -
18351835
(1.339 * log10(m_Mass / 30.0)) -
18361836
(1.601 * log10(v / 2.0)) +
1837-
(0.85 * LogMetallicityXi()) +
1837+
(0.85 * LogMetallicityXiAnders()) +
18381838
(1.07 * log10(teff / 20000.0));
18391839

18401840
rate = PPOW(10.0, logMdotOB);
@@ -1843,14 +1843,14 @@ double BaseStar::CalculateMassLossRateOBVink2001() const {
18431843
else if (utils::Compare(teff, VINK_MASS_LOSS_BISTABILITY_TEMP) > 0) {
18441844
SHOW_WARN_IF(utils::Compare(teff, VINK_MASS_LOSS_MAXIMUM_TEMP) > 0, ERROR::HIGH_TEFF_WINDS); // show warning if winds being used outside comfort zone
18451845

1846-
double v = 2.6; // v_inf/v_esc
1847-
v = v * PPOW(m_Metallicity / ZSOL, OPTIONS->ScaleTerminalWindVelocityWithMetallicityPower()); // Scale Vinf with metallicity
1846+
double v = 2.6; // v_inf/v_esc
1847+
v = v * PPOW(m_Metallicity / ZSOL_ANDERS, OPTIONS->ScaleTerminalWindVelocityWithMetallicityPower()); // Scale Vinf with metallicity
18481848

18491849
double logMdotOB = -6.697 +
18501850
(2.194 * log10(m_Luminosity / 1.0E5)) -
18511851
(1.313 * log10(m_Mass / 30.0)) -
18521852
(1.226 * log10(v / 2.0)) +
1853-
(0.85 * LogMetallicityXi()) +
1853+
(0.85 * LogMetallicityXiAnders()) +
18541854
(0.933 * log10(teff / 40000.0)) -
18551855
(10.92 * log10(teff / 40000.0) * log10(teff/40000.0));
18561856

@@ -1885,7 +1885,7 @@ double BaseStar::CalculateMassLossRateOBVinkSander2021() const {
18851885

18861886
double teff = m_Temperature * TSOL;
18871887
double Gamma = EDDINGTON_PARAMETER_FACTOR * m_Luminosity / m_Mass;
1888-
double charrho = -14.94 + (3.1857 * Gamma) + (zExp * LogMetallicityXi());
1888+
double charrho = -14.94 + (3.1857 * Gamma) + (zExp * LogMetallicityXiAnders());
18891889
double T2 = ( 61.2 + (2.59 * charrho) ) * 1000.0; // typically around 25000.0, higher jump first as in Vink python recipe
18901890
double T1 = ( 100.0 + (6.0 * charrho) ) * 1000.0; // typically around 20000.0, has similar behavior when fixed
18911891

@@ -1901,7 +1901,7 @@ double BaseStar::CalculateMassLossRateOBVinkSander2021() const {
19011901
(2.210 * logL5) -
19021902
(1.339 * logM30) -
19031903
(1.601 * log10(V / 2.0)) +
1904-
(zExp2001 * LogMetallicityXi()) +
1904+
(zExp2001 * LogMetallicityXiAnders()) +
19051905
(1.07 * logT20);
19061906

19071907
rate = PPOW(10.0, logMdotOB);
@@ -1914,7 +1914,7 @@ double BaseStar::CalculateMassLossRateOBVinkSander2021() const {
19141914
(2.210 * logL5) -
19151915
(1.339 * logM30) -
19161916
(1.601 * log10(V / 2.0)) +
1917-
(zExp2001 * LogMetallicityXi()) +
1917+
(zExp2001 * LogMetallicityXiAnders()) +
19181918
(1.07 * logT20);
19191919

19201920
rate = PPOW(10.0, logMdotOB);
@@ -1927,7 +1927,7 @@ double BaseStar::CalculateMassLossRateOBVinkSander2021() const {
19271927
(2.194 * logL5) -
19281928
(1.313 * logM30) -
19291929
(1.226 * log10(V / 2.0)) +
1930-
(zExp * LogMetallicityXi()) +
1930+
(zExp * LogMetallicityXiAnders()) +
19311931
(0.933 * logT40) -
19321932
(10.92 * logT40 * logT40);
19331933

@@ -1952,8 +1952,8 @@ double BaseStar::CalculateMassLossRateOBVinkSander2021() const {
19521952
* @return Mass loss rate for hot OB stars in Msol yr^-1
19531953
*/
19541954
double BaseStar::CalculateMassLossRateOBKrticka2018() const {
1955-
1956-
double logMdot = -5.70 + 0.50 * LogMetallicityXi() + (1.61 - 0.12 * LogMetallicityXi()) * log10(m_Luminosity / 1.0E6);
1955+
1956+
double logMdot = -5.70 + 0.50 * LogMetallicityXiAsplund() + (1.61 - 0.12 * LogMetallicityXiAsplund()) * log10(m_Luminosity / 1.0E6);
19571957

19581958
return PPOW(10.0, logMdot);
19591959
}
@@ -2299,7 +2299,7 @@ double BaseStar::CalculateMassLossRateWolfRayetSanderVink2020(const double p_Mu)
22992299
if (utils::Compare(p_Mu, 1.0) < 0) {
23002300

23012301
double logL = log10(m_Luminosity);
2302-
double logZ = LogMetallicityXi();
2302+
double logZ = LogMetallicityXiAnders();
23032303

23042304
// Calculate alpha, L0 and Mdot10
23052305
double alpha = 0.32 * logZ + 1.4; // Equation 18 in Sander & Vink 2020
@@ -2367,7 +2367,7 @@ double BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(c
23672367
*/
23682368
double BaseStar::CalculateMassLossRateHeliumStarVink2017() const {
23692369

2370-
double logMdot = -13.3 + (1.36 * log10(m_Luminosity)) + (0.61 * LogMetallicityXi()); // Vink 2017 Eq. 1.
2370+
double logMdot = -13.3 + (1.36 * log10(m_Luminosity)) + (0.61 * LogMetallicityXiAnders()); // Vink 2017 Eq. 1.
23712371

23722372
return PPOW(10.0, logMdot);
23732373
}

‎src/BaseStar.h‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -119,9 +119,11 @@ class BaseStar {
119119
bool IsSNIA() const { return (m_SupernovaDetails.events.current & SN_EVENT::SNIA) == SN_EVENT::SNIA; }
120120
bool IsUSSN() const { return (m_SupernovaDetails.events.current & SN_EVENT::USSN) == SN_EVENT::USSN; }
121121
bool LBV_PhaseFlag() const { return m_LBVphaseFlag; }
122-
double LogMetallicityRho() const { return LogMetallicityXi() + 1.0; } // rho in Hurley+ 2000
123-
double LogMetallicitySigma() const { return m_Log10Metallicity; } // sigma in Hurley+ 2000
124-
double LogMetallicityXi() const { return m_Log10Metallicity - LOG10_ZSOL; } // xi in Hurley+ 2000
122+
double LogMetallicityRho() const { return LogMetallicityXiHurley() + 1.0; } // rho in Hurley+ 2000
123+
double LogMetallicitySigma() const { return m_Log10Metallicity; } // sigma in Hurley+ 2000
124+
double LogMetallicityXiHurley() const { return m_Log10Metallicity - LOG10_ZSOL_HURLEY; } // xi in Hurley+ 2000
125+
double LogMetallicityXiAnders() const { return m_Log10Metallicity - LOG10_ZSOL_ANDERS; } // log10(Z / ZSOL_ANDERS)
126+
double LogMetallicityXiAsplund() const { return m_Log10Metallicity - LOG10_ZSOL_ASPLUND; } // log10(Z / ZSOL_ASPLUND)
125127
double Luminosity() const { return m_Luminosity; }
126128
double MainSequenceCoreMass() const { return m_MainSequenceCoreMass; }
127129
double Mass() const { return m_Mass; }

‎src/EAGB.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1046,7 +1046,7 @@ STELLAR_TYPE EAGB::ResolveEnvelopeLoss(bool p_Force) {
10461046

10471047
m_Age = HeGB::CalculateAgeOnPhase_Static(m_Mass, m_COCoreMass, timescales(tHeMS), m_GBParams);
10481048

1049-
HeHG::CalculateGBParams_Static(m_Mass0, m_Mass, LogMetallicityXi(), m_MassCutoffs, m_AnCoefficients, m_BnCoefficients, m_GBParams);
1049+
HeHG::CalculateGBParams_Static(m_Mass0, m_Mass, LogMetallicityXiHurley(), m_MassCutoffs, m_AnCoefficients, m_BnCoefficients, m_GBParams);
10501050
m_Luminosity = HeGB::CalculateLuminosityOnPhase_Static(m_COCoreMass, gbParams(B), gbParams(D));
10511051

10521052
double R1, R2;

‎src/GiantBranch.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -265,7 +265,7 @@ void GiantBranch::CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams)
265265
gbParams(AHe) = CalculateHeRateConstant_Static();
266266

267267
gbParams(B) = CalculateCoreMass_Luminosity_B_Static(p_Mass);
268-
gbParams(D) = CalculateCoreMass_Luminosity_D_Static(p_Mass, LogMetallicityXi(), m_MassCutoffs);
268+
gbParams(D) = CalculateCoreMass_Luminosity_D_Static(p_Mass, LogMetallicityXiHurley(), m_MassCutoffs);
269269

270270
gbParams(p) = CalculateCoreMass_Luminosity_p_Static(p_Mass, m_MassCutoffs);
271271
gbParams(q) = CalculateCoreMass_Luminosity_q_Static(p_Mass, m_MassCutoffs);
@@ -1299,7 +1299,7 @@ double GiantBranch::CalculateRemnantMassByMaltsev2024(const double p_COCoreMass,
12991299

13001300
ST_VECTOR mtHist = MassTransferDonorHistory(); // mass transfer history vector
13011301
MT_CASE massTransferCase = MT_CASE::OTHER;
1302-
double log10Z = m_Log10Metallicity - LOG10_ZSOL; // log_{10} (Z/Zsol), for convenience
1302+
double log10Z = m_Log10Metallicity - LOG10_ZSOL_ASPLUND; // log_{10} (Z/Zsol), for convenience
13031303
double M1, M2, M3;
13041304

13051305
if (utils::Compare(p_COCoreMass, MALTSEV2024_MMIN) < 0) // NS formation regardless of metallicity and MT history

‎src/MainSequence.cpp‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -848,7 +848,7 @@ double MainSequence::CalculateInitialMainSequenceCoreMass(const double p_Mass, c
848848
// After full mixing not at ZAMS, use the approach from Brcek+ (2025)
849849
else {
850850
double h = PPOW(10.0, p_HeliumAbundanceCore * (p_HeliumAbundanceCore + 2.0) / 4.0);
851-
fmix = BRCEK_FMIX_COEFFICIENTS[0] + BRCEK_FMIX_COEFFICIENTS[1] * std::exp(-p_Mass * h / BRCEK_FMIX_COEFFICIENTS[2]) * PPOW(1.0 - BRCEK_FMIX_COEFFICIENTS[4] / (p_Mass * h), BRCEK_FMIX_COEFFICIENTS[3]);
851+
fmix = (BRCEK_FMIX_COEFFICIENTS[0] + BRCEK_FMIX_COEFFICIENTS[1] * std::exp(-p_Mass * h / BRCEK_FMIX_COEFFICIENTS[2])) * PPOW(1.0 - BRCEK_FMIX_COEFFICIENTS[4] / (p_Mass * h), BRCEK_FMIX_COEFFICIENTS[3]);
852852
}
853853
return fmix * p_Mass;
854854
}
@@ -967,7 +967,7 @@ double MainSequence::CalculateLifetimeOnPhase(const double p_Mass, const double
967967
double tHook = mu * p_TBGB;
968968

969969
// For mass < Mhook, x > mu (i.e. for stars without a hook)
970-
double x = std::max(0.95, std::min((0.95 - (0.03 * (LogMetallicityXi() + 0.30103))), 0.99));
970+
double x = std::max(0.95, std::min((0.95 - (0.03 * (LogMetallicityXiHurley() + 0.30103))), 0.99));
971971

972972
return std::max(tHook, (x * p_TBGB));
973973

@@ -1326,8 +1326,8 @@ double MainSequence::InterpolateGeEtAlQCrit(const QCRIT_PRESCRIPTION p_qCritPres
13261326
double interpolatedQCritForZ = p_massTransferEfficiencyBeta * interpolatedQCritUpperEff + (1.0 - p_massTransferEfficiencyBeta) * interpolatedQCritLowerEff; // Don't need to use nearest neighbor for this, beta is always between 0 and 1
13271327
qCritPerMetallicity[ii] = interpolatedQCritForZ;
13281328
}
1329-
double logZlo = -3; // log10(0.001)
1330-
double logZhi = LOG10_ZSOL; // log10(0.02)
1329+
double logZlo = -3; // log10(0.001)
1330+
double logZhi = LOG10_ZSOL_HURLEY; // log10(0.02)
13311331

13321332
return qCritPerMetallicity[1] + (m_Log10Metallicity - logZhi)*(qCritPerMetallicity[1] - qCritPerMetallicity[0])/(logZhi - logZlo);
13331333
}
@@ -1354,9 +1354,9 @@ std::tuple <DBL_VECTOR, DBL_VECTOR, DBL_VECTOR> MainSequence::InterpolateShikauc
13541354
double logZ = std::log10(p_Metallicity);
13551355

13561356
// Coefficients are given for these metallicities
1357-
double low = std::log10(0.1 * ZSOL_ASPLUND);
1358-
double middle = std::log10(1.0 / 3.0 * ZSOL_ASPLUND);
1359-
double high = std::log10(ZSOL_ASPLUND);
1357+
double low = std::log10(0.1 * ZSOL_HURLEY);
1358+
double middle = std::log10(1.0 / 3.0 * ZSOL_HURLEY);
1359+
double high = std::log10(ZSOL_HURLEY);
13601360

13611361
// common factors
13621362
double middle_logZ = middle - logZ;

‎src/changelog.h‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1593,8 +1593,12 @@
15931593
// - Implemented a LOGNORMAL NS CCSN kick magnitude distribution based on Disberg & Mandel, 2025
15941594
// 03.20.07 IM - June 25, 2025 - Enhancement:
15951595
// - Added a maximum threshold of 1000 km/s for Disberg & Mandel (2025) LOGNORMAL kicks, matching paper
1596+
// 03.20.08 AB - Jun 26, 2025 - Defect repair:
1597+
// - Fix for issue #400; correct Zsol values are now used in stellar wind prescriptions
1598+
// - To avoid ambiguous ZSOL, we now use ZSOL_HURLEY = 0.02, ZSOL_ANDERS = 0.019, and ZSOL_ASPLUND = 0.0142
1599+
// - Fixed error in MainSequence::CalculateInitialMainSequenceCoreMass()
15961600

1597-
const std::string VERSION_STRING = "03.20.06";
1601+
const std::string VERSION_STRING = "03.20.08";
15981602

15991603

16001604
# endif // __changelog_h__

0 commit comments

Comments
 (0)