如何修改天体粒子模拟动画：仅逐帧显示粒子云并导出视频

需求说明 #

• 当前代码生成的粒子轨道为锯齿状线条，视觉效果差
• 显示100个粒子在每个时间点的动态位置，仅展示粒子云形态，不绘制轨迹
• 将动画保存为MP4格式

修改后的模拟代码 #

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation

# Constants
M_Sun = 1.989e30  # Solar Mass
M_bh = 4.3e6 * M_Sun  # Mass of Sgr A*
G = 6.67430e-11  # m^3 kg^(-1) s^(-2)
yr = 365 * 24 * 60 * 60  # 1 year in seconds
R = 6.371e3  # Radius of Sphere

# Number of particles
num_particles = 100

# Uniformly distributed particles in Sphere of radius R (注：此处代码未在初始条件中使用，若需替换初始位置可启用)
# phi = np.random.uniform(0, 2*np.pi, size=num_particles)
# costheta = np.random.uniform(-1, 1, size=num_particles)
# u = np.random.uniform(0, 1, size=num_particles)
# theta = np.arccos(costheta)
# r = R * (u**(1/3))
# x = r * np.sin(theta) * np.cos(phi)
# y = r * np.sin(theta) * np.sin(phi)
# z = r * np.cos(theta)

# Initial conditions for the particles (m and m/s)
initial_pos = np.random.uniform(0.8e13, 1.1e14, (num_particles, 3))
initial_vel = np.random.uniform(550e3, 730e3, (num_particles, 3))

# Steps
t_end = 1 * yr  # Total time of integration
dt_constant = 0.1
intervals = 1000000  # seconds after which to store pos & vel
next_interval_time = intervals

# Arrays to store pos and vel
pos = np.zeros((num_particles, 1, 3))  # shape: (#particles, timesteps taken, #coordinates)
vel = np.zeros((num_particles, 1, 3))  # shape: (#particles, timesteps taken, #coordinates)

# Initial conditions
time = np.zeros(1)  # time array
pos[:, 0] = initial_pos
vel[:, 0] = initial_vel
pos_output = []
vel_output = []
t_output = []

while time[-1] <= t_end:  # NOTE: We end up doing one more timestep after t_end
    r = np.linalg.norm(pos[:, -1], axis=1)
    acc = -G * M_bh / r[:, np.newaxis]**3 * pos[:, -1]  # np.newaxis for broadcasting with pos[:, i-1]

    # Calculate the time step for the current particle
    current_dt = dt_constant * np.sqrt(np.linalg.norm(pos[:, -1], axis=1)**3 / (G * M_bh))
    min_dt = np.min(current_dt)  # Use the minimum time step for all particles

    # leap frog integration
    # calculate velocity at half timestep
    half_vel = vel[:, -1] + 0.5 * acc * min_dt
    # current position only used in this timestep
    _pos_t = pos[:, -1] + half_vel * min_dt  # shape: (#particles, #coordinates)

    # Recalculate acceleration with the new position
    r = np.linalg.norm(_pos_t, axis=1)
    # Acceleration at timestep
    acc = -G * M_bh / r[:, np.newaxis]**3 * _pos_t  # np.newaxis for broadcasting with pos[:, i-1]
    # current velocity only used in this timestep
    _vel_t = half_vel + 0.5 * acc * min_dt  # shape: (#particles, #coordinates)
    
    # time at timestep t
    _time_t = time[-1] + min_dt

    # Check if the elapsed time has surpassed the next interval
    if _time_t >= next_interval_time:
        pos_output.append(_pos_t.copy())
        vel_output.append(_vel_t.copy())
        t_output.append(_time_t.copy())
        next_interval_time += intervals

    # add axis at position 1 to allow concatenation
    pos = np.concatenate((pos, _pos_t[:, np.newaxis]), axis=1)
    vel = np.concatenate((vel, _vel_t[:, np.newaxis]), axis=1)
    time = np.append(time, _time_t)

    # show current status by printing timestep number (-1 because initial conditions)
    print(f'timestep: {time.size -1} [progress: {_time_t/t_end*100:.3f}%]')

pos_output = np.array(pos_output).transpose((1, 0, 2))  # shape: (#objects, #stored timesteps, #coordinates)
np.save('pos_output.npy', pos_output)
vel_output = np.array(vel_output).transpose((1, 0, 2))  # shape: (#objects, #stored timesteps, #coordinates)
np.save('vel_output.npy', vel_output)
t_output = np.array(t_output)
np.save('t_output.npy', t_output)

print(pos_output.shape)
print(vel_output.shape)
print(t_output.shape)

# ANIMATION
from orbit_animation import animate_orbits
animate_orbits(pos_output)

修改后的动画代码（orbit_animation.py） #

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from matplotlib.animation import FFMpegWriter

def animate_orbits(pos, interval=50):
    fig = plt.figure(figsize=(8, 8))
    ax = fig.add_subplot(111, projection='3d')

    # Scatter plot for Sgr A*
    sgr_a_plot = ax.scatter([0], [0], [0], color='black', marker='o', s=100, label='Sgr A*')

    # Initialize scatter plot for particles (替换原线条为散点)
    cloud_plot = ax.scatter([], [], [], color='blue', marker='.', s=30, label='Cloud Particles')

    # Set plot labels and title
    ax.set_xlabel('X (m)')
    ax.set_ylabel('Y (m)')
    ax.set_zlabel('Z (m)')
    ax.legend(loc='upper right', bbox_to_anchor=(1.1, 1.1))
    ax.set_title('Cloud Particles orbiting Sgr A*')

    # 固定坐标轴范围（避免动画过程中抖动）
    x_min, x_max = np.min(pos[:, :, 0]) * 1.1, np.max(pos[:, :, 0]) * 1.1
    y_min, y_max = np.min(pos[:, :, 1]) * 1.1, np.max(pos[:, :, 1]) * 1.1
    z_min, z_max = np.min(pos[:, :, 2]) * 1.1, np.max(pos[:, :, 2]) * 1.1
    ax.set_xlim(x_min, x_max)
    ax.set_ylim(y_min, y_max)
    ax.set_zlim(z_min, z_max)

    # Animation update function
    def update(frame):
        # 更新粒子位置
        cloud_plot._offsets3d = (pos[:, frame, 0], pos[:, frame, 1], pos[:, frame, 2])
        return sgr_a_plot, cloud_plot

    # Create the animation
    animation = FuncAnimation(fig, update, frames=pos.shape[1], interval=interval, blit=True)

    # 保存为MP4格式（需提前安装ffmpeg：pip install ffmpeg-python）
    writer = FFMpegWriter(fps=1000//interval, metadata={'title': 'Particle Cloud Animation'}, bitrate=1800)
    animation.save('particle_cloud_animation.mp4', writer=writer)

    plt.show()

关键修改说明 #

1. 替换粒子展示方式：将原有的plot线条改为scatter散点，实现每个时间点仅显示当前粒子位置，不绘制轨迹
2. 固定坐标轴范围：提前计算所有时间点的粒子位置极值并扩大10%作为坐标轴范围，避免动画过程中坐标轴频繁缩放导致的抖动
3. 添加MP4保存功能：引入FFMpegWriter，将动画保存为MP4文件，需提前安装依赖pip install ffmpeg-python
4. 优化参数设置：调整粒子大小、黑洞标记大小，提升视觉效果；缩短动画间隔，让动画更流畅
5. 修正单位标注：原代码中坐标轴标注为km，但实际计算单位是m，已修正标注保持一致

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