Fix docs and add docstrings

compas_python_utils/preprocessing/compasConfigDefault.yaml

@@ -1,5 +1,5 @@
 ##~!!~## COMPAS option values
-##~!!~## File Created Tue May 14 19:39:51 2024 by COMPAS v02.46.00
+##~!!~## File Created Mon May 27 13:31:43 2024 by COMPAS v02.47.00
 ##~!!~##
 ##~!!~## The default COMPAS YAML file (``compasConfigDefault.yaml``), as distributed, has
 ##~!!~## all COMPAS option entries commented so that the COMPAS default value for the
@@ -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

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|

src/BaseBinaryStar.cpp

@@ -2371,70 +2371,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<double, double> 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<double, double> 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;
@@ -2466,12 +2417,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
 
@@ -2810,7 +2772,7 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
 
             // emit gravitational radiation if selected by user
             if (OPTIONS->EmitGravitationalRadiation()) {
-                EmitGW2(dt);
+                CalculateGravitationalRadiation(dt);
             }
 
             // check for problems

src/BaseBinaryStar.h

@@ -419,9 +419,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<double, double>    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);