Gravity simulation of the solar system

planets.py

import pygame


# units are kg, km and s


class Vector:
    def __init__(self, x, y, z):
        self.x = x
        self.y = y
        self.z = z

    def __add__(self, other):
        x = self.x + other.x
        y = self.y + other.y
        z = self.z + other.z
        return Vector(x, y, z)

    def __sub__(self, other):
        x = self.x - other.x
        y = self.y - other.y
        z = self.z - other.z
        return Vector(x, y, z)

    def __mul__(self, scalar):
        x = self.x * scalar
        y = self.y * scalar
        z = self.z * scalar
        return Vector(x, y, z)

    def __rmul__(self, scalar):
        return self * scalar

    def __truediv__(self, scalar):
        x = self.x / scalar
        y = self.y / scalar
        z = self.z / scalar
        return Vector(x, y, z)

    def __abs__(self):
        x = self.x
        y = self.y
        z = self.z
        return (x ** 2 + y ** 2 + z ** 2) ** 0.5


names = []
m = []  # masses
p = []  # positions
v = []  # velocities
object_count = 0


def add_object(name, mass, x, y, z, vx, vy, vz):
    global object_count
    names.append(name)
    m.append(mass)
    p.append(Vector(x, y, z))
    v.append(Vector(vx, vy, vz))
    object_count += 1


add_object(
    "Sonne",
    1.9885e30,
    -1.067706805380953E+06,
    -4.182752718194473E+05,
    3.086181725476820E+04,
    9.312571926520472E-03,
    -1.282475570794162E-02,
    -1.633507186350417E-04
)

add_object(
    "Merkur",
    3.3011e23,
    -2.052943316123468E+07,
    -6.733155053534345E+07,
    -3.648992526494771E+06,
    3.700430442920571E+01,
    -1.117724068132644E+01,
    -4.307791469376854E+00
)

add_object(
    "Venus",
    4.8675e24,
    -1.085242008575715E+08,
    -5.303290247691983E+06,
    6.166496116973171E+06,
    1.391218601189967E+00,
    -3.515311993215464E+01,
    -5.602056890007159E-01
)

add_object(
    "Erde",
    5.97237e24,
    -2.756674048281145E+07,
    1.442790215207299E+08,
    3.025066782881320E+04,
    -2.978494749851088E+01,
    -5.482119695478543E+00,
    1.843295986780902E-05
)

add_object(
    "Mars",
    6.4171e23,
    2.069805421180782E+08,
    -2.425286839967330E+06,
    -5.125451306151839E+06,
    1.171984953383225E+00,
    2.628323979722633E+01,
    5.221336703953376E-01
)

add_object(
    "Jupiter",
    1.8982e27,
    5.974999178516835E+08,
    4.391864046763535E+08,
    -1.519599985573271E+07,
    -7.900547720245487E+00,
    1.114339277065934E+01,
    1.307023308637314E-01
)

add_object(
    "Saturn",
    5.6834e26,
    9.573177603675805E+08,
    9.824380490886426E+08,
    -5.518214225187147E+07,
    -7.421900391331049E+00,
    6.723931000090287E+00,
    1.775749562901709E-01
)

add_object(
    "Uranus",
    8.6810e25,
    2.157907331374278E+09,
    -2.055043522898531E+09,
    -3.559460196107912E+07,
    4.646586371441413E+00,
    4.614774391283166E+00,
    -4.308124104538469E-02
)

add_object(
    "Neptun",
    1.02413e26,
    2.513978764682338E+09,
    -3.739132814439580E+09,
    1.906307923109627E+07,
    4.474587749877043E+00,
    3.063155425183056E+00,
    -1.664119935880864E-01
)


G = 6.67430e-20


def get_accelerations(p):
    a = [Vector(0, 0, 0) for i in range(object_count)]
    for i in range(object_count):
        for j in range(i + 1, object_count):
            r = p[j] - p[i]
            F = G * m[i] * m[j] * r / abs(r) ** 3
            a[i] += F
            a[j] -= F
        a[i] /= m[i]
    return a


Delta_t = 3600


def advance():
    # first set of intermediate values
    p1 = p
    v1 = v
    a1 = get_accelerations(p1)
    # second set of intermediate values
    p2 = object_count * [None]
    v2 = object_count * [None]
    for i in range(object_count):
        p2[i] = p[i] + v1[i] * Delta_t / 2
        v2[i] = v[i] + a1[i] * Delta_t / 2
    a2 = get_accelerations(p2)
    # third set of intermediate values
    p3 = object_count * [None]
    v3 = object_count * [None]
    for i in range(object_count):
        p3[i] = p[i] + v2[i] * Delta_t / 2
        v3[i] = v[i] + a2[i] * Delta_t / 2
    a3 = get_accelerations(p3)
    # fourth set of intermediate values
    p4 = object_count * [None]
    v4 = object_count * [None]
    for i in range(object_count):
        p4[i] = p[i] + v3[i] * Delta_t
        v4[i] = v[i] + a3[i] * Delta_t
    a4 = get_accelerations(p4)
    # add weighted averages
    for i in range(object_count):
        p[i] += (v1[i] + 2 * v2[i] + 2 * v3[i] + v4[i]) / 6 * Delta_t
        v[i] += (a1[i] + 2 * a2[i] + 2 * a3[i] + a4[i]) / 6 * Delta_t


pygame.init()
screen = pygame.display.set_mode((800, 800))
font = pygame.font.SysFont("DejaVu Sans", 10)
history = []
while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            exit()
    history.append(tuple(p))
    if len(history) > 256:
        history.pop(0)
    screen.fill((0, 0, 0))
    for i in range(object_count):
        for j in range(len(history)):
            x = int(round(400 + history[j][i].x / 2e7))
            y = int(round(400 - history[j][i].y / 2e7))
            color = (j / 2, 0, j)
            pygame.draw.circle(screen, color, (rounx, y), 1)
        text = font.render(names[i], True, (255, 255, 255))
        screen.blit(text, (x, y))
    pygame.display.flip()
    for i in range(24):
        advance()