From 0a5183c6b683aea067042d2bd6bb1feff8877e70 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 28 Nov 2024 08:02:51 +1100 Subject: [PATCH 01/12] Using core radii to determine CE survival --- src/BaseBinaryStar.cpp | 8 ++--- src/BaseStar.h | 3 ++ src/GiantBranch.h | 3 -- src/HG.cpp | 17 +++++++++ src/HG.h | 80 +++++++++++++++++++++--------------------- src/HeHG.cpp | 16 +++++++++ src/HeHG.h | 1 + src/HeMS.h | 2 ++ src/Remnants.h | 3 ++ src/Star.h | 2 ++ 10 files changed, 88 insertions(+), 47 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 4ef7a6487..43f3a11cf 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1640,6 +1640,10 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { } } + if (utils::Compare(m_SemiMajorAxis, 0.0) <= 0 || utils::Compare(m_Star1->CalculateRemnantRadius() + m_Star2->CalculateRemnantRadius(), m_SemiMajorAxis * AU_TO_RSOL) > 0) { // catch merger in CE here, do not update stars + m_Flags.stellarMerger = true; + } + if (!m_Flags.stellarMerger) { STELLAR_TYPE stellarType1 = m_Star1->StellarType(); // star 1 stellar type before resolving envelope loss @@ -1668,10 +1672,6 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { } } - if (utils::Compare(m_SemiMajorAxis, 0.0) <= 0 || utils::Compare(m_Star1->Radius() + m_Star2->Radius(), m_SemiMajorAxis * AU_TO_RSOL) > 0) { - m_Flags.stellarMerger = true; - } - // if stars are evolving as CHE stars, update their rotational frequency under the assumption of tidal locking if tides are not enabled if (OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::NONE) { double omega = OrbitalAngularVelocity(); // orbital angular velocity diff --git a/src/BaseStar.h b/src/BaseStar.h index 1a29906d9..bb8bff218 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -281,6 +281,9 @@ class BaseStar { virtual double CalculateRadialExtentConvectiveEnvelope() const { return 0.0; } // Default for stars with no convective envelope virtual double CalculateRadiusOnPhaseTau(const double p_Mass, const double p_Tau) const { return 0.0; } // Only defined for MS stars + + double CalculateRemnantRadius() const { return Radius(); } // relevant for MS stars, over-written for GB stars + void CalculateSNAnomalies(const double p_Eccentricity); diff --git a/src/GiantBranch.h b/src/GiantBranch.h index 93c432029..094e4063e 100755 --- a/src/GiantBranch.h +++ b/src/GiantBranch.h @@ -21,9 +21,6 @@ class GiantBranch: virtual public BaseStar, public MainSequence { GiantBranch(const BaseStar &p_BaseStar) : BaseStar(p_BaseStar), MainSequence(p_BaseStar) {} - virtual double CalculateRemnantRadius() const; - - protected: diff --git a/src/HG.cpp b/src/HG.cpp index 861de1f7f..41a45adf3 100755 --- a/src/HG.cpp +++ b/src/HG.cpp @@ -2,6 +2,7 @@ #include "HeMS.h" #include "HeWD.h" + /////////////////////////////////////////////////////////////////////////////////////// // // // COEFFICIENT AND CONSTANT CALCULATIONS ETC. // @@ -861,6 +862,22 @@ double HG::CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const } +/* + * Calculate radius of the remnant the star would become if it lost all of its + * envelope immediately (i.e. M = Mc, coreMass) + * + * Hurley et al. 2000, end of section 6 + * + * + * double CalculateRemnantRadius() + * + * @return Radius of remnant core in Rsol + */ +double HG::CalculateRemnantRadius() const { + return HeMS::CalculateRadiusOnPhase_Static(m_CoreMass, 0.0); +} + + /////////////////////////////////////////////////////////////////////////////////////// // // // MASS CALCULATIONS // diff --git a/src/HG.h b/src/HG.h index 1b7ea3e90..8f1657fac 100755 --- a/src/HG.h +++ b/src/HG.h @@ -9,54 +9,52 @@ #include #include "GiantBranch.h" +#include "HeMS.h" -// JR: todo: revisit this one day - sometimes HG works better as GiantBranch, sometimes not... -// Right now it is GiantBranch - figure out which is more appropriate - class BaseStar; class GiantBranch; class HG: virtual public BaseStar, public GiantBranch { - + public: - + HG() { m_StellarType = STELLAR_TYPE::HERTZSPRUNG_GAP; }; HG(const BaseStar &p_BaseStar, const bool p_Initialise = true) : BaseStar(p_BaseStar), GiantBranch(p_BaseStar) { - m_StellarType = STELLAR_TYPE::HERTZSPRUNG_GAP; // Set stellar type + m_StellarType = STELLAR_TYPE::HERTZSPRUNG_GAP; // Set stellar type if (p_Initialise) Initialise(); // Initialise if required } - + HG* Clone(const OBJECT_PERSISTENCE p_Persistence, const bool p_Initialise = true) { - HG* clone = new HG(*this, p_Initialise); - clone->SetPersistence(p_Persistence); - return clone; + HG* clone = new HG(*this, p_Initialise); + clone->SetPersistence(p_Persistence); + return clone; } - + static HG* Clone(HG& p_Star, const OBJECT_PERSISTENCE p_Persistence, const bool p_Initialise = true) { - HG* clone = new HG(p_Star, p_Initialise); - clone->SetPersistence(p_Persistence); - return clone; + HG* clone = new HG(p_Star, p_Initialise); + clone->SetPersistence(p_Persistence); + return clone; } - + MT_CASE DetermineMassTransferTypeAsDonor() const { return MT_CASE::B; } // Always case B - - + + protected: - + void Initialise() { - + m_Tau = 0.0; // Start of phase - + // update stellar properties at start of HG phase (since core definition changes) CalculateGBParams(); CalculateTimescales(); - // Initialise timescales + // Initialise timescales m_Age = m_Timescales[static_cast(TIMESCALE::tMS)]; // Set age appropriately - + // update effective "initial" mass (m_Mass0) so that the core mass is at least equal to the minimum core mass but no more than total mass - // (only relevant if RetainCoreMassDuringCaseAMassTransfer()) + // (only relevant if RetainCoreMassDuringCaseAMassTransfer()) if (utils::Compare(CalculateCoreMassOnPhase(m_Mass0, m_Age), std::min(m_Mass, MinimumCoreMass())) < 0) { double desiredCoreMass = std::min(m_Mass, MinimumCoreMass()); // desired core mass m_Mass0 = Mass0ToMatchDesiredCoreMass(this, desiredCoreMass); // use root finder to find new core mass estimate @@ -69,55 +67,57 @@ class HG: virtual public BaseStar, public GiantBranch { } EvolveOnPhase(0.0); } - - + + // member functions - alphabetically double CalculateCOCoreMassAtPhaseEnd() const { return 0.0; } // McCO(HG) = 0.0 double CalculateCOCoreMassOnPhase() const { return 0.0; } // McCO(HG) = 0.0 - + double CalculateCoreMassAt2ndDredgeUp(const DBL_VECTOR &p_GBParams) { return p_GBParams[static_cast(GBP::McDU)]; } // NO-OP double CalculateCoreMassAtPhaseEnd(const double p_Mass) const; double CalculateCoreMassAtPhaseEnd() const { return CalculateCoreMassAtPhaseEnd(m_Mass0); } // Use class member variables double CalculateCoreMassOnPhase(const double p_Mass, const double p_Time) const; double CalculateCoreMassOnPhase() const { return CalculateCoreMassOnPhase(m_Mass0, m_Age); } // Use class member variables double CalculateCoreMassOnPhaseIgnoringPreviousCoreMass(const double p_Mass, const double p_Time) const; // Ignore previous core mass constraint when computing expected core mass - + double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const; double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 0.25; } // As coded in BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details. - + double CalculateHeCoreMassAtPhaseEnd() const { return m_CoreMass; } // McHe(HG) = Core Mass double CalculateHeCoreMassOnPhase() const { return m_CoreMass; } // McHe(HG) = Core Mass - double CalculateHeliumAbundanceCoreAtPhaseEnd() const { return 1.0 - m_Metallicity; } - double CalculateHeliumAbundanceCoreOnPhase() const { return 1.0 - m_Metallicity; } // Use class member variables + double CalculateHeliumAbundanceCoreAtPhaseEnd() const { return 1.0 - m_Metallicity; } + double CalculateHeliumAbundanceCoreOnPhase() const { return 1.0 - m_Metallicity; } // Use class member variables double CalculateHeliumAbundanceSurfaceAtPhaseEnd() const { return CalculateHeliumAbundanceSurfaceOnPhase(); } - double CalculateHeliumAbundanceSurfaceOnPhase() const { return m_InitialHeliumAbundance; } // Use class member variables + double CalculateHeliumAbundanceSurfaceOnPhase() const { return m_InitialHeliumAbundance; } // Use class member variables - double CalculateHydrogenAbundanceCoreAtPhaseEnd() const { return CalculateHydrogenAbundanceCoreOnPhase(); } - double CalculateHydrogenAbundanceCoreOnPhase(const double p_Tau) const; - double CalculateHydrogenAbundanceCoreOnPhase() const { return 0.0; } // Star has exhausted hydrogen in its core + double CalculateHydrogenAbundanceCoreAtPhaseEnd() const { return CalculateHydrogenAbundanceCoreOnPhase(); } + double CalculateHydrogenAbundanceCoreOnPhase(const double p_Tau) const; + double CalculateHydrogenAbundanceCoreOnPhase() const { return 0.0; } // Star has exhausted hydrogen in its core - double CalculateHydrogenAbundanceSurfaceAtPhaseEnd() const { return CalculateHydrogenAbundanceSurfaceOnPhase(); } + double CalculateHydrogenAbundanceSurfaceAtPhaseEnd() const { return CalculateHydrogenAbundanceSurfaceOnPhase(); } double CalculateHydrogenAbundanceSurfaceOnPhase() const { return m_InitialHydrogenAbundance; } // Use class member variables - - + + double CalculateLambdaDewi() const; double CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) const; double CalculateLambdaNanjingEnhanced(const int p_MassIndex, const STELLAR_POPULATION p_StellarPop) const; - + double CalculateLuminosityAtPhaseEnd(const double p_Mass) const; double CalculateLuminosityAtPhaseEnd() const { return CalculateLuminosityAtPhaseEnd(m_Mass0);} // Use class member variables double CalculateLuminosityOnPhase(const double p_Age, const double p_Mass) const; double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase(m_Age, m_Mass0); } // Use class member variables - + double CalculateMassTransferRejuvenationFactor() { return 1.0; } - + double CalculateRadiusAtPhaseEnd(const double p_Mass) const; double CalculateRadiusAtPhaseEnd() const { return CalculateRadiusAtPhaseEnd(m_Mass); } // Use class member variables double CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const double p_RZAMS) const; double CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase(m_Mass0, m_Tau, m_RZAMS0); } // Use class member variables + + double CalculateRemnantRadius() const; double CalculateRho(const double p_Mass) const; diff --git a/src/HeHG.cpp b/src/HeHG.cpp index 56d411af7..34bfbd485 100755 --- a/src/HeHG.cpp +++ b/src/HeHG.cpp @@ -1,6 +1,7 @@ #include "HeHG.h" #include "HeGB.h" #include "COWD.h" +#include "HeWD.h" /* @@ -234,6 +235,21 @@ double HeHG::CalculatePerturbationMu() const { return std::max(5.0 * ((McMax - m_CoreMass) / McMax), 0.0); //return non-negative value to avoid round-off issues } +/* + * Calculate radius of the remnant the star would become if it lost all of its + * envelope immediately (i.e. M = Mc, coreMass) + * + * Hurley et al. 2000, end of section 6 + * + * + * double CalculateRemnantRadius() + * + * @return Radius of remnant core in Rsol + */ +double HeHG::CalculateRemnantRadius() const { + return HeWD::CalculateRadiusOnPhase_Static(m_CoreMass); +} + /////////////////////////////////////////////////////////////////////////////////////// // // diff --git a/src/HeHG.h b/src/HeHG.h index 038dea565..7956ef76e 100755 --- a/src/HeHG.h +++ b/src/HeHG.h @@ -40,6 +40,7 @@ class HeHG: virtual public BaseStar, public HeMS { double CalculateConvectiveCoreRadius () const { return std::min(5.0 * CalculateRemnantRadius(), m_Radius); } // Last paragraph of section 6 of Hurley+ 2000 static void CalculateGBParams_Static(const double p_Mass0, const double p_Mass, const double p_LogMetallicityXi, const DBL_VECTOR &p_MassCutoffs, const DBL_VECTOR &p_AnCoefficients, const DBL_VECTOR &p_BnCoefficients, DBL_VECTOR &p_GBParams); + double CalculateRemnantRadius() const; protected: diff --git a/src/HeMS.h b/src/HeMS.h index 00770d6fa..5b38a32e1 100644 --- a/src/HeMS.h +++ b/src/HeMS.h @@ -52,6 +52,8 @@ class HeMS: virtual public BaseStar, public TPAGB { static DBL_DBL CalculateRadiusAtPhaseEnd_Static(const double p_Mass, const double p_Luminosity); static double CalculateRadiusAtZAMS_Static(const double p_Mass); static double CalculateRadiusOnPhase_Static(const double p_Mass, const double p_Tau); + + double CalculateRemnantRadius() const { return Radius(); } MT_CASE DetermineMassTransferTypeAsDonor() const { return MT_CASE::OTHER; } // Not A, B, C, or NONE diff --git a/src/Remnants.h b/src/Remnants.h index 762c154bf..587616564 100644 --- a/src/Remnants.h +++ b/src/Remnants.h @@ -24,6 +24,9 @@ class Remnants: virtual public BaseStar, public HeGB { // member functions static double CalculateRemnantMass_Static(const double p_COCoreMass) { return utils::Compare(p_COCoreMass, 7.0)< 0 ? 1.17 + (0.09 * p_COCoreMass) : p_COCoreMass; } // Hurley et al., eq 92; Hurley says that a >7.0 CO core leads to a collapse into a BH, but is ambiguous about the BH mass in this case -- we will assume it's just the CO core + + double CalculateRemnantRadius() const { return Radius(); } + protected: diff --git a/src/Star.h b/src/Star.h index 2b34e090e..16ef92323 100755 --- a/src/Star.h +++ b/src/Star.h @@ -197,6 +197,8 @@ class Star { double CalculateRadiusOnPhaseTau(const double p_Mass, const double p_Tau) const { return m_Star->CalculateRadiusOnPhaseTau(p_Mass, p_Tau); } + double CalculateRemnantRadius() { return m_Star->CalculateRemnantRadius(); } + void CalculateSNAnomalies(const double p_Eccentricity) { m_Star->CalculateSNAnomalies(p_Eccentricity); } double CalculateSNKickMagnitude(const double p_RemnantMass, From fe208c95c987d22a7650bcb65b3d81e00cc29564 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 28 Nov 2024 08:36:45 +1100 Subject: [PATCH 02/12] Fix compilation --- src/BaseBinaryStar.cpp | 8 +++++++- src/GiantBranch.h | 2 ++ src/HG.cpp | 15 --------------- src/HG.h | 3 --- src/changelog.h | 4 +++- 5 files changed, 12 insertions(+), 20 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 43f3a11cf..bdff22918 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1640,8 +1640,12 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { } } + std::cout<<"some of radii"<CalculateRemnantRadius() + m_Star2->CalculateRemnantRadius()<<"sep"<CalculateRemnantRadius() + m_Star2->CalculateRemnantRadius(), m_SemiMajorAxis * AU_TO_RSOL) > 0) { // catch merger in CE here, do not update stars + m_MassTransferTrackerHistory = MT_TRACKING::MERGER; m_Flags.stellarMerger = true; + std::cout<<"Yes, merger!"<TidesPrescription() == TIDES_PRESCRIPTION::NONE) { + if (!m_Flags.stellarMerger && OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::NONE) { double omega = OrbitalAngularVelocity(); // orbital angular velocity if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star1->SetOmega(omega); if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star2->SetOmega(omega); @@ -1690,6 +1694,8 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { } (void)PrintCommonEnvelope(); // print (log) common envelope details + + std::cout<<"Flag"<StellarType()<<"TypeB"<<(int)m_Star2->StellarType()< #include "GiantBranch.h" -#include "HeMS.h" class BaseStar; @@ -117,8 +116,6 @@ class HG: virtual public BaseStar, public GiantBranch { double CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const double p_RZAMS) const; double CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase(m_Mass0, m_Tau, m_RZAMS0); } // Use class member variables - double CalculateRemnantRadius() const; - double CalculateRho(const double p_Mass) const; double CalculateTauAtPhaseEnd() const { return 1.0; } // tau = 1.0 at end of HG diff --git a/src/changelog.h b/src/changelog.h index a658ad9b8..bf410470c 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1401,6 +1401,8 @@ // - Fix issue (likely introduced in 03.08.00) with the accretor not gaining mass appropriately // 03.09.01 RTW - Nov 27, 2024 - Enhancement: // - Added systemic velocity components x, y, and z to the output -const std::string VERSION_STRING = "03.09.01"; +// 03.09.03 IM - Nov 28, 2024 - Enhancement, defect repair: +// - Delay changing stellar types until after checking for whether remnant cores would touch in a common enevelope, use core radii instead (fix to #1286) +const std::string VERSION_STRING = "03.09.03"; # endif // __changelog_h__ From 6f99f8af7d60340128bdb7d5c185bf7a27c8347e Mon Sep 17 00:00:00 2001 From: Jeff Riley Date: Thu, 28 Nov 2024 09:35:09 +1100 Subject: [PATCH 03/12] fixes forward evolution after stellar merger --- src/BaseBinaryStar.cpp | 128 +++++++++++++++++++++-------------------- 1 file changed, 66 insertions(+), 62 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 387f6a981..70d8ee41a 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1644,16 +1644,13 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { } } - std::cout<<"some of radii"<CalculateRemnantRadius() + m_Star2->CalculateRemnantRadius()<<"sep"<CalculateRemnantRadius() + m_Star2->CalculateRemnantRadius(), m_SemiMajorAxis * AU_TO_RSOL) > 0) { // catch merger in CE here, do not update stars m_MassTransferTrackerHistory = MT_TRACKING::MERGER; m_Flags.stellarMerger = true; - std::cout<<"Yes, merger!"<StellarType(); // star 1 stellar type before resolving envelope loss STELLAR_TYPE stellarType2 = m_Star2->StellarType(); // star 2 stellar type before resolving envelope loss @@ -1678,28 +1675,29 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { if (m_Star1->StellarType() != stellarType1 || m_Star2->StellarType() != stellarType2) { // stellar type change? (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::STELLAR_TYPE_CHANGE_DURING_CEE); // yes - print (log) detailed output } - } - // if stars are evolving as CHE stars, update their rotational frequency under the assumption of tidal locking if tides are not enabled - if (!m_Flags.stellarMerger && OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::NONE) { - double omega = OrbitalAngularVelocity(); // orbital angular velocity - if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star1->SetOmega(omega); - if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star2->SetOmega(omega); - } + // if stars are evolving as CHE stars, update their rotational frequency under the assumption of tidal locking if tides are not enabled + if (!m_Flags.stellarMerger && OPTIONS->TidesPrescription() == TIDES_PRESCRIPTION::NONE) { + double omega = OrbitalAngularVelocity(); // orbital angular velocity + if (m_Star1->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star1->SetOmega(omega); + if (m_Star2->StellarType() == STELLAR_TYPE::CHEMICALLY_HOMOGENEOUS) m_Star2->SetOmega(omega); + } - - m_Star1->SetPostCEEValues(); // squirrel away post CEE stellar values for star 1 - m_Star2->SetPostCEEValues(); // squirrel away post CEE stellar values for star 2 - SetPostCEEValues(m_SemiMajorAxis * AU_TO_RSOL, m_Eccentricity, rRLdfin1Rsol, rRLdfin2Rsol); // squirrel away post CEE binary values (checks for post-CE RLOF, so should be done at end) + m_Star1->SetPostCEEValues(); // squirrel away post CEE stellar values for star 1 + m_Star2->SetPostCEEValues(); // squirrel away post CEE stellar values for star 2 + SetPostCEEValues(m_SemiMajorAxis * AU_TO_RSOL, m_Eccentricity, rRLdfin1Rsol, rRLdfin2Rsol); // squirrel away post CEE binary values (checks for post-CE RLOF, so should be done at end) - if (m_RLOFDetails.immediateRLOFPostCEE == true && !OPTIONS->AllowImmediateRLOFpostCEToSurviveCommonEnvelope()) { // is there immediate post-CE RLOF which is not allowed? + if (m_RLOFDetails.immediateRLOFPostCEE == true && !OPTIONS->AllowImmediateRLOFpostCEToSurviveCommonEnvelope()) {// is there immediate post-CE RLOF which is not allowed? m_MassTransferTrackerHistory = MT_TRACKING::MERGER; m_Flags.stellarMerger = true; + } } +/* Ilya - should we execute the following statement if a stellar merger has occurred - we won't have stored post CEE stellar values for the stars... */ +/* I suppose we could still execute line 1688 above sto set the binary values */ +/* note we could have set the stellar merger flag as late as line 1692 above */ (void)PrintCommonEnvelope(); // print (log) common envelope details - std::cout<<"Flag"<StellarType()<<"TypeB"<<(int)m_Star2->StellarType()<IsSNevent() || m_Star2->IsSNevent())) { - EvaluateSupernovae(); // evaluate supernovae (both stars) if mass changes are responsible for a supernova - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_SN); // print (log) detailed output - if (HasOneOf({ STELLAR_TYPE::NEUTRON_STAR })) { - (void)PrintPulsarEvolutionParameters(PULSAR_RECORD_TYPE::POST_SN); // print (log) pulsar evolution parameters + if (!StellarMerger()) { // stellar merger? + // no - continue evolution + if ((m_Star1->IsSNevent() || m_Star2->IsSNevent())) { + EvaluateSupernovae(); // evaluate supernovae (both stars) if mass changes are responsible for a supernova + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_SN); // print (log) detailed output + if (HasOneOf({ STELLAR_TYPE::NEUTRON_STAR })) { + (void)PrintPulsarEvolutionParameters(PULSAR_RECORD_TYPE::POST_SN); // print (log) pulsar evolution parameters + } } - } - CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations + CalculateEnergyAndAngularMomentum(); // perform energy and angular momentum calculations - ProcessTides(p_Dt); // process tides if required + ProcessTides(p_Dt); // process tides if required - // assign new values to "previous" values, for following timestep - m_EccentricityPrev = m_Eccentricity; - m_SemiMajorAxisPrev = m_SemiMajorAxis; + // assign new values to "previous" values, for following timestep + m_EccentricityPrev = m_Eccentricity; + m_SemiMajorAxisPrev = m_SemiMajorAxis; - m_Star1->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star1 - m_Star2->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star2 + m_Star1->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star1 + m_Star2->UpdateMagneticFieldAndSpin(m_CEDetails.CEEnow, m_Dt * MYR_TO_YEAR * SECONDS_IN_YEAR, EPSILON_PULSAR); // update pulsar parameters for star2 + } } @@ -3146,45 +3147,48 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::PRE_STELLAR_TIMESTEP); // print (log) detailed output - error = EvolveOneTimestep(dt); // evolve the binary system one timestep - if (error != ERROR::NONE) { // SSE error for either constituent star? - evolutionStatus = EVOLUTION_STATUS::SSE_ERROR; // yes - stop evolution - } - else { // continue evolution + if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? + // yes + error = EvolveOneTimestep(dt); // evolve the binary system one timestep + if (error != ERROR::NONE) { // SSE error for either constituent star? + evolutionStatus = EVOLUTION_STATUS::SSE_ERROR; // yes - stop evolution + } + else { // continue evolution - (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::TIMESTEP_COMPLETED); // print (log) detailed output: this is after all changes made in the timestep + (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::TIMESTEP_COMPLETED); // print (log) detailed output: this is after all changes made in the timestep - if (stepNum >= OPTIONS->MaxNumberOfTimestepIterations()) evolutionStatus = EVOLUTION_STATUS::STEPS_UP; // number of timesteps for evolution exceeds maximum - else if (evolutionStatus == EVOLUTION_STATUS::CONTINUE && usingProvidedTimesteps && stepNum >= timesteps.size()) { // using user-provided timesteps and all consumed - evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_EXHAUSTED; // yes - set status - SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED); // show warning + if (stepNum >= OPTIONS->MaxNumberOfTimestepIterations()) evolutionStatus = EVOLUTION_STATUS::STEPS_UP; // number of timesteps for evolution exceeds maximum + else if (evolutionStatus == EVOLUTION_STATUS::CONTINUE && usingProvidedTimesteps && stepNum >= timesteps.size()) { // using user-provided timesteps and all consumed + evolutionStatus = EVOLUTION_STATUS::TIMESTEPS_EXHAUSTED; // yes - set status + SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED); // show warning + } } - } - if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? + if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) { // continue evolution? // yes - // if user selects to emit GWs, calculate the effects of radiation - // - note that this is placed before the call to ChooseTimestep() because when - // emitting GWs the timestep is a function of gravitational radiation - if (OPTIONS->EmitGravitationalRadiation()) CalculateGravitationalRadiation(); + // if user selects to emit GWs, calculate the effects of radiation + // - note that this is placed before the call to ChooseTimestep() because when + // emitting GWs the timestep is a function of gravitational radiation + if (OPTIONS->EmitGravitationalRadiation()) CalculateGravitationalRadiation(); - m_Star2->UpdatePreviousTimestepDuration(); // update stellar property for star2 - m_Star1->UpdatePreviousTimestepDuration(); // update stellar property for star1 + m_Star2->UpdatePreviousTimestepDuration(); // update stellar property for star2 + m_Star1->UpdatePreviousTimestepDuration(); // update stellar property for star1 - if (usingProvidedTimesteps) { // user-provided timesteps? - // select a timestep - // - don't quantise - // - don't apply timestep multiplier - // (we assume user wants the timesteps in the file) - // - // Open question: should we clamp this to NUCLEAR_MINIMUM_TIMESTEP? - dt = timesteps[stepNum]; - } - else { // no - not using user-provided timesteps - dt = ChooseTimestep(OPTIONS->TimestepMultiplier()); - } + if (usingProvidedTimesteps) { // user-provided timesteps? + // select a timestep + // - don't quantise + // - don't apply timestep multiplier + // (we assume user wants the timesteps in the file) + // + // Open question: should we clamp this to NUCLEAR_MINIMUM_TIMESTEP? + dt = timesteps[stepNum]; + } + else { // no - not using user-provided timesteps + dt = ChooseTimestep(OPTIONS->TimestepMultiplier()); + } - stepNum++; // increment stepNum + stepNum++; // increment stepNum + } } } From 1a412ce284858bc5ea62c21114b74f3ceef4fe61 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 28 Nov 2024 16:25:33 +1100 Subject: [PATCH 04/12] Switch both stars to Massless remnants during a CE merger, resolve #1265 --- src/BaseBinaryStar.cpp | 16 ++++++------- src/BaseStar.h | 2 ++ src/HG.h | 1 + src/MainSequence.cpp | 54 ++++++++++++++++++++---------------------- src/MainSequence.h | 2 ++ src/Star.h | 2 ++ src/changelog.h | 4 +++- 7 files changed, 44 insertions(+), 37 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 70d8ee41a..951daddd0 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1692,10 +1692,7 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { m_Flags.stellarMerger = true; } } - -/* Ilya - should we execute the following statement if a stellar merger has occurred - we won't have stored post CEE stellar values for the stars... */ -/* I suppose we could still execute line 1688 above sto set the binary values */ -/* note we could have set the stellar merger flag as late as line 1692 above */ + (void)PrintCommonEnvelope(); // print (log) common envelope details } @@ -1721,9 +1718,8 @@ void BaseBinaryStar::ResolveMainSequenceMerger() { double tau1 = m_Star1->Tau(); double tau2 = m_Star2->Tau(); - // /*ILYA*/ temporary solution, should use TAMS core mass - double TAMSCoreMass1 = 0.3 * mass1; - double TAMSCoreMass2 = 0.3 * mass2; + double TAMSCoreMass1 = m_Star1->TAMSCoreMass(); + double TAMSCoreMass2 = m_Star2->TAMSCoreMass(); double q = std::min(mass1 / mass2, mass2 / mass1); double phi = 0.3 * q / (1.0 + q) / (1.0 + q); // fraction of mass lost in merger, Wang+ 2022, https://www.nature.com/articles/s41550-021-01597-5 @@ -3088,8 +3084,12 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { if (StellarMerger() && !HasOneOf({ STELLAR_TYPE::MASSLESS_REMNANT })) { // have stars merged without merger already being resolved? if (m_Star1->IsOneOf(MAIN_SEQUENCE) && m_Star2->IsOneOf(MAIN_SEQUENCE) && OPTIONS->EvolveMainSequenceMergers()) // yes - both MS and evolving MS merger products? ResolveMainSequenceMerger(); // yes - handle main sequence mergers gracefully; no need to change evolution status - else + else { + // make both stars massless remnants if merging during CE, so this is recorded in the Switch log; eventually, will want to implement a more careful prescription for the merger product, perhaps allowing further evolution of the merger product + m_Star1->SwitchTo(STELLAR_TYPE::MASSLESS_REMNANT); + m_Star2->SwitchTo(STELLAR_TYPE::MASSLESS_REMNANT); evolutionStatus = EVOLUTION_STATUS::STELLAR_MERGER; // no - for now, stop evolution + } } else if (HasStarsTouching()) { // binary components touching? (should usually be avoided as MT or CE or merger should happen prior to this) evolutionStatus = EVOLUTION_STATUS::STARS_TOUCHING; // yes - stop evolution diff --git a/src/BaseStar.h b/src/BaseStar.h index bb8bff218..7830fe9f4 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -343,6 +343,8 @@ class BaseStar { void StashSupernovaDetails(const STELLAR_TYPE p_StellarType, const SSE_SN_RECORD_TYPE p_RecordType = SSE_SN_RECORD_TYPE::DEFAULT) { LOGGING->StashSSESupernovaDetails(this, p_StellarType, p_RecordType); } + virtual double TAMSCoreMass() const { return 0.0; } // except MS stars + virtual void UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { } // Default is NO-OP virtual void UpdateAgeAfterMassLoss() { } // Default is NO-OP diff --git a/src/HG.h b/src/HG.h index 5c2a0ea9e..cf23a5fc6 100755 --- a/src/HG.h +++ b/src/HG.h @@ -138,6 +138,7 @@ class HG: virtual public BaseStar, public GiantBranch { bool ShouldEvolveOnPhase() const { return (utils::Compare(m_Age, m_Timescales[static_cast(TIMESCALE::tBGB)]) < 0); } // Evolve on HG phase if age < Base Giant Branch timescale bool ShouldSkipPhase() const { return false; } // Never skip HG phase + double TAMSCoreMass() const { return 0.0; } void UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { } // Nothing to do for stars beyond the Main Sequence for now void UpdateAgeAfterMassLoss(); // Per Hurley et al. 2000, section 7.1 diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index fe69f6f3f..0836ece16 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -636,20 +636,7 @@ double MainSequence::CalculateConvectiveCoreRadius() const { * @return Mass of convective core in Msol */ double MainSequence::CalculateConvectiveCoreMass() const { - - // We need TAMSCoreMass, which is just the core mass at the start of the HG phase. - // Since we are on the main sequence here, we can clone this object as an HG object - // and, as long as it is initialised (to correctly set Tau to 0.0 on the HG phase), - // we can query the cloned object for its core mass. - // - // The clone should not evolve, and so should not log anything, but to be sure the - // clone does not participate in logging, we set its persistence to EPHEMERAL. - - HG *clone = HG::Clone(static_cast(const_cast(*this)), OBJECT_PERSISTENCE::EPHEMERAL); - double TAMSCoreMass = clone->CoreMass(); // get core mass from clone - delete clone; clone = nullptr; // return the memory allocated for the clone - - double finalConvectiveCoreMass = TAMSCoreMass; + double finalConvectiveCoreMass = TAMSCoreMass(); // core mass at TAMS double initialConvectiveCoreMass = finalConvectiveCoreMass / 0.6; return (initialConvectiveCoreMass - m_Tau * (initialConvectiveCoreMass - finalConvectiveCoreMass)); } @@ -849,23 +836,34 @@ STELLAR_TYPE MainSequence::ResolveEnvelopeLoss(bool p_Force) { */ void MainSequence::UpdateMinimumCoreMass() { if (OPTIONS->RetainCoreMassDuringCaseAMassTransfer()) { - - // We need TAMSCoreMass, which is just the core mass at the start of the HG phase. - // Since we are on the main sequence here, we can clone this object as an HG object - // and, as long as it is initialised (to correctly set Tau to 0.0 on the HG phase), - // we can query the cloned object for its core mass. - // - // The clone should not evolve, and so should not log anything, but to be sure the - // clone does not participate in logging, we set its persistence to EPHEMERAL. - - HG *clone = HG::Clone(static_cast(const_cast(*this)), OBJECT_PERSISTENCE::EPHEMERAL); - double TAMSCoreMass = clone->CoreMass(); // get core mass from clone - delete clone; clone = nullptr; // return the memory allocated for the clone - - m_MinimumCoreMass = std::max(m_MinimumCoreMass, CalculateTauOnPhase() * TAMSCoreMass); // update minimum core mass + m_MinimumCoreMass = std::max(m_MinimumCoreMass, CalculateTauOnPhase() * TAMSCoreMass()); // update minimum core mass } } +/* + * Return the expected core mass at terminal age main sequence, i.e., at the start of the HG phase + * + * double TAMSCoreMass() const + * + * + * @return TAMS core Mass (Msol) + * + */ +double MainSequence::TAMSCoreMass() const { + // Since we are on the main sequence here, we can clone this object as an HG object + // and, as long as it is initialised (to correctly set Tau to 0.0 on the HG phase), + // we can query the cloned object for its core mass. + // + // The clone should not evolve, and so should not log anything, but to be sure the + // clone does not participate in logging, we set its persistence to EPHEMERAL. + + HG *clone = HG::Clone(static_cast(const_cast(*this)), OBJECT_PERSISTENCE::EPHEMERAL); + double TAMSCoreMass = clone->CoreMass(); // get core mass from clone + delete clone; clone = nullptr; // return the memory allocated for the clone + + return TAMSCoreMass; +} + /* * Sets the mass and age of a merge product of two main sequence stars diff --git a/src/MainSequence.h b/src/MainSequence.h index a83b4cbe7..36aa11b2c 100644 --- a/src/MainSequence.h +++ b/src/MainSequence.h @@ -107,6 +107,8 @@ class MainSequence: virtual public BaseStar { STELLAR_TYPE ResolveEnvelopeLoss(bool p_Force = false); bool ShouldEvolveOnPhase() const { return (m_Age < m_Timescales[static_cast(TIMESCALE::tMS)]); } // Evolve on MS phase if age in MS timescale + + double TAMSCoreMass() const; void UpdateInitialMass() { m_Mass0 = m_Mass; } // Per Hurley et al. 2000, section 7.1 diff --git a/src/Star.h b/src/Star.h index 16ef92323..c486cc1a9 100755 --- a/src/Star.h +++ b/src/Star.h @@ -262,6 +262,8 @@ class Star { STELLAR_TYPE SwitchTo(const STELLAR_TYPE p_StellarType, bool p_SetInitialType = false); + + double TAMSCoreMass() const { return m_Star->TAMSCoreMass(); } void UpdateAfterMerger(double p_Mass, double p_HydrogenMass) { m_Star->UpdateAfterMerger(p_Mass, p_HydrogenMass); } void UpdateAgeAfterMassLoss() { m_Star->UpdateAgeAfterMassLoss(); } diff --git a/src/changelog.h b/src/changelog.h index deb5f6fb3..975cb42be 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1405,7 +1405,9 @@ // - Fixed bugs in vector3d related to indexing and rotation // - Added tweak for circular systems at first SN, to fix the x-axis along the separation vector // 03.09.03 IM - Nov 28, 2024 - Enhancement, defect repair: -// - Delay changing stellar types until after checking for whether remnant cores would touch in a common enevelope, use core radii instead (fix to #1286) +// - Delay changing stellar types until after checking for whether remnant cores would touch in a common enevelope, use core radii instead (partial fix to #1286) +// - Define a new function MainSequence::TAMSCoreMass(); use it for determining the amount of He in a star during MS mergers +// - Switch both stars to Massless remnants during a CE merger, resolve #1265 const std::string VERSION_STRING = "03.09.03"; # endif // __changelog_h__ From 50f6028612fc7ffcea53794a66eb405fde5016a6 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 28 Nov 2024 16:52:04 +1100 Subject: [PATCH 05/12] Minor cleanup --- src/changelog.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/changelog.h b/src/changelog.h index 975cb42be..060e61c1a 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1344,7 +1344,7 @@ // - Changed the prescription for Tonset in the Picker+ models to take advantage of improved metallicity-dependent fits // 03.05.02 IM - Oct 10, 2024 - Enhancement, defect repair: // - Reverted IsCCSN() to include USSN following a change in 3.00.00 that inadvertently led to no binary orbit updates following USSNe -// - Include a call to EvolveBinary(0.0) on initialisation of evolved stellar types; this ensures that Switch logs include consistent data for the new stellar type +// - Include a call to EvolveOnPhase(0.0) on initialisation of evolved stellar types; this ensures that Switch logs include consistent data for the new stellar type // 03.06.00 IM - Oct 14, 2024 - Enhancement, code cleanup: // - Incorporating the Maltsev+ (2024) prescription for supernova remnant masses // - Minor fixes, including in fallback fraction for Schneider SN prescription, documentation From da0c89738b9b8e2b632c76ded58fd38826647f72 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 28 Nov 2024 16:55:32 +1100 Subject: [PATCH 06/12] Documentation update --- online-docs/pages/whats-new.rst | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index dfaef4514..a5fd6420a 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -3,6 +3,13 @@ What's new Following is a brief list of important updates to the COMPAS code. A complete record of changes can be found in the file ``changelog.h``. +**03.09.00 Nov 25, 2024** + +Improved nuclear timescale mass transfer: the nuclear timescale mass transfer rate is now set by the requirement that the star +ends the time step just filling its Roche lobe. +Fixed several significant mass-transfer issues, such as accretors not gaining mass appropriately and failures +in the root solver for fitting the star into the Roche lobe that were leading to artificial common envelopes and mergers. + **03.08.02 Nov 18, 2024** Updated implementation of the mass transfer stability critical mass ratio tables from the team of Hongwei Ge. From 8e44b2f08e0b102e5b426e333902eb83251819c7 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 28 Nov 2024 17:02:56 +1100 Subject: [PATCH 07/12] Documentation update --- .../Program options/program-options-list-defaults.rst | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst index d362c6d05..a58e637b3 100644 --- a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst +++ b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst @@ -847,7 +847,7 @@ Options: { JEANS, ISOTROPIC, CIRCUMBINARY, MACLEOD_LINEAR, ARBITRARY } |br| Default = ISOTROPIC **--mass-transfer-fa** |br| -Mass Transfer fraction accreted. |br| +Mass Transfer fraction accreted (beta). |br| Used when ``--mass-transfer-accretion-efficiency-prescription = FIXED_FRACTION``. |br| Default = 0.5 @@ -1122,10 +1122,11 @@ Enable mass loss due to pulsational-pair-instability (PPI). |br| Default = TRUE **--pulsational-pair-instability-prescription** |br| -Pulsational pair instability prescription. |br| +Pulsational pair instability prescription (only relevant when using ``--pulsational-pair-instability``. |br| Options: { HENDRIKS, COMPAS, STARTRACK, MARCHANT, FARMER } |br| ``HENDRIKS`` implements the prescription from Hendriks et al. 2023 |br| -``COMPAS``, ``STARTRACK`` and ``MARCHANT`` follow Woosley 2017, Belczynski et al. 2016, and Marchant et al. 2018, all as implemented in Stevenson et al. 2019. |br| +``COMPAS``, ``STARTRACK`` and ``MARCHANT`` follow Woosley 2017, Belczynski et al. 2016, and Marchant et al. 2018, +all as implemented in Stevenson et al. 2019. |br| ``FARMER`` follows Farmer et al. 2019 |br| Default = MARCHANT From 760fc47b7417b37dab056dfbc98c05378f82c86b Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 28 Nov 2024 17:49:08 +1100 Subject: [PATCH 08/12] Minor fixes, including to #1255, #1258 --- src/BaseStar.cpp | 4 ++++ src/MainSequence.cpp | 7 ++++--- src/changelog.h | 1 + 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index a550f7750..500fb76d3 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -2478,6 +2478,8 @@ double BaseStar::CalculateMassLossRateWolfRayetSanderVink2020(const double p_Mu) */ double BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(const double p_Mdot) const { + if (p_Mdot <= 0.0) return 0.0; // nothing to adjust + const double teffRef = 141.0E3; // reference effective temperature in Kelvin const double teffMin = 100.0E3; // minimum effective temperature in Kelvin to apply correction @@ -2493,6 +2495,8 @@ double BaseStar::CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(c logMdotCorrected = logMdotUncorrected; } + std::cout< Date: Thu, 28 Nov 2024 17:52:06 +1100 Subject: [PATCH 09/12] Added closing parenthesis :) --- .../Program options/program-options-list-defaults.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst index a58e637b3..5a380bd6d 100644 --- a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst +++ b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst @@ -1122,7 +1122,7 @@ Enable mass loss due to pulsational-pair-instability (PPI). |br| Default = TRUE **--pulsational-pair-instability-prescription** |br| -Pulsational pair instability prescription (only relevant when using ``--pulsational-pair-instability``. |br| +Pulsational pair instability prescription (only relevant when using ``--pulsational-pair-instability``). |br| Options: { HENDRIKS, COMPAS, STARTRACK, MARCHANT, FARMER } |br| ``HENDRIKS`` implements the prescription from Hendriks et al. 2023 |br| ``COMPAS``, ``STARTRACK`` and ``MARCHANT`` follow Woosley 2017, Belczynski et al. 2016, and Marchant et al. 2018, From 6bfcd7c16de2812ed206f6138615c2065560da9a Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 28 Nov 2024 20:45:49 +1100 Subject: [PATCH 10/12] Forgot _virtual_ ! :) --- src/BaseStar.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BaseStar.h b/src/BaseStar.h index 7830fe9f4..ad4262c37 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -282,7 +282,7 @@ class BaseStar { virtual double CalculateRadiusOnPhaseTau(const double p_Mass, const double p_Tau) const { return 0.0; } // Only defined for MS stars - double CalculateRemnantRadius() const { return Radius(); } // relevant for MS stars, over-written for GB stars + virtual double CalculateRemnantRadius() const { return Radius(); } // relevant for MS stars, over-written for GB stars void CalculateSNAnomalies(const double p_Eccentricity); From 8fcd7eaab9fc4ab17233fdf7625629b52055dccd Mon Sep 17 00:00:00 2001 From: jeffriley Date: Fri, 29 Nov 2024 06:50:41 +1100 Subject: [PATCH 11/12] Update MainSequence.cpp - formatting --- src/MainSequence.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index b807f9bc2..e9d9d20d2 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -494,8 +494,8 @@ double MainSequence::CalculateRadiusOnPhase(const double p_Mass, const double p_ // pow() is slow - use multiplication where it makes sense double tau_3 = tau * tau * tau; - double tau_10 = tau Date: Fri, 29 Nov 2024 06:51:35 +1100 Subject: [PATCH 12/12] Update MainSequence.cpp - more formatting --- src/MainSequence.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index e9d9d20d2..04b7657fe 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -837,7 +837,7 @@ STELLAR_TYPE MainSequence::ResolveEnvelopeLoss(bool p_Force) { */ void MainSequence::UpdateMinimumCoreMass() { if (OPTIONS->RetainCoreMassDuringCaseAMassTransfer()) { - m_MinimumCoreMass = std::max(m_MinimumCoreMass, CalculateTauOnPhase() * TAMSCoreMass()); // update minimum core mass + m_MinimumCoreMass = std::max(m_MinimumCoreMass, CalculateTauOnPhase() * TAMSCoreMass()); // update minimum core mass } }