Skip to content

Commit

Permalink
Defect repair:
Browse files Browse the repository at this point in the history 
Fix for issue TeamCOMPAS#1034
  • Loading branch information
@jeffriley
jeffriley committed Dec 27, 2023
1 parent acf3283 commit a3b1fac
Show file tree
Hide file tree
Showing 5 changed files with 66 additions and 120 deletions.
5 changes: 5 additions & 0 deletions online-docs/pages/whats-new.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
112 changes: 36 additions & 76 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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));
}


Expand All @@ -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;
}


Expand Down
54 changes: 13 additions & 41 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();

Expand Down Expand Up @@ -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);

Expand All @@ -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,
Expand Down
2 changes: 1 addition & 1 deletion src/BaseStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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; }
Expand Down
13 changes: 11 additions & 2 deletions src/changelog.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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__

0 comments on commit a3b1fac

Please sign in to comment.