我正在尝试模拟一个有点现实的程序,地球和月球可以在引力作用下相互作用。现在的问题是月球一直在朝着地球旋转,我不明白为什么。
这是我的代码:
from math import sin,cos,sqrt,atan2,pi
import pygame
pygame.init()
class Planet:
dt = 1/100
G = 6.67428e-11 #G constant
scale = 1/(1409466.667) #1 m = 1/1409466.667 pixels
def __init__(self,x=0,y=0,radius=0,color=(0,0,0),mass=0,vx=0,vy=0):
self.x = x #x-coordinate pygame-window
self.y = y #y-coordinate pygame-window
self.radius = radius
self.color = color
self.mass = mass
self.vx = vx #velocity in the x axis
self.vy = vy #velocity in the y axis
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 #Calculates acceleration in x- and y-axis for body 1.
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)) #changes position considering each body's velocity.
self.y += ((self.vy * Planet.dt))
def move(self,body):
dx = (self.x - body.x) #Calculates difference in x- and y-axis between the bodies
dy = (self.y - body.y)
r = (sqrt((dy**2)+(dx**2))) #Calculates the distance between the bodies
angle = atan2(dy, dx) #Calculates the angle between the bodies with atan2!
if r < self.radius: #Checks if the distance between the bodies is less than the radius of the bodies. Uses then Gauss gravitational law to calculate force.
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) #Newtons gravitational formula.
Fx = cos(angle) * F
Fy = sin(angle) * F
return Fx,Fy
def motion():
for i in range(0,len(bodies)):
Fnx = 0 #net force
Fny = 0
for j in range(0,len(bodies)):
if bodies[i] != bodies[j]:
Fnx += (bodies[i].move(bodies[j]))[0]
Fny += (bodies[i].move(bodies[j]))[1]
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]) #width - height
trace = pygame.Surface((900, 650))
pygame.display.set_caption("Moon simulation")
FPS = 150 #how quickly/frames per second our game should update.
earth = Planet(450,325,30,(0,0,255),5.97219*10**(24),-24.947719394204714/2) #450= xpos,325=ypos,30=radius
luna = Planet(450,(575/11),10,(128,128,128),7.349*10**(22),1023)
bodies = [earth,luna]
running = True
clock = pygame.time.Clock()
while running: #if user clicks close window
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()
一旦我让地月系统开始工作,我想扩展它并尝试拥有三个物体(这就是为什么有这么多否则会是“不必要的”代码的原因)
我愿意接受建议或/和建议!谢谢
您应该一次性更新系统状态;这意味着,您还需要一次计算力,而不是逐个移动身体。
而不是计算力during更新循环:
def motion():
for i in range(0,len(bodies)):
Fnx = 0 #net force
Fny = 0
for j in range(0,len(bodies)):
if bodies[i] != bodies[j]:
Fnx += (bodies[i].move(bodies[j]))[0]
Fny += (bodies[i].move(bodies[j]))[1]
elif bodies[i] == bodies[j]:
continue
bodies[i].update_vel(Fnx,Fny)
bodies[i].draw(screen)
bodies[i].orbit(trace)
Fnx,Fny=0,0
安顿部队预先:
def motion():
force = [ (
sum([(bodies[i].move(bodies[j]))[0] for j in range(0, len(bodies)) if i != j ]),
sum([(bodies[i].move(bodies[j]))[1] for j in range(0, len(bodies)) if i != j ])
) for i in range(0,len(bodies)) ]
for i in range(0,len(bodies)):
Fnx = force[i][0]
Fny = force[i][1]
bodies[i].update_state(Fnx,Fny)
bodies[i].draw(screen)
bodies[i].orbit(trace)
Fnx,Fny=0,0
(我平时不用Python写,所以风格不是很完美。)
以下文字来自之前的回答。它可能有帮助,但不是解决问题所必需的;你可以在这里停止阅读。
您可以使用更精细的方法(如 Runge-Kutta)进一步减少数字截断错误。为此:
update_vel
和update_pos
,相反,尝试写一个update_state
同时结合两者的方法;重要的是,方程的左边要么是 delta 要么是新状态,方程的右边是旧状态(高阶 Runge-Kutta 将有一些中间状态,Planet.dt
的分数)如果 Runge-Kutta 的起点太重,请考虑 MacCormack 或 Lax-Wendroff。
对于 Lax-Wendroff'ish 方式,而不是:
def update_vel(self,Fnx,Fny):
ax = Fnx/self.mass
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 update_state(self,Fnx,Fny):
ax = Fnx/self.mass
ay = Fny/self.mass
self.x += (self.vx * Planet.dt) - (ax/Planet.scale) * Planet.dt**2
self.y += (self.vy * Planet.dt) - (ay/Planet.scale) * Planet.dt**2
self.vx -= (ax/Planet.scale) * Planet.dt
self.vy -= (ay/Planet.scale) * Planet.dt