Интеграция задачи двух тел в Python

Я не знаю, исправит ли это все, но вам, возможно, стоит обратить внимание на то, как вы назначаете dVector в PointMass.acceleration_due_to_gravity.

У вас это определено как self.position - other.position. Это противоположно тому, что говорится в вашем комментарии; обычно, если вам нужен вектор от self до other, вы берете other.position - self.position (см. эту диаграмму).

Это потому, что если бы вам нужен вектор от self до other, вы бы хотели иметь возможность добраться до other, добавив vector к self. То есть вы хотите:

self + vector = other

Который вы можете использовать свойство вычитания, чтобы переписать как vector = other - self.

Изменение порядка означает, что ваш вектор ускорения будет указывать в противоположном направлении, в котором вы хотите, что может испортить движение ваших точечных масс.

TL;DR: переключите назначение dVector на dVector = other.position - self.position и посмотрите, изменится ли что-нибудь.

В вашем коде много отдельных проблем, некоторые из которых устраняются! Самая большая проблема с точки зрения переполнения стека заключается в том, что вы не минимизировали пример: в вашем коде много постороннего материала, который не имеет отношения к делу.

Ваш RK4 просто неправильный. В обновлениях позиций отсутствуют необходимые dt*. Более незначительным является то, что k4_v должно зависеть только от k3_r и т. д., а не от 0.5*k3_r. Аналогично для k4_r.

Как отметил Эндрю Йим, вектор смещения направлен в направлении, противоположном предполагаемому. Однако эта ошибка как раз и компенсируется ошибками в процедурах ускорения и, собственно, тем, какой объект здесь притягивается (см. return planet.acceleration_due_to_gravity(temp_mass), а не наоборот).

Для сферически-симметричных тел «расстояние» должно быть между центрами и не должно корректироваться по радиусу.

Использовать км для определения местоположения, а м для скорости и ускорения — это настоящий кошмар!

Ваша начальная скорость будет определять, будет ли конечная орбита эллипсом или кругом. То, что у вас есть на данный момент, дает эллипс (это нормально). Если вам нужна правильная скорость для круговой орбиты, сбалансируйте центростремительную силу с гравитационной силой.

Вот сокращенная версия вашего кода в виде одного файла. Он опускает ненужную информацию. Последняя орбита закрыта.

import pygame
import sys
import time
import numpy as np
import math


class Constants:
    gravitational_constant = 6.6740e-11
    mu_earth = 3.986004418e14  # in m^3 s^-2

class PointMass:
    def __init__(self, position, mass, radius, colour, velocity, acceleration):
        self.position = np.array(position)          # in KM (!OUCH!)
        self.mass = mass                            # in kg
        self.radius = radius                        # in meters (bizarre, given that position is in km!)
        self.colour = colour
        self.velocity = np.array(velocity)          # in m/s
        self.acceleration = np.array(acceleration)  # in m/s^2

    def acceleration_due_to_gravity( self, other ):
        dVector = self.position - other.position
        distance_km = np.linalg.norm(dVector)
        distance_m = distance_km * 1000
        unit_vector = (dVector / distance_km)
        acceleration_magnitude = Constants.mu_earth / distance_m**2
        return acceleration_magnitude * unit_vector


# Initialize Pygame
pygame.init()
window_size = (1200, 800)
window = pygame.display.set_mode(window_size)
pygame.display.set_caption("Two Body Simulation")
font = pygame.font.SysFont(None, 24)

planet = PointMass(
    position=[600.0,400.0,0.0],
    mass=5.9722e24,
    radius=6.371e6,
    velocity=[0.0,0.0,0.0],    
    acceleration=[0.0,0.0,0.0],
    colour=[125,100,100]
)

satellite = PointMass(
    position=[280.0,400.0,0.0],
    mass=100,
    radius=1,
    velocity=[0.0,7718.0*2,0.0],
    acceleration=[0.0,0.0,0.0],
    colour=[255,255,255]
)


grid_spacing = 50
# Function to draw the grid
def draw_grid():
    for x in range(0, window_size[0], grid_spacing):
        pygame.draw.line(window, [255,255,255], (x, 0), (x, window_size[1]))
    for y in range(0, window_size[1], grid_spacing):
        pygame.draw.line(window, [255,255,255], (0, y), (window_size[0], y))

def rk4_integration(satellite, planet, dt):
    def get_acceleration(position, velocity):
        temp_mass = PointMass(position,satellite.mass,satellite.radius,satellite.colour,(velocity), np.zeros_like(satellite.acceleration))
        return planet.acceleration_due_to_gravity(temp_mass)
    
    k1_v = dt * get_acceleration(satellite.position, (satellite.velocity))
    k1_r = dt * (satellite.velocity / 1000)
    
    k2_v = dt * get_acceleration(satellite.position + 0.5 * k1_r, (satellite.velocity) + 0.5 * k1_v)
    k2_r = dt * (satellite.velocity + 0.5 * k1_v) / 1000
    
    k3_v = dt * get_acceleration(satellite.position + 0.5 * k2_r, (satellite.velocity) + 0.5 * k2_v)
    k3_r = dt * (satellite.velocity + 0.5 * k2_v) / 1000
    
    k4_v = dt * get_acceleration(satellite.position + k3_r, (satellite.velocity) + k3_v)
    k4_r = dt * (satellite.velocity + k3_v) / 1000
    
    satellite.position +=(k1_r + 2*k2_r + 2*k3_r + k4_r) / 6 
    satellite.velocity +=(k1_v + 2*k2_v + 2*k3_v + k4_v) / 6
        

#actual drawing of objects on screen
circle_radius = 10
c1_pos = (planet.position[0],planet.position[1])
c2_pos = (int(satellite.position[0]), int(satellite.position[1]))
   
clock = pygame.time.Clock()
fps = 50
positions = []


# Main loop
running = True
while running:
    # Event handling
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False
            
    dt = clock.tick(fps) / 1000.0
    
    rk4_integration(satellite, planet, dt)
    positions.append((int(satellite.position[0]), satellite.position[1]))

    #update the position of the drawn object on screen
    c1_pos = (planet.position[0],planet.position[1])
    c2_pos = (int(satellite.position[0]), int(satellite.position[1]))
   
    #Draw
    window.fill((0,0,0))
    pygame.draw.circle(window, planet.colour, c1_pos, 25)
    pygame.draw.circle(window, satellite.colour, c2_pos, 5)

    if len(positions) > 1:
        for i in range(len(positions) - 1):
            pygame.draw.line(window,(205,205,205),positions[i], positions[i + 1], 1)

    draw_grid()
    
    #Refresh display
    pygame.display.flip()
    clock.tick(fps)

pygame.quit()
sys.exit()

Две орбиты совпадают: