I am trying to simulate a somewhat realistic program where Earth and the moon can interact gravitationally with each other. Now the problem is that the moon keeps on spiraling towards Earth and I don't understand why.
This is my code:
from math import sin,cos,sqrt,atan2,pi
import pygame
pygame.init()
class Planet:
dt = 1/100
G = 6.67428e-11 #G constant
scale = 1/(1409466.667) #1 m = 1/1409466.667 pixels
def __init__(self,x=0,y=0,radius=0,color=(0,0,0),mass=0,vx=0,vy=0):
self.x = x #x-coordinate pygame-window
self.y = y #y-coordinate pygame-window
self.radius = radius
self.color = color
self.mass = mass
self.vx = vx #velocity in the x axis
self.vy = vy #velocity in the y axis
def draw(self,screen):
pygame.draw.circle(screen, self.color, (self.x, self.y), self.radius)
def orbit(self,trace):
pygame.draw.rect(trace, self.color, (self.x, self.y, 2, 2))
def update_vel(self,Fnx,Fny):
ax = Fnx/self.mass #Calculates acceleration in x- and y-axis for body 1.
ay = Fny/self.mass
self.vx -= ((ax * Planet.dt)/Planet.scale)
self.vy -= ((ay * Planet.dt)/Planet.scale)
self.update_pos()
def update_pos(self):
self.x += ((self.vx * Planet.dt)) #changes position considering each body's velocity.
self.y += ((self.vy * Planet.dt))
def move(self,body):
dx = (self.x - body.x) #Calculates difference in x- and y-axis between the bodies
dy = (self.y - body.y)
r = (sqrt((dy**2)+(dx**2))) #Calculates the distance between the bodies
angle = atan2(dy, dx) #Calculates the angle between the bodies with atan2!
if r < self.radius: #Checks if the distance between the bodies is less than the radius of the bodies. Uses then Gauss gravitational law to calculate force.
F = 4/3 * pi * r
Fx = cos(angle) * F
Fy = sin(angle) * F
else:
F = (Planet.G*self.mass*body.mass)/((r/Planet.scale)**2) #Newtons gravitational formula.
Fx = cos(angle) * F
Fy = sin(angle) * F
return Fx,Fy
def motion():
for i in range(0,len(bodies)):
Fnx = 0 #net force
Fny = 0
for j in range(0,len(bodies)):
if bodies[i] != bodies[j]:
Fnx += (bodies[i].move(bodies[j]))[0]
Fny += (bodies[i].move(bodies[j]))[1]
elif bodies[i] == bodies[j]:
continue
bodies[i].update_vel(Fnx,Fny)
bodies[i].draw(screen)
bodies[i].orbit(trace)
Fnx,Fny=0,0
screen = pygame.display.set_mode([900,650]) #width - height
trace = pygame.Surface((900, 650))
pygame.display.set_caption("Moon simulation")
FPS = 150 #how quickly/frames per second our game should update.
earth = Planet(450,325,30,(0,0,255),5.97219*10**(24),-24.947719394204714/2) #450= xpos,325=ypos,30=radius
luna = Planet(450,(575/11),10,(128,128,128),7.349*10**(22),1023)
bodies = [earth,luna]
running = True
clock = pygame.time.Clock()
while running: #if user clicks close window
clock.tick(FPS)
for event in pygame.event.get():
if event.type == pygame.QUIT:
running = False
screen.fill((0,0,0))
pygame.Surface.blit(screen, trace, (0, 0))
motion()
pygame.display.flip()
pygame.quit()
Once I get the earth-moon system to work I would like to expand on it and try having three bodies (the reason why there is so much of what otherwise would be "unnecessary" code)
I am open to suggestions or/and advices! Thanks
You need to compute the forces in one go, not body-move by body-move.
Instead of computing the forces during the update loop, while shifting positions:
def motion():
for i in range(0,len(bodies)):
Fnx = 0 #net force
Fny = 0
for j in range(0,len(bodies)):
if bodies[i] != bodies[j]:
Fnx += (bodies[i].move(bodies[j]))[0]
Fny += (bodies[i].move(bodies[j]))[1]
elif bodies[i] == bodies[j]:
continue
bodies[i].update_vel(Fnx,Fny)
bodies[i].draw(screen)
bodies[i].orbit(trace)
Fnx,Fny=0,0
settle the forces upfront:
def motion():
force = [ (
sum([(bodies[i].move(bodies[j]))[0] for j in range(0, len(bodies)) if i != j ]),
sum([(bodies[i].move(bodies[j]))[1] for j in range(0, len(bodies)) if i != j ])
) for i in range(0,len(bodies)) ]
for i in range(0,len(bodies)):
Fnx = force[i][0]
Fny = force[i][1]
bodies[i].update_vel(Fnx,Fny)
bodies[i].draw(screen)
bodies[i].orbit(trace)
Fnx,Fny=0,0
(I don't write in Python usually, so the style is not perfect.)
The following text is from a previous answer. It may be helpful, but it is not required to solve the problem asked; you may stop reading here.
You can reduce numeric truncation errors further with more elaborated methods, like Runge-Kutta. For that:
update_vel and update_pos one-after the other, instead, try to write an update_state method which combines both simultaneously; important is that the left hand side of the equations is either the delta or the new state, the right hand side of the equations is the old state, exclusively (higher-order Runge-Kutta will have some intermediate states, fractions of Planet.dt)If Runge-Kutta is too heavy to start with, consider MacCormack or Lax-Wendroff.
For a Lax-Wendroff'ish way, instead of:
def update_vel(self,Fnx,Fny):
ax = Fnx/self.mass
ay = Fny/self.mass
self.vx -= ((ax * Planet.dt)/Planet.scale)
self.vy -= ((ay * Planet.dt)/Planet.scale)
self.update_pos()
def update_pos(self):
self.x += ((self.vx * Planet.dt))
self.y += ((self.vy * Planet.dt))
try this:
def update_state(self,Fnx,Fny):
ax = Fnx/self.mass
ay = Fny/self.mass
self.x += (self.vx * Planet.dt) - (ax/Planet.scale) * Planet.dt**2
self.y += (self.vy * Planet.dt) - (ay/Planet.scale) * Planet.dt**2
self.vx -= (ax/Planet.scale) * Planet.dt
self.vy -= (ay/Planet.scale) * Planet.dt