From 3a9426799cb95fbe1bb77296a104fe7589c6471c Mon Sep 17 00:00:00 2001 From: Jeff Riley Date: Thu, 28 Dec 2023 22:02:43 +1100 Subject: [PATCH 1/3] Changed CalculateMomentOfInertia() to properly implement Hurley et al., 2000 eq 109 --- src/BH.h | 4 ++-- src/BaseStar.h | 4 ++-- src/GiantBranch.cpp | 6 +++--- src/GiantBranch.h | 5 ++++- src/MR.h | 8 ++++---- src/MainSequence.h | 3 +++ src/NS.h | 4 ++-- src/Star.h | 4 ++-- src/changelog.h | 3 ++- 9 files changed, 24 insertions(+), 17 deletions(-) diff --git a/src/BH.h b/src/BH.h index c7ea512b1..0fdb19d0b 100755 --- a/src/BH.h +++ b/src/BH.h @@ -54,8 +54,8 @@ class BH: virtual public BaseStar, public Remnants { double CalculateMassLossRate() { return 0.0; } // Ensure that BHs don't lose mass in winds - double CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return (2.0 / 5.0) * m_Mass * m_Radius * m_Radius; } - double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return CalculateMomentOfInertia(p_RemnantRadius * RSOL_TO_AU) * RSOL_TO_AU * RSOL_TO_AU; } + double CalculateMomentOfInertia() const { return (2.0 / 5.0) * m_Mass * m_Radius * m_Radius; } + double CalculateMomentOfInertiaAU() const { return CalculateMomentOfInertia() * RSOL_TO_AU * RSOL_TO_AU; } double CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase_Static(m_Mass); } // Use class member variables - returns radius in Rsol }; diff --git a/src/BaseStar.h b/src/BaseStar.h index ee608578e..665462be8 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -186,8 +186,8 @@ class BaseStar { double CalculateMassLossValues(const bool p_UpdateMDot = false, const bool p_UpdateMDt = false); // JR: todo: better name? - virtual double CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return (0.1 * (m_Mass - m_CoreMass) * m_Radius * m_Radius) + (0.21 * m_CoreMass * p_RemnantRadius * p_RemnantRadius); } // k2 = 0.1 and k3 = 0.21 as defined in Hurley et al. 2000, after eq 109 - virtual double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return CalculateMomentOfInertia(p_RemnantRadius * RSOL_TO_AU) * RSOL_TO_AU * RSOL_TO_AU; } + virtual double CalculateMomentOfInertia() const { return (0.1 * (m_Mass) * m_Radius * m_Radius); } // Defaults to MS. k2 = 0.1 as defined in Hurley et al. 2000, after eq 109 + virtual double CalculateMomentOfInertiaAU() const { return CalculateMomentOfInertia() * RSOL_TO_AU * RSOL_TO_AU; } double CalculateNuclearTimescale() const { return CalculateNuclearTimescale_Static(m_Mass, m_Luminosity); } // Use class member variables diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index 6bb424a0b..0f4b3deff 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -1151,10 +1151,10 @@ STELLAR_TYPE GiantBranch::CalculateRemnantTypeByMuller2016(const double p_COCore * STELLAR_TYPE CalculateRemnantTypeBySchneider2020(const double p_COCoreMass) * * @param [IN] p_COCoreMass COCoreMass in Msol - * @param [IN] useSchneiderAlt Whether to use the Schneider alt prescription + * @param [IN] p_UseSchneiderAlt Whether to use the Schneider alt prescription * @return Remnant mass in Msol */ -double GiantBranch::CalculateRemnantMassBySchneider2020(const double p_COCoreMass, const bool p_useSchneiderAlt) { +double GiantBranch::CalculateRemnantMassBySchneider2020(const double p_COCoreMass, const bool p_UseSchneiderAlt) { double logRemnantMass; STYPE_VECTOR mtHist = MassTransferDonorHistory(); @@ -1188,7 +1188,7 @@ double GiantBranch::CalculateRemnantMassBySchneider2020(const double p_COCoreMas case MT_CASE::NONE: // No history of MT - if (!p_useSchneiderAlt) { // Use standard or alternative remnant mass prescription for effectively single stars? + if (!p_UseSchneiderAlt) { // Use standard or alternative remnant mass prescription for effectively single stars? // standard prescription if (utils::Compare(p_COCoreMass, 6.357) < 0) { logRemnantMass = log10(0.03357*p_COCoreMass + 1.31780); } else if (utils::Compare(p_COCoreMass, 7.311) < 0) { logRemnantMass = -0.02466*p_COCoreMass + 1.28070; } diff --git a/src/GiantBranch.h b/src/GiantBranch.h index 9714d86b4..bc366ae63 100755 --- a/src/GiantBranch.h +++ b/src/GiantBranch.h @@ -84,9 +84,12 @@ class GiantBranch: virtual public BaseStar, public MainSequence { DBL_DBL CalculateRemnantMassByFryer2022(const double p_Mass, const double p_COCoreMass); double CalculateRemnantMassByMuller2016(const double p_Mass, const double p_COCoreMass); double CalculateRemnantMassByMullerMandel(const double p_COCoreMass, const double p_HeCoreMass); - double CalculateRemnantMassBySchneider2020(const double p_COCoreMass, const bool p_useSchneiderAlt = false); + double CalculateRemnantMassBySchneider2020(const double p_COCoreMass, const bool p_UseSchneiderAlt = false); double CalculateRemnantMassBySchneider2020Alt(const double p_COCoreMass) { return CalculateRemnantMassBySchneider2020(p_COCoreMass, true); } + double CalculateMomentOfInertia() const { double Rc = CalculateRemnantRadius(); return (0.1 * (m_Mass - m_CoreMass) * m_Radius * m_Radius) + (0.21 * m_CoreMass * Rc * Rc); } // k2 = 0.1 and k3 = 0.21 as defined in Hurley et al. 2000, after eq 109 + double CalculateMomentOfInertiaAU() const { return CalculateMomentOfInertia() * RSOL_TO_AU * RSOL_TO_AU; } + double CalculatePerturbationMu() const; double CalculateRadialExtentConvectiveEnvelope() const; diff --git a/src/MR.h b/src/MR.h index 58dbc83e9..67490b91e 100755 --- a/src/MR.h +++ b/src/MR.h @@ -43,13 +43,13 @@ class MR: virtual public BaseStar, public Remnants { // member functions - double CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return 0.0; } // No moment of inertia for massless remnants - use 0.0 - double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return 0.0; } // No moment of inertia for massless remnants - use 0.0 + double CalculateMomentOfInertia() const { return 0.0; } // No moment of inertia for massless remnants - use 0.0 + double CalculateMomentOfInertiaAU() const { return 0.0; } // No moment of inertia for massless remnants - use 0.0 void SetPulsarParameters() const { } // NO-OP - bool ShouldEvolveOnPhase() const { return true; } // Always - bool ShouldSkipPhase() const { return false; } // Don't skip + bool ShouldEvolveOnPhase() const { return true; } // Always + bool ShouldSkipPhase() const { return false; } // Don't skip }; #endif // __MR_h__ diff --git a/src/MainSequence.h b/src/MainSequence.h index 7760f7041..c38e62b05 100644 --- a/src/MainSequence.h +++ b/src/MainSequence.h @@ -52,6 +52,9 @@ class MainSequence: virtual public BaseStar { double CalculateLuminosityOnPhase(const double p_Time, const double p_Mass, const double p_LZAMS) const; double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase(m_Age, m_Mass0, m_LZAMS0); } // Use class member variables + double CalculateMomentOfInertia() const { return (0.1 * (m_Mass) * m_Radius * m_Radius); } // k2 = 0.1 as defined in Hurley et al. 2000, after eq 109 + double CalculateMomentOfInertiaAU() const { return CalculateMomentOfInertia() * RSOL_TO_AU * RSOL_TO_AU; } + double CalculatePerturbationMu() const { return 5.0; } // mu(MS) = 5.0 (Hurley et al. 2000, eqs 97 & 98) double CalculateRadialExtentConvectiveEnvelope() const; diff --git a/src/NS.h b/src/NS.h index aad0adb40..8fe5c1423 100755 --- a/src/NS.h +++ b/src/NS.h @@ -72,8 +72,8 @@ class NS: virtual public BaseStar, public Remnants { double CalculateMassLossRate() { return 0.0; } // Ensure that NSs don't lose mass in winds double CalculateMomentOfInertiaCGS() const; // MoI in CGS - double CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return CalculateMomentOfInertiaCGS() / MSOL_TO_G / RSOL_TO_CM / RSOL_TO_CM; } // MoI (default is solar units) - double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return CalculateMomentOfInertia() * RSOL_TO_AU * RSOL_TO_AU; } + double CalculateMomentOfInertia() const { return CalculateMomentOfInertiaCGS() / MSOL_TO_G / RSOL_TO_CM / RSOL_TO_CM; } // MoI (default is solar units) + double CalculateMomentOfInertiaAU() const { return CalculateMomentOfInertia() * RSOL_TO_AU * RSOL_TO_AU; } double CalculatePulsarBirthMagneticField(); diff --git a/src/Star.h b/src/Star.h index a18b7579a..674aa2c91 100755 --- a/src/Star.h +++ b/src/Star.h @@ -171,8 +171,8 @@ class Star { double CalculateMassLossValues(const bool p_UpdateMDot = false, const bool p_UpdateMDt = false) { return m_Star->CalculateMassLossValues(p_UpdateMDot, p_UpdateMDt); } - double CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return m_Star->CalculateMomentOfInertia(p_RemnantRadius); } - double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return m_Star->CalculateMomentOfInertiaAU(p_RemnantRadius); } + double CalculateMomentOfInertia() const { return m_Star->CalculateMomentOfInertia(); } + double CalculateMomentOfInertiaAU() const { return m_Star->CalculateMomentOfInertiaAU(); } void CalculateSNAnomalies(const double p_Eccentricity) { m_Star->CalculateSNAnomalies(p_Eccentricity); } diff --git a/src/changelog.h b/src/changelog.h index bb2ff3911..b610e63da 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1081,7 +1081,8 @@ // . BaseBinaryStar::CalculateTotalEnergy(), and // . BaseStar::AngularMomentum() // to use moment of inertia rather than gyration radius. -// It also removes CalculateGyrationRadius() from all classes, and changes code that called CalculateGyrationRadius(). +// This fix changes CalculateMomentOfInertia to properly implement Hurley et al., 2000 eq 109 +// This fix also removes CalculateGyrationRadius() from all classes, and changes code that called CalculateGyrationRadius(). // These changes have wider implications than just issue #1034 and may change DCO yields slightly. // - Removed some unused functions. // - Change to functionality (noted above) noted in 'What's New' online documentation page From 1354bbfa193112bab4e82828af66e5b20288e437 Mon Sep 17 00:00:00 2001 From: Jeff Riley Date: Thu, 28 Dec 2023 22:05:00 +1100 Subject: [PATCH 2/3] Update What's New --- online-docs/pages/whats-new.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 78e12c3b9..1834652bd 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -9,6 +9,7 @@ Following is a brief list of important updates to the COMPAS code. A complete r **02.41.03 Dec 28, 2023** * The functions ``BaseBinaryStar::CalculateAngularMomentum()``, ``BaseBinaryStar::CalculateTotalEnergy()``, and ``BaseStar::AngularMomentum()`` changed to use moment of inertia instead of gyration radius. +* CalculateMomentOfInertia() to properly implement Hurley et al., 2000 eq 109. * This change may change DCO yields slightly when compared to previous versions of the code. **02.41.00 Nov 02, 2023** From 18217b20627f636a720ed6cd5e2470b3609b1b94 Mon Sep 17 00:00:00 2001 From: Jeff Riley Date: Thu, 28 Dec 2023 22:06:21 +1100 Subject: [PATCH 3/3] Update What's New --- online-docs/pages/whats-new.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 1834652bd..ee4360849 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -9,7 +9,7 @@ Following is a brief list of important updates to the COMPAS code. A complete r **02.41.03 Dec 28, 2023** * The functions ``BaseBinaryStar::CalculateAngularMomentum()``, ``BaseBinaryStar::CalculateTotalEnergy()``, and ``BaseStar::AngularMomentum()`` changed to use moment of inertia instead of gyration radius. -* CalculateMomentOfInertia() to properly implement Hurley et al., 2000 eq 109. +* Changed CalculateMomentOfInertia() to properly implement Hurley et al., 2000 eq 109. * This change may change DCO yields slightly when compared to previous versions of the code. **02.41.00 Nov 02, 2023**