Integrations for Two Body Problem in python

I am trying to create as the first stage of a larger project goal a two body Sim in python with pygame to display on screen the objects.

My main problem at the moment is that the orbiting satellite is oscaliting around the target planet when it should be in a stable 320km circular orbit.

I have made four different functions for four different integrations. Euler, Leapfrog, Verlet and RK4.I'd expect Euler & leapfrog to have some inaccuracy but not RK4 or Verlet. My calculus knowledge is limited so i'd like some extra eyes to check.

They are as follows:

def leapfrog_integration(satellite, planet, dt): #most accurate under 5 orbits
    # Update velocity by half-step
    satellite.velocity += 0.5 * satellite.acceleration * dt
    # Update position
    satellite.position += (satellite.velocity / 1000)
    # Calculate new acceleration
    satellite.acceleration = satellite.acceleration_due_to_gravity(planet)
    # Update velocity by another half-step
    satellite.velocity += 0.5 * satellite.acceleration * dt

def euler_integration(satellite, planet, dt):
   # satellite.accleration = satellite.calculate_gravity(planet)
    satellite.acceleration = satellite.acceleration_due_to_gravity(planet)
    satellite.velocity += satellite.acceleration * dt
    satellite.position += (satellite.velocity / 1000) #convert into kilometers
    
def verlet_integration(satellite, planet, dt):
    acc_c = (satellite.acceleration_due_to_gravity(planet) / 1000)#convert to km/s
    satellite.velocity = (satellite.position - satellite.previous_position)
    new_pos = 2 * satellite.position - satellite.previous_position + (acc_c * dt) 
    satellite.previous_position = satellite.position #km
    satellite.position = new_pos #km
    satellite.velocity = (satellite.position - satellite.previous_position)

def rk4_integration(satellite, planet, dt):# need to resolve the conversion to km for position. If i remove the DT from the kx_r then it's the excat same as Verlet and Euler
    def get_acceleration(position, velocity):
        temp_mass = PointMass(position,satellite.mass,satellite.radius,satellite.colour,(velocity), np.zeros_like(satellite.acceleration))
        return temp_mass.acceleration_due_to_gravity(planet)
    
    k1_v = dt * get_acceleration(satellite.position, (satellite.velocity))
    k1_r = (satellite.velocity / 1000)
    
    k2_v = dt * get_acceleration(satellite.position + 0.5 * k1_r, (satellite.velocity) + 0.5 * k1_v)
    k2_r = (satellite.velocity + 0.5 * k1_v) / 1000
    
    k3_v = dt * get_acceleration(satellite.position + 0.5 * k2_r, (satellite.velocity) + 0.5 * k2_v)
    k3_r = (satellite.velocity + 0.5 * k2_v) / 1000
    
    k4_v = dt * get_acceleration(satellite.position + 0.5 * k3_r, (satellite.velocity) + 0.5 * k3_v)
    k4_r = (satellite.velocity + 0.5 * k3_v) / 1000
    
    satellite.position +=(k1_r + 2*k2_r + 2*k3_r + k4_r) / 6 
    satellite.velocity +=(k1_v + 2*k2_v + 2*k3_v + k4_v) / 6

To give you the class for the pointmasses:

class PointMass:
    def __init__(self, position, mass, radius, colour, velocity, acceleration):
        self.position = np.array(position) #in KM
        self.mass = mass #in Kilograms
        self.radius = radius #in meters
        self.colour = colour
        self.velocity = np.array(velocity) #in m/s
        self.acceleration = np.array(acceleration) #in m/s per second
        self.gForce = None #This is in Newtons
        self.previous_position = self.position - (self.velocity / 1000)  # Initialize previous position for Verlet integration

def acceleration_due_to_gravity(self,other):
        dVector = self.position - other.position # distance vector from self to the other point mass in pixels(km)
        distance_km = np.linalg.norm(dVector)  #Compute the Euclidean distance in km

        distance_m = distance_km * 1000 + other.radius #the distance including the radius to the centre in meters
        unit_vector = (dVector / distance_km) #the unit vector for the direction of the force
        acceleration_magnitude = -Constants.mu_earth / distance_m**2
        return acceleration_magnitude * (unit_vector * 1000) #Return the acceleration vector by multiplying the magnitude with the unit vector(converted to meters)
        #the returned acceleration vector is in m/s

With the intial conditions being:

planet = PointMass(
    position=[600.0,400.0,0.0],
    mass=5.9722e24,
    radius=6.371e6,
    velocity=[0.0,0.0,0.0], #in m/s
    acceleration=[0.0,0.0,0.0], #in m/s^2
    colour=[125,100,100]
)

satellite = PointMass(
    position=[280.0,400.0,0.0],
    mass=100,
    radius=1,
    velocity=[0.0,7718.0,0.0], #need an intial velocity or else it'll jsut fall to the central mass
    acceleration=[1.0,0.0,0.0],#added an non-zero acceleration jsut to make sure there's no issues with the integrations.
    colour=[255,255,255]
)

The code is a bit messy since I'm constantly trying new things and haven't really formatted it.

Now I'm not 100% the problem lays with the integrations but I'd thought i would start there and go form that. It could easily be my gravity force function, the delta time being applied wrong, unit conversion causing problems or how it's being placed on screen that is the problem.

Github Link

Above is the repo link just to keep things in here tidy, since it could be anything. Screen shots of tests etc.

Now it could be just expect behaviour from python with loosing accuracy with floating point numbers but I'd hit X on that as it's more likely my poor calculus knowledge.

After extensive testing, I got a DT of 0.021 as the best result for delta time. Leapfrog seems the most accurate for either low orbit counts or high orbit counts.

Euler and Verlet seem about even under 100 orbits where Euler goes a little buts and RK4 seems unstable in the fact it slowly slows down and thus the orbits get larger and larger.

I have the conversion from Meters to KM in different places for each integration since they wouldn't work right depending on where i put it.

I had to remove the DT from some of the integration parts since if i left them in the object would just spiral slightly and fall to the central mass.

There are a lot of separate problems with your code - some of them cancelling out! The biggest problem from the perspective of Stack Overflow is that you didn't minimise the example: there is a lot of extraneous material in your code which isn't relevant.

Your RK4 is simply wrong. It is missing the necessary dt* in the position updates. More minor is that k4_v should depend on just k3_r etc rather than 0.5*k3_r. Similarly for k4_r.

As Andrew Yim has pointed out, the displacement vector is in the opposite direction to that intended. However, this error is precisely cancelled by errors in the acceleration routines and, indeed, on which object is doing the gravitating here (see return planet.acceleration_due_to_gravity(temp_mass) and not the other way round).

For spherically-symmetric bodies, the "distance" should be centre-to-centre and should not be adjusted for radius.

Your use of km for position but m for velocity and acceleration is a nightmare to sort out!

Your initial velocity will determine whether the final orbit is an ellipse or a circle. What you have at the moment gives an ellipse (which is OK). Balance centripetal force with gravitational force if you want the correct velocity for a circular orbit.

Here is a condensed version of your code as a single file. It omits the unnecessary information. The final orbit is closed.

import pygame
import sys
import time
import numpy as np
import math


class Constants:
    gravitational_constant = 6.6740e-11
    mu_earth = 3.986004418e14  # in m^3 s^-2

class PointMass:
    def __init__(self, position, mass, radius, colour, velocity, acceleration):
        self.position = np.array(position)          # in KM (!OUCH!)
        self.mass = mass                            # in kg
        self.radius = radius                        # in meters (bizarre, given that position is in km!)
        self.colour = colour
        self.velocity = np.array(velocity)          # in m/s
        self.acceleration = np.array(acceleration)  # in m/s^2

    def acceleration_due_to_gravity( self, other ):
        dVector = self.position - other.position
        distance_km = np.linalg.norm(dVector)
        distance_m = distance_km * 1000
        unit_vector = (dVector / distance_km)
        acceleration_magnitude = Constants.mu_earth / distance_m**2
        return acceleration_magnitude * unit_vector


# Initialize Pygame
pygame.init()
window_size = (1200, 800)
window = pygame.display.set_mode(window_size)
pygame.display.set_caption("Two Body Simulation")
font = pygame.font.SysFont(None, 24)

planet = PointMass(
    position=[600.0,400.0,0.0],
    mass=5.9722e24,
    radius=6.371e6,
    velocity=[0.0,0.0,0.0],    
    acceleration=[0.0,0.0,0.0],
    colour=[125,100,100]
)

satellite = PointMass(
    position=[280.0,400.0,0.0],
    mass=100,
    radius=1,
    velocity=[0.0,7718.0*2,0.0],
    acceleration=[0.0,0.0,0.0],
    colour=[255,255,255]
)


grid_spacing = 50
# Function to draw the grid
def draw_grid():
    for x in range(0, window_size[0], grid_spacing):
        pygame.draw.line(window, [255,255,255], (x, 0), (x, window_size[1]))
    for y in range(0, window_size[1], grid_spacing):
        pygame.draw.line(window, [255,255,255], (0, y), (window_size[0], y))

def rk4_integration(satellite, planet, dt):
    def get_acceleration(position, velocity):
        temp_mass = PointMass(position,satellite.mass,satellite.radius,satellite.colour,(velocity), np.zeros_like(satellite.acceleration))
        return planet.acceleration_due_to_gravity(temp_mass)
    
    k1_v = dt * get_acceleration(satellite.position, (satellite.velocity))
    k1_r = dt * (satellite.velocity / 1000)
    
    k2_v = dt * get_acceleration(satellite.position + 0.5 * k1_r, (satellite.velocity) + 0.5 * k1_v)
    k2_r = dt * (satellite.velocity + 0.5 * k1_v) / 1000
    
    k3_v = dt * get_acceleration(satellite.position + 0.5 * k2_r, (satellite.velocity) + 0.5 * k2_v)
    k3_r = dt * (satellite.velocity + 0.5 * k2_v) / 1000
    
    k4_v = dt * get_acceleration(satellite.position + k3_r, (satellite.velocity) + k3_v)
    k4_r = dt * (satellite.velocity + k3_v) / 1000
    
    satellite.position +=(k1_r + 2*k2_r + 2*k3_r + k4_r) / 6 
    satellite.velocity +=(k1_v + 2*k2_v + 2*k3_v + k4_v) / 6
        

#actual drawing of objects on screen
circle_radius = 10
c1_pos = (planet.position[0],planet.position[1])
c2_pos = (int(satellite.position[0]), int(satellite.position[1]))
   
clock = pygame.time.Clock()
fps = 50
positions = []


# Main loop
running = True
while running:
    # Event handling
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False
            
    dt = clock.tick(fps) / 1000.0
    
    rk4_integration(satellite, planet, dt)
    positions.append((int(satellite.position[0]), satellite.position[1]))

    #update the position of the drawn object on screen
    c1_pos = (planet.position[0],planet.position[1])
    c2_pos = (int(satellite.position[0]), int(satellite.position[1]))
   
    #Draw
    window.fill((0,0,0))
    pygame.draw.circle(window, planet.colour, c1_pos, 25)
    pygame.draw.circle(window, satellite.colour, c2_pos, 5)

    if len(positions) > 1:
        for i in range(len(positions) - 1):
            pygame.draw.line(window,(205,205,205),positions[i], positions[i + 1], 1)

    draw_grid()
    
    #Refresh display
    pygame.display.flip()
    clock.tick(fps)

pygame.quit()
sys.exit()

Two orbits coincide: