<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Additional forces
REBOUND is a gravitational N-body integrator. But you can also use it to integrate systems with additional, non-gravitational forces.
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 be adding two particles, the Sun and an Earth-like planet to REBOUND.
Step1:  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 acting on the planet and is pointing 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
Step2:  Next, we need to tell REBOUND about this function.
Step3:  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.
Step4:  And let's plot the result.
Step5:  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 -&gt; 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.
Step6:  But we change the additional force to be
Step7:  We need to let REBOUND know that our force is velocity dependent. Otherwise, REBOUND will not update the velocities of the particles.
Step8:  Now, we integrate as before. But this time we keep track of the semi-major axis instead of the eccentricity.
<ASSISTANT_TASK:>
Python Code:
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
Explanation: Additional forces
REBOUND is a gravitational N-body integrator. But you can also use it to integrate systems with additional, non-gravitational forces.
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 be adding two particles, the Sun and an Earth-like planet to REBOUND.
End of explanation
ps = sim.particles
c = 0.01
def starkForce(reb_sim):
    ps[1].ax += c
Explanation: 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 acting on the planet and is pointing 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
End of explanation
sim.additional_forces = starkForce
Explanation: Next, we need to tell REBOUND about this function.
End of explanation
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 
Explanation: 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.
End of explanation
%matplotlib inline
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(15,5))
ax = plt.subplot(111)
plt.plot(times, es);
Explanation: And let's plot the result.
End of explanation
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
Explanation: 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 -&gt; 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.
End of explanation
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
Explanation: But we change the additional force to be
End of explanation
sim.additional_forces = migrationForce
sim.force_is_velocity_dependent = 1
Explanation: We need to let REBOUND know that our force is velocity dependent. Otherwise, REBOUND will not update the velocities of the particles.
End of explanation
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);
Explanation: Now, we integrate as before. But this time we keep track of the semi-major axis instead of the eccentricity.
End of explanation
<END_TASK>