# -*- coding:utf-8 -*- # title :二维天体运行模拟 # description :二维天体运行模拟，适合初学者学习。 # author :Python超人 # date :2025-03-23 # link :https://gitcode.net/pythoncr/ # python_version :3.8 # ============================================================================== import math import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation # 设置 matplotlib 支持中文显示 plt.rcParams['font.sans-serif'] = ['SimHei'] # 使用黑体字体 plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题 plt.rcParams['font.family'] = 'Times New Roman' # https://miktex.org/download # 安装：https://mirrors.ustc.edu.cn/CTAN/systems/win32/miktex/setup/windows-x64/basic-miktex-24.1-x64.exe # plt.rcParams['text.usetex'] = 'True' # 定义物理常量 # 万有引力常量，单位：m^3 kg^-1 s^-2 万有引力常量 = 6.67430e-11 # 定义星球类 class 星球: def __init__(self, 质量, 位置, 速度, 颜色, 大小=None, 名称=None, 名称偏移=[2, 0], 轨迹最大长度=2200): """ 初始化星球对象。 :param 质量: 星球的质量，单位：kg :param 位置: 星球的初始位置，是一个二维列表或数组，单位：m :param 速度: 星球的初始速度，是一个二维列表或数组，单位：m/s :param 颜色: 星球的颜色，用于绘图显示 :param 大小: 星球在图中的显示大小，若无指定则根据质量计算 :param 名称: 星球的名称，用于在图中显示 :param 名称偏移: 星球名称显示位置相对于星球位置的偏移量，是一个二维列表或数组，单位：m :param 轨迹最大长度: 星球轨迹的最大记录长度，无单位 """ if 大小 is None: # 大小 = np.sqrt(质量) / 1000 大小 = np.power(质量, 0.16) / 10 self.质量 = 质量 self.位置 = np.array(位置, dtype=float) self.速度 = np.array(速度, dtype=float) self.颜色 = 颜色 self.大小 = 大小 self.名称 = 名称 self.名称偏移 = 名称偏移 self.轨迹 = [self.位置.copy()] self.轨迹最大长度 = 轨迹最大长度 self.加速度 = np.zeros(2) def 更新(self, 星球列表, 时间步长): """ 更新星球的位置和速度。 :param 星球列表: 包含所有星球对象的列表 :param 时间步长: 每次更新的时间间隔，单位：s """ 总力 = np.zeros(2) for 星球 in 星球列表: if 星球 is not self: 距离向量 = 星球.位置 - self.位置 # 距离指的是两个物体之间的长度，是一个标量，仅具有大小 距离 = np.linalg.norm(距离向量) # 距离向量:[30,40] 距离:50 # 计算万有引力，力的单位：N（kg m s^-2） # 万有引力是一个矢量，它的方向是沿着两个物体质心的连线方向。 # 所以，要得到万有引力的矢量形式，就需要在大小的基础上乘以一个单位向量，该单位向量的方向与距离向量的方向相同。 力 = 万有引力常量 * self.质量 * 星球.质量 * 距离向量 / (距离 ** 3) 总力 += 力 # 计算加速度，单位：m s^-2 加速度 = 总力 / self.质量 self.加速度 = 加速度 # 更新速度，单位：m/s self.速度 += 加速度 * 时间步长 # 更新位置，单位：m self.位置 += self.速度 * 时间步长 self.轨迹.append(self.位置.copy()) if len(self.轨迹) > self.轨迹最大长度: self.轨迹.pop(0) # 初始化图形 图形, 坐标轴 = plt.subplots() 星球列表 = [] # 时间步长，单位：s 时间步长 = 0.01 def 场景1(): global 星球列表, 时间步长 时间步长 = 0.03 # 创建星球列表 星球列表 = [ 星球(质量=5e12, 位置=[0, 0], 速度=[0, 0.15], 颜色='red', 名称偏移=[1, 0]), 星球(质量=2e11, 位置=[25, 0], 速度=[0, 2.1], 颜色='blue'), 星球(质量=2e11, 位置=[-25, 0], 速度=[0, -1.8], 颜色='green') ] 坐标轴.set_xlim(-50, 50) 坐标轴.set_ylim(-20, 80) def 场景2(): global 星球列表, 时间步长 时间步长 = 0.01 r1, r2 = 20, 10 # A、B 质量 mA = 1e13 mB = 2 * mA pxA, pyA = 0, r1 pxB, pyB = r1 + r2, r1 # A、B 坐标 pA, pB = [pxA, pyA], [pxB, pyB] vxA, vxB = 0, 0 # https://latex.codecogs.com/svg.image?\frac{G.mA.mB}{(r1+r2)^{2}}%20=%20mA.\frac{vA^{2}}{r1}=%20mB.\frac{vB^{2}}{r1} # https://latex.codecogs.com/png.image?\dpi{300}\frac{G.m_{A}.m_{B}}{(r_{1}+r_{2})^{2}}%20=%20m_{A}.\frac{v{_{A}}^{2}}{r_{1}}=%20m_{B}.\frac{v{_{B}}^{2}}{r_{1}} # 万有引力常量 *mA * mB/ math.pow(r1+r2,2) = mA*math.pow(vyA,2)/r1 # https://latex.codecogs.com/svg.image?vyA%20=%20\sqrt{\frac{G.mA.mB}{(r1+r2)^{2}}*\frac{r1}{mA}} # https://latex.codecogs.com/svg.image?v_{Ay}%20=%20\sqrt{\frac{G.m_{A}.m_{B}}{(r_{1}+r_{2})^{2}}*\frac{r_{1}}{m_{A}}} vyA = -math.sqrt((万有引力常量 * mA * mB) / math.pow(r1 + r2, 2) * r1 / mA) # 万有引力常量 *mA * mB/ math.pow(r1+r2,2) = mB*math.pow(vyB,2)/r2 # https://latex.codecogs.com/svg.image?vyB%20=%20\sqrt{\frac{G.mA.mB}{(r1+r2)^{2}}*\frac{r2}{mB}} # https://latex.codecogs.com/svg.image?v_{By}%20=%20\sqrt{\frac{G.m_{A}.m_{B}}{(r_{1}+r_{2})^{2}}*\frac{r_{2}}{m_{B}}} vyB = math.sqrt((万有引力常量 * mA * mB) / math.pow(r1 + r2, 2) * r2 / mB) vA, vB = [vxA, vyA], [vxB, vyB] 星球列表 = [ 星球(质量=mA, 位置=pA, 速度=vA, 颜色='blue', 名称="A"), 星球(质量=mB, 位置=pB, 速度=vB, 颜色='green', 名称="B") ] 坐标轴.set_xlim(-5, 60) 坐标轴.set_ylim(-10, 50) 场景1() # 场景2() # 设置纵横比为固定值，使得 X、Y 刻度比例不变 坐标轴.set_aspect('equal') # 调整边距 plt.subplots_adjust(left=0.05, bottom=0.04, right=0.97, top=0.97, wspace=0.4, hspace=0.4) 点列表 = [坐标轴.plot([], [], 'o', color=星球.颜色, markersize=星球.大小)[0] for 星球 in 星球列表] 轨迹列表 = [坐标轴.plot([], [], '--', color=星球.颜色)[0] for 星球 in 星球列表] 名称列表 = [坐标轴.text(0, 0, 星球.名称, color="white", size=30, fontname='SimHei') for 星球 in 星球列表] 速度线列表 = [plt.Arrow(0, 0, 0, 0, color=星球.颜色) for 星球 in 星球列表] 加速度线列表 = [plt.Arrow(0, 0, 0, 0, color=星球.颜色) for 星球 in 星球列表] 速度文本列表 = [坐标轴.text(0, 0, "", size=15, fontname='Consolas', color=星球.颜色) for 星球 in 星球列表] 加速度文本列表 = [坐标轴.text(0, 0, "", size=15, fontname='Consolas', color="black") for 星球 in 星球列表] for 速度线 in 速度线列表: 坐标轴.add_patch(速度线) for 加速度线 in 加速度线列表: 坐标轴.add_patch(加速度线) # 初始化函数 def 初始化(): """ 初始化动画的状态。 """ for 点 in 点列表: 点.set_data([], []) for 轨迹 in 轨迹列表: 轨迹.set_data([], []) for 名称 in 名称列表: 名称.set_position((0, 0)) # for 速度线 in 速度线列表: # 速度线.set_data([0, 0], [0, 0]) # for 加速度线 in 加速度线列表: # 加速度线.set_data([0, 0], [0, 0]) # for 速度文本 in 速度文本列表: # 速度文本.set_text("") # for 加速度文本 in 加速度文本列表: # 加速度文本.set_text("") return 点列表 + 轨迹列表 + 名称列表 + 速度线列表 + 加速度线列表 + 速度文本列表 + 加速度文本列表 # 更新函数 def 更新(帧): """ 更新每一帧动画中星球的位置和轨迹。 :param 帧: 当前动画的帧数，无单位 """ global 时间步长 for 星球 in 星球列表: 星球.更新(星球列表, 时间步长) 速度线长率 = 1.5 加速度线长率 = 4 for i, 星球 in enumerate(星球列表): x, y = 星球.位置[0], 星球.位置[1] 点列表[i].set_data(x, y) 轨迹 = np.array(星球.轨迹) 轨迹列表[i].set_data(轨迹[:, 0], 轨迹[:, 1]) 偏移量_x, 偏移量_y = 0, 0 # 星球.名称偏移 星球.位置[0] + 偏移量_x + 0.5, 星球.位置[1] + 偏移量_y 名称列表[i].set_position((x - 0.8, y - 1)) min_length, max_length = 0, 100 # 更新速度箭头 速度线列表[i].remove() dx, dy = 星球.速度[0] * 速度线长率, 星球.速度[1] * 速度线长率 distance = (dx ** 2 + dy ** 2) ** 0.5 if distance > max_length: # 重新设置 dx,dy ，使得 (dx ** 2 + dy ** 2) ** 0.5 = 50 dx, dy = dx / distance * max_length, dy / distance * max_length elif distance < min_length: # 重新设置 dx,dy ，使得 (dx ** 2 + dy ** 2) ** 0.5 = 5 dx, dy = dx / distance * min_length, dy / distance * min_length 速度线 = plt.arrow(x, y, dx, dy, width=0.02, head_width=0.8, color=星球.颜色) 速度线列表[i] = 速度线 坐标轴.add_patch(速度线) # 更新加速度箭头 加速度线列表[i].remove() dx, dy = 星球.加速度[0] * 加速度线长率, 星球.加速度[1] * 加速度线长率 distance = (dx ** 2 + dy ** 2) ** 0.5 if distance > max_length: # 重新设置 dx,dy ，使得 (dx ** 2 + dy ** 2) ** 0.5 = 50 dx, dy = dx / distance * max_length, dy / distance * max_length elif distance < min_length: # 重新设置 dx,dy ，使得 (dx ** 2 + dy ** 2) ** 0.5 = 5 dx, dy = dx / distance * min_length, dy / distance * min_length 加速度线 = plt.arrow(x, y, dx, dy, width=0.02, head_width=0.8, color="black") 加速度线列表[i] = 加速度线 坐标轴.add_patch(加速度线) # 更新速度文本 速度文本列表[i].set_position((x + 偏移量_x - 2, y + 偏移量_y + 3)) # 速度文本列表[i].set_text(f"v:{星球.速度[0]:+.2f},{星球.速度[1]:+.2f}(m/s)") 速度文本列表[i].set_text(f"v:{np.linalg.norm(星球.速度):.2f}(m/s)") # 更新加速度文本 加速度文本列表[i].set_position((x + 偏移量_x - 2, y + 偏移量_y + 5)) # 加速度文本列表[i].set_text(f"a:{星球.加速度[0]:+.2f},{星球.加速度[1]:+.2f}(m/s²)") 加速度文本列表[i].set_text(f"a:{np.linalg.norm(星球.加速度):.2f}(m/s²)") return 点列表 + 轨迹列表 + 名称列表 + 速度线列表 + 加速度线列表 + 速度文本列表 + 加速度文本列表 def 设置窗口(): # 设置窗口大小（单位：英寸） width = 9 height = 11 图形.set_size_inches(width, height) # 设置窗口位置（x, y, width, height） x_pos = 1020 y_pos = 0 try: manager = plt.get_current_fig_manager() if hasattr(manager, 'window'): manager.window.setGeometry(x_pos, y_pos, int(width * 100), int(height * 100)) except AttributeError: print("无法设置窗口位置，可能是因为使用的后端不支持此操作。") # 创建动画 动画 = FuncAnimation(图形, 更新, frames=range(1000), init_func=初始化, blit=True, interval=2) # 显示动画 plt.show()