Sign-error in force calculation for Jupiter

[heinzs57]

Hi there, I've been playing around with your very nice solver--thanks for putting those into the public domain!

I noticed that there's a sign error for the Jupiter force calculation. It probably doesn't matter for Jupiter mass planets, but when integrating planets closer to each other in mass, then that error will become important. The force function should be

def force(r,planet,ro,vo):
if planet == 'earth':
return force_es(r) + force_ej(r,ro)
if planet == 'jupiter':
return force_js(r) + force_ej(r,ro)

Also, it would be beneficial to forward-integrate the motions of both planets in the Runge-Kutta solver symmetrically for cases where the planets are closer to each other and/or similar in mass. Here's a way to do that:

def RK4Solver(t,re,ve,rj,vj,h):
k11e = dr_dt(t,re,ve,'earth',rj,vj)
k21e = dv_dt(t,re,ve,'earth',rj,vj)
k11j = dr_dt(t,rj,vj,'jupiter',re,ve)
k21j = dv_dt(t,rj,vj,'jupiter',re,ve)

r1e=re+k11e*h*0.5
v1e=ve+k21e*h*0.5
r1j=rj+k11j*h*0.5
v1j=vj+k21j*h*0.5

k12e = dr_dt(t + 0.5*h,r1e,v1e,'earth',r1j,v1j)
k22e = dv_dt(t + 0.5*h,r1e,v1e,'earth',r1j,v1j)
k12j = dr_dt(t + 0.5*h,r1j,v1j,'jupiter',r1e,v1e)
k22j = dv_dt(t + 0.5*h,r1j,v1j,'jupiter',r1e,v1e)

r2e=re+k12e*h*0.5
v2e=ve+k22e*h*0.5
r2j=rj+k12j*h*0.5
v2j=vj+k22j*h*0.5

k13e = dr_dt(t + 0.5*h,r2e,v2e,'earth',r2j,v2j)
k23e = dv_dt(t + 0.5*h,r2e,v2e,'earth',r2j,v2j)
k13j = dr_dt(t + 0.5*h,r2j,v2j,'jupiter',r2e,v2e)
k23j = dv_dt(t + 0.5*h,r2j,v2j,'jupiter',r2e,v2e)

r3e=re+k13e*h
v3e=ve+k23e*h
r3j=rj+k13j*h
v3j=vj+k23j*h

k14e = dr_dt(t + h,r3e,v3e,'earth',r3j,v3j)
k24e = dv_dt(t + h,r3e,v3e,'earth',r3j,v3j)
k14j = dr_dt(t + h,r3j,v3j,'earth',r3e,v3e)
k24j = dv_dt(t + h,r3j,v3j,'earth',r3e,v3e)

y0e = re + h * (k11e + 2.*k12e + 2.*k13e + k14e) / 6.
y1e = ve + h * (k21e + 2.*k22e + 2.*k23e + k24e) / 6.
y0j = rj + h * (k11j + 2.*k12j + 2.*k13j + k14j) / 6.
y1j = vj + h * (k21j + 2.*k22j + 2.*k23j + k24j) / 6.

ze = np.zeros([2,2])
ze = [y0e, y1e]
zj = np.zeros([2,2])
zj = [y0j, y1j]

return ze,zj:

Finally, you can avoid the case distinction in the gravity term if you write the force terms as (e.g.):

def force_es(r):
return = -rGGMe*Ms/(r[0]**2+r[1]2)(3./2.)