Gravity Problems

Question

I've been trying to make a code with pygame to simulate simple gravity. At the moment, there is only one object (HOM) which is orbiting the sun. However, for reasons unknown to me, whenever I run the code, HOM travels round the sun in an orbit at the start, but then accelerates away from the sun when it reaches ~135 degrees from vertical.

Does anyone know why this is happening and how I can fix it? I have been printing some variables to try and source the problem, but have had no luck so far.

Code:

import pygame,sys,time
from math import *

screen=pygame.display.set_mode((800,600))

G = 5

class Object: #Just an object, like a moon or planet
    def __init__(self,mass,init_cds,init_vel,orbit_obj='Sun'):
        self.mass = mass
        self.cds = init_cds
        self.velocity = init_vel
        self.accel = [0,0]
        self.angle = 0
        self.orb_obj = orbit_obj

    def display(self):
        int_cds = (round(self.cds[0]),round(self.cds[1]))#Stores its co-ordinates as floats, has to convert to integers for draw function
        pygame.draw.circle(screen,(255,0,0),int_cds,10)

    def calc_gravity(self):
        if self.orb_obj == 'Sun':
            c_x,c_y = 400,300
            c_mass = 10000
        else:
            c_x,c_y = self.orb_obj.cds
            c_mass = self.orb_obj.mass
        d_x = self.cds[0]-c_x
        d_y = self.cds[1]-c_y
        dist = sqrt(d_x**2+d_y**2) #Find direct distance
        angle = atan(d_x/d_y) #Find angle
        print(d_x,d_y)
        print(dist,degrees(angle))
        if dist == 0:
            acc = 0
        else:
            acc = G*c_mass/(dist**2) #F=G(Mm)/r^2, a=F/m -> a=GM/r^2
        print(acc)
        acc_x = acc*sin(angle) #Convert acceleration from magnitude+angle -> x and y components
        acc_y = acc*cos(angle)
        self.accel = [acc_x,acc_y]
        print(self.accel)
        self.velocity = [self.velocity[0]+self.accel[0],self.velocity[1]+self.accel[1]] #Add acceleration to velocity
        print(self.velocity)
        self.cds = (self.cds[0]+self.velocity[0],self.cds[1]+self.velocity[1]) #Change co-ordinates by velocity
        print(self.cds)
        print('-------------------') #For seperating each run of the function when printing variables

HOM = Object(1000000,(400,100),[10,0]) #The problem planet

clock = pygame.time.Clock()

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            sys.exit()

    screen.fill((0,0,0))
    pygame.draw.circle(screen,(255,255,0),(400,300),25)
    HOM.display()
    HOM.calc_gravity()

    clock.tick(30)

    pygame.display.flip()

Answer 1

Your main issue has to do with this line:

angle = atan(d_x/d_y) #Find angle

The atan function is very limited in its ability to compute angles because it can't tell the signs of the coordinates you combined in your division. For instance, it will give the same result for atan(1/1) and atan(-1/-1), since both divisions compute the same slope (1).

Instead you should use atan2, and pass the coordinates separately. Since this will let the code see both coordinates, it can pick an angle in the right quadrant of the circle every time.

But there's an even better fix. Instead of computing an angle and then immediately converting it back to a unit vector (by calling sin and cos on it), why not compute the unit vector directly? You already have the original vector's length! Instead of:

acc_x = acc*sin(angle) #Convert acceleration from magnitude+angle -> x and y components
acc_y = acc*cos(angle)

Use:

acc_x = acc * d_x / distance
acc_y = acc * d_y / distance

The d_x / distance and d_y / distance values are the same as the sin and cos values you were getting before (for the angles when they were working correctly), but there's no need for the trigonometry. You can get rid of the line I quoted up top completely!

Note that you might need to reverse the way you're computing d_x and d_y, so that you get a vector that points from the orbiting object towards the object it's orbiting around (instead of pointing the other way, from the center of the orbit towards the orbiting object). I'm not sure if I'm reading your code correctly, but it looks to me like you have it the other way around right now. That means that you were actually getting the wrong results from atan in the cases where your current code was working the way you expected, and the bad behavior (flying off into nowhere) is the code working "correctly" (from a mathematical point of view). Alternatively, you could compute acc to be negative, rather than positive.

As several commenters mentioned, you may have other issues related to your choice of integration algorithm, but those errors are not going to be as large as the main issue with the acceleration angle. They'll crop up as you run your simulation over longer time periods (and try to use larger time steps to make the simulation go faster). Your current algorithm is good enough for an orbit or two, but if you're simulating dozens or hundreds of orbits, you'll start seeing errors accumulate and so you should pick a better integrator.