Question: 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

Answer: 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:  don't do 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