Example Notebooks and code cleanup

rebound/integrators/trace.py

@@ -41,11 +41,11 @@ class IntegratorTRACE(ctypes.Structure):
 
     """
     def __repr__(self):
-        return '<{0}.{1} object at {2}, r_crit_hill={3}, peri_mode=={4}, peri_crit_eta=={5},peri_crit_fdot=={6},peri_crit_distance={7}>'.format(self.__module__, type(self).__name__, hex(id(self)), self.r_crit_hill, self.peri_crit_fdot, self.peri_crit_distance)
+        return '<{0}.{1} object at {2}, r_crit_hill={3}, peri_mode=={4}, peri_crit_eta=={5},peri_crit_fdot=={6},peri_crit_distance={7}>'.format(self.__module__, type(self).__name__, hex(id(self)), self.r_crit_hill, self.peri_mode, self.peri_crit_eta, self.peri_crit_fdot, self.peri_crit_distance)
 
     _fields_ = [("_S", ctypes.CFUNCTYPE(ctypes.c_int, ctypes.POINTER(Simulation), ctypes.c_uint, ctypes.c_uint)),
                 ("_S_peri", ctypes.CFUNCTYPE(ctypes.c_int, ctypes.POINTER(Simulation), ctypes.c_uint)),
-                ("peri_mode", ctypes.c_uint),
+                ("_peri_mode", ctypes.c_uint),
                 ("r_crit_hill", ctypes.c_double),
                 ("peri_crit_eta", ctypes.c_double),
                 ("peri_crit_fdot", ctypes.c_double),
@@ -107,19 +107,18 @@ def peri_mode(self):
         - ``'PARTIAL_BS'`` (Integrate only the Kepler step with BS)
         - ``'FULL_IAS15'`` (Integrate entire system with IAS15)
         """
-        i = self.peri_mode
+        i = self._peri_mode
         for name, _i in TRACE_PERI_MODES.items():
             if i==_i:
                 return name
         return i
     @peri_mode.setter
     def peri_mode(self, value):
         if isinstance(value, int):
-            self.peri_mode = ctypes.c_uint(value)
+            self._peri_mode = ctypes.c_uint(value)
         elif isinstance(value, basestring):
-            value = value.lower().replace(" ", "")
             if value in TRACE_PERI_MODES: 
-                self.peri_mode = TRACE_PERI_MODES[value]
+                self._peri_mode = TRACE_PERI_MODES[value]
             else:
                 raise ValueError("Warning. Pericenter switching mode not found.")
 

src/integrator_bs.c

@@ -283,7 +283,7 @@ static int tryStep(struct reb_simulation* r, const int Ns, const int k, const in
         double* yDot = odes[s]->yDot;
         const int length = odes[s]->length;
         for (int i = 0; i < length; ++i) {
-            y1[i] = 0.5 * (yTmp[i] + y1[i] + subStep * yDot[i]); // = 0.25*(y_(2n-1) + 2*y_n(2) + y_(2n+1))     Eq (9.13c)
+            y1[i] = 0.5 * (yTmp[i] + y1[i] + subStep * yDot[i]); 
         }
     }
 
@@ -292,7 +292,7 @@ static int tryStep(struct reb_simulation* r, const int Ns, const int k, const in
 
 static void extrapolate(const struct reb_ode* ode, double * const coeff, const int k) {
     double* const y1 = ode->y1;
-    double* const C = ode->C;  // C and D values follow Numerical Recipes 
+    double* const C = ode->C;
     double** const D =  ode->D;
     double const length = ode->length;
     for (int j = 0; j < k; ++j) {
@@ -302,8 +302,8 @@ static void extrapolate(const struct reb_ode* ode, double * const coeff, const i
         double facD = xim1/(xi-xim1);
         for (int i = 0; i < length; ++i) {
             double CD = C[i] - D[k - j -1][i];
-            C[i] = facC * CD; // Only need to keep one C value
-            D[k - j - 1][i] = facD * CD; // Keep all D values for recursion
+            C[i] = facC * CD;
+            D[k - j - 1][i] = facD * CD;
         }
     }
     for (int i = 0; i < length; ++i) {
@@ -339,7 +339,7 @@ static void nbody_derivatives(struct reb_ode* ode, double* const yDot, const dou
 
 void reb_integrator_bs_part1(struct reb_simulation* r){
     struct reb_ode** odes = r->odes;
-    int Ns = r->N_odes; // Number of ode sets
+    int Ns = r->N_odes;
     for (int s=0; s < Ns; s++){
         const int length = odes[s]->length;
         double* y0 = odes[s]->y;

src/integrator_trace.c

@@ -93,13 +93,11 @@ int reb_integrator_trace_switch_peri_distance(struct reb_simulation* const r, co
 int reb_integrator_trace_switch_peri_fdot(struct reb_simulation* const r, const unsigned int j){
     const struct reb_integrator_trace* const ri_trace = &(r->ri_trace);
     const double pfdot = ri_trace->peri_crit_fdot;
-    const double pdist = ri_trace->peri_crit_distance;
 
     const double dx = r->particles[j].x;
     const double dy = r->particles[j].y;
     const double dz = r->particles[j].z;
     const double d2 = dx*dx + dy*dy + dz*dz;
-    if (d2 < pdist*pdist) return 1;
 
     const double dvx = r->particles[j].vx;
     const double dvy = r->particles[j].vy;
@@ -114,8 +112,6 @@ int reb_integrator_trace_switch_peri_fdot(struct reb_simulation* const r, const
     const double fdot2 = h2 / (d2*d2);
     const double peff2 = (4 * M_PI * M_PI) / fdot2; // effective period squared
 
-
-    // Failsafe: use pericenter pericenter distance
     return peff2 < pfdot*pfdot * r->dt*r->dt;
 }
 
@@ -325,7 +321,6 @@ void reb_integrator_trace_nbody_derivatives(struct reb_ode* ode, double* const y
     reb_integrator_trace_update_particles(r, y);
     reb_simulation_update_acceleration(r);
 
-    // TLu Levison & Duncan 22, 23 EoMs
     double px=0., py=0., pz=0.;
     int* map = r->ri_trace.encounter_map;
     int N = r->ri_trace.encounter_N;
@@ -722,7 +717,6 @@ static void reb_integrator_trace_step(struct reb_simulation* const r){
         const double old_t = r->t;
         r->gravity = REB_GRAVITY_BASIC;
         r->ri_trace.mode = REB_TRACE_MODE_FULL; // for collision search
-        //reb_integrator_trace_dh_to_inertial(r);
         switch (r->ri_trace.peri_mode){
             case REB_TRACE_PERI_FULL_IAS15:
                 // Run default IAS15 integration
@@ -773,7 +767,6 @@ static void reb_integrator_trace_step(struct reb_simulation* const r){
 
                         reb_integrator_bs_update_particles(r, nbody_ode->y);
 
-			// Does this need the star velocity changes too? Or is it OK because it's in inertial
                         reb_collision_search(r);
                     }
                     reb_ode_free(nbody_ode);
@@ -786,7 +779,6 @@ static void reb_integrator_trace_step(struct reb_simulation* const r){
         r->gravity = REB_GRAVITY_TRACE;
         r->t = old_t; // final time will be set later
         r->dt = old_dt;
-        //reb_integrator_trace_inertial_to_dh(r); // TODO: This should be optimized.
     }
 }
 

src/rebound.h

@@ -257,7 +257,7 @@ struct reb_integrator_sei {
     double tandtz;      // Cached tan(), z axis
 };
 
-// TRACE (Lu et al. 2023)
+// TRACE (Lu Hernandez & Rein 2024)
 struct reb_integrator_trace {
     int (*S) (struct reb_simulation* const r, const unsigned int i, const unsigned int j);
     int (*S_peri) (struct reb_simulation* const r, const unsigned int j);
@@ -288,16 +288,16 @@ struct reb_integrator_trace {
     unsigned int N_allocated_additional_forces;
     unsigned int tponly_encounter; // 0 if any encounters are between two massive bodies. 1 if encounters only involve test particles
 
-    struct reb_particle* REB_RESTRICT particles_backup; //  TLu contains coordinates before the entire step
-    struct reb_particle* REB_RESTRICT particles_backup_kepler; //  TLu contains coordinates before kepler step
+    struct reb_particle* REB_RESTRICT particles_backup; //  Contains coordinates before the entire step
+    struct reb_particle* REB_RESTRICT particles_backup_kepler; //  Contains coordinates before kepler step
     struct reb_particle* REB_RESTRICT particles_backup_additional_forces; // For additional forces
 
     int* encounter_map;             // Map to represent which particles are integrated with BS
     struct reb_vec3d com_pos;       // Used to keep track of the centre of mass during the timestep
     struct reb_vec3d com_vel;
 
-    int* current_Ks; // TLu tracking K for the entire timestep
-    unsigned int current_C; // TLu tracking C for the entire timestep
+    int* current_Ks; // Tracking K_ij for the entire timestep
+    unsigned int current_C; // Tracking C for the entire timestep
     unsigned int force_accept; // Force accept for irreversible steps: collisions and adding particles
 };