# Additional forces (iPython)¶

This tutorial gives you a very quick overview of how that works. Implementing additional forces in python as below will typically be a factor of a few slower than a C implementation. For a library that has C implementations for several commonly used additional effects (with everything callable from Python), see REBOUNDx.

**Stark problem**

We’ll start by adding two particles, the Sun and an Earth-like planet, to REBOUND.

```
import rebound
sim = rebound.Simulation()
sim.integrator = "whfast"
sim.add(m=1.)
sim.add(m=1e-6,a=1.)
sim.move_to_com() # Moves to the center of momentum frame
```

We could integrate this system, and the planet would go around the star at a fixed orbit with \(a=1\) forever. Let’s add an additional constant force that acts on the planet and points in one direction \(F_x = m\cdot c\), where \(m\) is the planet’s mass and \(c\) a constant. This is called the Stark problem. In python we can describe this with the following function

```
ps = sim.particles
c = 0.01
def starkForce(reb_sim):
ps[1].ax += c
```

Next, we need to tell REBOUND about this function.

```
sim.additional_forces = starkForce
```

Now we can just integrate as usual. Let’s keep track of the eccentricity as we integrate as it will change due to the additional force.

```
import numpy as np
Nout = 1000
es = np.zeros(Nout)
times = np.linspace(0.,100.*2.*np.pi,Nout)
for i, time in enumerate(times):
sim.integrate(time, exact_finish_time=0) # integrate to the nearest timestep so WHFast's timestep stays constant
es[i] = sim.particles[1].e
```

And let’s plot the result.

```
%matplotlib inline
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(15,5))
ax = plt.subplot(111)
plt.plot(times, es);
```

You can see that the eccentricity is oscillating between 0 and almost 1.

Note that the function `starkForce(reb_sim)`

above receives the
argument `reb_sim`

when it is called. This is a pointer to the
simulation structure. Instead of using the global `ps`

variable to
access particle data, one could also use `reb_sim.contents.particles`

.
This could be useful when one is running multiple simulations in
parallel or when the particles get added and removed (in those cases
`particles`

might change). The `contents`

has the same meaning a
`->`

in C, i.e. follow the memory address. To find out more about
pointers, check out the ctypes
documentation.

**Non-conservative forces**

The previous example assumed a conservative force, i.e. we could describe it as a potential as it is velocity independent. Now, let’s assume we have a velocity dependent force. This could be a migration force in a protoplanetary disk or PR drag. We’ll start from scratch and add the same two particles as before.

```
sim = rebound.Simulation()
sim.integrator = "ias15"
sim.add(m=1.)
sim.add(m=1e-6,a=1.)
sim.move_to_com() # Moves to the center of momentum frame
```

But we change the additional force to be

```
ps = sim.particles
tau = 1000.
def migrationForce(reb_sim):
ps[1].ax -= ps[1].vx/tau
ps[1].ay -= ps[1].vy/tau
ps[1].az -= ps[1].vz/tau
```

We need to let REBOUND know that our force is velocity dependent. Otherwise, REBOUND will not update the velocities of the particles.

```
sim.additional_forces = migrationForce
sim.force_is_velocity_dependent = 1
```

Now, we integrate as before. But this time we keep track of the semi-major axis instead of the eccentricity.

```
Nout = 1000
a_s = np.zeros(Nout)
times = np.linspace(0.,100.*2.*np.pi,Nout)
for i, time in enumerate(times):
sim.integrate(time)
a_s[i] = sim.particles[1].a
fig = plt.figure(figsize=(15,5))
ax = plt.subplot(111)
ax.set_xlabel("time")
ax.set_ylabel("semi-major axis")
plt.plot(times, a_s);
```

The semi-major axis decays exponentially on a timescale `tau`

.

In the above example, REBOUND is calling a python function at every
timestep. This can be slow. Note that you can also set
`rebound.additional_forces`

to a c function pointer. This let’s you
speed up the simulation significantly. However, you need to write you
own c function/library that knows how to calculate the forces. Or, you
use Dan Tamayo’s new migration library (in preparation).