Small bug fix for MEGNO calculation when adaptive timesteps are used.

hannorein · Jan 25, 2025 · 6c2578b · 6c2578b
1 parent 43ef972
commit 6c2578b
Show file tree

Hide file tree

Showing 5 changed files with 7 additions and 7 deletions.
diff --git a/src/integrator_eos.c b/src/integrator_eos.c
@@ -679,7 +679,7 @@ void reb_integrator_eos_part2(struct reb_simulation* const r){
         r->gravity_ignore_terms = 0;
         reb_calculate_acceleration_var(r);
         double dY = r->dt * 2. * r->t * reb_tools_megno_deltad_delta(r);
-        reb_tools_megno_update(r, dY);
+        reb_tools_megno_update(r, dY, r->dt);
     }
 
 }

diff --git a/src/integrator_ias15.c b/src/integrator_ias15.c
@@ -674,7 +674,7 @@ static int reb_integrator_ias15_step(struct reb_simulation* r) {
 
     if (r->calculate_megno){
         double dY = dt_done*integrator_megno_thisdt;
-        reb_tools_megno_update(r, dY);
+        reb_tools_megno_update(r, dY, dt_done);
     }
 
     // Swap particle buffers

diff --git a/src/integrator_whfast.c b/src/integrator_whfast.c
@@ -1190,7 +1190,7 @@ void reb_integrator_whfast_part2(struct reb_simulation* const r){
         // Update MEGNO in middle of timestep as we need synchonized x/v/a.
         if (r->calculate_megno){
             double dY = r->dt * 2. * r->t * reb_tools_megno_deltad_delta(r);
-            reb_tools_megno_update(r, dY);
+            reb_tools_megno_update(r, dY, dt);
         }
         if (ri_whfast->keep_unsynchronized){
             memcpy(p_j,sync_pj,r->N*sizeof(struct reb_particle));

diff --git a/src/tools.c b/src/tools.c
@@ -1500,12 +1500,12 @@ double reb_tools_megno_deltad_delta(struct reb_simulation* const r){
     return deltad/delta2;
 }
 
-void reb_tools_megno_update(struct reb_simulation* r, double dY){
+void reb_tools_megno_update(struct reb_simulation* r, double dY, double dt_done){
 	// Calculate running Y(t)
 	r->megno_Ys += dY;
 	double Y = r->megno_Ys/r->t;
 	// Calculate averge <Y> 
-	r->megno_Yss += Y * r->dt;
+	r->megno_Yss += Y * dt_done;
 	// Update covariance of (Y,t) and variance of t
 	r->megno_n++;
 	double _d_t = r->t - r->megno_mean_t;

diff --git a/src/tools.h b/src/tools.h
@@ -38,11 +38,11 @@ struct reb_particles;
 double reb_tools_megno_deltad_delta(struct reb_simulation* const r);
 
 /**
- * @brief Update MEGNO after a successful timestep by adding dY (=ddelta/delta*dt)
+ * @brief Update MEGNO after a successful timestep by adding dY (=ddelta/delta*dt_done)
  * @param r REBOUND simulation to be considered.
  * @param dY Increment for MEGNO Y
  */
-void reb_tools_megno_update(struct reb_simulation* r, double dY);
+void reb_tools_megno_update(struct reb_simulation* r, double dY, double dt_done);
 
 /**
  * @brief Init random number generator based on time and process id.