Adding custom force

[lenunez13]

Hello!

I am using REBOUNDx to (1) add additional forces on a massless planet in a 2-body system and (2) track the planet's eccentricity as a function of time. As per the latest version of Custom_Effects.ipynb, I use rebx.add_force(myforce) to call my Python-based orbital evolution function. This function is a line-by-line translation of some Mercury subroutine MFO_USER.FOR (which was motivated by the appendix in Wolff et al. 2012). It forces additional acceleration and velocity terms with an arbitrarily assigned (cosine) function of e-dot and updates the terms via lines like ps[1].ax += ax and ps[1].vx += vx. I then compare REBOUND's results to Mercury's and to the corresponding eccentricity model.

My issue is that REBOUND does not agree with the model, whereas Mercury does, despite the implementation of similar code structures. See figure below. My suspicion is that the issue is rooted in the user-defined velocities ps[1].vx += vx, since the results look much better (albeit not perfect) when I turn the velocities off (i.e., do not update them).

I would like to get a better sense of what REBOUND is doing with these velocities. Where in the source code can I find the velocity calculations?

Thank you, Luis

[dtamayo]

Hi,

You should only be updating the accelerations in your custom routine. The REBOUND integration routines will then use the accelerations to update the positions and velocities. If you still don't get the correct behavior, I'm happy to take a look at it if you want to post your code!

[lenunez13]

Hi,

The Mercury subroutine and the aforementioned paper require that additional velocity terms also be updated. Can REBOUND update velocities by ps[1].vx += as is done with acceleration?

Sure! If you need to see additional details, see my code here.

Thanks!

[dtamayo]

So unless I'm misinterpreting something, this looks like a prescription for updating the positions and velocities of all the particles, x += dx/dtdt, vx += dxdot/dt dt etc. You would add this as an operator rather than an additional force. See our REBOUNDx paper on arxiv.

A different version of this is already implemented in REBOUNDx (modify_orbits_forces as a force and modify_orbits_direct as an operator). Do you need something about this implementation in particular? You can look at the migration and eccentricity damping examples in ipython_examples for how they work.

[dtamayo]

I'm going to close this for now, but feel free to reopen if you have any follow-up questions!

[lenunez13]

Thanks for the help! Using the operator proved a better approach.

A follow-up question. I successfully added my own effect by copying a version of modify_orbits_direct.c. I would like to implement a time-dependent eccentricity, so I added const double t to the static struct and changed the operator to (for example) o.e += 0.5*cos(t/(*tau_e)). In general, is this a correct approach for implementing time-dependence?

[dtamayo]

Hi Luis,

Are you trying to enforce a specific e(t), or add a specific e(t) on top of whatever the eccentricity is doing due to other planets?

If the former, you can just put in the e(t) you want e.g., o.e = sin(omega * sim->t)

If the latter, you can do what's in modify_orbits_direct, which is just a simple euler step. o.e += omegacos(omegasim->t)*dt

In both cases you should probably use the simulation time sim->t rather than adding a const double t, and you should probably add a separate variable like omega rather than reusing tau_e, so that you (or others) could in principle use both effects at the same time. If it's just a one-off thing I guess it doesn't matter.

Note that in both cases it will not work with an adaptive timestep (IAS15) or Mercurius (which only allows forces). For IAS15 you can set sim.ri_ias15.epsilon = 0 to make it use a constant timestep sim.dt.