I'm trying to plot a stellar orbit in a given potential. First, I initialize the position and velocity, and derive the acceleration from the position according to the given potential.
I then advance the time with a defined timestep and calculate the orbit. The problem arises when I try to store the calculated positions in an empty list. Here is my code:
## Initial position, velocity, and acceleration
r = np.array([20., 20., 0.])
v = np.array([0., 1., 0.])
g = acc(*r) #calculates acceleration from a function
## empty list to store position data
posdata = []
## Orbit integration
dt = 0.1
for t in np.arange(0, 1000, dt):
v += 0.5 * g * dt
r += v * dt
if t%100 == 0:
print(r) #check if r actually changes
g = acc(*r)
v += 0.5 * g * dt
posdata.append(r)
This is what I'm expecting to get:
posdata
>>> [array([19.999875, 20.099875, 0.]), array([19.99950125, 20.19950001, 0.]), array([19.99887999, 20.29887502, 0.]), ...]
But I actually get this:
>>> [array([-17.57080611, -34.03696644, 0.]), array([-17.57080611, -34.03696644, 0.]), array([-17.57080611, -34.03696644, 0.])]
All the elements are identical to the last element calculated. As you can see, I checked if r actually changes, and it does. I think it has to do with the fact that r is an array, but I don't know how to correct this.
rstays the same object during your loop, so you are appending the same object toposdataat every iteration. Quick fix:r = r + v*dt,