From a07ae2a4a8b190fdfda3bdcf3e372f4dea0a7d0d Mon Sep 17 00:00:00 2001 From: Adam Brcek Date: Fri, 20 Jun 2025 17:00:27 +1000 Subject: [PATCH 1/6] Fixes and upgrades to CHE and MS mergers --- src/BaseBinaryStar.cpp | 11 +++++------ src/BaseStar.cpp | 1 + src/BaseStar.h | 2 ++ src/CH.cpp | 4 ++++ src/CH.h | 2 ++ src/MS_gt_07.h | 2 +- src/MainSequence.cpp | 41 ++++++++++++++++++++++++++++------------- src/MainSequence.h | 4 +--- src/Star.h | 1 + src/changelog.h | 2 ++ src/constants.h | 4 ++++ 11 files changed, 51 insertions(+), 23 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index dea202841..eed6da0c7 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1778,13 +1778,12 @@ void BaseBinaryStar::ResolveMainSequenceMerger() { (utils::Compare(m_Star1->MZAMS(), BRCEK_LOWER_MASS_LIMIT) >= 0) && (utils::Compare(m_Star2->MZAMS(), BRCEK_LOWER_MASS_LIMIT) >= 0)) { - double coreMass1 = m_Star1->MainSequenceCoreMass(); - double coreMass2 = m_Star2->MainSequenceCoreMass(); + // total hydrogen masses (from the core, envelope, and the intermediate region between the core and envelope) + double hydrogenMass1 = m_Star1->HydrogenAbundanceCore() * m_Star1->MainSequenceCoreMass() + (m_Star1->HydrogenAbundanceCore() + initialHydrogenFraction) * (m_Star1->InitialMainSequenceCoreMass() - m_Star1->MainSequenceCoreMass()) / 2.0 + initialHydrogenFraction * (m_Star1->Mass() - m_Star1->InitialMainSequenceCoreMass()); + double hydrogenMass2 = m_Star2->HydrogenAbundanceCore() * m_Star2->MainSequenceCoreMass() + (m_Star2->HydrogenAbundanceCore() + initialHydrogenFraction) * (m_Star2->InitialMainSequenceCoreMass() - m_Star2->MainSequenceCoreMass()) / 2.0 + initialHydrogenFraction * (m_Star2->Mass() - m_Star2->InitialMainSequenceCoreMass()); - double coreHeliumMass1 = m_Star1->HeliumAbundanceCore() * coreMass1; - double coreHeliumMass2 = m_Star2->HeliumAbundanceCore() * coreMass2; - - finalHydrogenMass = finalMass * initialHydrogenFraction - coreHeliumMass1 - coreHeliumMass2; + // final hydrogen mass in the merger product, assuming that the lost mass comes from the envelope + finalHydrogenMass = hydrogenMass1 + hydrogenMass2 - (m_Star1->Mass() + m_Star2->Mass() - finalMass) * initialHydrogenFraction; } else finalHydrogenMass = finalMass * initialHydrogenFraction - tau1 * TAMSCoreMass1 * initialHydrogenFraction - tau2 * TAMSCoreMass2 * initialHydrogenFraction; diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index 0e0f303a5..0cf90bd29 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -107,6 +107,7 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed, // initialise remaining member variables // Zero age main sequence parameters + m_InitialMainSequenceCoreMass = DEFAULT_INITIAL_DOUBLE_VALUE; // initialised in MS_gt_07 class if BRCEK core mass prescription is used m_RZAMS = CalculateRadiusAtZAMS(m_MZAMS); m_LZAMS = CalculateLuminosityAtZAMS(m_MZAMS); m_TZAMS = CalculateTemperatureOnPhase_Static(m_LZAMS, m_RZAMS); diff --git a/src/BaseStar.h b/src/BaseStar.h index 2108dea48..b563f5a78 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -106,6 +106,7 @@ class BaseStar { double HydrogenAbundanceSurface() const { return m_HydrogenAbundanceSurface; } double InitialHeliumAbundance() const { return m_InitialHeliumAbundance; } double InitialHydrogenAbundance() const { return m_InitialHydrogenAbundance; } + double InitialMainSequenceCoreMass() const { return m_InitialMainSequenceCoreMass; } bool IsAIC() const { return (m_SupernovaDetails.events.current & SN_EVENT::AIC) == SN_EVENT::AIC; } bool IsCCSN() const { return (m_SupernovaDetails.events.current & SN_EVENT::CCSN) == SN_EVENT::CCSN; } bool IsHeSD() const { return (m_SupernovaDetails.events.current & SN_EVENT::HeSD) == SN_EVENT::HeSD; } @@ -404,6 +405,7 @@ class BaseStar { // Zero Age Main Sequence double m_InitialHeliumAbundance; // Initial helium abundance (Y) double m_InitialHydrogenAbundance; // Initial hydrogen abundance (X) + double m_InitialMainSequenceCoreMass; // Initial main sequence core mass (used in BRCEK core mass prescription) double m_LZAMS; // ZAMS Luminosity double m_MZAMS; // ZAMS Mass double m_OmegaZAMS; // ZAMS Angular Frequency diff --git a/src/CH.cpp b/src/CH.cpp index 0eacf0679..7506b25ca 100755 --- a/src/CH.cpp +++ b/src/CH.cpp @@ -423,6 +423,10 @@ STELLAR_TYPE CH::EvolveToNextPhase() { if (m_Age < m_Timescales[static_cast(TIMESCALE::tMS)]) { // evolving off because of age? stellarType = STELLAR_TYPE::MS_GT_07; // no - must have spun down - evolve as MS star now m_CHE = false; // evolved CH->MS + + // if BRCEK core mass calculations enabled, initialise the core mass based on current mass and central helium fraction + if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0)) + m_MainSequenceCoreMass = MainSequence::CalculateInitialMainSequenceCoreMass(m_Mass, m_HeliumAbundanceCore); } else { // yes stellarType = STELLAR_TYPE::NAKED_HELIUM_STAR_MS; // evolve as HeMS star now diff --git a/src/CH.h b/src/CH.h index 77741f7f7..ebf8ffa54 100755 --- a/src/CH.h +++ b/src/CH.h @@ -89,6 +89,8 @@ class CH: virtual public BaseStar, public MS_gt_07 { bool ShouldEvolveOnPhase() const { return m_Age < m_Timescales[static_cast(TIMESCALE::tMS)] && (OPTIONS->OptimisticCHE() || Omega() >= m_OmegaCHE); } // Evolve on CHE phase if age in MS timescale and spinning at least as fast as CHE threshold void UpdateAgeAfterMassLoss(); + + void UpdateMainSequenceCoreMass(const double p_Dt, const double p_TotalMassLossRate) { }; // Do not use core mass calculations during CHE phase }; diff --git a/src/MS_gt_07.h b/src/MS_gt_07.h index 5714d9bbb..e4bb2ffa4 100755 --- a/src/MS_gt_07.h +++ b/src/MS_gt_07.h @@ -48,7 +48,7 @@ class MS_gt_07: virtual public BaseStar, public MainSequence { utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0 && // ZAMS mass >= BRCEK_LOWER_MASS_LIMIT? m_Time <= 0.0) { // star not yet aged past creation? // yes - initialise - m_InitialMainSequenceCoreMass = MainSequence::CalculateInitialMainSequenceCoreMass(m_MZAMS); + m_InitialMainSequenceCoreMass = MainSequence::CalculateInitialMainSequenceCoreMass(m_MZAMS, m_InitialHeliumAbundance); m_MainSequenceCoreMass = m_InitialMainSequenceCoreMass; m_Luminosity = MainSequence::CalculateLuminosityOnPhase(m_Age, m_Mass0, m_LZAMS0); m_Radius = MainSequence::CalculateRadiusOnPhase(m_Mass, m_Tau, m_RZAMS0); diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index 7f1e6d177..755c7fd49 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -289,8 +289,9 @@ double MainSequence::CalculateLuminosityOnPhase(const double p_Time, const doubl #define timescales(x) m_Timescales[static_cast(TIMESCALE::x)] // for convenience and readability - undefined at end of function // If BRCEK core prescription is used, return luminosity from Shikauchi et al. (2024) during core hydrogen burning (valid for MZAMS >= 15 Msol) or - // luminosity that smoothly connects MS and HG during MS hook (valid for MZAMS >= BRCEK_LOWER_MASS_LIMIT) - if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0)) { + // luminosity that smoothly connects MS and HG during MS hook (valid for MZAMS >= BRCEK_LOWER_MASS_LIMIT); do not use Shikauchi luminosity + // prescription during CHE + if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0) && !m_CHE) { if (utils::Compare(p_Time, 0.99 * timescales(tMS)) > 0) // star in MS hook? return CalculateLuminosityTransitionToHG(p_Mass, p_Time, p_LZAMS); else { @@ -826,18 +827,30 @@ DBL_DBL MainSequence::CalculateMainSequenceCoreMassBrcek(const double p_Dt, cons /* - * Calculate the initial convective core mass of a main sequence star using Equation (A3) from Shikauchi et al. (2024), - * also used for calculating core mass after MS merger + * Calculate the initial convective core mass of a main sequence star at ZAMS using Equation (A3) from Shikauchi+ (2024), + * or after full mixing (due to merger or CHE) for an arbitrary central helium fraction using the approach + * described in Brcek+ (2025) * - * double CalculateInitialMainSequenceCoreMass(const double p_MZAMS) + * double CalculateInitialMainSequenceCoreMass(const double p_Mass, const double p_HeliumAbundanceCore) * - * @param [IN] p_MZAMS Mass at ZAMS or after merger in Msol + * @param [IN] p_Mass Mass at ZAMS, after merger or after spin down of CH star in Msol + * @param [IN] p_HeliumAbundanceCore Central helium fraction * @return Mass of the convective core at ZAMS or after merger in Msol */ -double MainSequence::CalculateInitialMainSequenceCoreMass(const double p_MZAMS) const { - DBL_VECTOR fmixCoefficients = std::get<1>(SHIKAUCHI_COEFFICIENTS); - double fmix = fmixCoefficients[0] + fmixCoefficients[1] * std::exp(-p_MZAMS / fmixCoefficients[2]); - return fmix * p_MZAMS; +double MainSequence::CalculateInitialMainSequenceCoreMass(const double p_Mass, const double p_HeliumAbundanceCore) const { + + double fmix = 0.0; + // At ZAMS, use the approach from Shikauchi+ (2024) + if (utils::Compare(p_HeliumAbundanceCore, m_InitialHeliumAbundance) == 0) { + DBL_VECTOR fmixCoefficients = std::get<1>(SHIKAUCHI_COEFFICIENTS); + fmix = fmixCoefficients[0] + fmixCoefficients[1] * std::exp(-p_Mass / fmixCoefficients[2]); + } + // After full mixing not at ZAMS, use the apporoach from Brcek+ (2025) + else { + double h = PPOW(10.0, p_HeliumAbundanceCore * (p_HeliumAbundanceCore + 2.0) / 4.0); + 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]); + } + return fmix * p_Mass; } @@ -1142,12 +1155,14 @@ void MainSequence::UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { m_Age = m_Tau * timescales(tMS); m_HeliumAbundanceCore = 1.0 - m_Metallicity - p_HydrogenMass / p_Mass; - m_HydrogenAbundanceCore = 1.0 - m_Metallicity - m_HeliumAbundanceCore; + m_HeliumAbundanceSurface = m_HeliumAbundanceCore; // abundances are the same throughout the star, assuming uniform mixing after merger + m_HydrogenAbundanceSurface = m_HydrogenAbundanceCore; + if ((OPTIONS->MainSequenceCoreMassPrescription() == CORE_MASS_PRESCRIPTION::BRCEK) && (utils::Compare(m_MZAMS, BRCEK_LOWER_MASS_LIMIT) >= 0)) { - m_InitialMainSequenceCoreMass = CalculateInitialMainSequenceCoreMass(p_Mass); // update initial mixing core mass - m_MainSequenceCoreMass = m_InitialMainSequenceCoreMass; // update core mass + m_InitialMainSequenceCoreMass = CalculateInitialMainSequenceCoreMass(p_Mass, m_HeliumAbundanceCore); // update initial mixing core mass + m_MainSequenceCoreMass = m_InitialMainSequenceCoreMass; // update core mass } UpdateAttributesAndAgeOneTimestep(0.0, 0.0, 0.0, true); diff --git a/src/MainSequence.h b/src/MainSequence.h index e567d515c..d378b1559 100644 --- a/src/MainSequence.h +++ b/src/MainSequence.h @@ -25,9 +25,7 @@ class MainSequence: virtual public BaseStar { protected: // member variables - double m_HeliumAbundanceCoreOut = m_InitialHeliumAbundance; // Helium abundance just outside the core, used for rejuvenation calculations - double m_InitialMainSequenceCoreMass = 0.0; // Initial mass of the mixing core is initialised in MS_gt_07 class // member functions - alphabetically double CalculateAlphaL(const double p_Mass) const; @@ -78,7 +76,7 @@ class MainSequence: virtual public BaseStar { double CalculateLuminosityShikauchi(const double p_CoreMass, const double p_HeliumAbundanceCore) const; double CalculateLuminosityTransitionToHG(const double p_Mass, const double p_Age, double const p_LZAMS) const; DBL_DBL CalculateMainSequenceCoreMassBrcek(const double p_Dt, const double p_MassLossRate); - double CalculateInitialMainSequenceCoreMass(const double p_MZAMS) const; + double CalculateInitialMainSequenceCoreMass(const double p_Mass, const double p_HeliumAbundanceCore) const; double CalculateMomentOfInertia() const { return (0.1 * (m_Mass) * m_Radius * m_Radius); } // k2 = 0.1 as defined in Hurley et al. 2000, after eq 109 double CalculatePerturbationMu() const { return 5.0; } // mu(MS) = 5.0 (Hurley et al. 2000, eqs 97 & 98) diff --git a/src/Star.h b/src/Star.h index be87e0d1c..57cc7e388 100755 --- a/src/Star.h +++ b/src/Star.h @@ -97,6 +97,7 @@ class Star { double HydrogenAbundanceSurface() const { return m_Star->HydrogenAbundanceSurface(); } double InitialHeliumAbundance() const { return m_Star->InitialHeliumAbundance(); } double InitialHydrogenAbundance() const { return m_Star->InitialHydrogenAbundance(); } + double InitialMainSequenceCoreMass() const { return m_Star->InitialMainSequenceCoreMass(); } bool IsAIC() const { return m_Star->IsAIC(); } bool IsCCSN() const { return m_Star->IsCCSN(); } bool IsDegenerate() const { return m_Star->IsDegenerate(); } diff --git a/src/changelog.h b/src/changelog.h index 8940c8ab8..4682f8287 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1583,6 +1583,8 @@ // - TPAGB stars should no longer experience supernovae if SN conditions are not satisfied, rather than defaulting to CCSN (corrects the partial fix in 03.10.02) // - Added new parameter (threshold mass, generally expected to be MCH or MECS) to CalculateCoreMassAtSupernova_Static() // - Removed McSN from GBParams, instead computed on the fly when needed +// 03.20.04 AB - Jun 20, 2025 - Defect repair: +// - Fixed behavior for MS mergers and CH stars if the BRCEK core mass prescription is used const std::string VERSION_STRING = "03.20.03"; diff --git a/src/constants.h b/src/constants.h index f7da2f0d1..7a4899c8b 100755 --- a/src/constants.h +++ b/src/constants.h @@ -3811,5 +3811,9 @@ const std::vector SHIKAUCHI_L_COEFFICIENTS = { {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 {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 }; +// Coefficients used to determine the initial convective core mass of MS star after full mixing (due to a merger or CHE) +// from Brcek et al. (2025) +const DBL_VECTOR BRCEK_FMIX_COEFFICIENTS = {0.897487516925644, -0.570499075062176, 56.9817940103868, 0.54483109433479, 1.72842835324164}; + #endif // __constants_h__ From 880ca510f90ab2d48af569d05d5662bc2ce1114a Mon Sep 17 00:00:00 2001 From: Adam Brcek Date: Mon, 23 Jun 2025 10:55:39 +1000 Subject: [PATCH 2/6] Update changelog.h --- src/changelog.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/changelog.h b/src/changelog.h index 4682f8287..19ce2eed8 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1583,10 +1583,11 @@ // - TPAGB stars should no longer experience supernovae if SN conditions are not satisfied, rather than defaulting to CCSN (corrects the partial fix in 03.10.02) // - Added new parameter (threshold mass, generally expected to be MCH or MECS) to CalculateCoreMassAtSupernova_Static() // - Removed McSN from GBParams, instead computed on the fly when needed -// 03.20.04 AB - Jun 20, 2025 - Defect repair: -// - Fixed behavior for MS mergers and CH stars if the BRCEK core mass prescription is used +// 03.20.04 AB - Jun 23, 2025 - Defect repair, enhancement: +// - Fixes to MS mergers and CHE when BRCEK core mass prescription is used -- MS core mass is now correctly initialised after full mixing in MS +// mergers and CH stars that spun down -const std::string VERSION_STRING = "03.20.03"; +const std::string VERSION_STRING = "03.20.04"; # endif // __changelog_h__ From f8b1be07c83306b32652f55eedbb59dce8d6168c Mon Sep 17 00:00:00 2001 From: Adam Brcek Date: Mon, 23 Jun 2025 11:02:14 +1000 Subject: [PATCH 3/6] Typo fix --- src/MainSequence.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index 755c7fd49..72890f76c 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -843,9 +843,9 @@ double MainSequence::CalculateInitialMainSequenceCoreMass(const double p_Mass, c // At ZAMS, use the approach from Shikauchi+ (2024) if (utils::Compare(p_HeliumAbundanceCore, m_InitialHeliumAbundance) == 0) { DBL_VECTOR fmixCoefficients = std::get<1>(SHIKAUCHI_COEFFICIENTS); - fmix = fmixCoefficients[0] + fmixCoefficients[1] * std::exp(-p_Mass / fmixCoefficients[2]); + fmix = fmixCoefficients[0] + fmixCoefficients[1] * std::exp(-p_Mass / fmixCoefficients[2]); } - // After full mixing not at ZAMS, use the apporoach from Brcek+ (2025) + // After full mixing not at ZAMS, use the approach from Brcek+ (2025) else { double h = PPOW(10.0, p_HeliumAbundanceCore * (p_HeliumAbundanceCore + 2.0) / 4.0); 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]); From c23bc48f53059575cc5718f8661a127bd6b6ed0e Mon Sep 17 00:00:00 2001 From: Adam Brcek Date: Mon, 23 Jun 2025 14:45:03 +1000 Subject: [PATCH 4/6] Fix to hydrogen mass calculations --- src/BaseBinaryStar.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index eed6da0c7..ba3a122ce 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1779,13 +1779,13 @@ void BaseBinaryStar::ResolveMainSequenceMerger() { (utils::Compare(m_Star2->MZAMS(), BRCEK_LOWER_MASS_LIMIT) >= 0)) { // total hydrogen masses (from the core, envelope, and the intermediate region between the core and envelope) - double hydrogenMass1 = m_Star1->HydrogenAbundanceCore() * m_Star1->MainSequenceCoreMass() + (m_Star1->HydrogenAbundanceCore() + initialHydrogenFraction) * (m_Star1->InitialMainSequenceCoreMass() - m_Star1->MainSequenceCoreMass()) / 2.0 + initialHydrogenFraction * (m_Star1->Mass() - m_Star1->InitialMainSequenceCoreMass()); - double hydrogenMass2 = m_Star2->HydrogenAbundanceCore() * m_Star2->MainSequenceCoreMass() + (m_Star2->HydrogenAbundanceCore() + initialHydrogenFraction) * (m_Star2->InitialMainSequenceCoreMass() - m_Star2->MainSequenceCoreMass()) / 2.0 + initialHydrogenFraction * (m_Star2->Mass() - m_Star2->InitialMainSequenceCoreMass()); + double hydrogenMass1 = m_Star1->HydrogenAbundanceCore() * m_Star1->MainSequenceCoreMass() + (m_Star1->HydrogenAbundanceCore() + m_Star1->HydrogenAbundanceSurface()) * (std::min(m_Star1->InitialMainSequenceCoreMass(), m_Star1->Mass()) - m_Star1->MainSequenceCoreMass()) / 2.0 + initialHydrogenFraction * std::max(m_Star1->Mass() - m_Star1->InitialMainSequenceCoreMass(), 0.0); + double hydrogenMass2 = m_Star2->HydrogenAbundanceCore() * m_Star2->MainSequenceCoreMass() + (m_Star2->HydrogenAbundanceCore() + m_Star2->HydrogenAbundanceSurface()) * (std::min(m_Star2->InitialMainSequenceCoreMass(), m_Star2->Mass()) - m_Star2->MainSequenceCoreMass()) / 2.0 + initialHydrogenFraction * std::max(m_Star2->Mass() - m_Star1->InitialMainSequenceCoreMass(), 0.0); - // final hydrogen mass in the merger product, assuming that the lost mass comes from the envelope - finalHydrogenMass = hydrogenMass1 + hydrogenMass2 - (m_Star1->Mass() + m_Star2->Mass() - finalMass) * initialHydrogenFraction; + // final hydrogen mass in the merger product, assuming that the lost mass comes from the envelope of the star with higher surface hydrogen abundance + finalHydrogenMass = hydrogenMass1 + hydrogenMass2 - (m_Star1->Mass() + m_Star2->Mass() - finalMass) * std::max(m_Star1->HydrogenAbundanceSurface(), m_Star2->HydrogenAbundanceSurface()); } - else finalHydrogenMass = finalMass * initialHydrogenFraction - tau1 * TAMSCoreMass1 * initialHydrogenFraction - tau2 * TAMSCoreMass2 * initialHydrogenFraction; + else finalHydrogenMass = finalMass * initialHydrogenFraction - tau1 * TAMSCoreMass1 * initialHydrogenFraction - tau2 * TAMSCoreMass2 * initialHydrogenFraction; m_SemiMajorAxis = std::numeric_limits::infinity(); // set separation to infinity to avoid subsequent fake interactions with a massless companion (RLOF, CE, etc.) From 02468cd94c0f1685885fc298794ac6fa64dacb40 Mon Sep 17 00:00:00 2001 From: Adam Brcek Date: Mon, 23 Jun 2025 14:55:13 +1000 Subject: [PATCH 5/6] Typo fix --- src/BaseBinaryStar.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index ba3a122ce..dc0290715 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1780,7 +1780,7 @@ void BaseBinaryStar::ResolveMainSequenceMerger() { // total hydrogen masses (from the core, envelope, and the intermediate region between the core and envelope) double hydrogenMass1 = m_Star1->HydrogenAbundanceCore() * m_Star1->MainSequenceCoreMass() + (m_Star1->HydrogenAbundanceCore() + m_Star1->HydrogenAbundanceSurface()) * (std::min(m_Star1->InitialMainSequenceCoreMass(), m_Star1->Mass()) - m_Star1->MainSequenceCoreMass()) / 2.0 + initialHydrogenFraction * std::max(m_Star1->Mass() - m_Star1->InitialMainSequenceCoreMass(), 0.0); - double hydrogenMass2 = m_Star2->HydrogenAbundanceCore() * m_Star2->MainSequenceCoreMass() + (m_Star2->HydrogenAbundanceCore() + m_Star2->HydrogenAbundanceSurface()) * (std::min(m_Star2->InitialMainSequenceCoreMass(), m_Star2->Mass()) - m_Star2->MainSequenceCoreMass()) / 2.0 + initialHydrogenFraction * std::max(m_Star2->Mass() - m_Star1->InitialMainSequenceCoreMass(), 0.0); + double hydrogenMass2 = m_Star2->HydrogenAbundanceCore() * m_Star2->MainSequenceCoreMass() + (m_Star2->HydrogenAbundanceCore() + m_Star2->HydrogenAbundanceSurface()) * (std::min(m_Star2->InitialMainSequenceCoreMass(), m_Star2->Mass()) - m_Star2->MainSequenceCoreMass()) / 2.0 + initialHydrogenFraction * std::max(m_Star2->Mass() - m_Star2->InitialMainSequenceCoreMass(), 0.0); // final hydrogen mass in the merger product, assuming that the lost mass comes from the envelope of the star with higher surface hydrogen abundance finalHydrogenMass = hydrogenMass1 + hydrogenMass2 - (m_Star1->Mass() + m_Star2->Mass() - finalMass) * std::max(m_Star1->HydrogenAbundanceSurface(), m_Star2->HydrogenAbundanceSurface()); From 5ab3409192390b587db7a61604f637432bf55d44 Mon Sep 17 00:00:00 2001 From: Adam Brcek Date: Mon, 23 Jun 2025 16:28:14 +1000 Subject: [PATCH 6/6] Fix to hydrogen mass calculation --- src/BaseBinaryStar.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index dc0290715..7a70594e4 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1779,8 +1779,8 @@ void BaseBinaryStar::ResolveMainSequenceMerger() { (utils::Compare(m_Star2->MZAMS(), BRCEK_LOWER_MASS_LIMIT) >= 0)) { // total hydrogen masses (from the core, envelope, and the intermediate region between the core and envelope) - double hydrogenMass1 = m_Star1->HydrogenAbundanceCore() * m_Star1->MainSequenceCoreMass() + (m_Star1->HydrogenAbundanceCore() + m_Star1->HydrogenAbundanceSurface()) * (std::min(m_Star1->InitialMainSequenceCoreMass(), m_Star1->Mass()) - m_Star1->MainSequenceCoreMass()) / 2.0 + initialHydrogenFraction * std::max(m_Star1->Mass() - m_Star1->InitialMainSequenceCoreMass(), 0.0); - double hydrogenMass2 = m_Star2->HydrogenAbundanceCore() * m_Star2->MainSequenceCoreMass() + (m_Star2->HydrogenAbundanceCore() + m_Star2->HydrogenAbundanceSurface()) * (std::min(m_Star2->InitialMainSequenceCoreMass(), m_Star2->Mass()) - m_Star2->MainSequenceCoreMass()) / 2.0 + initialHydrogenFraction * std::max(m_Star2->Mass() - m_Star2->InitialMainSequenceCoreMass(), 0.0); + double hydrogenMass1 = m_Star1->HydrogenAbundanceCore() * m_Star1->MainSequenceCoreMass() + (m_Star1->HydrogenAbundanceCore() + m_Star1->HydrogenAbundanceSurface()) * (std::min(m_Star1->InitialMainSequenceCoreMass(), m_Star1->Mass()) - m_Star1->MainSequenceCoreMass()) / 2.0 + m_Star1->HydrogenAbundanceSurface() * std::max(m_Star1->Mass() - m_Star1->InitialMainSequenceCoreMass(), 0.0); + double hydrogenMass2 = m_Star2->HydrogenAbundanceCore() * m_Star2->MainSequenceCoreMass() + (m_Star2->HydrogenAbundanceCore() + m_Star2->HydrogenAbundanceSurface()) * (std::min(m_Star2->InitialMainSequenceCoreMass(), m_Star2->Mass()) - m_Star2->MainSequenceCoreMass()) / 2.0 + m_Star2->HydrogenAbundanceSurface() * std::max(m_Star2->Mass() - m_Star2->InitialMainSequenceCoreMass(), 0.0); // final hydrogen mass in the merger product, assuming that the lost mass comes from the envelope of the star with higher surface hydrogen abundance finalHydrogenMass = hydrogenMass1 + hydrogenMass2 - (m_Star1->Mass() + m_Star2->Mass() - finalMass) * std::max(m_Star1->HydrogenAbundanceSurface(), m_Star2->HydrogenAbundanceSurface());