Adding custom force about reboundx

[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/dt*dt, 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.