Update BaseBinaryStar.cpp

Matching Jeff's PR, #1

src/BaseBinaryStar.cpp

@@ -2554,27 +2554,26 @@ void BaseBinaryStar::EmitGravitationalWave(const 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.
+ * Returns a timestep based on the minimal timesteps of the component stars, 
+ * adjusted if relevant by the orbital evolution due to GW radiation
  * 
  *
- * double ChooseTimestep(const double p_Dt)
+ * double ChooseTimestep(const double p_Multiplier)
  * 
- * @param   [IN]    p_Dt                        previous timestep in Myr
+ * @param   [IN]    p_Multiplier                timestep multiplier
  * @return                                      new timestep in Myr
  */
-double BaseBinaryStar::ChooseTimestep(const double p_Dt) {
+double BaseBinaryStar::ChooseTimestep(const double p_Multiplier) {
 
-    double newDt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep());        			// new timestep
+    double dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep());                   // timestep based on orbital timescale
 
-    if (OPTIONS->EmitGravitationalRadiation()) {                                                                        // emitting GWs?
-        newDt = std::min(newDt, -1.0E-2 * m_SemiMajorAxis / m_DaDtGW);                                                  // reduce timestep if necessary to ensure that the orbital separation does not change by more than ~1% per timestep due to GW emission
+    if (OPTIONS->EmitGravitationalRadiation()) {                                                        // emitting GWs?
+        dt = std::min(dt, -1.0E-2 * m_SemiMajorAxis / m_DaDtGW);                                        // yes - reduce timestep if necessary to ensure that the orbital separation does not change by more than ~1% per timestep due to GW emission
     }
 
-    newDt *= OPTIONS->TimestepMultiplier();	
+    dt *= p_Multiplier;	
 
-    return std::max(std::round(newDt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP);                // quantised and not less than minimum
+    return std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP);    // quantised and not less than minimum
 }
 
 
@@ -2780,7 +2779,28 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
             // evolve the current binary up to the maximum evolution time (and number of steps)
 
             double dt;                                                                                                                  // timestep
-            unsigned long int stepNum = 0;                                                                                              // initialise step number
+            if (usingProvidedTimesteps) {                                                                                               // user-provided timesteps?
+                // get new timestep
+                //   - don't quantise
+                //   - don't apply timestep multiplier
+                // (we assume user wants the timesteps in the file)
+                // 
+                // Open question: should we clamp this to NUCLEAR_MINIMUM_TIMESTEP?
+                dt = timesteps[0];
+            }
+            else {                                                                                                                      // no - not using user-provided timesteps
+                // if user selects to emit GWs, calculate the effects of radiation
+                //     - note that this is placed before the call to ChooseTimestep() because when
+                //       emitting GWs the timestep is a function of graviational radiation
+                if (OPTIONS->EmitGravitationalRadiation()) {
+                    CalculateGravitationalRadiation();
+                }
+
+                // we want the first timestep to be small - calculate timestep and divide by 1000.0
+                dt = ChooseTimestep(OPTIONS->TimestepMultiplier() / 1000.0);                                                            // calculate timestep - make first step small
+            }
+
+            unsigned long int stepNum = 1; 
 
             while (evolutionStatus == EVOLUTION_STATUS::CONTINUE) {                                                                     // perform binary evolution - iterate over timesteps until told to stop
 
@@ -2818,6 +2838,9 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
                     evolutionStatus = EVOLUTION_STATUS::SSE_ERROR;                                                                      // yes - stop evolution
                 }
                 else {                                                                                                                  // continue evolution
+		    if (OPTIONS->EmitGravitationalRadiation()) {                                                                        // emitting GWs?
+                        EmitGravitationalWave(dt);                                                                                      // yes - emit graviataional wave
+                    }
 
                     (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_STELLAR_TIMESTEP);                                   // print (log) detailed output
 
@@ -2900,6 +2923,33 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
                     SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED);                                                                              // show warning
                 }
 
+                if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) {                                                                    // continue evolution?
+
+                    // if user selects to emit GWs, calculate the effects of radiation
+                    //     - note that this is placed before the call to ChooseTimestep() because when
+                    //       emitting GWs the timestep is a function of graviational radiation                    
+                    if (OPTIONS->EmitGravitationalRadiation()) {
+                        CalculateGravitationalRadiation();
+                    }
+
+                    m_Star2->UpdatePreviousTimestepDuration();
+                    m_Star1->UpdatePreviousTimestepDuration();
+
+                    if (usingProvidedTimesteps) {                                                                                       // user-provided timesteps?
+                        // select a timestep
+                        //   - don't quantise
+                        //   - don't apply timestep multiplier
+                        // (we assume user wants the timesteps in the file)
+                        // 
+                        // Open question: should we clamp this to NUCLEAR_MINIMUM_TIMESTEP?
+                        dt = timesteps[stepNum];
+                    }
+                    else {                                                                                                              // no - not using user-provided timesteps
+                        dt = ChooseTimestep(OPTIONS->TimestepMultiplier());
+                    }
+
+                    stepNum++;                                                                                                          // increment stepNum
+                }
             }
 
             if (usingProvidedTimesteps && timesteps.size() > stepNum) {                                                                 // all user-defined timesteps consumed?