天体引力弹弓效应

天体引力弹弓效应
模拟航天器利用行星引力改变轨道的过程。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib.patches as patches

# 设置中文字体 - 改为系统常用中文字体
plt.rcParams["font.family"] = ["SimHei", "Microsoft YaHei", "Arial Unicode MS", "sans-serif"]
plt.rcParams["axes.unicode_minus"] = False  # 解决负号显示问题

# 物理常量
G = 6.67430e-11  # 引力常量 (m^3 kg^-1 s^-2)
AU = 1.496e11  # 天文单位 (米)
DAY = 86400  # 一天的秒数
YEAR = 365 * DAY  # 一年的秒数


# 天体参数
class CelestialBody:
    def __init__(self, name, mass, radius, color, pos, vel):
        self.name = name
        self.mass = mass  # kg
        self.radius = radius  # 可视化半径 (AU)
        self.color = color
        self.pos = np.array(pos, dtype=float)  # 位置 (m)
        self.vel = np.array(vel, dtype=float)  # 速度 (m/s)
        self.trail = []  # 轨迹点列表
        self.patch = None  # 图形补丁
        self.trail_line = None  # 轨迹线

    def update_position(self, dt):
        self.pos += self.vel * dt
        self.trail.append(self.pos.copy())

    def draw(self, ax):
        # 绘制天体
        if self.patch:
            self.patch.remove()

        self.patch = patches.Circle(
            (self.pos[0] / AU, self.pos[1] / AU),
            radius=self.radius,
            color=self.color,
            alpha=0.8
        )
        ax.add_patch(self.patch)

        # 绘制轨迹
        if len(self.trail) > 1:
            if self.trail_line:
                self.trail_line.remove()

            trail_x = [pos[0] / AU for pos in self.trail]
            trail_y = [pos[1] / AU for pos in self.trail]
            self.trail_line, = ax.plot(trail_x, trail_y, color=self.color, alpha=0.3, linewidth=1)

        return self.patch


# 创建天体
sun = CelestialBody(
    "太阳",
    mass=1.989e30,
    radius=0.2,
    color='yellow',
    pos=[0, 0],
    vel=[0, 0]
)

earth = CelestialBody(
    "地球",
    mass=5.972e24,
    radius=0.05,
    color='blue',
    pos=[1 * AU, 0],
    vel=[0, 29780]  # 地球轨道速度 (m/s)
)

jupiter = CelestialBody(
    "木星",
    mass=1.898e27,
    radius=0.1,
    color='orange',
    pos=[5.2 * AU, 0],
    vel=[0, 13070]  # 木星轨道速度 (m/s)
)

spacecraft = CelestialBody(
    "航天器",
    mass=1000,
    radius=0.02,
    color='red',
    pos=[1.1 * AU, 0],
    vel=[0, 35000]  # 初始速度 (m/s)
)


# 计算引力
def calculate_gravitational_force(body1, body2):
    """计算两个天体之间的引力"""
    r = body2.pos - body1.pos
    distance = np.linalg.norm(r)
    if distance == 0:  # 避免除以零
        return np.zeros(2)
    force_magnitude = G * body1.mass * body2.mass / (distance ** 2)
    force_direction = r / distance
    return force_direction * force_magnitude


# 更新所有天体的位置和速度
def update_all_bodies(dt, bodies):
    # 计算所有引力
    forces = {}
    for body in bodies:
        total_force = np.zeros(2)
        for other_body in bodies:
            if body is not other_body:
                total_force += calculate_gravitational_force(body, other_body)
        forces[body] = total_force

    # 更新速度和位置
    for body in bodies:
        acceleration = forces[body] / body.mass
        body.vel += acceleration * dt
        body.update_position(dt)


# 查找最近距离
def find_closest_approach(bodies, target_body, steps, dt):
    """计算航天器与目标天体的最近距离"""
    min_distance = float('inf')
    min_step = 0

    # 保存当前状态以便恢复
    initial_states = {body: (body.pos.copy(), body.vel.copy(), body.trail.copy()) for body in bodies}

    for step in range(steps):
        update_all_bodies(dt, bodies)

        # 计算与目标的距离
        distance = np.linalg.norm(spacecraft.pos - target_body.pos)
        if distance < min_distance:
            min_distance = distance
            min_step = step

    # 恢复初始状态
    for body, (pos, vel, trail) in initial_states.items():
        body.pos = pos
        body.vel = vel
        body.trail = trail

    return min_distance, min_step


# 创建图表
fig, ax = plt.subplots(figsize=(10, 8))
ax.set_xlim(-8, 8)  # AU
ax.set_ylim(-8, 8)  # AU
ax.set_xlabel('X (AU)')
ax.set_ylabel('Y (AU)')
ax.set_title('航天器引力弹弓效应模拟')
ax.grid(True, linestyle='--', alpha=0.7)

# 绘制图例
legend_elements = [
    patches.Patch(color=sun.color, label=sun.name),
    patches.Patch(color=earth.color, label=earth.name),
    patches.Patch(color=jupiter.color, label=jupiter.name),
    patches.Patch(color=spacecraft.color, label=spacecraft.name)
]
ax.legend(handles=legend_elements, loc='upper right')

# 初始化速度向量
velocity_vector = None

# 创建信息文本
info_text = ax.text(0.02, 0.95, '', transform=ax.transAxes,
                    bbox=dict(facecolor='white', alpha=0.8))

# 计算最近距离
simulation_steps = 2 * 365  # 模拟2年
time_step = 1 * DAY  # 每天一步
closest_distance, closest_step = find_closest_approach(
    [sun, earth, jupiter, spacecraft],
    jupiter,
    simulation_steps,
    time_step
)


# 动画更新函数
def update(frame):
    global velocity_vector

    # 更新所有天体
    update_all_bodies(time_step, [sun, earth, jupiter, spacecraft])

    # 绘制天体和轨迹
    sun.draw(ax)
    earth.draw(ax)
    jupiter.draw(ax)
    spacecraft.draw(ax)

    # 更新速度向量
    if velocity_vector:
        try:
            velocity_vector.remove()
        except ValueError:
            pass  # 如果向量已被移除，则忽略错误

    velocity_vector = ax.arrow(
        spacecraft.pos[0] / AU, spacecraft.pos[1] / AU,
        spacecraft.vel[0] / 500000, spacecraft.vel[1] / 500000,
        head_width=0.1, head_length=0.2, fc='green', ec='green', alpha=0.7,
        label='航天器速度'
    )

    # 更新信息文本
    days = frame * (time_step / DAY)
    speed = np.linalg.norm(spacecraft.vel)
    distance_to_jupiter = np.linalg.norm(spacecraft.pos - jupiter.pos) / AU

    info = (f"时间: {days:.1f} 天\n"
            f"航天器速度: {speed:.1f} m/s\n"
            f"与木星距离: {distance_to_jupiter:.2f} AU")

    if frame == closest_step:
        info += f"\n\n最近距离: {closest_distance / AU:.3f} AU"

    info_text.set_text(info)

    # 更新标题
    ax.set_title(f'航天器引力弹弓效应模拟 - 第 {days:.0f} 天')

    return (sun.patch, earth.patch, jupiter.patch, spacecraft.patch,
            velocity_vector, info_text)


# 创建动画
ani = FuncAnimation(
    fig,
    update,
    frames=simulation_steps,
    interval=50,
    blit=False,  # 设为False以避免更新问题
    repeat=False
)

# 显示动画
plt.tight_layout()
plt.show()