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.
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:
- Is it possible to call setter as a method?
- Accessing a Python traceback from the C API
- Scraping hierarchical website in a specific category
- How to offset a marker in matplotlib?
- How to access current user in Django class based view
- Is there a way to define list(obj) method on a user defined class in python?
- Deserializing Json list into class objects not able to extract attrributes
- What is the pythonic way to detect the last element in a 'for' loop?
- How to get integer values from a string in Python?
- Optimize this Python code that involves matrix inversion
- how to read a file with dictionary?
- How to extract only audio from downloading video? Python yt-dlp
- How to solve Invalid element(s) received for the 'color' property?
- Plotly Value error - Invalid property for colour
- no module named distutils....but distutils installed?
- how can i change the venv of my current shell from within a python script
- What is the difference between "cast(...)" and "... # type: ..."?
- getting the file size with pycurl
- export conda environment without prefix variable which shows local path to executable
- How can I fix datetime error when trying to use the now syntax?
- Python 3 asyncio with aioboto3 seems sequential
- How can I set OpenCV to use origin as bottom-left corner of image?
- Can I use sleep in the main thread without blocking other threads?
- Getting number of elements in an iterator in Python
- Merging tuples if they have one common element
- How to specify the response_model in FastAPI on a non-default return?
- PyCharm doesn’t see files in project
- How to tune PID controller on non-linear model (here: diving of a body)
- Create list of factorials using list comprehension
- Rewriting URL request function to satisfy GitHub CodeQL server side request forgery (SSRF) warning