Python行星模拟故障求助：月球无法正常环绕地球

问题概述 #

我正在用Python编写天体受力模拟程序，计算地球、月球及假想月球间的引力，根据速度和受力变化实现符合现实的天体运动，通过缩放因子转换米与像素单位。初始代码无法实现月球环绕地球，修改update_vel和update_pos方法后，月球虽能环绕，但地球出现向左移动且呈抛物线轨迹的异常。调整时间步长dt和缩放因子位置均未解决问题，在update_vel中除以缩放因子虽得到非线性轨迹，但因单位混合导致轨迹非圆形。

原始代码 #

from math import sin,cos,sqrt,atan2,pi
import pygame
pygame.init()

class Planet:
    dt = 1/100
    G = 6.67428e-11 #引力常数G
    scale = 1/(1409466.667) #1米对应1/1409466.667像素
    def __init__(self,x=0,y=0,radius=0,color=(0,0,0),mass=0,vx=0,vy=0):
        self.x = x #Pygame窗口x坐标
        self.y = y #Pygame窗口y坐标
        self.radius = radius
        self.color = color
        self.mass = mass
        self.vx = vx #x轴速度
        self.vy = vy #y轴速度
        
    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 #计算天体x、y轴加速度
        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))
        
    def move(self,body):
        dx = (self.x - body.x) #计算两天体x、y轴差值
        dy = (self.y - body.y)
        r = (sqrt((dy**2)+(dx**2))) #计算两天体间距离
        angle = atan2(dy, dx) #用atan2计算两天体间夹角
        if r < self.radius: #距离小于天体半径时用高斯引力定律算力
            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) #牛顿万有引力公式
            Fx = cos(angle) * F
            Fy = sin(angle) * F
        return Fx,Fy

def motion():
    for i in range(0,len(bodies)):
        Fnx = 0 #合外力x分量
        Fny = 0 #合外力y分量
        for j in range(0,len(bodies)):
            if bodies[i] != bodies[j]:
                Fnx += (bodies[i].move(bodies[j]))[0] #天体j对i的x方向合力
                Fny += (bodies[i].move(bodies[j]))[1] #天体j对i的y方向合力
            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]) #宽-高
trace = pygame.Surface((900, 650))
pygame.display.set_caption("月球模拟")
FPS = 60 #每秒帧数

earth = Planet(450,325,30,(0,0,255),5.97219*10**(24)) #450=x坐标,325=y坐标,30=半径
luna = Planet(450,(575/11),10,(128,128,128),7.349*10**(22),1023)
moon = Planet() #第二个月球
bodies = [earth,luna]

running = True
clock = pygame.time.Clock()

while running: #用户点击关闭窗口时退出
    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()

修改后的方法 #

def update_vel(self,Fnx,Fny):
    ax = Fnx/self.mass #计算天体x、y轴加速度
    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))

解决方案 #

核心问题分析 #

1. 单位混合混乱：代码中同时用像素存储位置、物理单位存储速度，导致加速度与速度的单位不匹配，计算逻辑错误。
2. 质心未静止：初始状态下地球速度为0，月球有速度，系统质心会移动，导致地球出现不合理的抛物线运动。
3. 内部引力公式错误：距离小于天体半径时的引力计算不符合高斯引力定律，公式完全错误。
4. 力的方向错误：move方法返回的力方向与实际引力方向相反，导致加速度方向反转。

修复步骤与完整代码 #

1. 统一物理单位与绘制单位

内部计算全部使用物理单位（米、米/秒），仅在绘制时转换为像素，彻底避免单位混合。

2. 修正内部引力公式

均匀球体内部的引力与距离成正比，公式应为 F = (G * self.mass * body.mass * r) / (self.radius_m ** 3)，其中radius_m是天体的物理半径（米）。

3. 设置质心静止的初始速度

根据动量守恒，地球的初始速度应为 vx = -(luna.mass * luna.vx) / earth.mass，抵消月球动量，保持系统质心静止。

4. 修正力的方向

确保返回的力指向对方天体，符合引力的吸引特性。

from math import sin, cos, sqrt, atan2, pi
import pygame
pygame.init()

class Planet:
    dt = 1/100
    G = 6.67428e-11  # 引力常数G
    scale = 1/(1409466.667)  # 1米对应1/1409466.667像素
    screen_center_x = 450
    screen_center_y = 325

    def __init__(self, x_m=0, y_m=0, radius_pix=0, color=(0,0,0), mass=0, vx=0, vy=0):
        # 物理单位（米）存储位置
        self.x_m = x_m
        self.y_m = y_m
        self.radius_pix = radius_pix
        self.color = color
        self.mass = mass
        self.vx = vx  # 物理单位：米/秒
        self.vy = vy  # 物理单位：米/秒
        # 物理半径（米）：根据像素半径和缩放因子转换
        self.radius_m = radius_pix / Planet.scale

    def get_screen_pos(self):
        # 转换物理位置到屏幕像素位置，以屏幕中心为原点
        screen_x = Planet.screen_center_x + self.x_m * Planet.scale
        screen_y = Planet.screen_center_y - self.y_m * Planet.scale
        return (int(screen_x), int(screen_y))

    def draw(self, screen):
        pygame.draw.circle(screen, self.color, self.get_screen_pos(), self.radius_pix)
    
    def orbit(self, trace):
        screen_x, screen_y = self.get_screen_pos()
        pygame.draw.rect(trace, self.color, (screen_x, screen_y, 2, 2))
        
    def update_vel(self, Fnx, Fny):
        # 加速度计算（物理单位：米/秒²）
        ax = Fnx / self.mass
        ay = Fny / self.mass
        # 更新速度（物理单位）
        self.vx += ax * Planet.dt
        self.vy += ay * Planet.dt
        self.update_pos()
     
    def update_pos(self):
        # 更新物理位置（米）
        self.x_m += self.vx * Planet.dt
        self.y_m += self.vy * Planet.dt
        
    def move(self, body):
        # 计算物理距离（米）
        dx_m = self.x_m - body.x_m
        dy_m = self.y_m - body.y_m
        r_m = sqrt(dy_m**2 + dx_m**2)
        angle = atan2(dy_m, dx_m)

        if r_m < self.radius_m:
            # 均匀球体内部引力：与距离成正比
            F = (Planet.G * self.mass * body.mass * r_m) / (self.radius_m ** 3)
        else:
            # 牛顿万有引力公式（物理单位）
            F = (Planet.G * self.mass * body.mass) / (r_m ** 2)
        
        Fx = cos(angle) * F
        Fy = sin(angle) * F
        # 返回的力是body对self的引力，方向指向body，所以取负
        return -Fx, -Fy

def motion():
    for i in range(len(bodies)):
        Fnx = 0
        Fny = 0
        for j in range(len(bodies)):
            if i != j:
                fx, fy = bodies[i].move(bodies[j])
                Fnx += fx
                Fny += fy
        bodies[i].update_vel(Fnx, Fny)
        bodies[i].draw(screen)
        bodies[i].orbit(trace)

# 初始化窗口
screen = pygame.display.set_mode([900, 650])
trace = pygame.Surface((900, 650))
pygame.display.set_caption("月球模拟")
FPS = 60

# 物理参数：地月距离约3.844e8米
earth_mass = 5.97219e24
luna_mass = 7.349e22
moon_distance = 3.844e8
# 月球环绕速度约1023米/秒
luna_vx = 0
luna_vy = 1023
# 地球初始速度：抵消月球动量，保持质心静止
earth_vx = -(luna_mass * luna_vx) / earth_mass
earth_vy = -(luna_mass * luna_vy) / earth_mass

# 创建天体：地球在物理原点，月球在y轴正方向
earth = Planet(x_m=0, y_m=0, radius_pix=30, color=(0,0,255), mass=earth_mass, vx=earth_vx, vy=earth_vy)
luna = Planet(x_m=0, y_m=moon_distance, radius_pix=10, color=(128,128,128), mass=luna_mass, vx=luna_vx, vy=luna_vy)
bodies = [earth, luna]

running = True
clock = pygame.time.Clock()

while running:
    clock.tick(FPS)    
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False
            
    screen.fill((0,0,0))
    screen.blit(trace, (0, 0))
    motion()

    pygame.display.flip()

pygame.quit()

修复说明 #

• 物理计算与绘制完全分离，内部所有位置、速度、力的计算都使用国际单位制，仅在绘制时转换为像素。
• 修正了内部引力的计算逻辑，符合高斯引力定律。
• 通过动量守恒设置地球初始速度，确保地月系统质心静止，地球只会做微小摆动而非抛物线运动。
• 调整了力的方向，保证引力指向对方天体，加速度方向符合物理规律。

内容的提问来源于stack exchange，提问作者Hale