Gravity calculation using pygame

Question

import pygame  
from pygame import *  
import math  
class orb:  
    def __init__(self, mass=1, s0=[0,0], v0=[0,0], a=[0,0], r=10, color=(255, 255, 255), fix = False, F=0, Fx=0, Fy=0):  
        self.a= a  
        self.mass=mass  
        self.r=r  
        self.t = 0  
        self.v0=v0  
        self.v=[self.v0[0], self.v0[1]]  
        self.s0=s0  
        self.s=[self.s0[0], self.s0[1]]  
        self.color=color  
        self.fix = fix  
        self.F=F  
        self.Fx = Fx  
        self.Fy = Fy  

calculating the position of the earth

    def comPos(self, dt = 0.01):
        if not self.fix:  
            self.v[0] = self.v0[0] + self.a[0] * dt
            self.v[1] = self.v0[1] + self.a[1] * dt
            ds_x = self.v0[0] * dt + (1/2) * self.a[0] * dt**2
            ds_y = self.v0[1] * dt + (1/2) * self.a[1] * dt**2
            self.s[0] = self.s0[0] + ds_x
            self.s[1] = self.s0[1] + ds_y
            self.v0[0] = self.v[0]
            self.v0[1] = self.v[1]
            self.s0[0] = self.s[0]
            self.s0[1] = self.s[1]
            self.t += dt

calculating the gravity of the earth

    def comForce(self, other, G = 10000.0, F = 0.0, dist_cr = 0.0, Fx = 0.0, Fy = 0.0):
        self.F = G * self.mass * other.mass / ( (self.s[0]- other.s[0])**2 + (self.s[1] - other.s[1])**2 )
        self.Fx = (self.F * math.fabs(math.cos(math.atan((other.s[1] - self.s[1]) / (other.s[0] - self.s[0])))))
        self.Fy = (self.F * math.fabs(math.sin(math.atan((other.s[1] - self.s[1]) / (other.s[0] - self.s[0])))))
        if ((other.s[0] - self.s[0])**2 + (other.s[1] - self.s[1])**2)**(1/2) >= dist_cr:
            if other.s[0] - self.s[0] == 0:
                self.a[0] = 0
            else:
                self.a[0] = (( (other.s[0] - self.s[0]) / math.fabs(other.s[0] - self.s[0])) * (self.Fx/self.mass))
            if other.s[1] - self.s[1] == 0:
                self.a[1] = 0
            else:
                self.a[1] = (( (other.s[1] - self.s[1]) / math.fabs(other.s[1] - self.s[1])) * self.Fy / self.mass)

Displaying on the screen

screen = pygame.display.set_mode((800, 600))
clock = pygame.time.Clock()
earth = orb(s0=[150, 300], v0=[0,60], r=10, mass= 3)
sun = orb(s0=[400, 300], r=70, color = [255, 0, 0], fix=True, mass=100)
orb = [sun, earth]

running = True
while running:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False
    screen.fill((0, 0, 0))
    pygame.draw.circle(screen, earth.color, earth.s, earth.r, 0)
    pygame.draw.circle(screen, sun.color, sun.s, sun.r, 0)
    earth.comForce(sun)
    earth.comPos()
    pygame.display.update()
    clock.tick(720)  

I'm going to simulate the Earth that goes around the sun. But as time goes on, the orbit changes little by little. I tried to find the reason, but I don't know the cause no matter how hard I try. Can you tell me why the trajectory changes?

I want the trajectory to remain constant without changing over time.

Answer 1

• You can use float64 for higher accuracies.
• Note that using the math fabs, sin, cos, etc., reduces the accuracies.
• Also, you can use small dt.
• There are several other ways to make it accurate.

Code:

import pygame
import numpy as np


class Orb:
    def __init__(self, mass=1, s0=[0, 0], v0=[0, 0], r=10, color=(255, 255, 255), fix=False):
        self.mass = np.float64(mass)
        self.r = r
        self.s = np.array(s0, dtype=np.float64)
        self.v = np.array(v0, dtype=np.float64)
        self.a = np.array([0.0, 0.0], dtype=np.float64)
        self.color = color
        self.fix = fix

    def comForce(self, other, G=np.float64(10000.0)):
        dx = other.s[0] - self.s[0]
        dy = other.s[1] - self.s[1]
        dist_sq = dx ** 2 + dy ** 2
        dist = np.sqrt(dist_sq)
        F = G * self.mass * other.mass / dist_sq
        Fx = F * dx / dist
        Fy = F * dy / dist
        return np.array([Fx, Fy], dtype=np.float64)

    def update(self, other, dt=np.float64(0.001)):
        if not self.fix:
            force = self.comForce(other)
            self.a = force / self.mass
            self.s += self.v * dt + 0.5 * self.a * dt**2
            new_force = self.comForce(other)
            new_a = new_force / self.mass
            self.v += 0.5 * (self.a + new_a) * dt


pygame.init()
screen = pygame.display.set_mode((900, 600))
clock = pygame.time.Clock()

earth = Orb(mass=np.float64(2), s0=[150, 300], v0=[0, 60], r=10, color=(0, 0, 255))
sun = Orb(mass=np.float64(100), s0=[420, 300], v0=[0, 0], r=50, color=(255, 0, 0), fix=True)
orbs = [sun, earth]

running = True
frames = 1000
while running:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False

    for _ in range(frames):
        earth.update(sun, dt=np.float64(0.001))

    screen.fill((0, 0, 0))
    for orb in orbs:
        pygame.draw.circle(screen, orb.color, (int(orb.s[0]), int(orb.s[1])), orb.r)

    pygame.display.flip()
    clock.tick(60)

pygame.quit()