From 5aaae537ecf886078f37b8cff6963c531e0546d2 Mon Sep 17 00:00:00 2001 From: Adam Boesky Date: Mon, 27 May 2024 13:55:39 -0700 Subject: [PATCH] Fix docs and add docstrings --- .../preprocessing/compasConfigDefault.yaml | 2 +- .../program-options-list-defaults.rst | 2 +- src/BaseBinaryStar.cpp | 82 +++++-------------- src/BaseBinaryStar.h | 4 +- 4 files changed, 25 insertions(+), 65 deletions(-) diff --git a/compas_python_utils/preprocessing/compasConfigDefault.yaml b/compas_python_utils/preprocessing/compasConfigDefault.yaml index 168c5d8ab..2c99a4810 100644 --- a/compas_python_utils/preprocessing/compasConfigDefault.yaml +++ b/compas_python_utils/preprocessing/compasConfigDefault.yaml @@ -14,7 +14,7 @@ booleanChoices: # --enable-warnings: False # Default: False # option to enable/disable warning messages # --errors-to-file: False # Default: False # --evolve-unbound-systems: True # Default: True -# --emit-gravitational-radiation: True # Default: True +# --emit-gravitational-radiation: False # Default: False # --population-data-printing: False # Default: False # --print-bool-as-string: False # Default: False # --quiet: False # Default: False diff --git a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst index f8a242d5b..47a161750 100644 --- a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst +++ b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst @@ -351,7 +351,7 @@ Default = 1.0 **--emit-gravitational-radiation** |br| Emit gravitational radiation at each timestep of binary evolution according to Peters 1964. |br| -Default = TRUE +Default = FALSE **--enable-warnings** |br| Display warning messages to stdout. |br| diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index ee86be2e4..4d8162423 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2388,70 +2388,21 @@ void BaseBinaryStar::ResolveMassChanges() { /* - * Update the GW effects on the binary + * Calculate and emit gravitational radiation. * - * Use Peters 1964's approximations of the effects of emitting GWs by: - * - Updating the binary's semi-major axis (eq 5.6) - * - Updating the binary's eccentricity (eq 5.7) + * This function uses Peters 1964 to approximate the effects of GW emission with two steps: + * - Update m_SemiMajorAxis using the change in semi-major axis per time given by eq 5.6. + * - Update m_Eccentricity using the change in eccentricity per time given by eq 5.7. * - * std::tuple CalculateGravitationalRadiation() + * Notably, m_DaDtGW and m_DeDtGW are updated so that they can be used to calculate the timestep dynamically. + * m_SemiMajorAxisPrev and m_SemiMajorAxisPrev are also modified for later use. * - * @return A tuple of the rate of change in semimajor axis (AU/Myr) and - * eccentricity (1/Myr) due to GW emission. -*/ -std::tuple BaseBinaryStar::CalculateGravitationalRadiation() { - - // Useful values - double eccentricitySquared = m_Eccentricity * m_Eccentricity; - double oneMinusESq = 1.0 - eccentricitySquared; - double oneMinusESq_5 = oneMinusESq * oneMinusESq * oneMinusESq * oneMinusESq * oneMinusESq; - double G_AU_Msol_yr_3 = G_AU_Msol_yr * G_AU_Msol_yr * G_AU_Msol_yr; - double C_AU_Yr_5 = C_AU_yr * C_AU_yr * C_AU_yr * C_AU_yr * C_AU_yr; - double m_SemiMajorAxis_3 = m_SemiMajorAxis * m_SemiMajorAxis * m_SemiMajorAxis; - double massAndGAndCTerm = G_AU_Msol_yr_3 * m_Star1->Mass() * m_Star2->Mass() * (m_Star1->Mass() + m_Star2->Mass()) / C_AU_Yr_5; // G^3 * m1 * m2(m1 + m2) / c^5 in units of Msol, AU and yr - - // Approximate rate of change in semimajor axis - double numeratorA = -64.0 * massAndGAndCTerm; - double denominatorA = 5.0 * m_SemiMajorAxis_3 * std::sqrt(oneMinusESq_5 * oneMinusESq * oneMinusESq); - double DaDtGW = (numeratorA / denominatorA) * (1.0 + (73.0 / 24.0) * eccentricitySquared + (37.0 / 96.0) * eccentricitySquared * eccentricitySquared) * MYR_TO_YEAR; // units of AU Myr^-1 - - // Approximate rate of change in eccentricity - double numeratorE = -304.0 * m_Eccentricity * massAndGAndCTerm; - double denominatorE = 15.0 * m_SemiMajorAxis_3 * m_SemiMajorAxis * std::sqrt(oneMinusESq_5); - double DeDtGW = (numeratorE / denominatorE) * (1.0 + (121.0 / 304.0) * eccentricitySquared) * YEAR_TO_MYR; // units of Myr^-1 - - return std::make_tuple(DaDtGW, DeDtGW); -} - - -/* - * Emit a GW based on the effects calculated by BaseBinaryStar::CalculateGravitationalRadiation(). - * - * This function updates the semi-major axis, eccentricity, and previous eccentricity values - * (m_SemiMajorAxis, m_Eccentricity, and m_EccentricityPrev) as a result of emitting the GW. + * void CalculateGravitationalRadiation() * - * void EmitGravitationalRadiation(const double p_Dt) - * * @param [IN] p_Dt timestep in Myr - * @param [IN] DaDtGW change in semimajor axis per time in AU/Myr due to GW radiation - * @param [IN] DeDtGW change in eccentricity per time in Myr^-1 due to GW radiation + * */ -void BaseBinaryStar::EmitGravitationalWave(const double p_Dt, const double DaDtGW, const double DeDtGW) { - - // Update semimajor axis - double aNew = m_SemiMajorAxis + (DaDtGW * p_Dt); - m_SemiMajorAxis = utils::Compare(aNew, 0.0) > 0 ? aNew : 1E-20; // if <0, set to arbitrarily small number - - // Update the eccentricity - m_Eccentricity += DeDtGW * p_Dt; - - // Save values as previous timestep - m_SemiMajorAxisPrev = m_SemiMajorAxis; - m_EccentricityPrev = m_Eccentricity; -} - - -void BaseBinaryStar::EmitGW2(double p_Dt) { +void BaseBinaryStar::CalculateGravitationalRadiation(double p_Dt) { // Useful values double eccentricitySquared = m_Eccentricity * m_Eccentricity; @@ -2483,12 +2434,23 @@ void BaseBinaryStar::EmitGW2(double p_Dt) { } +/* + * Choose a timestep based on the parameters of the binary. + * + * This function will return the minimum of (i) a timestep based on the + * orbital timescale of the binary and (ii) (if configured to emit GWs) + * a timestep based on the magnitude of gravitational radiation. + * + * double ChooseTimestep(double p_Dt) + * + * @param [IN] p_Dt previous timestep in Myr + * @return new timestep in Myr +*/ double BaseBinaryStar::ChooseTimestep(double p_Dt) { m_Star1->UpdatePreviousTimestepDuration(); m_Star2->UpdatePreviousTimestepDuration(); - p_Dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep()) * OPTIONS->TimestepMultiplier(); // update dt for orbital timescale p_Dt = std::round(p_Dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM; // quantised @@ -2827,7 +2789,7 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() { // emit gravitational radiation if selected by user if (OPTIONS->EmitGravitationalRadiation()) { - EmitGW2(dt); + CalculateGravitationalRadiation(dt); } // check for problems diff --git a/src/BaseBinaryStar.h b/src/BaseBinaryStar.h index f4930d4d1..dda57c0fe 100644 --- a/src/BaseBinaryStar.h +++ b/src/BaseBinaryStar.h @@ -423,9 +423,7 @@ class BaseBinaryStar { 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()); } - std::tuple CalculateGravitationalRadiation(); - void EmitGravitationalWave(const double p_Dt, const double DaDtGW, const double DeDtGW); - void EmitGW2(double p_Dt); + void CalculateGravitationalRadiation(const double p_Dt); double ChooseTimestep(double p_Dt);