猫史档案馆

import turtle
import random
import time


turtle.bgcolor('#555555') #(55)₁₆ = 85 = 255/3

t = turtle.Pen()
t.speed(0)
t.pencolor('#FF0000')
t.hideturtle()

u = turtle.Pen()
u.speed(0)
u.pencolor('#00FF00')
u.hideturtle()

v = turtle.Pen()
v.speed(0)
v.pencolor('#0000FF')
v.hideturtle()

Time = 0 #时间

power = 2 #万有引力与距离的power次方成反比

G = 1000 #万有引力常数

if input('{0>>default, 1>>random} = ') == '0':
    m = [1,1,1] #质量
    ts = [100,-50] #坐标
    us = [0,150]
    tv = [0,1.5] #速度
    uv = [-1,-0.5]
else:
    m = []
    ts = []
    us = []
    tv = []
    uv = []
    
    for i in range(2):
        m.append(random.uniform(0,2))
        ts.append(random.uniform(-200,200))
        us.append(random.uniform(-200,200))
        tv.append(random.uniform(-2,2))
        uv.append(random.uniform(-2,2))
    m.append(random.uniform(max(m[0], m[1]), 2))

print()
print('m = '+str(m))
print('ts = '+str(ts))
print('us = '+str(us))
print('tv = '+str(tv))
print('uv = '+str(uv))
print()
    
time.sleep(1)


if m[2] == 0:
    vs = [0,0]
else:
    vs = [-m[0]/m[2]*ts[0]-m[1]/m[2]*us[0],-m[0]/m[2]*ts[1]-m[1]/m[2]*us[1]]

pt = 30000
p = pt/50
n = pt//30
if n == 0:
    n = 1

t.penup()
u.penup()
v.penup()

t.goto(ts)
u.goto(us)
v.goto(vs)

t.pendown()
u.pendown()
v.pendown()

def atu(xy):
    l = (ts[0]-us[0])**2+(ts[1]-us[1])**2
    return m[1]*(us[xy]-ts[xy])/l/l**((power-1)/2)

def atv(xy):
    l = (ts[0]-vs[0])**2+(ts[1]-vs[1])**2
    return m[2]*(vs[xy]-ts[xy])/l/l**((power-1)/2)

def aut(xy):
    l = (us[0]-ts[0])**2+(us[1]-ts[1])**2
    return m[0]*(ts[xy]-us[xy])/l/l**((power-1)/2)

def auv(xy):
    l = (us[0]-vs[0])**2+(us[1]-vs[1])**2
    return m[2]*(vs[xy]-us[xy])/l/l**((power-1)/2)

while True:
    for i in range(abs(n)):
        p = pt*(((ts[0]-us[0])**2+(ts[1]-us[1])**2)**(-1/2)+((us[0]-vs[0])**2+(us[1]-vs[1])**2)**(-1/2)+((vs[0]-ts[0])**2+(vs[1]-ts[1])**2)**(-1/2)) #实时刷新精度。能量不守恒是误差的体现，若时时计算总能量，会发现误差主要是在恒星的距离小、速度大时产生的，所以在此时增大精度是最优的选择，既节省时间，又减小误差；同时还避免了“飞星”现象的发生
        ta = [G*(atu(0)+atv(0)),G*(atu(1)+atv(1))]
        ua = [G*(aut(0)+auv(0)),G*(aut(1)+auv(1))]
        
        tv = [tv[0]+ta[0]/p,tv[1]+ta[1]/p]
        uv = [uv[0]+ua[0]/p,uv[1]+ua[1]/p]
        
        ts = [ts[0]+tv[0]/p,ts[1]+tv[1]/p]
        us = [us[0]+uv[0]/p,us[1]+uv[1]/p]
        if m[2] != 0:
            vs = [-m[0]/m[2]*ts[0]-m[1]/m[2]*us[0],-m[0]/m[2]*ts[1]-m[1]/m[2]*us[1]]

        Time += 1/p

    t.goto(ts)
    u.goto(us)
    v.goto(vs)

    print(Time)

#近似模拟三颗恒星（近似为质点）在互相的万有引力作用下的运动轨迹。
#初始参数（万有引力常数G，质量m，坐标s，速度v，模拟精度pt）可在代码中自行调整，终端区显示时间。