I cannot claim very deep computational expertise, so in case it might help you, I can just post some code that contains two versions of a fixed-step straight-forward symplectic algorithms (leap-frog / Verlet's type). I however am not so proficient in C++, so I can offer a python version, so maybe you can translate it into C++. On the bright side, the python is more transparent in terms of methodological ideas and organization, so that could be beneficial. The code also does not feature the realistic constants, but scaled versions.
Also, the code is set up to two mass-points in a plane, because you specified that you want to start with a two-body simulation. I hope that's ok to get you started.
'''
Two body leap-forg algorithm simulation
'''
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
'''
Calculation of the gravitational acceleration
r is a 2 x 2 matrix, each row is the position of one of the two mass-points
'''
def g_accel(r, mass1, mass2):
r_12 = r[1,:] - r[0,:]
g_force = mass1 * mass2 * r_12 / np.linalg.norm( r_12 )**3
return np.array([g_force / mass1, - g_force / mass2])
'''
propagation/numerical integration using leap-from / verlet's type time-reversible symplectic algorithm of convergence rate dt**2
'''
def propagate_system(r1_in, r2_in, v1_in, v2_in, mass1, mass2, n_steps, dt):
r = np.empty((n_steps, 2, 2),dtype=float)
v = np.empty((n_steps, 2, 2),dtype=float)
r[0,0,:] = r1_in
r[0,1,:] = r2_in
v[0,0,:] = v1_in
v[0,1,:] = v2_in
for n in range(n_steps-1):
r[n+1,:,:] = r[n,:,:] + dt * v[n,:,:] / 2
v[n+1,:,:] = v[n,:,:] + dt * g_accel(r[n+1,:], mass1, mass2)
r[n+1,:,:] = r[n+1,:,:] + dt * v[n+1,:,:] / 2
return r, v
'''
another propagation/numerical integration using leap-from / verlet's type time-reversible symplectic algorithm of convergence rate dt**4
'''
def propagate_system_1(r1_in, r2_in, v1_in, v2_in, mass1, mass2, n_steps, dt):
w0 = pow(2, 1/3)
w1 = 1 / (2 - w0)
w0 = - w0*w1
c = np.array([w1, w0+w1 , w0+w1, w1]) / 2
d = np.array([w1, w0, w1])
r = np.empty((n_steps, 2, 2),dtype=float)
v = np.empty((n_steps, 2, 2),dtype=float)
r[0,0,:] = r1_in
r[0,1,:] = r2_in
v[0,0,:] = v1_in
v[0,1,:] = v2_in
for n in range(n_steps-1):
r[n+1,:,:] = r[n,:,:]
v[n+1,:,:] = v[n,:,:]
for i in range(3):
r[n+1,:,:] = r[n+1,:,:] + c[i]*dt * v[n+1,:,:]
v[n+1,:,:] = v[n+1,:,:] + d[i]*dt * g_accel(r[n+1,:,:], mass1, mass2)
r[n+1,:,:] = r[n+1,:,:] + c[3]*dt * v[n+1,:,:]
return r, v
'''
Test simulation:
Initial conditions:
'''
mass1 = 2
mass2 = 1
initial_position1 = np.array([-1,0])
initial_position2 = np.array([1.3,0])
initial_velocity1 = np.array([0, 0.5])
initial_velocity2 = - mass1 * initial_velocity1 / mass2
'''
Integration parameters:
'''
t_step = 0.05
N = 20000
'''
System time-propagation (numerical integration):
'''
r, v = propagate_system_1(initial_position1,
initial_position2,
initial_velocity1,
initial_velocity2,
mass1, mass2,
N, t_step)
'''
reducing the frequency of time-measurements for faster simulation speed
'''
r = r[np.arange(1, N, 10), :, :]
'''
animation plot of time-evolution:
'''
plt.style.use('seaborn-whitegrid')
fig = plt.figure()
ax = plt.axes()
ax.set_aspect('equal')
ax.set_xlim(-12, 12)
ax.set_ylim(-12, 12)
line = np.empty(2, dtype=type(ax.plot(r[0, 0, 0], r[0, 0, 1])))
for point in range(2):
line[point], = ax.plot( r[0, point, 0], r[0, point, 1] )
def animate(i):
'''
update plot
'''
for p in range(2):
line[p].set_xdata(r[0:i, p, 0])
line[p].set_ydata(r[0:i, p, 1])
return line
intervals = 50
frames = N
frames = int(frames)
anim = FuncAnimation(fig, animate, frames=frames, interval=intervals)
plt.show()
```
