diff --git a/docs/index.md b/docs/index.md
index 845ff3406..5f7dec029 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -110,3 +110,4 @@ REBOUND is free software: you can redistribute it and/or modify it under the ter
REBOUND is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with REBOUND. If not, see .
+
diff --git a/docs/integrators.md b/docs/integrators.md
index 06c3f1ab0..bde2bd59a 100644
--- a/docs/integrators.md
+++ b/docs/integrators.md
@@ -404,7 +404,7 @@ The `reb_integrator_mercurius` structure contains the configuration and data str
TRACE is a hybrid time-reversible integrator, based on the algorithm described in [Hernandez & Dehnen 2023](https://ui.adsabs.harvard.edu/abs/2023MNRAS.522.4639H/abstract).
It uses WHFast for long term integrations but switches time-reversibly to BS or IAS15 for all close encounters. TRACE is appropriate for systems with a dominant central mass that will occasionally have close encounters.
-The TRACE implementation is described in UPDATE LINK [Lu et al in prep](https://ui.adsabs.harvard.edu/abs/2023MNRAS.522.4639H/abstract)
+A paper describing the TRACE implementation is in preparation.
The following code enables TRACE and sets the critical radius to 4 Hill radii
@@ -429,7 +429,8 @@ The `reb_integrator_trace` structure contains the configuration and data structu
If NULL (the default), the default switching function will be used.
The arguments `i` and `j` are the indices of the two particles considered.
The return value is a scalar.
- If this function returns a negative value, it means a close encounter has been flagged.
+ A negative return value means a close encounter has been flagged.
+ If the return values of both this function and the central switching function below are always positive, then the integrator effectively becomes the standard Wisdom-Holman integrator.
- Default switching function
@@ -458,14 +459,14 @@ The `reb_integrator_trace` structure contains the configuration and data structu
If NULL (the default), the default switching function will be used.
The argument `j` is the index of the non-central particle considered.
The return value is a scalar.
- If this function returns a negative value, it means a close encounter has been flagged. If the return values of both this function and the non-central switching function are always positive, then the integrator effectively becomes the standard Wisdom-Holman integrator.
+ A negative value means a close encounter has been flagged.
- Default switching function
This switching function checks if a body is close to its pericenter by considering a timescale derived from high-order derivatives of the particle's herliocentric position, inspired by [Pham, Rein, and Spiegel 2024](https://ui.adsabs.harvard.edu/abs/2024OJAp....7E...1P/abstract).
```c
- double reb_integrator_trace_switch_peri_pham(const struct reb_simulation* const r, const unsigned int j);
+ double reb_integrator_trace_switch_peri_default(const struct reb_simulation* const r, const unsigned int j);
```
- Fdot switching function
@@ -488,14 +489,19 @@ The `reb_integrator_trace` structure contains the configuration and data structu
=== "C"
```c
struct reb_simulation* r = reb_create_simulation();
+ r->ri_trace.S_peri = reb_integrator_trace_switch_peri_default; // default
r->ri_trace.S_peri = reb_integrator_trace_switch_peri_fdot;
r->ri_trace.S_peri = reb_integrator_trace_switch_peri_distance;
+ r->ri_trace.S_peri = reb_integrator_trace_switch_peri_none; // Turn off pericenter switching
```
=== "Python"
```python
sim = rebound.Simulation()
+ sim.ri_trace.S_peri = "default" # Following Pham et al 2024
+ sim.ri_trace.S_peri = "fdot"
sim.ri_trace.S_peri = "distance"
+ sim.ri_trace.S_peri = "none" # Turn off pericenter switching
```
`double r_crit_hill`
diff --git a/rebound/integrators/trace.py b/rebound/integrators/trace.py
index 72ba658a2..6f979be7d 100644
--- a/rebound/integrators/trace.py
+++ b/rebound/integrators/trace.py
@@ -85,7 +85,7 @@ def S_peri(self):
@S_peri.setter
def S_peri(self, func):
if func == "default":
- self._S_peri = cast(clibrebound.reb_integrator_trace_switch_peri_pham,TRACECF)
+ self._S_peri = cast(clibrebound.reb_integrator_trace_switch_peri_default,TRACECF)
elif func == "fdot":
self._S_peri = cast(clibrebound.reb_integrator_trace_switch_peri_fdot,TRACECF)
elif func == "distance":
diff --git a/src/integrator_trace.c b/src/integrator_trace.c
index b18650f20..cdf27dd60 100644
--- a/src/integrator_trace.c
+++ b/src/integrator_trace.c
@@ -115,8 +115,8 @@ int reb_integrator_trace_switch_peri_fdot(struct reb_simulation* const r, const
return peff2 < pfdot*pfdot * r->dt*r->dt;
}
-int reb_integrator_trace_switch_peri_pham(struct reb_simulation* const r, const unsigned int j){
- // Many square roots, can this be fixed?
+int reb_integrator_trace_switch_peri_default(struct reb_simulation* const r, const unsigned int j){
+ // Following Pham et al (2024)
const struct reb_integrator_trace* const ri_trace = &(r->ri_trace);
double GM = r->G*r->particles[0].m; // Not sure if this is the right mass to use.
@@ -544,7 +544,7 @@ void reb_integrator_trace_pre_ts_check(struct reb_simulation* const r){
const int N = r->N;
const int Nactive = r->N_active==-1?r->N:r->N_active;
int (*_switch) (struct reb_simulation* const r, const unsigned int i, const unsigned int j) = ri_trace->S ? ri_trace->S : reb_integrator_trace_switch_default;
- int (*_switch_peri) (struct reb_simulation* const r, const unsigned int j) = ri_trace->S_peri ? ri_trace->S_peri : reb_integrator_trace_switch_peri_pham;
+ int (*_switch_peri) (struct reb_simulation* const r, const unsigned int j) = ri_trace->S_peri ? ri_trace->S_peri : reb_integrator_trace_switch_peri_default;
// Clear encounter map
for (unsigned int i=1; iN; i++){
@@ -621,7 +621,7 @@ double reb_integrator_trace_post_ts_check(struct reb_simulation* const r){
const int N = r->N;
const int Nactive = r->N_active==-1?r->N:r->N_active;
int (*_switch) (struct reb_simulation* const r, const unsigned int i, const unsigned int j) = ri_trace->S ? ri_trace->S : reb_integrator_trace_switch_default;
- int (*_switch_peri) (struct reb_simulation* const r, const unsigned int j) = ri_trace->S_peri ? ri_trace->S_peri : reb_integrator_trace_switch_peri_pham;
+ int (*_switch_peri) (struct reb_simulation* const r, const unsigned int j) = ri_trace->S_peri ? ri_trace->S_peri : reb_integrator_trace_switch_peri_default;
int new_close_encounter = 0; // New CEs
// Clear encounter maps
diff --git a/src/rebound.h b/src/rebound.h
index fdf4585b6..f3431e2a0 100644
--- a/src/rebound.h
+++ b/src/rebound.h
@@ -849,7 +849,7 @@ DLLEXPORT double reb_integrator_mercurius_L_C5(const struct reb_simulation* cons
// Built in trace switching functions
-DLLEXPORT int reb_integrator_trace_switch_peri_pham(struct reb_simulation* const r, const unsigned int j);
+DLLEXPORT int reb_integrator_trace_switch_peri_default(struct reb_simulation* const r, const unsigned int j);
DLLEXPORT int reb_integrator_trace_switch_peri_fdot(struct reb_simulation* const r, const unsigned int j);
DLLEXPORT int reb_integrator_trace_switch_peri_distance(struct reb_simulation* const r, const unsigned int j);
DLLEXPORT int reb_integrator_trace_switch_peri_none(struct reb_simulation* const r, const unsigned int j);