diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 439a2c37e..ee4360849 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -6,6 +6,12 @@ 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. +* 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** * Added a naive tides implementation. diff --git a/src/BH.h b/src/BH.h index b597c966a..f5e62165f 100755 --- a/src/BH.h +++ b/src/BH.h @@ -50,14 +50,11 @@ class BH: virtual public BaseStar, public Remnants { double CalculateEddingtonCriticalRate() const { return 2.6E-8 * m_Mass * MYR_TO_YEAR; } // E.g., Marchant+, 2017, Eq. 3, assuming accretion efficiency of 10% - double CalculateGyrationRadius() const { return 0.0; } // No tidal coupling to a BH - double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase_Static(); } 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 CalculateRadiusOnPhase() const { return CalculateRadiusOnPhase_Static(m_Mass); } // Use class member variables - returns radius in Rsol }; diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 518c587f7..8c75a118d 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->CalculateMomentOfInertiaAU(); + double momentOfInertia2 = m_Star2->CalculateMomentOfInertiaAU(); - 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..63cd97bc8 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->CalculateMomentOfInertiaAU(), m_Star2->CalculateMomentOfInertiaAU()); } 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->CalculateMomentOfInertiaAU(), m_Star2->CalculateMomentOfInertiaAU()); } 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..665462be8 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 CalculateMomentOfInertiaAU() * 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; } @@ -175,8 +175,6 @@ class BaseStar { virtual void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { } // Default is NO-OP virtual void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables - virtual double CalculateGyrationRadius() const { return 0.0; } // Default is 0.0 - void CalculateLambdas() { CalculateLambdas(m_Mass - m_CoreMass); } // Use class member variables void CalculateLambdas(const double p_EnvMass); @@ -188,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.0; } // Use inheritance hierarchy - virtual double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return 0.0; } // Use inheritance hierarchy + 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/BinaryConstituentStar.cpp b/src/BinaryConstituentStar.cpp index 1017c3ddf..ec7982b88 100644 --- a/src/BinaryConstituentStar.cpp +++ b/src/BinaryConstituentStar.cpp @@ -384,7 +384,8 @@ double BinaryConstituentStar::CalculateCircularisationTimescale(const double p_S */ double BinaryConstituentStar::CalculateSynchronisationTimescale(const double p_SemiMajorAxis) { - double gyrationRadiusSquared_1 = 1.0 / CalculateGyrationRadius(); + double gyrationRadiusSquared = CalculateMomentOfInertia() / Mass() / Radius() / Radius(); + double gyrationRadiusSquared_1 = 1.0 / gyrationRadiusSquared; double rOverA = Radius() / p_SemiMajorAxis; double rOverA_6 = rOverA * rOverA * rOverA * rOverA * rOverA * rOverA; double q2 = m_Companion->Mass() / Mass(); diff --git a/src/CHeB.h b/src/CHeB.h index fe0a7b0f1..187cd5fa7 100755 --- a/src/CHeB.h +++ b/src/CHeB.h @@ -71,8 +71,6 @@ class CHeB: virtual public BaseStar, public FGB { double CalculateHeCoreMassAtPhaseEnd() const { return m_CoreMass; } - double CalculateGyrationRadius() const { return 0.21; } // Hurley et al., 2000, after eq 109 for n=3/2 polytrope or dense convective core. Single number approximation. - double CalculateLambdaDewi() const; double CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) const; double CalculateLambdaNanjingEnhanced(const int p_MassInd, const int p_Zind) const; diff --git a/src/EAGB.h b/src/EAGB.h index 3212ed4a2..85dc59373 100755 --- a/src/EAGB.h +++ b/src/EAGB.h @@ -53,8 +53,6 @@ class EAGB: virtual public BaseStar, public CHeB { double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return GiantBranch::CalculateCriticalMassRatioHurleyHjellmingWebbink(); } - double CalculateGyrationRadius() const { return 0.1; } // Hurley et al., 2000, after eq 109 for giants. Single number approximation. JR: todo: should this be in constants.h? - double CalculateHeCoreMassAtPhaseEnd() const { return CalculateHeCoreMassOnPhase(); } // Same as on phase double CalculateHeCoreMassOnPhase() const { return m_HeCoreMass; } // NO-OP diff --git a/src/FGB.h b/src/FGB.h index 2a08f0e70..d189b6c95 100755 --- a/src/FGB.h +++ b/src/FGB.h @@ -47,8 +47,6 @@ class FGB: virtual public BaseStar, public HG { double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const { return GiantBranch::CalculateCriticalMassRatioClaeys14(p_AccretorIsDegenerate); } // Skip HG double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return GiantBranch::CalculateCriticalMassRatioHurleyHjellmingWebbink(); } - - double CalculateGyrationRadius() const { return 0.1; } // Hurley et al., 2000, after eq 109 for giants. Single number approximation. double CalculateHeCoreMassAtPhaseEnd() const { return CalculateHeCoreMassOnPhase(); } // Same as on phase double CalculateHeCoreMassOnPhase() const { return m_CoreMass; } // McHe(FGB) = Core Mass diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index 6bb424a0b..05878c085 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -692,7 +692,7 @@ double GiantBranch::CalculateRemnantRadius() const { * * @return Radial extent of the star's convective envelope in Rsol */ -double GiantBranch::CalculateRadialExtentConvectiveEnvelope() const{ +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 @@ -1108,6 +1108,29 @@ double GiantBranch::CalculateLifetimeToHeIgnition(const double p_Mass, const dou } +/////////////////////////////////////////////////////////////////////////////////////// +// // +// ROTATION CALCULATIONS // +// // +/////////////////////////////////////////////////////////////////////////////////////// + +/* + * Calculate moment of inertia + * + * Hurley et al., 2000, paragraph immediately following eq 109 + * + * + * double GiantBranch::CalculateMomentOfInertia() + * + * @return Moment of inertia (Msol AU^2) + */ +double GiantBranch::CalculateMomentOfInertia() const { + double Rc = CalculateRemnantRadius(); + + return (0.1 * (m_Mass - m_CoreMass) * m_Radius * m_Radius) + (0.21 * m_CoreMass * Rc * Rc); +} + + /////////////////////////////////////////////////////////////////////////////////////// // // // SUPERNOVA CALCULATIONS // @@ -1151,10 +1174,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 +1211,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; } @@ -1673,12 +1696,6 @@ double GiantBranch::CalculateRemnantMassByBelczynski2002(const double p_Mass, co } - - - - - - /* * Driver function for Core Collapse Supernovas * diff --git a/src/GiantBranch.h b/src/GiantBranch.h index 04232c322..3630cd043 100755 --- a/src/GiantBranch.h +++ b/src/GiantBranch.h @@ -84,10 +84,10 @@ 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 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 - 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; double CalculatePerturbationMu() const; diff --git a/src/HG.h b/src/HG.h index 53a92a201..676125d63 100755 --- a/src/HG.h +++ b/src/HG.h @@ -71,8 +71,6 @@ class HG: virtual public BaseStar, public GiantBranch { double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const; double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 0.25; } // As coded in BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details. - double CalculateGyrationRadius() const { return 0.21; } // Hurley et al., 2000, after eq 109 for n=3/2 polytrope or dense convective core. Single number approximation. - double CalculateHeCoreMassAtPhaseEnd() const { return m_CoreMass; } // McHe(HG) = Core Mass double CalculateHeCoreMassOnPhase() const { return m_CoreMass; } // McHe(HG) = Core Mass diff --git a/src/HeGB.h b/src/HeGB.h index 950dde742..15c619c2b 100755 --- a/src/HeGB.h +++ b/src/HeGB.h @@ -52,8 +52,6 @@ class HeGB: virtual public BaseStar, public HeHG { 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. - double CalculateGyrationRadius() const { return 0.1; } // Hurley et al., 2000, after eq 109 for giants. Single number approximation. - double CalculateLuminosityOnPhase(const double p_CoreMass, const double p_GBPB, const double p_GBPD) const { return CalculateLuminosityOnPhase_Static(p_CoreMass, p_GBPB, p_GBPD); } double CalculateLuminosityOnPhase() const { return CalculateLuminosityOnPhase(m_CoreMass, m_GBParams[static_cast(GBP::B)], m_GBParams[static_cast(GBP::D)]); } diff --git a/src/HeHG.h b/src/HeHG.h index 389e007a3..a198e8983 100755 --- a/src/HeHG.h +++ b/src/HeHG.h @@ -66,8 +66,6 @@ class HeHG: virtual public BaseStar, public HeMS { void CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams); void CalculateGBParams() { CalculateGBParams(m_Mass0, m_GBParams); } // Use class member variables - double CalculateGyrationRadius() const { return 0.21; } // Hurley et al., 2000, after eq 109 for n=3/2 polytrope or dense convective core. Single number approximation. - double CalculateHeCoreMassAtPhaseEnd() const { return CalculateHeCoreMassOnPhase(); } // Same as on phase double CalculateHeCoreMassOnPhase() const { return m_Mass; } // NO-OP @@ -79,9 +77,6 @@ class HeHG: virtual public BaseStar, public HeMS { double CalculateMassTransferRejuvenationFactor() const; - double CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return GiantBranch::CalculateMomentOfInertia(p_RemnantRadius); } // Skip HeMS - double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return GiantBranch::CalculateMomentOfInertiaAU(p_RemnantRadius); } // Skip HeMS - double CalculatePerturbationMu() const; double CalculatePerturbationMuAtPhaseEnd() const { return m_Mu; } // NO-OP diff --git a/src/HeMS.h b/src/HeMS.h index 890cfee8a..6fb67fb4f 100644 --- a/src/HeMS.h +++ b/src/HeMS.h @@ -64,8 +64,6 @@ class HeMS: virtual public BaseStar, public TPAGB { double CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const; double CalculateCriticalMassRatioHurleyHjellmingWebbink() const { return 0.33; } // As coded in BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details. - double CalculateGyrationRadius() const { return 0.1; } - double CalculateHeCoreMassOnPhase() const { return m_Mass; } // McHe(HeMS) = Mass double CalculateHeCoreMassAtPhaseEnd() const { return CalculateHeCoreMassOnPhase(); } // Same as on phase @@ -86,9 +84,6 @@ class HeMS: virtual public BaseStar, public TPAGB { double CalculateMassTransferRejuvenationFactor() const; - double CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return MainSequence::CalculateMomentOfInertia(p_RemnantRadius); } // Use MainSequence - double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return MainSequence::CalculateMomentOfInertiaAU(p_RemnantRadius); } // Use MainSequence - 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 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.cpp b/src/MainSequence.cpp index 4d5a04a10..d5b9f2f22 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -562,50 +562,6 @@ void MainSequence::UpdateAgeAfterMassLoss() { } -/////////////////////////////////////////////////////////////////////////////////////// -// // -// ROTATIONAL / GYRATION / FREQUENCY CALCULATIONS // -// // -/////////////////////////////////////////////////////////////////////////////////////// - -/* - * Calculate gyration radius - * - * Define gyration radius 'k=r_g^2' using fit from de Mink et al. 2013, calling k_definition function - * Original created by Alejandro Vigna-Gomez on 11/2015. Rewritten June 2019, JR. - * - * The original fits from de Mink+2013 were made for MS stars a Z=0.02. - * - * Uses class member variables instead of passing in parameters - * - * - * double CalculateGyrationRadius() - * - * @return Gyration radius in Rsol - * - */ -double MainSequence::CalculateGyrationRadius() const { - - double log10M = log10(m_Mass); - - double cLower = 0.0; // correction factor 'c' (lowercase 'c') in de Mink et al., 2013 eq A1 - if ((utils::Compare(log10M, 1.3) > 0)) { // log10(M) > 1.3 (de Mink doesn't include '=' - we assume it not to be here)) - double log10M_13 = log10M - 1.3; - cLower = -0.055 * log10M_13 * log10M_13; - } - - double CUpper = -2.5; // exponent 'C' (uppercase 'C') in de Mink et al., 2013 eq A2 - if ((utils::Compare(log10M, 0.2) > 0)) CUpper = -1.5; // log10(M) > 0.2 - else if ((utils::Compare(log10M, 0.0) > 0)) CUpper = -2.5 + (5.0 * log10M); // 0.2 <= log10(M) > 0.0 (de Mink doesn't include '=' - we assume it to be here (and for log10(M) <= 0.0)) - - double k0 = cLower + std::min(0.21, std::max(0.09 - (0.27 * log10M), 0.037 + (0.033 * log10M)));// gyration radius squared for ZAMS stars - - double radiusRatio = m_Radius / m_RZAMS; - - return ((k0 - 0.025) * PPOW(radiusRatio, CUpper)) + (0.025 * PPOW(radiusRatio, -0.1)); // gyration radius -} - - /////////////////////////////////////////////////////////////////////////////////////// // // // MISCELLANEOUS FUNCTIONS / CONTROL FUNCTIONS // diff --git a/src/MainSequence.h b/src/MainSequence.h index 42f85cffe..a5ce35fa2 100644 --- a/src/MainSequence.h +++ b/src/MainSequence.h @@ -42,8 +42,6 @@ class MainSequence: virtual public BaseStar { double CalculateCoreMassAtPhaseEnd() const { return OPTIONS->RetainCoreMassDuringCaseAMassTransfer() ? MinimumCoreMass() : 0.0; } // Accounts for minimal core mass built up prior to mass loss through mass transfer double CalculateCoreMassOnPhase() const { return 0.0; } // Mc(MS) = 0.0 (Hurley et al. 2000, just before eq 28) - double CalculateGyrationRadius() const; - double CalculateHeCoreMassAtPhaseEnd() const { return CalculateCoreMassAtPhaseEnd(); } // Same as He core mass double CalculateHeCoreMassOnPhase() const { return 0.0; } // McHe(MS) = 0.0 @@ -54,8 +52,7 @@ 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 double p_RemnantRadius = 0.0) const { return CalculateGyrationRadius() * 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 (0.1 * (m_Mass) * m_Radius * m_Radius); } // k2 = 0.1 as defined in Hurley et al. 2000, after eq 109 double CalculatePerturbationMu() const { return 5.0; } // mu(MS) = 5.0 (Hurley et al. 2000, eqs 97 & 98) diff --git a/src/NS.h b/src/NS.h index aad0adb40..6a37bc8e9 100755 --- a/src/NS.h +++ b/src/NS.h @@ -72,8 +72,7 @@ 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 CalculatePulsarBirthMagneticField(); diff --git a/src/Remnants.h b/src/Remnants.h index c3c9d7003..24cdb94aa 100644 --- a/src/Remnants.h +++ b/src/Remnants.h @@ -40,8 +40,6 @@ class Remnants: virtual public BaseStar, public HeGB { 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 CalculateGyrationRadius() const { return 0.21; } // Hurley et al., 2000, after eq 109 for n=3/2 polytrope or dense convective core. Single number approximation. - double CalculateHeCoreMassOnPhase() const { return m_Mass; } // Return m_Mass double CalculateInitialSupernovaMass() const { return GiantBranch::CalculateInitialSupernovaMass(); } // Use GiantBranch @@ -58,9 +56,6 @@ class Remnants: virtual public BaseStar, public HeGB { double CalculateMassLossRateHurley() { return 0.0; } double CalculateMassLossRateBelczynski2010() { return 0.0; } - double CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return GiantBranch::CalculateMomentOfInertia(p_RemnantRadius); } // Default to GiantBranch - double CalculateMomentOfInertiaAU(const double p_RemnantRadius = 0.0) const { return GiantBranch::CalculateMomentOfInertiaAU(p_RemnantRadius); } // Default to GiantBranch - double CalculatePerturbationMuOnPhase() const { return m_Mu; } // NO-OP double CalculateRadialExtentConvectiveEnvelope() const { return BaseStar::CalculateRadialExtentConvectiveEnvelope(); } // WD stars don't have a convective envelope diff --git a/src/Star.h b/src/Star.h index a6b1a41eb..674aa2c91 100755 --- a/src/Star.h +++ b/src/Star.h @@ -161,7 +161,6 @@ class Star { double CalculateConvectiveEnvelopeMass() { return m_Star->CalculateConvectiveEnvelopeMass(); } double CalculateEddyTurnoverTimescale() { return m_Star->CalculateEddyTurnoverTimescale(); } - double CalculateGyrationRadius() const { return m_Star->CalculateGyrationRadius(); } void CalculateLambdas() { m_Star->CalculateLambdas(); } void CalculateLambdas(const double p_EnvMass) { m_Star->CalculateLambdas(p_EnvMass); } @@ -172,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/TPAGB.h b/src/TPAGB.h index 37b6dd363..06cd85cc3 100755 --- a/src/TPAGB.h +++ b/src/TPAGB.h @@ -44,8 +44,6 @@ class TPAGB: virtual public BaseStar, public EAGB { 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 CalculateGyrationRadius() const { return 0.1; } // Hurley et al., 2000, after eq 109 for giants. Single number approximation. - double CalculateHeCoreMassAtPhaseEnd() const { return m_HeCoreMass; } // NO-OP double CalculateHeCoreMassOnPhase() const { return m_CoreMass; } // NO-OP diff --git a/src/changelog.h b/src/changelog.h index 40a299b8e..b610e63da 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1072,10 +1072,22 @@ // - 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 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 -const std::string VERSION_STRING = "02.41.02"; +const std::string VERSION_STRING = "02.41.03"; # endif // __changelog_h__