From a3b1fac313616876e6feb35843f3d89a77659fa5 Mon Sep 17 00:00:00 2001 From: Jeff Riley Date: Thu, 28 Dec 2023 08:25:23 +1100 Subject: [PATCH] Defect repair: Fix for issue #1034 --- online-docs/pages/whats-new.rst | 5 ++ src/BaseBinaryStar.cpp | 112 ++++++++++---------------------- src/BaseBinaryStar.h | 54 ++++----------- src/BaseStar.h | 2 +- src/changelog.h | 13 +++- 5 files changed, 66 insertions(+), 120 deletions(-) diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 439a2c37e..78e12c3b9 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -6,6 +6,11 @@ Following is a brief list of important updates to the COMPAS code. A complete r **LATEST RELEASE** |br| +**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. +* This change may change DCO yields slightly when compared to previous versions of the code. + **02.41.00 Nov 02, 2023** * Added a naive tides implementation. diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 518c587f7..a00ab1768 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -272,12 +272,12 @@ void BaseBinaryStar::SetRemainingValues() { m_SemiMajorAxisAtDCOFormation = DEFAULT_INITIAL_DOUBLE_VALUE; m_EccentricityAtDCOFormation = DEFAULT_INITIAL_DOUBLE_VALUE; - double gyrationRadius1 = m_Star1->CalculateGyrationRadius(); - double gyrationRadius2 = m_Star2->CalculateGyrationRadius(); + double momentOfInertia1 = m_Star1->CalculateMomentOfInertia(); + double momentOfInertia2 = m_Star2->CalculateMomentOfInertia(); - m_TotalEnergy = CalculateTotalEnergy(m_SemiMajorAxis, m_Star1->Mass(), m_Star2->Mass(), m_Star1->RZAMS(), m_Star2->RZAMS(), m_Star1->Omega(), m_Star2->Omega(), gyrationRadius1, gyrationRadius2); + m_TotalEnergy = CalculateTotalEnergy(m_SemiMajorAxis, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Omega(), m_Star2->Omega(), momentOfInertia1, momentOfInertia2); - m_TotalAngularMomentum = CalculateAngularMomentum(m_SemiMajorAxis, m_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->RZAMS(), m_Star2->RZAMS(), m_Star1->Omega(), m_Star2->Omega(), gyrationRadius1, gyrationRadius2); + m_TotalAngularMomentum = CalculateAngularMomentum(m_SemiMajorAxis, m_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Omega(), m_Star2->Omega(), momentOfInertia1, momentOfInertia2); m_TotalAngularMomentumPrev = m_TotalAngularMomentum; m_Omega = 0.0; @@ -2046,51 +2046,32 @@ void BaseBinaryStar::InitialiseMassTransfer() { * double CalculateTotalEnergy(const double p_SemiMajorAxis, * const double p_Star1Mass, * const double p_Star2Mass, - * const double p_Star1Radius, - * const double p_Star2Radius, - * const double p_Star1_SpinAngularVelocity, - * const double p_Star1_SpinAngularVelocity, - * const double p_Star1_GyrationRadius, - * const double p_Star2_GyrationRadius) + * const double p_Star1SpinAngularVelocity, + * const double p_Star2SpinAngularVelocity, + * const double p_Star1MomentOfInertia, + * const double p_Star2MomentOfInertia) * * @param [IN] p_SemiMajorAxis Semi-major axis of the binary * @param [IN] p_Star1Mass Mass of star 1 * @param [IN] p_Star2Mass Mass of star 2 - * @param [IN] p_Star1Radius Radius of star 1 - * @param [IN] p_Star2Radius Radius of star 2 - * @param [IN] p_Star1_SpinAngularVelocity Spin angular velocity of star 1 - * @param [IN] p_Star1_SpinAngularVelocity Spin angular velocity of star 1 - * @param [IN] p_Star1_GyrationRadius Gyration radius of star 1 - * @param [IN] p_Star2_GyrationRadius Gyration radius of star 2 + * @param [IN] p_Star1SpinAngularVelocity Spin angular velocity of star 1 + * @param [IN] p_Star2SpinAngularVelocity Spin angular velocity of star 2 + * @param [IN] p_Star1MomentOfInertia Moment of inertia of star 1 + * @param [IN] p_Star2MomentOfInertia Moment of inertia of star 2 * @return Total energy of the binary */ double BaseBinaryStar::CalculateTotalEnergy(const double p_SemiMajorAxis, const double p_Star1Mass, const double p_Star2Mass, - const double p_Star1Radius, - const double p_Star2Radius, - const double p_Star1_SpinAngularVelocity, - const double p_Star2_SpinAngularVelocity, - const double p_Star1_GyrationRadius, - const double p_Star2_GyrationRadius) const { - double m1 = p_Star1Mass; - double m2 = p_Star2Mass; + const double p_Star1SpinAngularVelocity, + const double p_Star2SpinAngularVelocity, + const double p_Star1MomentOfInertia, + const double p_Star2MomentOfInertia) const { - double R1 = p_Star1Radius; - double R2 = p_Star2Radius; + double w1_2 = p_Star1SpinAngularVelocity * p_Star1SpinAngularVelocity; + double w2_2 = p_Star2SpinAngularVelocity * p_Star2SpinAngularVelocity; - double w1 = p_Star1_SpinAngularVelocity; - double w2 = p_Star2_SpinAngularVelocity; - - double ks1 = p_Star1_GyrationRadius; - double ks2 = p_Star2_GyrationRadius; - - constexpr double RSOL_TO_AU_2 = RSOL_TO_AU * RSOL_TO_AU; - - double Is1 = ks1 * m1 * R1 * R1 * RSOL_TO_AU_2; - double Is2 = ks2 * m2 * R2 * R2 * RSOL_TO_AU_2; - - return (0.5 * Is1 * w1 * w1) + (0.5 * Is2 * w2 * w2) - (0.5 * G_AU_Msol_yr * m1 * m2 / p_SemiMajorAxis); + return 0.5 * ((p_Star1MomentOfInertia * w1_2) + (p_Star2MomentOfInertia * w2_2) - (G_AU_Msol_yr * p_Star1Mass * p_Star2Mass / p_SemiMajorAxis)); } @@ -2101,57 +2082,36 @@ double BaseBinaryStar::CalculateTotalEnergy(const double p_SemiMajorAxis, * * * double CalculateAngularMomentum(const double p_SemiMajorAxis, - * const double p_Eccentricity, - * const double p_Star1Mass, - * const double p_Star2Mass, - * const double p_Star1Radius, - * const double p_Star2Radius, - * const double p_Star1_SpinAngularVelocity, - * const double p_Star1_SpinAngularVelocity, - * const double p_Star1_GyrationRadius, - * const double p_Star2_GyrationRadius) + * const double p_Eccentricity, + * const double p_Star1Mass, + * const double p_Star2Mass, + * const double p_Star1SpinAngularVelocity, + * const double p_Star1SpinAngularVelocity, + * const double p_Star1MomentOfInertia, + * const double p_Star2MomentOfInertia) * * @param [IN] p_SemiMajorAxis Semi-major axis of the binary * @param [IN] p_Eccentricity Eccentricity of the binary * @param [IN] p_Star1Mass Mass of the primary * @param [IN] p_Star2Mass Mass of the secondary - * @param [IN] p_Star1Radius Radius of the primary - * @param [IN] p_Star2Radius Radius of the secondary - * @param [IN] p_Star1_SpinAngularVelocity Orbital frequency of the primary - * @param [IN] p_Star1_SpinAngularVelocity Orbital frequency of the secondary - * @param [IN] p_Star1_GyrationRadius Gyration radius of the primary - * @param [IN] p_Star2_GyrationRadius Gyration radius of the secondary + * @param [IN] p_Star1SpinAngularVelocity Orbital frequency of the primary + * @param [IN] p_Star1SpinAngularVelocity Orbital frequency of the secondary + * @param [IN] p_Star1MomentOfInertia Moment of inertia of the primary + * @param [IN] p_Star2MomentOfInertia Moment of inertia of the secondary * @return Angular momentum of the binary */ double BaseBinaryStar::CalculateAngularMomentum(const double p_SemiMajorAxis, const double p_Eccentricity, const double p_Star1Mass, const double p_Star2Mass, - const double p_Star1Radius, - const double p_Star2Radius, - const double p_Star1_SpinAngularVelocity, - const double p_Star2_SpinAngularVelocity, - const double p_Star1_GyrationRadius, - const double p_Star2_GyrationRadius) const { - - double m1 = p_Star1Mass; - double m2 = p_Star2Mass; - - double R1 = p_Star1Radius * RSOL_TO_AU; - double R2 = p_Star2Radius * RSOL_TO_AU; - - double w1 = p_Star1_SpinAngularVelocity; - double w2 = p_Star2_SpinAngularVelocity; - - double ks1 = p_Star1_GyrationRadius; - double ks2 = p_Star2_GyrationRadius; - - double Is1 = ks1 * m1 * R1 * R1; - double Is2 = ks2 * m2 * R2 * R2; + const double p_Star1SpinAngularVelocity, + const double p_Star2SpinAngularVelocity, + const double p_Star1MomentOfInertia, + const double p_Star2MomentOfInertia) const { - double Jorb = CalculateOrbitalAngularMomentum(m1, m2, p_SemiMajorAxis, p_Eccentricity); + double Jorb = CalculateOrbitalAngularMomentum(p_Star1Mass, p_Star2Mass, p_SemiMajorAxis, p_Eccentricity); - return (Is1 * w1) + (Is2 * w2) + Jorb; + return (p_Star1MomentOfInertia * p_Star1SpinAngularVelocity) + (p_Star2MomentOfInertia * p_Star2SpinAngularVelocity) + Jorb; } diff --git a/src/BaseBinaryStar.h b/src/BaseBinaryStar.h index f6671dd8b..977e46f8b 100644 --- a/src/BaseBinaryStar.h +++ b/src/BaseBinaryStar.h @@ -407,35 +407,16 @@ class BaseBinaryStar { // member functions - alphabetical in groups (sort of...) - // CalculateAngularMomentum - the actual function takes 10 parameters because of the various calling permutations - // - various signatures are defined here - they just assemble the parameters as required - // and call the actual function - // JR: todo: note in the original code the binary orbital velicity was passed in as a parameter but never used - I removed it - - void SetInitialValues(const unsigned long int p_Seed, const long int p_Id); - void SetRemainingValues(); - double CalculateAngularMomentum(const double p_SemiMajorAxis, const double p_Eccentricity, const double p_Star1Mass, const double p_Star2Mass, - const double p_Star1Radius, - const double p_Star2Radius, - const double p_Star1_SpinAngularVelocity, - const double p_Star2_SpinAngularVelocity, - const double p_Star1_GyrationRadius, - const double p_Star2_GyrationRadius) const; - - double CalculateAngularMomentum() const { return CalculateAngularMomentum(m_SemiMajorAxis, m_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Radius(), m_Star2->Radius(), m_Star1->Omega(), m_Star2->Omega(), m_Star1->CalculateGyrationRadius(), m_Star2->CalculateGyrationRadius()); } - - double CalculateAngularMomentum(const double p_SemiMajorAxis, - const double p_Eccentricity, - const double p_Star1_SpinAngularVelocity, - const double p_Star2_SpinAngularVelocity, - const double p_Star1_GyrationRadius, - const double p_Star2_GyrationRadius) const { return CalculateAngularMomentum(p_SemiMajorAxis, p_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Radius(), m_Star2->Radius(), p_Star1_SpinAngularVelocity, p_Star2_SpinAngularVelocity, p_Star1_GyrationRadius, p_Star2_GyrationRadius); } + const double p_Star1SpinAngularVelocity, + const double p_Star2SpinAngularVelocity, + const double p_Star1MomentOfInertia, + const double p_Star2MomentOfInertia) const; - double CalculateAngularMomentumPrev() const { return CalculateAngularMomentum(m_SemiMajorAxisPrev, m_EccentricityPrev, m_Star1->MassPrev(), m_Star2->MassPrev(), m_Star1->RadiusPrev(), m_Star2->RadiusPrev(), m_Star1->OmegaPrev(), m_Star2->OmegaPrev(), m_Star1->CalculateGyrationRadius(), m_Star2->CalculateGyrationRadius()); } + double CalculateAngularMomentum() const { return CalculateAngularMomentum(m_SemiMajorAxis, m_Eccentricity, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Omega(), m_Star2->Omega(), m_Star1->CalculateMomentOfInertia(), m_Star2->CalculateMomentOfInertia()); } void CalculateEnergyAndAngularMomentum(); @@ -466,27 +447,15 @@ class BaseBinaryStar { double CalculateTimeToCoalescence(double a0, double e0, double m1, double m2) const; - // CalculateTotalEnergy - the actual function takes 9 parameters because of the various calling permutations - // - various signatures are defined here - they just assemble the parameters as required - // and call the actual function double CalculateTotalEnergy(const double p_SemiMajorAxis, const double p_Star1Mass, const double p_Star2Mass, - const double p_Star1Radius, - const double p_Star2Radius, - const double p_Star1_SpinAngularVelocity, - const double p_Star2_SpinAngularVelocity, - const double p_Star1GyrationRadius, - const double p_Star2GyrationRadius) const; - - double CalculateTotalEnergy() const { return CalculateTotalEnergy(m_SemiMajorAxis, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Radius(), m_Star2->Radius(), m_Star1->Omega(), m_Star2->Omega(), m_Star1->CalculateGyrationRadius(), m_Star2->CalculateGyrationRadius()); } - - double CalculateTotalEnergy(const double p_SemiMajorAxis, - const double p_Star1_SpinAngularVelocity, - const double p_Star2_SpinAngularVelocity, - const double p_Star1_GyrationRadius, - const double p_Star2_GyrationRadius) const { return CalculateTotalEnergy(p_SemiMajorAxis, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Radius(), m_Star2->Radius(), p_Star1_SpinAngularVelocity, p_Star2_SpinAngularVelocity, p_Star1_GyrationRadius, p_Star2_GyrationRadius); } + const double p_Star1SpinAngularVelocity, + const double p_Star2SpinAngularVelocity, + const double p_Star1MomentOfInertia, + const double p_Star2MomentOfInertia) const; + double CalculateTotalEnergy() const { return CalculateTotalEnergy(m_SemiMajorAxis, m_Star1->Mass(), m_Star2->Mass(), m_Star1->Omega(), m_Star2->Omega(), m_Star1->CalculateMomentOfInertia(), m_Star2->CalculateMomentOfInertia()); } void EvaluateBinary(const double p_Dt); @@ -500,6 +469,9 @@ class BaseBinaryStar { void ResolveMassChanges(); bool ResolveSupernova(); + void SetInitialValues(const unsigned long int p_Seed, const long int p_Id); + void SetRemainingValues(); + void SetPostCEEValues(const double p_SemiMajorAxis, const double p_Eccentricity, const double p_RocheLobe1to2, diff --git a/src/BaseStar.h b/src/BaseStar.h index 517bee446..34672250a 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -40,7 +40,7 @@ class BaseStar { // getters - alphabetically double Age() const { return m_Age; } - double AngularMomentum() const { return CalculateGyrationRadius() * m_Radius * RSOL_TO_AU * m_Radius * RSOL_TO_AU * m_Omega; } + double AngularMomentum() const { return CalculateMomentOfInertia() * m_Omega; } double BindingEnergy_Fixed() const { return m_BindingEnergies.fixed; } double BindingEnergy_Nanjing() const { return m_BindingEnergies.nanjing; } double BindingEnergy_Loveridge() const { return m_BindingEnergies.loveridge; } diff --git a/src/changelog.h b/src/changelog.h index 40a299b8e..92ab962c2 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1072,10 +1072,19 @@ // - Fix for issue #1022 - incorrect index used for last array entry. // - A little code cleanup // 02.41.02 JR - Dec 15, 2023 - Defect repair: -// - 2.41.00 backed-ou the changes made in 2.40.00 - this puts them back +// - 2.41.00 backed-out the changes made in 2.40.00 - this puts them back // - Calling it a defect repair so we get a new version number - just in case we need it... +// 02.41.03 JR - Dec 28, 2023 - Defect repair: +// - Fix for issue #1034 +// - This fix changes the functions +// . BaseBinaryStar::CalculateAngularMomentum(), +// . BaseBinaryStar::CalculateTotalEnergy(), and +// . BaseStar::AngularMomentum() +// to use moment of inertia rather than gyration radius. This has wider implications than just issue #1034 +// and may change DCO yields slightly. +// - Change to functionality (noted above) noted in 'What's New' online documentation page -const std::string VERSION_STRING = "02.41.02"; +const std::string VERSION_STRING = "02.41.03"; # endif // __changelog_h__