From 920853ae48164f57814d76f71ecedad09fd1ca77 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Mon, 15 Apr 2024 11:20:58 +1000 Subject: [PATCH 01/10] Initial commit --- src/BaseBinaryStar.cpp | 15 ++- src/BaseStar.cpp | 36 +++++- src/BaseStar.h | 14 ++- src/BinaryConstituentStar.cpp | 20 +-- src/CHeB.h | 3 +- src/FGB.h | 2 +- src/GiantBranch.cpp | 36 ++---- src/GiantBranch.h | 5 +- src/HG.cpp | 10 +- src/HeGB.h | 1 + src/HeHG.h | 7 +- src/HeMS.h | 6 +- src/MainSequence.cpp | 47 ++++++- src/MainSequence.h | 4 +- src/Makefile | 16 +-- src/Remnants.h | 230 +++++++++++++++++----------------- src/Star.h | 12 +- src/TPAGB.h | 2 + 18 files changed, 271 insertions(+), 195 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 01e93a127..ce770f81e 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1583,16 +1583,21 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { // Two-stage common envelope, Hirai & Mandel (2022) else if ( OPTIONS->CommonEnvelopeFormalism() == CE_FORMALISM::TWO_STAGE ) { - double convectiveEnvelopeMass1 = m_Star1->CalculateConvectiveEnvelopeMass(); + double convectiveEnvelopeMass1, maxConvectiveEnvelopeMass1; + std::tie(convectiveEnvelopeMass1, maxConvectiveEnvelopeMass1) = m_Star1->CalculateConvectiveEnvelopeMass(); double radiativeIntershellMass1 = m_MassEnv1 - convectiveEnvelopeMass1; double endOfFirstStageMass1 = m_Mass1Final + radiativeIntershellMass1; - double convectiveEnvelopeMass2 = m_Star2->CalculateConvectiveEnvelopeMass(); + double convectiveEnvelopeMass2, maxConvectiveEnvelopeMass2; + std::tie(convectiveEnvelopeMass2, maxConvectiveEnvelopeMass2) = m_Star2->CalculateConvectiveEnvelopeMass(); double radiativeIntershellMass2 = m_MassEnv2 - convectiveEnvelopeMass2; double endOfFirstStageMass2 = m_Mass2Final + radiativeIntershellMass2; - // Stage 1: convective envelope removal on a dynamical timescale; assumes lambda = 1.5, motivated by bottom panel of Figure 3 of Hirai & Mandel 2022, including internal energy - double k1 = m_Star1->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (1.5 * alphaCE)) * m_Star1->Mass() * convectiveEnvelopeMass1 / m_Star1->Radius(); - double k2 = m_Star2->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (1.5 * alphaCE)) * m_Star2->Mass() * convectiveEnvelopeMass2 / m_Star2->Radius(); + // Stage 1: convective envelope removal on a dynamical timescale; assumes lambda = lambda_He + double lambda1 = m_Star1->CalculateConvectiveEnvelopeLambdaPicker(convectiveEnvelopeMass1, maxConvectiveEnvelopeMass1); + double lambda2 = m_Star1->CalculateConvectiveEnvelopeLambdaPicker(convectiveEnvelopeMass2, maxConvectiveEnvelopeMass2); + + double k1 = m_Star1->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (lambda1 * alphaCE)) * m_Star1->Mass() * convectiveEnvelopeMass1 / m_Star1->Radius(); + double k2 = m_Star2->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (lambda2 * alphaCE)) * m_Star2->Mass() * convectiveEnvelopeMass2 / m_Star2->Radius(); double k3 = m_Star1->Mass() * m_Star2->Mass() / periastronRsol; // assumes immediate circularisation at periastron at start of CE double k4 = (endOfFirstStageMass1 * endOfFirstStageMass2); double aFinalRsol = k4 / (k1 + k2 + k3); diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index 3ffb1a8ec..f4fdab1f6 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -2208,7 +2208,7 @@ double BaseStar::CalculateMassLossRateRSGVinkSabhahit2023() const { if (utils::Compare(logL, logLkink) < 0) { logMdot = -8.0 + 0.7 * logL - 0.7 * logM; } - else if (utils::Compare(logL, logLkink) >= 0) { + else { logMdot = -24.0 + 4.77 * logL - 3.99 * logM; } @@ -3804,19 +3804,43 @@ void BaseStar::CalculateBindingEnergies(const double p_CoreMass, const double p_ * Calculate convective envelope binding energy for the two-stage Hirai & Mandel (2022) common envelope formalism * * - * double CalculateConvectiveEnvelopeBindingEnergy(const double p_CoreMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_lambda) + * double CalculateConvectiveEnvelopeBindingEnergy(const double p_TotalMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_lambda) * - * @param [IN] p_CoreMass Core mass of the star (Msol) + * @param [IN] p_TotalMass Total mass of the star (Msol) * @param [IN] p_ConvectiveEnvelopeMass Mass of the convective outer envelope (Msol) * @param [IN] p_Radius Radius of the star (Rsol) - * @param [IN] p_Lambda Dimensionless parameter defining the binding energy + * @param [IN] p_Lambda Lambda parameter for the convective envelope * @return Binding energy (erg) */ -double BaseStar::CalculateConvectiveEnvelopeBindingEnergy(const double p_CoreMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_Lambda) { - return CalculateBindingEnergy(p_CoreMass, p_ConvectiveEnvelopeMass, p_Radius, p_Lambda); +double BaseStar::CalculateConvectiveEnvelopeBindingEnergy(const double p_TotalMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_lambda) { + return CalculateBindingEnergy(p_TotalMass - p_ConvectiveEnvelopeMass, p_ConvectiveEnvelopeMass, p_Radius, p_lambda); } +/* + * Approximates the lambda binding energy parameter of the outer convective envelope + * + * This is needed for the Hirai & Mandel (2022) two-stage CE formalism. + * Follows the fits of Picker, Hirai, Mandel (2024), arXiv:2402.13180 for lambda_He + * + * + * double BaseStar::CalculateConvectiveEnvelopeLambdaPicker(double p_convectiveEnvelopeMass, double p_maxConvectiveEnvelopeMass) + * + * @param [IN] p_convectiveEnvelopeMass Mass of the outer convective envelope shell + * @param [IN] p_maxConvectiveEnvelopeMass Maximum mass of the outher convective envelope shell at the onset of carbon burning + * @return Lambda binding energy parameter for the outer convective envelope + */ +double BaseStar::CalculateConvectiveEnvelopeLambdaPicker(double p_convectiveEnvelopeMass, double p_maxConvectiveEnvelopeMass) { + + double log10Z = log10 (m_Metallicity); + double m2 = 0.0023 * log10Z * log10Z + 0.0088 * log10Z + 0.013; // Eq. (12) and Table 1 of Picker, Hirai, Mandel (2024) + double b1 = m2 * m_Mass - 0.23; // Eq. (11) of Picker+ (2024) + double logLambda = p_convectiveEnvelopeMass / p_maxConvectiveEnvelopeMass > 0.3 ? + 0.42 * p_convectiveEnvelopeMass / p_maxConvectiveEnvelopeMass + b1 : 0.3 * 0.42 + b1; + + return exp(logLambda); +} + /////////////////////////////////////////////////////////////////////////////////////// // // // MISCELLANEOUS FUNCTIONS / CONTROL FUNCTIONS // diff --git a/src/BaseStar.h b/src/BaseStar.h index 551e22da2..95996fa6c 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -163,9 +163,13 @@ class BaseStar { void CalculateBindingEnergies(const double p_CoreMass, const double p_EnvMass, const double p_Radius); - double CalculateConvectiveEnvelopeBindingEnergy(const double p_CoreMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_Lambda); + virtual double CalculateConvectiveCoreMass () const {return 0.0;} + virtual double CalculateConvectiveCoreRadius () const {return 0.0;} - virtual double CalculateConvectiveEnvelopeMass() const { return 0.0; } + double CalculateConvectiveEnvelopeBindingEnergy(const double p_TotalMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_Lambda); + + double CalculateConvectiveEnvelopeLambdaPicker(double p_convectiveEnvelopeMass, double p_maxConvectiveEnvelopeMass); + virtual DBL_DBL CalculateConvectiveEnvelopeMass() const { return std::tuple (0.0, 0.0); } virtual double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate); virtual double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const { return 0.0; } // Default is 0.0 @@ -201,6 +205,8 @@ class BaseStar { double CalculateRadialExpansionTimescale() const { return CalculateRadialExpansionTimescale_Static(m_StellarType, m_StellarTypePrev, m_Radius, m_RadiusPrev, m_DtPrev); } // Use class member variables + virtual double CalculateRadialExtentConvectiveEnvelope() const { return 0.0; } // default for stars with no convective envelope + void CalculateSNAnomalies(const double p_Eccentricity); double CalculateSNKickMagnitude(const double p_RemnantMass, const double p_EjectaMass, const STELLAR_TYPE p_StellarType); @@ -426,7 +432,7 @@ class BaseStar { virtual double CalculateCOCoreMassAtPhaseEnd() const { return m_COCoreMass; } // Default is NO-OP virtual double CalculateCOCoreMassOnPhase() const { return m_COCoreMass; } // Default is NO-OP - + virtual double CalculateCoreMassAtPhaseEnd() const { return m_CoreMass; } // Default is NO-OP static double CalculateCoreMassGivenLuminosity_Static(const double p_Luminosity, const DBL_VECTOR &p_GBParams); virtual double CalculateCoreMassOnPhase() const { return m_CoreMass; } // Default is NO-OP @@ -532,8 +538,6 @@ class BaseStar { const double p_RadiusPrev, const double p_DtPrev); - virtual double CalculateRadialExtentConvectiveEnvelope() const { return m_Radius; } // default for stars with no convective envelope - virtual double CalculateRadiusAtPhaseEnd() const { return m_Radius; } // Default is NO-OP double CalculateRadiusAtZAMS(const double p_MZAMS) const; virtual double CalculateRadiusOnPhase() const { return m_Radius; } // Default is NO-OP diff --git a/src/BinaryConstituentStar.cpp b/src/BinaryConstituentStar.cpp index ec7982b88..c85b47280 100644 --- a/src/BinaryConstituentStar.cpp +++ b/src/BinaryConstituentStar.cpp @@ -236,7 +236,9 @@ void BinaryConstituentStar::SetPostCEEValues() { * m_CEDetails.CoreMass * m_CEDetails.bindingEnergy * m_CEDetails.lambda - * + * m_CEDetails.convectiveEnvelopeMass + * m_CEDetails.radiativeIntershellMass + * m_CEDetails.convectiveEnvelopeBindingEnergy * * void CalculateCommonEnvelopeValues() */ @@ -247,10 +249,6 @@ void BinaryConstituentStar::CalculateCommonEnvelopeValues() { m_CEDetails.CoreMass = CoreMass(); m_CEDetails.lambda = 0.0; // default - - m_CEDetails.convectiveEnvelopeMass = CalculateConvectiveEnvelopeMass(); - m_CEDetails.radiativeIntershellMass = Mass() - CoreMass() - m_CEDetails.convectiveEnvelopeMass; - m_CEDetails.convectiveEnvelopeBindingEnergy = 0.0; switch (OPTIONS->CommonEnvelopeLambdaPrescription()) { // which common envelope lambda prescription? @@ -278,7 +276,7 @@ void BinaryConstituentStar::CalculateCommonEnvelopeValues() { m_CEDetails.lambda = Lambda_Dewi(); m_CEDetails.bindingEnergy = BindingEnergy_Dewi(); break; - + default: // unknown prescription SHOW_WARN(ERROR::UNKNOWN_CE_LAMBDA_PRESCRIPTION, "Lambda = 0.0"); // show warning } @@ -286,7 +284,15 @@ void BinaryConstituentStar::CalculateCommonEnvelopeValues() { if (utils::Compare(m_CEDetails.lambda, 0.0) <= 0) m_CEDetails.lambda = 0.0; // force non-positive lambda to 0 m_CEDetails.lambda *= OPTIONS->CommonEnvelopeLambdaMultiplier(); // multiply by constant (program option, default = 1.0) - m_CEDetails.convectiveEnvelopeBindingEnergy = CalculateConvectiveEnvelopeBindingEnergy(CoreMass(), m_CEDetails.convectiveEnvelopeMass, Radius(), m_CEDetails.lambda); + + // Properties relevant for the Hirai & Mandel (2022) formalism + double maxConvectiveEnvelopeMass; + std::tie(m_CEDetails.convectiveEnvelopeMass, maxConvectiveEnvelopeMass) = CalculateConvectiveEnvelopeMass(); + m_CEDetails.radiativeIntershellMass = Mass() - CoreMass() - m_CEDetails.convectiveEnvelopeMass; + if ( OPTIONS->CommonEnvelopeFormalism() == CE_FORMALISM::TWO_STAGE ) + m_CEDetails.lambda = CalculateConvectiveEnvelopeLambdaPicker(m_CEDetails.convectiveEnvelopeMass, maxConvectiveEnvelopeMass); + + m_CEDetails.convectiveEnvelopeBindingEnergy = CalculateConvectiveEnvelopeBindingEnergy(Mass(), m_CEDetails.convectiveEnvelopeMass, Radius(), m_CEDetails.lambda); } diff --git a/src/CHeB.h b/src/CHeB.h index 187cd5fa7..390b57bc7 100755 --- a/src/CHeB.h +++ b/src/CHeB.h @@ -63,6 +63,7 @@ class CHeB: virtual public BaseStar, public FGB { double CalculateCOCoreMassOnPhase() const { return 0.0; } // McCO(CHeB) = 0.0 + double CalculateConvectiveCoreRadius () const { return CalculateRemnantRadius (); } // Last paragraph of section 6 of Hurley+ 2000 double CalculateCoreMassAtPhaseEnd() const { return CalculateCoreMassAtBAGB(m_Mass0); } // Use class member variables double CalculateCoreMassOnPhase(const double p_Mass, const double p_Tau) const; double CalculateCoreMassOnPhase() const { return CalculateCoreMassOnPhase(m_Mass0, m_Tau); } // Use class member variables @@ -85,8 +86,6 @@ class CHeB: virtual public BaseStar, public FGB { double CalculateLuminosityOnPhase(const double p_Mass, const double p_Tau) const; double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase(m_Mass0, m_Tau); } - double CalculateRadialExtentConvectiveEnvelope() const { return BaseStar::CalculateRadialExtentConvectiveEnvelope(); } // CHeB stars don't have a convective envelope - double CalculateRadiusAtBluePhaseEnd(const double p_Mass) const; double CalculateRadiusAtBluePhaseStart(const double p_Mass) const; diff --git a/src/FGB.h b/src/FGB.h index d189b6c95..dcec8efd7 100755 --- a/src/FGB.h +++ b/src/FGB.h @@ -56,7 +56,7 @@ class FGB: virtual public BaseStar, public HG { double CalculateLuminosityOnPhase(const double p_Time) const; double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase(m_Age); } // Use class member variables - double CalculateRadialExtentConvectiveEnvelope() const { return GiantBranch::CalculateRadialExtentConvectiveEnvelope(); } // Skip HG + double CalculateRadialExtentConvectiveEnvelope() const { return Radius() - CalculateConvectiveCoreRadius(); } // Skip HG double CalculateRadiusAtPhaseEnd(const double p_Mass, const double p_Luminosity) const { return GiantBranch::CalculateRadiusOnPhase(p_Mass, p_Luminosity); } // Skip HG - same as on phase double CalculateRadiusAtPhaseEnd() const { return CalculateRadiusAtPhaseEnd(m_Mass, m_Luminosity); } // Use class member variables diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index 842c73378..f896a8984 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -683,21 +683,14 @@ double GiantBranch::CalculateRemnantRadius() const { /* * Calculate the radial extent of the star's convective envelope (if it has one) * - * Hurley et al. 2000, sec. 2.3, particularly subsec. 2.3.1, eqs 36-40 - * - * (Technically not a radius calculation I suppose, but "radial extent" is close enough to put it with the radius calculations...) - * + * Hurley et al. 2002, sec. 2.3, particularly subsec. 2.3.1, eqs 36-40 * * double CalculateRadialExtentConvectiveEnvelope() * * @return Radial extent of the star's convective envelope in Rsol */ double GiantBranch::CalculateRadialExtentConvectiveEnvelope() const { - - BaseStar clone = *this; // clone this star so can manipulate without changes persisiting - clone.ResolveEnvelopeLoss(true); // update clone's attributes after envelope is lost - - return m_Radius - clone.Radius(); + return m_Radius - CalculateConvectiveCoreRadius(); } @@ -1044,29 +1037,26 @@ double GiantBranch::CalculateZetaConstantsByEnvelope(ZETA_PRESCRIPTION p_ZetaPre * Approximates the mass of the outer convective envelope. * * This is needed for the Hirai & Mandel (2022) two-stage CE formalism. - * Follows the fits of Picker, Hirai, Mandel (2023). + * Follows the fits of Picker, Hirai, Mandel (2024), arXiv:2402.13180 * * * double GiantBranch::CalculateConvectiveEnvelopeMass() * - * @return Mass of the outer convective envelope + * @return Tuple containing the mass of the outer convective envelope and its maximum value */ -double GiantBranch::CalculateConvectiveEnvelopeMass() const { +DBL_DBL GiantBranch::CalculateConvectiveEnvelopeMass() const { double log10Z = log10 (m_Metallicity); - HG clone = *this; // Create an HG star clone to query its core mass just after TAMS - double log10Ltams = log10 (clone.Luminosity()); + double MinterfMcoref = -0.021 * log10Z + 0.0038; //Eq. (8) of Picker+ 2024 + double Tonset = -139.8 * log10Z * log10Z - 981.7 * log10Z + 2798.3; //Eq. (6) of Picker+ 2024 + EAGB clone = *this; // Create an HG star clone to query its core mass just after BAGB + double Tmin=this->Temperature(); double Mcorefinal = CalculateCoreMassAtBAGB(m_Mass); - double Mconvmax = m_Mass - 1.1 * Mcorefinal; - double b1 = 14.4 * log10Z * log10Z + 57.4 * log10Z + 95.7; - double a2 = -16.9 * log10Z * log10Z - 81.9 * log10Z - 47.9; - double b2 = 184.0 * log10Z * log10Z + 872.2 * log10Z + 370.0; - double c2 = -660.1 * log10Z * log10Z - 3482.0 * log10Z + 1489.0; - double Tnorm = a2 * log10Ltams * log10Ltams + b2 * log10Ltams + c2; - double convectiveEnvelopeMass = Mconvmax / (1+exp(b1*(m_Temperature*TSOL-Tnorm)/Tnorm)); - convectiveEnvelopeMass = std::max(std::min(convectiveEnvelopeMass, (m_Mass - m_CoreMass)), 0.0); // Ensure that convective envelope mass is limited to [0, envelope mass] + std::cout<<"temp"<CoreMass()< (convectiveEnvelopeMass, Mconvmax); } diff --git a/src/GiantBranch.h b/src/GiantBranch.h index 3630cd043..e1f742479 100755 --- a/src/GiantBranch.h +++ b/src/GiantBranch.h @@ -25,8 +25,9 @@ class GiantBranch: virtual public BaseStar, public MainSequence { // member functions - alphabetically (sort of - some are grouped by functionality) - double CalculateConvectiveEnvelopeMass() const; - + double CalculateConvectiveCoreMass() const { return m_CoreMass; } + double CalculateConvectiveEnvelopeLambdaPicker(double p_convectiveEnvelopeMass, double p_maxConvectiveEnvelopeMass); + DBL_DBL CalculateConvectiveEnvelopeMass() const; static double CalculateCoreMassAt2ndDredgeUp_Static(const double p_McBAGB); double CalculateCoreMassAtBAGB(const double p_Mass) const; static double CalculateCoreMassAtBAGB_Static(const double p_Mass, const DBL_VECTOR &p_BnCoefficients); diff --git a/src/HG.cpp b/src/HG.cpp index 8240b6b9f..60d88b8f6 100755 --- a/src/HG.cpp +++ b/src/HG.cpp @@ -796,9 +796,7 @@ double HG::CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const /* * Calculate the radial extent of the star's convective envelope (if it has one) * - * Hurley et al. 2000, sec. 2.3, particularly subsec. 2.3.1, eqs 36-40 - * - * (Technically not a radius calculation I suppose, but "radial extent" is close enough to put it with the radius calculations...) + * Hurley et al. 2002, sec. 2.3, particularly subsec. 2.3.1, eqs 39-40 * * * double CalculateRadialExtentConvectiveEnvelope() @@ -806,11 +804,7 @@ double HG::CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const * @return Radial extent of the star's convective envelope in Rsol */ double HG::CalculateRadialExtentConvectiveEnvelope() const { - - BaseStar clone = *this; // clone this star so can manipulate without changes persisiting - clone.ResolveEnvelopeLoss(true); // update clone's attributes after envelope is lost - - return std::sqrt(m_Tau) * (m_Radius - clone.Radius()); + return std::sqrt(m_Tau) * (m_Radius - CalculateConvectiveCoreRadius()); } diff --git a/src/HeGB.h b/src/HeGB.h index 15c619c2b..a8526a16d 100755 --- a/src/HeGB.h +++ b/src/HeGB.h @@ -34,6 +34,7 @@ class HeGB: virtual public BaseStar, public HeHG { static double CalculateLuminosityOnPhase_Static(const double p_CoreMass, const double p_GBPB, const double p_GBPD); + double CalculateRadialExtentConvectiveEnvelope() const { return FGB::CalculateRadialExtentConvectiveEnvelope(); } static DBL_DBL CalculateRadiusOnPhase_Static(const double p_Mass, const double p_Luminosity); diff --git a/src/HeHG.h b/src/HeHG.h index cc711e14f..e76d0b31d 100755 --- a/src/HeHG.h +++ b/src/HeHG.h @@ -50,8 +50,10 @@ class HeHG: virtual public BaseStar, public HeMS { // member functions - aphabetically - double CalculateCOCoreMassAtPhaseEnd() const { return m_COCoreMass; } // NO-OP + double CalculateCOCoreMassAtPhaseEnd() const { return m_COCoreMass; } //*ILYA* check // NO-OP double CalculateCOCoreMassOnPhase() const; + + double CalculateConvectiveCoreMass() const { return m_CoreMass; } double CalculateCoreMassAtBAGB() const { return m_Mass0; } // McBAGB = M0 (Hurely et al. 2000, discussion just before eq 89) double CalculateCoreMassAtPhaseEnd() const { return m_CoreMass; } // NO-OP @@ -60,6 +62,7 @@ class HeHG: virtual public BaseStar, public HeMS { static double CalculateCoreMass_Luminosity_B_Static() { return 4.1E4; } static double CalculateCoreMass_Luminosity_D_Static(const double p_Mass) { return 5.5E4 / (1.0 + (0.4 * p_Mass * p_Mass * p_Mass * p_Mass)); } // pow() is slow - use multiplication + double CalculateConvectiveCoreRadius () const { return 5.0 * CalculateRemnantRadius (); } // Last paragraph of section 6 of Hurley+ 2000 double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const; double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 1.28; } // From BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details. @@ -82,7 +85,7 @@ class HeHG: virtual public BaseStar, public HeMS { double CalculatePerturbationMu() const; double CalculatePerturbationMuAtPhaseEnd() const { return m_Mu; } // NO-OP - double CalculateRadialExtentConvectiveEnvelope() const { return GiantBranch::CalculateRadialExtentConvectiveEnvelope(); } // Skip HeMS + double CalculateRadialExtentConvectiveEnvelope() const { return HG::CalculateRadialExtentConvectiveEnvelope(); } double CalculateRadiusAtPhaseEnd() const { return m_Radius; } // NO-OP double CalculateRadiusOnPhase() const; diff --git a/src/HeMS.h b/src/HeMS.h index 806fdeb20..53e7a09f3 100644 --- a/src/HeMS.h +++ b/src/HeMS.h @@ -58,7 +58,9 @@ class HeMS: virtual public BaseStar, public TPAGB { double CalculateCOCoreMassAtPhaseEnd() const { return CalculateCOCoreMassOnPhase(); } // Same as on phase double CalculateCOCoreMassOnPhase() const { return 0.0; } // McCO(HeMS) = 0.0 - double CalculateCoreMassAtPhaseEnd() const { return CalculateHeCoreMassOnPhase(); } // Same as on phase + double CalculateConvectiveCoreMass() const { return MainSequence::CalculateConvectiveCoreMass(); } // Temporary solution, until we have tested the rate at which the convective core recedes in HeMS stars + double CalculateConvectiveCoreRadius () const { return 0.5 * m_Radius; } // Temporary solution, until we have tested the core radii of HeMS stars + double CalculateCoreMassAtPhaseEnd() const { return CalculateHeCoreMassOnPhase(); } // Same as on phase /*ILYA*/ To fix, not everything will become CO core double CalculateCoreMassOnPhase() const { return 0.0; } // Mc(HeMS) = 0.0 double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const; @@ -88,7 +90,7 @@ class HeMS: virtual public BaseStar, public TPAGB { double CalculatePerturbationMu() const { return 5.0; } // Hurley et al. 2000, eqs 97 & 98 - double CalculateRadialExtentConvectiveEnvelope() const { return BaseStar::CalculateRadialExtentConvectiveEnvelope(); } // HeMS stars don't have a convective envelope + double CalculateRadialExtentConvectiveEnvelope() const { return 0.0; } // HeMS stars don't have a convective envelope double CalculateRadiusAtPhaseEnd(const double p_Mass) const { return CalculateRadiusAtPhaseEnd_Static(p_Mass); } double CalculateRadiusAtPhaseEnd() const { return CalculateRadiusAtPhaseEnd(m_Mass); } // Use class member variables diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index 58b4b3659..665c90c4a 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -470,11 +470,7 @@ double MainSequence::CalculateRadiusOnPhase(const double p_Mass, const double p_ /* * Calculate the radial extent of the star's convective envelope (if it has one) * - * Hurley et al. 2000, sec. 2.3, particularly subsec. 2.3.1, eqs 36-40 - * - * (Technically not a radius calculation I suppose, but "radial extent" is close enough to put it with the radius calculations...) - * - * JR: todo: original code for MS is broken for mass < 1.25 - check this (see calculateRadialExtentConvectiveEnvelope()) + * Hurley et al. 2002, sec. 2.3, particularly subsec. 2.3.1, eqs 36-38 * * * double CalculateRadialExtentConvectiveEnvelope() @@ -482,11 +478,50 @@ double MainSequence::CalculateRadiusOnPhase(const double p_Mass, const double p_ * @return Radial extent of the star's convective envelope in Rsol */ double MainSequence::CalculateRadialExtentConvectiveEnvelope() const { - return utils::Compare(m_Mass, 0.35) <= 0 ? m_Radius * std::sqrt(std::sqrt(1.0 - m_Tau)) : 0.0; + double radiusEnvelope0 = m_Radius; + if ( utils::Compare(m_Mass, 1.25) >= 0) + radiusEnvelope0 = 0.0; + else if (utils::Compare(m_Mass, 0.35) > 0) { + double radiusM035 = m_Radius; // /*ILYA*/ to fix: radius of a 0.35 solar mass star of fractional age Tau + radiusEnvelope0 = radiusM035 * std::sqrt((1.25 - m_Mass) / 0.9); + } + return radiusEnvelope0 * std::sqrt(std::sqrt(1.0 - m_Tau)); +} + +double MainSequence::CalculateConvectiveCoreRadius() const { + if (utils::Compare(m_Mass, 1.25) < 0) // /*ILYA*/ To check + return 0; + return ( m_Mass * (0.06 + 0.05 * exp(-m_Mass / 61.57))); // Preliminary fit from Minori } /////////////////////////////////////////////////////////////////////////////////////// +// // +// MASS CALCULATIONS // +// // +/////////////////////////////////////////////////////////////////////////////////////// + + +/* + * Calculate the mass of the convective core + * + * Based on Shikauchi, Hirai, Mandel (2024), core mass shrinks to 60% of initial value over the course of the MS + * + * + * double CalculateConvectiveCoreMass() const + * + * @return Mass of convective core in Msol + */ +double MainSequence::CalculateConvectiveCoreMass() const { + HG clone = *this; //create an HG star clone to query its core mass just after TAMS + double TAMSCoreMass = clone.CoreMass(); + double finalConvectiveCoreMass = TAMSCoreMass; + double initialConvectiveCoreMass = finalConvectiveCoreMass / 0.6; + return ( initialConvectiveCoreMass - m_Tau * (initialConvectiveCoreMass - finalConvectiveCoreMass) ); +} + + + // // // LIFETIME / AGE CALCULATIONS // // // diff --git a/src/MainSequence.h b/src/MainSequence.h index a5ce35fa2..c87539ee5 100644 --- a/src/MainSequence.h +++ b/src/MainSequence.h @@ -25,7 +25,9 @@ class MainSequence: virtual public BaseStar { double CalculateAlphaL(const double p_Mass) const; double CalculateAlphaR(const double p_Mass) const; - double CalculateConvectiveEnvelopeMass() const { return 0.0; } + double CalculateConvectiveCoreMass() const; + double CalculateConvectiveCoreRadius() const; + DBL_DBL CalculateConvectiveEnvelopeMass() const { return std::tuple (0.0, 0.0); } double CalculateBetaL(const double p_Mass) const; double CalculateBetaR(const double p_Mass) const; diff --git a/src/Makefile b/src/Makefile index 30a76774b..6221f5b55 100644 --- a/src/Makefile +++ b/src/Makefile @@ -9,16 +9,16 @@ CPP := g++ # gsl directories -GSLINCDIR := /include -GSLLIBDIR := /lib +GSLINCDIR := /opt/homebrew/Cellar/gsl/2.7.1/include +GSLLIBDIR := /opt/homebrew/Cellar/gsl/2.7.1/lib # boost directories -BOOSTINCDIR := /include -BOOSTLIBDIR := /lib +BOOSTINCDIR := /opt/homebrew/Cellar/boost/1.84.0_1/include +BOOSTLIBDIR := /opt/homebrew/Cellar/boost/1.84.0_1/lib # hdf5 directories -HDF5INCDIR := /usr/include/hdf5/serial -HDF5LIBDIR := /usr/lib/x86_64-linux-gnu/hdf5/serial +HDF5INCDIR := /opt/homebrew/Cellar/hdf5/1.14.3_1/include +HDF5LIBDIR := /opt/homebrew/Cellar/hdf5/1.14.3_1/lib EXE := COMPAS @@ -40,7 +40,7 @@ ifneq ($(filter staticfast,$(MAKECMDGOALS)),) endif -CXXFLAGS := -std=c++11 -Wall -Woverloaded-virtual $(OPTFLAGS) +CXXFLAGS := -std=c++11 -Wall $(OPTFLAGS) ICFLAGS := -I$(GSLINCDIR) -I$(BOOSTINCDIR) -I$(HDF5INCDIR) -I. LIBS := -lm -lz -ldl -lpthread @@ -121,7 +121,7 @@ $(EXE)_STATIC: $(OBJI) @echo $(OBJI) $(CPP) $(OBJI) $(LFLAGS) -static -o $@ -%.o: %.cpp +.cpp.o: $(SOURCES) $(INCL) Makefile $(CPP) $(CXXFLAGS) $(ICFLAGS) -c $? .phony: clean static fast staticfast diff --git a/src/Remnants.h b/src/Remnants.h index 284315988..589ffd579 100644 --- a/src/Remnants.h +++ b/src/Remnants.h @@ -1,114 +1,116 @@ -#ifndef __Remnants_h__ -#define __Remnants_h__ - -#include "constants.h" -#include "typedefs.h" -#include "profiling.h" -#include "utils.h" - -#include "HeGB.h" - - -class BaseStar; -class HeGB; - -class Remnants: virtual public BaseStar, public HeGB { - -public: - - Remnants(const BaseStar &p_BaseStar, const bool p_Initialise = true) : BaseStar(p_BaseStar), HeGB(p_BaseStar, false) { - if (p_Initialise) Initialise(); - } - - Remnants& operator = (const BaseStar &p_BaseStar) { static_cast(*this) = p_BaseStar; return *this; } - - - // member functions - -protected: - - - // member functions - alphabetically - double CalculateCOCoreMassOnPhase() const { return m_Mass; } // Return m_Mass - - double CalculateConvectiveEnvelopeMass() const { return 0.0; } - - double CalculateCoreMassOnPhase() const { return m_Mass; } // Return m_Mass - - double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate) { return 0.0; } // Should never be called... - - void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { GiantBranch::CalculateGBParams(p_Mass, p_GBParams); } // Default to GiantBranch - void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables - - double CalculateHeCoreMassOnPhase() const { return m_Mass; } // Return m_Mass - - double CalculateInitialSupernovaMass() const { return GiantBranch::CalculateInitialSupernovaMass(); } // Use GiantBranch - - double CalculateLambdaDewi() const { return BaseStar::CalculateLambdaDewi(); } // Not supported - use BaseStar - double CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) const { return BaseStar::CalculateLambdaNanjingStarTrack(0.0, 0.0); } // Not supported - use BaseStar (0.0 are dummy values) JR: todo: check this (type 10 not mentioned as not supported in original code) - - DBL_DBL CalculateMassAcceptanceRate(const double p_DonorMassRate, - const double p_AccretorMassRate); - DBL_DBL CalculateMassAcceptanceRate(const double p_DonorMassRate, - const double p_AccretorMassRate, - const bool p_IsHeRich) { return CalculateMassAcceptanceRate(p_DonorMassRate, p_AccretorMassRate); } // Ignore the He content for non-WDs - - double CalculateMassLossRateHurley() { return 0.0; } - double CalculateMassLossRateBelczynski2010() { return 0.0; } - - double CalculatePerturbationMuOnPhase() const { return m_Mu; } // NO-OP - - double CalculateRadialExtentConvectiveEnvelope() const { return BaseStar::CalculateRadialExtentConvectiveEnvelope(); } // WD stars don't have a convective envelope - - std::tuple CalculateRadiusAndStellarTypeOnPhase() const { return BaseStar::CalculateRadiusAndStellarTypeOnPhase(); } - - double CalculateTauOnPhase() const { return m_Tau; } // NO-OP - - double CalculateThermalTimescale(const double p_Radius = 1.0) const { return CalculateDynamicalTimescale(); } // Parameter is ignored - double CalculateThermalTimescale() const { return CalculateThermalTimescale(m_Radius); } // Use inheritance hierarchy - - double CalculateThermalMassLossRate() const { return BaseStar::CalculateThermalMassLossRate(); } // Set thermal mass gain rate to be effectively infinite, using dynamical timescale (in practice, will be Eddington limited), avoid division by zero - - void CalculateTimescales(const double p_Mass, DBL_VECTOR &p_Timescales) { return TPAGB::CalculateTimescales(p_Mass, p_Timescales); } // Use TPAGB - void CalculateTimescales() { CalculateTimescales(m_Mass0, m_Timescales); } // Use class member variables - - double CalculateZetaConstantsByEnvelope(ZETA_PRESCRIPTION p_ZetaPrescription) { return 0.0; } // Should never be called... - - double ChooseTimestep(const double p_Time) const; - - ENVELOPE DetermineEnvelopeType() const { return ENVELOPE::REMNANT; } // Always REMNANT - - void EvolveOneTimestepPreamble() { BaseStar::EvolveOneTimestepPreamble(); } // Default to BaseStar - - STELLAR_TYPE EvolveToNextPhase() { return BaseStar::EvolveToNextPhase(); } // Default to BaseStar - - bool IsDegenerate() const { return true; } // White Dwarfs, NS and BH are degenerate - - bool IsEndOfPhase() const { return !ShouldEvolveOnPhase(); } // Phase ends when envelope loss or going supernova - - bool IsSupernova() const { return false; } // Default - - void PerturbLuminosityAndRadius() { } // NO-OP - - STELLAR_TYPE ResolveEnvelopeLoss(bool p_NoCheck = false) { return BaseStar::ResolveEnvelopeLoss(p_NoCheck); } // Default to BaseStar - - void ResolveEnvelopeMassAtPhaseEnd(const double p_Tau) const { ResolveEnvelopeMassOnPhase(p_Tau); } // Same as on phase - void ResolveEnvelopeMassOnPhase(const double p_Tau) const { } // NO-OP - - void ResolveMassLoss(const bool p_UpdateMDt = true) { } // NO-OP - - STELLAR_TYPE ResolveSkippedPhase() { return BaseStar::ResolveSkippedPhase(); } // Default to BaseStar - // - void ResolveShellChange(const double p_AccretedMass) { } // NO-OP - // - STELLAR_TYPE ResolveSupernova() { return BaseStar::ResolveSupernova(); } // Default to BaseStar - - void SetPulsarParameters() const { } // NO-OP - - bool ShouldEnvelopeBeExpelledByPulsations() const { return false; } // No envelope to lose by pulsations - bool ShouldEvolveOnPhase() const { return true; } // Default - bool ShouldSkipPhase() const { return false; } // Don't skip WD phase - -}; - -#endif // __Remnants_h__ +#ifndef __Remnants_h__ +#define __Remnants_h__ + +#include "constants.h" +#include "typedefs.h" +#include "profiling.h" +#include "utils.h" + +#include "HeGB.h" + + +class BaseStar; +class HeGB; + +class Remnants: virtual public BaseStar, public HeGB { + +public: + + Remnants(const BaseStar &p_BaseStar, const bool p_Initialise = true) : BaseStar(p_BaseStar), HeGB(p_BaseStar, false) { + if (p_Initialise) Initialise(); + } + + Remnants& operator = (const BaseStar &p_BaseStar) { static_cast(*this) = p_BaseStar; return *this; } + + + // member functions + +protected: + + + // member functions - alphabetically + double CalculateCOCoreMassOnPhase() const { return m_Mass; } // Return m_Mass + + double CalculateConvectiveCoreRadius () const { return m_Radius; } // All core + DBL_DBL CalculateConvectiveEnvelopeMass() const { return std::tuple (0.0, +0.0); } + + double CalculateCoreMassOnPhase() const { return m_Mass; } // Return m_Mass + + double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate) { return 0.0; } // Should never be called... + + void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { GiantBranch::CalculateGBParams(p_Mass, p_GBParams); } // Default to GiantBranch + void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables + + double CalculateHeCoreMassOnPhase() const { return m_Mass; } // Return m_Mass + + double CalculateInitialSupernovaMass() const { return GiantBranch::CalculateInitialSupernovaMass(); } // Use GiantBranch + + double CalculateLambdaDewi() const { return BaseStar::CalculateLambdaDewi(); } // Not supported - use BaseStar + double CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) const { return BaseStar::CalculateLambdaNanjingStarTrack(0.0, 0.0); } // Not supported - use BaseStar (0.0 are dummy values) JR: todo: check this (type 10 not mentioned as not supported in original code) + + DBL_DBL CalculateMassAcceptanceRate(const double p_DonorMassRate, + const double p_AccretorMassRate); + DBL_DBL CalculateMassAcceptanceRate(const double p_DonorMassRate, + const double p_AccretorMassRate, + const bool p_IsHeRich) { return CalculateMassAcceptanceRate(p_DonorMassRate, p_AccretorMassRate); } // Ignore the He content for non-WDs + + double CalculateMassLossRateHurley() { return 0.0; } + double CalculateMassLossRateBelczynski2010() { return 0.0; } + + double CalculatePerturbationMuOnPhase() const { return m_Mu; } // NO-OP + + double CalculateRadialExtentConvectiveEnvelope() const { return 0.0; } // WD stars don't have a convective envelope + + std::tuple CalculateRadiusAndStellarTypeOnPhase() const { return BaseStar::CalculateRadiusAndStellarTypeOnPhase(); } + + double CalculateTauOnPhase() const { return m_Tau; } // NO-OP + + double CalculateThermalTimescale(const double p_Radius = 1.0) const { return CalculateDynamicalTimescale(); } // Parameter is ignored + double CalculateThermalTimescale() const { return CalculateThermalTimescale(m_Radius); } // Use inheritance hierarchy + + double CalculateThermalMassLossRate() const { return BaseStar::CalculateThermalMassLossRate(); } // Set thermal mass gain rate to be effectively infinite, using dynamical timescale (in practice, will be Eddington limited), avoid division by zero + + void CalculateTimescales(const double p_Mass, DBL_VECTOR &p_Timescales) { return TPAGB::CalculateTimescales(p_Mass, p_Timescales); } // Use TPAGB + void CalculateTimescales() { CalculateTimescales(m_Mass0, m_Timescales); } // Use class member variables + + double CalculateZetaConstantsByEnvelope(ZETA_PRESCRIPTION p_ZetaPrescription) { return 0.0; } // Should never be called... + + double ChooseTimestep(const double p_Time) const; + + ENVELOPE DetermineEnvelopeType() const { return ENVELOPE::REMNANT; } // Always REMNANT + + void EvolveOneTimestepPreamble() { BaseStar::EvolveOneTimestepPreamble(); } // Default to BaseStar + + STELLAR_TYPE EvolveToNextPhase() { return BaseStar::EvolveToNextPhase(); } // Default to BaseStar + + bool IsDegenerate() const { return true; } // White Dwarfs, NS and BH are degenerate + + bool IsEndOfPhase() const { return !ShouldEvolveOnPhase(); } // Phase ends when envelope loss or going supernova + + bool IsSupernova() const { return false; } // Default + + void PerturbLuminosityAndRadius() { } // NO-OP + + STELLAR_TYPE ResolveEnvelopeLoss(bool p_NoCheck = false) { return BaseStar::ResolveEnvelopeLoss(p_NoCheck); } // Default to BaseStar + + void ResolveEnvelopeMassAtPhaseEnd(const double p_Tau) const { ResolveEnvelopeMassOnPhase(p_Tau); } // Same as on phase + void ResolveEnvelopeMassOnPhase(const double p_Tau) const { } // NO-OP + + void ResolveMassLoss(const bool p_UpdateMDt = true) { } // NO-OP + + STELLAR_TYPE ResolveSkippedPhase() { return BaseStar::ResolveSkippedPhase(); } // Default to BaseStar + // + void ResolveShellChange(const double p_AccretedMass) { } // NO-OP + // + STELLAR_TYPE ResolveSupernova() { return BaseStar::ResolveSupernova(); } // Default to BaseStar + + void SetPulsarParameters() const { } // NO-OP + + bool ShouldEnvelopeBeExpelledByPulsations() const { return false; } // No envelope to lose by pulsations + bool ShouldEvolveOnPhase() const { return true; } // Default + bool ShouldSkipPhase() const { return false; } // Don't skip WD phase + +}; + +#endif // __Remnants_h__ diff --git a/src/Star.h b/src/Star.h index 1670b807a..2474acf02 100755 --- a/src/Star.h +++ b/src/Star.h @@ -154,11 +154,15 @@ class Star { const double p_EnvMass, const double p_Radius) { m_Star->CalculateBindingEnergies(p_CoreMass, p_EnvMass, p_Radius); } - double CalculateConvectiveEnvelopeBindingEnergy(const double p_CoreMass, + double CalculateConvectiveCoreMass () { return m_Star->CalculateConvectiveCoreMass(); } + double CalculateConvectiveCoreRadius () { return m_Star->CalculateConvectiveCoreRadius(); } + + double CalculateConvectiveEnvelopeBindingEnergy(const double p_TotalMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, - const double p_Lambda) { return m_Star->CalculateConvectiveEnvelopeBindingEnergy(p_CoreMass, p_ConvectiveEnvelopeMass, p_Radius, p_Lambda); } - double CalculateConvectiveEnvelopeMass() { return m_Star->CalculateConvectiveEnvelopeMass(); } + const double p_Lambda) { return m_Star->CalculateConvectiveEnvelopeBindingEnergy(p_TotalMass, p_ConvectiveEnvelopeMass, p_Radius, p_Lambda); } + double CalculateConvectiveEnvelopeLambdaPicker( double p_convectiveEnvelopeMass, double p_maxConvectiveEnvelopeMass ) { return m_Star->CalculateConvectiveEnvelopeLambdaPicker(p_convectiveEnvelopeMass, p_maxConvectiveEnvelopeMass); } + DBL_DBL CalculateConvectiveEnvelopeMass() { return m_Star->CalculateConvectiveEnvelopeMass(); } double CalculateEddyTurnoverTimescale() { return m_Star->CalculateEddyTurnoverTimescale(); } @@ -173,6 +177,8 @@ class Star { double CalculateMomentOfInertia() const { return m_Star->CalculateMomentOfInertia(); } double CalculateMomentOfInertiaAU() const { return m_Star->CalculateMomentOfInertiaAU(); } + + double CalculateRadialExtentConvectiveEnvelope() { return m_Star->CalculateRadialExtentConvectiveEnvelope(); } void CalculateSNAnomalies(const double p_Eccentricity) { m_Star->CalculateSNAnomalies(p_Eccentricity); } diff --git a/src/TPAGB.h b/src/TPAGB.h index 06cd85cc3..51f328984 100755 --- a/src/TPAGB.h +++ b/src/TPAGB.h @@ -40,6 +40,8 @@ class TPAGB: virtual public BaseStar, public EAGB { double CalculateCOCoreMassAtPhaseEnd() const { return (utils::Compare(m_COCoreMass, m_GBParams[static_cast(GBP::McSN)]) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0) ? m_COCoreMass : m_Mass; } double CalculateCOCoreMassOnPhase() const { return CalculateCoreMassOnPhase(m_Mass0, m_Age); } // McCO(TPAGB) = Mc(TPAGB)Same as on phase + double CalculateConvectiveCoreRadius () const { return 5.0 * CalculateRemnantRadius (); } // Last paragraph of section 6 of Hurley+ 2000 + double CalculateCoreMassAtPhaseEnd() const { return m_CoreMass; } // NO-OP double CalculateCoreMassOnPhase(const double p_Mass, const double p_Time) const; double CalculateCoreMassOnPhase() const { return CalculateCoreMassOnPhase(m_Mass0, m_Age); } // Use class member variables From 2d139d997e43ca8d22fba55eaf99750571d6378b Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Mon, 15 Apr 2024 15:35:32 +1000 Subject: [PATCH 02/10] Removing debugging print message --- src/GiantBranch.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index f896a8984..9f36bbc8b 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -1050,11 +1050,11 @@ DBL_DBL GiantBranch::CalculateConvectiveEnvelopeMass() const { double MinterfMcoref = -0.021 * log10Z + 0.0038; //Eq. (8) of Picker+ 2024 double Tonset = -139.8 * log10Z * log10Z - 981.7 * log10Z + 2798.3; //Eq. (6) of Picker+ 2024 EAGB clone = *this; // Create an HG star clone to query its core mass just after BAGB - double Tmin=this->Temperature(); + clone.UpdateAttributesAndAgeOneTimestep(0.0, 0.0, 0.0, true); // Otherwise, temperature not updated + double Tmin=clone.Temperature(); double Mcorefinal = CalculateCoreMassAtBAGB(m_Mass); - std::cout<<"temp"<CoreMass()< (convectiveEnvelopeMass, Mconvmax); } From f9ee31fa16fd601211ca498d958d758715d4d410 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Mon, 15 Apr 2024 17:04:13 +1000 Subject: [PATCH 03/10] Updated changelog, fixed typo in comments --- src/BaseStar.cpp | 2 +- src/changelog.h | 9 +++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index f4fdab1f6..a790c8287 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -3827,7 +3827,7 @@ double BaseStar::CalculateConvectiveEnvelopeBindingEnergy(const double p_TotalMa * double BaseStar::CalculateConvectiveEnvelopeLambdaPicker(double p_convectiveEnvelopeMass, double p_maxConvectiveEnvelopeMass) * * @param [IN] p_convectiveEnvelopeMass Mass of the outer convective envelope shell - * @param [IN] p_maxConvectiveEnvelopeMass Maximum mass of the outher convective envelope shell at the onset of carbon burning + * @param [IN] p_maxConvectiveEnvelopeMass Maximum mass of the outer convective envelope shell at the onset of carbon burning * @return Lambda binding energy parameter for the outer convective envelope */ double BaseStar::CalculateConvectiveEnvelopeLambdaPicker(double p_convectiveEnvelopeMass, double p_maxConvectiveEnvelopeMass) { diff --git a/src/changelog.h b/src/changelog.h index 4334b5e37..55d0dc5fc 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1117,8 +1117,13 @@ // - Added Hirai pulsar rocket kick, and related options // 02.43.01 SS - Apr 8, 2024 - Defect repair // - Fix CalculateMassLossRateBjorklundEddingtonFactor to use LSOLW (in SI) rather than LSOL (in cgs) -// +// 02.43.02 IM - Apr 15, 2024 - Enhancement +// - Updated fits for the mass and binding energy of the outer convective envelope based on Picker, Hirai, Mandel (2024) +// - Added functionality for CalculateConvectiveEnvelopeMass(), CalculateConvectiveCoreMass(), CalculateConvectiveCoreRadius() +// - Defect repair +// - Fixes to CalculateRadialExtentConvectiveEnvelope(), comments + -const std::string VERSION_STRING = "02.43.01"; +const std::string VERSION_STRING = "02.43.02"; # endif // __changelog_h__ From 1f0f43bfab51aa2fea65f1cfe29d8e6cb2a13e7f Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Wed, 17 Apr 2024 11:47:26 +1000 Subject: [PATCH 04/10] Updated envelope masses for low-mass MS stars --- src/GiantBranch.cpp | 2 +- src/MainSequence.cpp | 18 ++++++++++++++++++ src/MainSequence.h | 2 +- 3 files changed, 20 insertions(+), 2 deletions(-) diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index 9f36bbc8b..65bcad7cf 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -1057,7 +1057,7 @@ DBL_DBL GiantBranch::CalculateConvectiveEnvelopeMass() const { double convectiveEnvelopeMass = Mconvmax / (1.0 + exp(4.6 * (Tmin + Tonset - 2.0 * m_Temperature)/(Tmin-Tonset))); //Eq. (7) of Picker+ 2024 return std::tuple (convectiveEnvelopeMass, Mconvmax); -} +} // /*ILYA*/ check consistency with HG convective envelope radii and masses from Hurley+ 2002, 2000 /////////////////////////////////////////////////////////////////////////////////////// diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index 665c90c4a..16624bd1b 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -520,6 +520,24 @@ double MainSequence::CalculateConvectiveCoreMass() const { return ( initialConvectiveCoreMass - m_Tau * (initialConvectiveCoreMass - finalConvectiveCoreMass) ); } +/* + * Calculate the mass of the convective envelope + * + * Based on section 7.2 (after Eq. 111) of Hurley, Pols, Tout (2000) + * + * + * double CalculateConvectiveEnvelopeMass() const + * + * @return Mass of convective envelope in Msol + */ +double MainSequence::CalculateConvectiveEnvelopeMass() const { + if (utils::Compare(m_Mass, 1.25) > 0) + return 0; + double massEnvelope0 = 0.35; + if(utils::Compare(m_Mass, 0.35) > 0) + massEnvelope0 = 0.35 * (1.25 - m_Mass) * (1.25 - m_Mass) / 0.81; + return massEnvelope0 * sqrt(sqrt(1.0 - m_Tau)); +} // // diff --git a/src/MainSequence.h b/src/MainSequence.h index c87539ee5..81a650cd1 100644 --- a/src/MainSequence.h +++ b/src/MainSequence.h @@ -27,7 +27,7 @@ class MainSequence: virtual public BaseStar { double CalculateConvectiveCoreMass() const; double CalculateConvectiveCoreRadius() const; - DBL_DBL CalculateConvectiveEnvelopeMass() const { return std::tuple (0.0, 0.0); } + DBL_DBL CalculateConvectiveEnvelopeMass() const; double CalculateBetaL(const double p_Mass) const; double CalculateBetaR(const double p_Mass) const; From adb3ccddd81cd1ab6f9b3cee510b0ae87776ebab Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 18 Apr 2024 11:01:08 +1000 Subject: [PATCH 05/10] Typo in MainSequence::CalculateConvectiveEnvelopeMass() return type --- src/MainSequence.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index 16624bd1b..0acbc3931 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -530,13 +530,14 @@ double MainSequence::CalculateConvectiveCoreMass() const { * * @return Mass of convective envelope in Msol */ -double MainSequence::CalculateConvectiveEnvelopeMass() const { +DBL_DBL MainSequence::CalculateConvectiveEnvelopeMass() const { if (utils::Compare(m_Mass, 1.25) > 0) - return 0; + return std::tuple (0.0, 0.0); double massEnvelope0 = 0.35; if(utils::Compare(m_Mass, 0.35) > 0) massEnvelope0 = 0.35 * (1.25 - m_Mass) * (1.25 - m_Mass) / 0.81; - return massEnvelope0 * sqrt(sqrt(1.0 - m_Tau)); + double massEnvelope = massEnvelope0 * sqrt(sqrt(1.0 - m_Tau)); + return std::tuple (massEnvelope, massEnvelope); } From 3dc90a674eb26e4610c473527601a533211a5495 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 18 Apr 2024 14:06:23 +1000 Subject: [PATCH 06/10] Changelog conflict resolution --- src/changelog.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/changelog.h b/src/changelog.h index 55d0dc5fc..bdfca3dc9 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1117,13 +1117,15 @@ // - Added Hirai pulsar rocket kick, and related options // 02.43.01 SS - Apr 8, 2024 - Defect repair // - Fix CalculateMassLossRateBjorklundEddingtonFactor to use LSOLW (in SI) rather than LSOL (in cgs) -// 02.43.02 IM - Apr 15, 2024 - Enhancement +// 02.43.02 JR - Apr 15, 2024 - Defect repair +// - Fix for issue #1074 - SSE Supernova records duplicated +// 02.43.03 IM - Apr 15, 2024 - Enhancement // - Updated fits for the mass and binding energy of the outer convective envelope based on Picker, Hirai, Mandel (2024) // - Added functionality for CalculateConvectiveEnvelopeMass(), CalculateConvectiveCoreMass(), CalculateConvectiveCoreRadius() // - Defect repair // - Fixes to CalculateRadialExtentConvectiveEnvelope(), comments -const std::string VERSION_STRING = "02.43.02"; +const std::string VERSION_STRING = "02.43.03"; # endif // __changelog_h__ From 1d1631d8c0489f4f53499bb460405bd6eedf17f9 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 18 Apr 2024 19:42:40 +1000 Subject: [PATCH 07/10] Makefile revert --- src/Makefile | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/Makefile b/src/Makefile index 6221f5b55..30a76774b 100644 --- a/src/Makefile +++ b/src/Makefile @@ -9,16 +9,16 @@ CPP := g++ # gsl directories -GSLINCDIR := /opt/homebrew/Cellar/gsl/2.7.1/include -GSLLIBDIR := /opt/homebrew/Cellar/gsl/2.7.1/lib +GSLINCDIR := /include +GSLLIBDIR := /lib # boost directories -BOOSTINCDIR := /opt/homebrew/Cellar/boost/1.84.0_1/include -BOOSTLIBDIR := /opt/homebrew/Cellar/boost/1.84.0_1/lib +BOOSTINCDIR := /include +BOOSTLIBDIR := /lib # hdf5 directories -HDF5INCDIR := /opt/homebrew/Cellar/hdf5/1.14.3_1/include -HDF5LIBDIR := /opt/homebrew/Cellar/hdf5/1.14.3_1/lib +HDF5INCDIR := /usr/include/hdf5/serial +HDF5LIBDIR := /usr/lib/x86_64-linux-gnu/hdf5/serial EXE := COMPAS @@ -40,7 +40,7 @@ ifneq ($(filter staticfast,$(MAKECMDGOALS)),) endif -CXXFLAGS := -std=c++11 -Wall $(OPTFLAGS) +CXXFLAGS := -std=c++11 -Wall -Woverloaded-virtual $(OPTFLAGS) ICFLAGS := -I$(GSLINCDIR) -I$(BOOSTINCDIR) -I$(HDF5INCDIR) -I. LIBS := -lm -lz -ldl -lpthread @@ -121,7 +121,7 @@ $(EXE)_STATIC: $(OBJI) @echo $(OBJI) $(CPP) $(OBJI) $(LFLAGS) -static -o $@ -.cpp.o: $(SOURCES) $(INCL) Makefile +%.o: %.cpp $(CPP) $(CXXFLAGS) $(ICFLAGS) -c $? .phony: clean static fast staticfast From 75bf3ee3b2cfe0d03623f3233121df481c285e14 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 18 Apr 2024 20:40:17 +1000 Subject: [PATCH 08/10] Better approximation for radial extent of low-mass convective envelopes --- src/MainSequence.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index 0acbc3931..540027a88 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -482,7 +482,7 @@ double MainSequence::CalculateRadialExtentConvectiveEnvelope() const { if ( utils::Compare(m_Mass, 1.25) >= 0) radiusEnvelope0 = 0.0; else if (utils::Compare(m_Mass, 0.35) > 0) { - double radiusM035 = m_Radius; // /*ILYA*/ to fix: radius of a 0.35 solar mass star of fractional age Tau + double radiusM035 = CalculateRadiusAtZAMS(0.35); // uses radius of a 0.35 solar mass star at ZAMS rather than at fractional age Tau, but such low-mass stars only grow by a maximum factor of 1.5 [just above Eq. (10) in Hurley, Pols, Tout (2000), so this is a reasonable approximation radiusEnvelope0 = radiusM035 * std::sqrt((1.25 - m_Mass) / 0.9); } return radiusEnvelope0 * std::sqrt(std::sqrt(1.0 - m_Tau)); From f5b43655aaf7906c6e32f74e153d28ef1c11d208 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Fri, 19 Apr 2024 08:52:18 +1000 Subject: [PATCH 09/10] Minor fixes in response to reviewer comments --- src/BaseBinaryStar.cpp | 2 +- src/BaseStar.cpp | 7 +++---- src/BaseStar.h | 2 +- src/FGB.h | 4 ++-- src/GiantBranch.cpp | 16 +--------------- src/GiantBranch.h | 3 +-- src/HG.cpp | 15 --------------- src/HG.h | 2 +- src/MainSequence.cpp | 4 ++-- src/Star.h | 2 +- 10 files changed, 13 insertions(+), 44 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index ce770f81e..ed6bde5bb 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1599,7 +1599,7 @@ void BaseBinaryStar::ResolveCommonEnvelopeEvent() { double k1 = m_Star1->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (lambda1 * alphaCE)) * m_Star1->Mass() * convectiveEnvelopeMass1 / m_Star1->Radius(); double k2 = m_Star2->IsOneOf(COMPACT_OBJECTS) ? 0.0 : (2.0 / (lambda2 * alphaCE)) * m_Star2->Mass() * convectiveEnvelopeMass2 / m_Star2->Radius(); double k3 = m_Star1->Mass() * m_Star2->Mass() / periastronRsol; // assumes immediate circularisation at periastron at start of CE - double k4 = (endOfFirstStageMass1 * endOfFirstStageMass2); + double k4 = endOfFirstStageMass1 * endOfFirstStageMass2; double aFinalRsol = k4 / (k1 + k2 + k3); m_SemiMajorAxis = aFinalRsol * RSOL_TO_AU; diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index bb3da744b..194cc6c76 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -3825,16 +3825,15 @@ double BaseStar::CalculateConvectiveEnvelopeBindingEnergy(const double p_TotalMa * Follows the fits of Picker, Hirai, Mandel (2024), arXiv:2402.13180 for lambda_He * * - * double BaseStar::CalculateConvectiveEnvelopeLambdaPicker(double p_convectiveEnvelopeMass, double p_maxConvectiveEnvelopeMass) + * double BaseStar::CalculateConvectiveEnvelopeLambdaPicker(const double p_convectiveEnvelopeMass, const double p_maxConvectiveEnvelopeMass) const * * @param [IN] p_convectiveEnvelopeMass Mass of the outer convective envelope shell * @param [IN] p_maxConvectiveEnvelopeMass Maximum mass of the outer convective envelope shell at the onset of carbon burning * @return Lambda binding energy parameter for the outer convective envelope */ -double BaseStar::CalculateConvectiveEnvelopeLambdaPicker(double p_convectiveEnvelopeMass, double p_maxConvectiveEnvelopeMass) { +double BaseStar::CalculateConvectiveEnvelopeLambdaPicker(const double p_convectiveEnvelopeMass, const double p_maxConvectiveEnvelopeMass) const { - double log10Z = log10 (m_Metallicity); - double m2 = 0.0023 * log10Z * log10Z + 0.0088 * log10Z + 0.013; // Eq. (12) and Table 1 of Picker, Hirai, Mandel (2024) + double m2 = 0.0023 * m_Log10Metallicity * m_Log10Metallicity + 0.0088 * m_Log10Metallicity + 0.013; // Eq. (12) and Table 1 of Picker, Hirai, Mandel (2024) double b1 = m2 * m_Mass - 0.23; // Eq. (11) of Picker+ (2024) double logLambda = p_convectiveEnvelopeMass / p_maxConvectiveEnvelopeMass > 0.3 ? 0.42 * p_convectiveEnvelopeMass / p_maxConvectiveEnvelopeMass + b1 : 0.3 * 0.42 + b1; diff --git a/src/BaseStar.h b/src/BaseStar.h index f73063b2b..5f2b73b03 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -169,7 +169,7 @@ class BaseStar { double CalculateConvectiveEnvelopeBindingEnergy(const double p_TotalMass, const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_Lambda); - double CalculateConvectiveEnvelopeLambdaPicker(double p_convectiveEnvelopeMass, double p_maxConvectiveEnvelopeMass); + double CalculateConvectiveEnvelopeLambdaPicker(const double p_convectiveEnvelopeMass, const double p_maxConvectiveEnvelopeMass) const; virtual DBL_DBL CalculateConvectiveEnvelopeMass() const { return std::tuple (0.0, 0.0); } virtual double CalculateCriticalMassRatio(const bool p_AccretorIsDegenerate); diff --git a/src/FGB.h b/src/FGB.h index dcec8efd7..55e7134c9 100755 --- a/src/FGB.h +++ b/src/FGB.h @@ -56,9 +56,9 @@ class FGB: virtual public BaseStar, public HG { double CalculateLuminosityOnPhase(const double p_Time) const; double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase(m_Age); } // Use class member variables - double CalculateRadialExtentConvectiveEnvelope() const { return Radius() - CalculateConvectiveCoreRadius(); } // Skip HG + double CalculateRadialExtentConvectiveEnvelope() const { return m_Radius - CalculateConvectiveCoreRadius(); } // Skip HG - double CalculateRadiusAtPhaseEnd(const double p_Mass, const double p_Luminosity) const { return GiantBranch::CalculateRadiusOnPhase(p_Mass, p_Luminosity); } // Skip HG - same as on phase + double CalculateRadiusAtPhaseEnd(const double p_Mass, const double p_Luminosity) const { return GiantBranch::CalculateRadiusOnPhase(p_Mass, p_Luminosity); } // Skip HG - same as on phase double CalculateRadiusAtPhaseEnd() const { return CalculateRadiusAtPhaseEnd(m_Mass, m_Luminosity); } // Use class member variables double CalculateRadiusOnPhase(const double p_Mass, const double p_Luminosity) const { return GiantBranch::CalculateRadiusOnPhase(p_Mass, p_Luminosity); } // Skip HG double CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase(m_Mass, m_Luminosity); } // Use class member variables diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index 1b6b6d017..deb31f216 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -680,20 +680,6 @@ double GiantBranch::CalculateRemnantRadius() const { } -/* - * Calculate the radial extent of the star's convective envelope (if it has one) - * - * Hurley et al. 2002, sec. 2.3, particularly subsec. 2.3.1, eqs 36-40 - * - * double CalculateRadialExtentConvectiveEnvelope() - * - * @return Radial extent of the star's convective envelope in Rsol - */ -double GiantBranch::CalculateRadialExtentConvectiveEnvelope() const { - return m_Radius - CalculateConvectiveCoreRadius(); -} - - /////////////////////////////////////////////////////////////////////////////////////// // // // MASS CALCULATIONS // @@ -1053,7 +1039,7 @@ DBL_DBL GiantBranch::CalculateConvectiveEnvelopeMass() const { clone.UpdateAttributesAndAgeOneTimestep(0.0, 0.0, 0.0, true); // Otherwise, temperature not updated double Tmin=clone.Temperature(); double Mcorefinal = CalculateCoreMassAtBAGB(m_Mass); - double Mconvmax = std::max(m_Mass - Mcorefinal * (1 + MinterfMcoref), 0.0); //Eq. (9) of Picker+ 2024 + double Mconvmax = std::max(m_Mass - Mcorefinal * (1.0 + MinterfMcoref), 0.0); //Eq. (9) of Picker+ 2024 double convectiveEnvelopeMass = Mconvmax / (1.0 + exp(4.6 * (Tmin + Tonset - 2.0 * m_Temperature)/(Tmin-Tonset))); //Eq. (7) of Picker+ 2024 return std::tuple (convectiveEnvelopeMass, Mconvmax); diff --git a/src/GiantBranch.h b/src/GiantBranch.h index e1f742479..e8744be00 100755 --- a/src/GiantBranch.h +++ b/src/GiantBranch.h @@ -26,7 +26,6 @@ class GiantBranch: virtual public BaseStar, public MainSequence { // member functions - alphabetically (sort of - some are grouped by functionality) double CalculateConvectiveCoreMass() const { return m_CoreMass; } - double CalculateConvectiveEnvelopeLambdaPicker(double p_convectiveEnvelopeMass, double p_maxConvectiveEnvelopeMass); DBL_DBL CalculateConvectiveEnvelopeMass() const; static double CalculateCoreMassAt2ndDredgeUp_Static(const double p_McBAGB); double CalculateCoreMassAtBAGB(const double p_Mass) const; @@ -92,7 +91,7 @@ class GiantBranch: virtual public BaseStar, public MainSequence { double CalculatePerturbationMu() const; - double CalculateRadialExtentConvectiveEnvelope() const; + double CalculateRadialExtentConvectiveEnvelope() const { return (m_Radius - CalculateConvectiveCoreRadius()); } //Hurley et al. 2002, sec. 2.3, particularly subsec. 2.3.1, eqs 36-40 double CalculateRadiusAtHeIgnition(const double p_Mass) const; double CalculateRadiusOnPhase(const double p_Mass, const double p_Luminosity) const { return CalculateRadiusOnPhase_Static(p_Mass, p_Luminosity, m_BnCoefficients); } diff --git a/src/HG.cpp b/src/HG.cpp index 60d88b8f6..9602e235c 100755 --- a/src/HG.cpp +++ b/src/HG.cpp @@ -793,21 +793,6 @@ double HG::CalculateRadiusOnPhase(const double p_Mass, const double p_Tau, const } -/* - * Calculate the radial extent of the star's convective envelope (if it has one) - * - * Hurley et al. 2002, sec. 2.3, particularly subsec. 2.3.1, eqs 39-40 - * - * - * double CalculateRadialExtentConvectiveEnvelope() - * - * @return Radial extent of the star's convective envelope in Rsol - */ -double HG::CalculateRadialExtentConvectiveEnvelope() const { - return std::sqrt(m_Tau) * (m_Radius - CalculateConvectiveCoreRadius()); -} - - /////////////////////////////////////////////////////////////////////////////////////// // // // MASS CALCULATIONS // diff --git a/src/HG.h b/src/HG.h index d631284f2..45c1e978d 100755 --- a/src/HG.h +++ b/src/HG.h @@ -92,7 +92,7 @@ class HG: virtual public BaseStar, public GiantBranch { double CalculateMassTransferRejuvenationFactor() const; - double CalculateRadialExtentConvectiveEnvelope() const; + double CalculateRadialExtentConvectiveEnvelope() const { return (std::sqrt(m_Tau) * (m_Radius - CalculateConvectiveCoreRadius())); } // Hurley et al. 2002, sec. 2.3, particularly subsec. 2.3.1, eqs 39-40 double CalculateRadiusAtPhaseEnd(const double p_Mass) const; double CalculateRadiusAtPhaseEnd() const { return CalculateRadiusAtPhaseEnd(m_Mass); } // Use class member variables diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index 540027a88..d369b31ab 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -490,8 +490,8 @@ double MainSequence::CalculateRadialExtentConvectiveEnvelope() const { double MainSequence::CalculateConvectiveCoreRadius() const { if (utils::Compare(m_Mass, 1.25) < 0) // /*ILYA*/ To check - return 0; - return ( m_Mass * (0.06 + 0.05 * exp(-m_Mass / 61.57))); // Preliminary fit from Minori + return 0.0; + return ( m_Mass * (0.06 + 0.05 * exp(-m_Mass / 61.57))); // Preliminary fit from Minori Shikauchi @ ZAMS, does not take evolution into account yet } diff --git a/src/Star.h b/src/Star.h index 2474acf02..32f636735 100755 --- a/src/Star.h +++ b/src/Star.h @@ -161,7 +161,7 @@ class Star { const double p_ConvectiveEnvelopeMass, const double p_Radius, const double p_Lambda) { return m_Star->CalculateConvectiveEnvelopeBindingEnergy(p_TotalMass, p_ConvectiveEnvelopeMass, p_Radius, p_Lambda); } - double CalculateConvectiveEnvelopeLambdaPicker( double p_convectiveEnvelopeMass, double p_maxConvectiveEnvelopeMass ) { return m_Star->CalculateConvectiveEnvelopeLambdaPicker(p_convectiveEnvelopeMass, p_maxConvectiveEnvelopeMass); } + double CalculateConvectiveEnvelopeLambdaPicker(const double p_convectiveEnvelopeMass, const double p_maxConvectiveEnvelopeMass ) const { return m_Star->CalculateConvectiveEnvelopeLambdaPicker(p_convectiveEnvelopeMass, p_maxConvectiveEnvelopeMass); } DBL_DBL CalculateConvectiveEnvelopeMass() { return m_Star->CalculateConvectiveEnvelopeMass(); } double CalculateEddyTurnoverTimescale() { return m_Star->CalculateEddyTurnoverTimescale(); } From b94f7aa86c510c76eb584ef5c271ab9d4625598a Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Fri, 19 Apr 2024 13:01:39 +1000 Subject: [PATCH 10/10] Fixes in response to Veome's review --- src/MainSequence.cpp | 4 ++-- src/TPAGB.h | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index d369b31ab..db3869415 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -533,11 +533,11 @@ double MainSequence::CalculateConvectiveCoreMass() const { DBL_DBL MainSequence::CalculateConvectiveEnvelopeMass() const { if (utils::Compare(m_Mass, 1.25) > 0) return std::tuple (0.0, 0.0); - double massEnvelope0 = 0.35; + double massEnvelope0 = m_Mass; if(utils::Compare(m_Mass, 0.35) > 0) massEnvelope0 = 0.35 * (1.25 - m_Mass) * (1.25 - m_Mass) / 0.81; double massEnvelope = massEnvelope0 * sqrt(sqrt(1.0 - m_Tau)); - return std::tuple (massEnvelope, massEnvelope); + return std::tuple (massEnvelope, massEnvelope0); } diff --git a/src/TPAGB.h b/src/TPAGB.h index 51f328984..ba29f442f 100755 --- a/src/TPAGB.h +++ b/src/TPAGB.h @@ -38,13 +38,13 @@ class TPAGB: virtual public BaseStar, public EAGB { // member functions - alphabetically double CalculateCOCoreMassAtPhaseEnd() const { return (utils::Compare(m_COCoreMass, m_GBParams[static_cast(GBP::McSN)]) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0) ? m_COCoreMass : m_Mass; } - double CalculateCOCoreMassOnPhase() const { return CalculateCoreMassOnPhase(m_Mass0, m_Age); } // McCO(TPAGB) = Mc(TPAGB)Same as on phase + double CalculateCOCoreMassOnPhase() const { return CalculateCoreMassOnPhase(m_Mass0, m_Age); } // McCO(TPAGB) = Mc(TPAGB)Same as on phase - double CalculateConvectiveCoreRadius () const { return 5.0 * CalculateRemnantRadius (); } // Last paragraph of section 6 of Hurley+ 2000 + double CalculateConvectiveCoreRadius () const { return std::min(5.0 * CalculateRemnantRadius(), m_Radius); } // Last paragraph of section 6 of Hurley+ 2000 double CalculateCoreMassAtPhaseEnd() const { return m_CoreMass; } // NO-OP 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 CalculateCoreMassOnPhase() const { return CalculateCoreMassOnPhase(m_Mass0, m_Age); } // Use class member variables double CalculateHeCoreMassAtPhaseEnd() const { return m_HeCoreMass; } // NO-OP double CalculateHeCoreMassOnPhase() const { return m_CoreMass; } // NO-OP