Skip to content

Commit

Permalink
Merge pull request #1046 from jeffriley/issue-1034
Browse files Browse the repository at this point in the history 
Issue 1034
  • Loading branch information
@jeffriley
jeffriley authored Dec 28, 2023
2 parents 133e4be + 18217b2 commit 2624083
Show file tree
Hide file tree
Showing 10 changed files with 25 additions and 17 deletions.
1 change: 1 addition & 0 deletions online-docs/pages/whats-new.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
* 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**
Expand Down
4 changes: 2 additions & 2 deletions src/BH.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
};
Expand Down
4 changes: 2 additions & 2 deletions src/BaseStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
6 changes: 3 additions & 3 deletions src/GiantBranch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -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; }
Expand Down
5 changes: 4 additions & 1 deletion src/GiantBranch.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
8 changes: 4 additions & 4 deletions src/MR.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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__
3 changes: 3 additions & 0 deletions src/MainSequence.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
4 changes: 2 additions & 2 deletions src/NS.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();

Expand Down
4 changes: 2 additions & 2 deletions src/Star.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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); }

Expand Down
3 changes: 2 additions & 1 deletion src/changelog.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 2624083

Please sign in to comment.