带电粒子在行星磁场中运动的3D动画及绘图问题求助

我通过求解运动方程，尝试绘制并动画展示行星磁场中带电粒子的运动，已得到2D和3D绘图，但3D绘图仅呈现出类似单张2D图的效果。目前已获得随时间（1e5秒）变化的r（径向距离）、theta（纬度）、phi（经度）及其时间导数（速度）的解，希望能得到3D动画实现的帮助，以直观呈现粒子绕球体（行星）的运动。

现有代码 #

import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi
from scipy.integrate import odeint

scales = np.array([1e2, 0.1, 1, 1e-15, 10, 0.1])

GM = 379312077e8    # m^3/s^2       
β = 9.67e7  

def odes(p, t):
    r,x,θ,y,ϕ,z = p*scales   # assigning each ODE to a vector element
    
    # constants
    R = 60268e3
    g_10 = 21141e-9
    Ω = 9.74e-3
    B_θ = (R/r)**3*g_10*sin(θ)
    B_r = 2*(R/r)**3*g_10*cos(θ)
    β = 9.67e7

    # defining the ODEs only Lorentz Force
    drdt = x
    dxdt = r*(y**2 +(z+Ω)**2*sin(θ)**2-β*z*sin(θ)*B_θ)
    dθdt = y
    dydt = (-2*x*y +r*(z+Ω)**2*sin(θ)*cos(θ)+β*r*z*sin(θ)*B_r)/r 
    dϕdt = z
    dzdt = (-2*(z+Ω)*(x*sin(θ)+r*y*cos(θ))+β*(x*B_θ-r*y*B_r))/(r*sin(θ))

    return np.array([drdt,dxdt,dθdt,dydt,dϕdt,dzdt])/scales

# initial conditions
r0 = 6.7e+07
x0 = 0.0
θ0 = 88.0
y0 = 0.0
ϕ0 = 0.0
z0 = 0.022

# time window
t = np.arange(0,3600*24,360)
p0 = np.array([r0,x0,θ0,y0,ϕ0,z0])
p = odeint(odes,p0,t, atol=1e-8, rtol=1e-8)

r,x,θ,y,ϕ,z = p.T*scales[:,None]

# 2D plot the results
fig,ax=plt.subplots(2,3,figsize=(8,4))
plt.ylabel('parameters')
for a,u in zip(ax.flatten(),[r,x,θ,y,ϕ,z]):
    a.plot(t,u); a.grid()

plt.tight_layout(); plt.show()

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from mpl_toolkits import mplot3d

# 3D plots
fig = plt.figure()
omega=50
x_line= r
y_line =θ
z_line =ϕ
ax = plt.axes(projection="3d")
ax.plot3D(x_line,y_line,z_line, 'red')
ax.view_init(120,120)
plt.show()

2D绘图为r、theta、phi、dr/dt、dtheta/dt、dphi/dt随时间变化的曲线，3D绘图为r、theta、phi的三维曲线。

解决方案 #

1. 球坐标转笛卡尔坐标 #

当前3D图直接使用球坐标参数作为坐标轴，无法直观展示粒子在空间中的位置。需要将球坐标(r, θ, φ)转换为笛卡尔坐标(x, y, z)：

# 注意：原theta是纬度（0-90°），转换为极角（从z轴正方向算起）需用90°-theta，再转弧度
theta_rad = np.deg2rad(90 - θ)
phi_rad = np.deg2rad(ϕ)

x = r * np.sin(theta_rad) * np.cos(phi_rad)
y = r * np.sin(theta_rad) * np.sin(phi_rad)
z = r * np.cos(theta_rad)

2. 添加行星球体模型 #

为了对比粒子与行星的位置关系，在3D空间中添加代表行星的球体：

R = 60268e3  # 行星半径，与ODE中一致
# 生成球体网格
u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
planet_x = R * np.cos(u) * np.sin(v)
planet_y = R * np.sin(u) * np.sin(v)
planet_z = R * np.cos(v)

3. 实现3D运动动画 #

使用matplotlib.animation.FuncAnimation制作逐帧动画，展示粒子轨迹和实时位置：

import matplotlib.animation as animation

# 创建3D画布
fig = plt.figure(figsize=(10, 10))
ax = plt.axes(projection='3d')

# 绘制行星
ax.plot_surface(planet_x, planet_y, planet_z, color='blue', alpha=0.3)

# 初始化轨迹和粒子标记
trajectory, = ax.plot3D([], [], [], 'red', linewidth=1)
particle, = ax.plot3D([], [], [], 'o', color='orange', markersize=6)

# 设置坐标轴范围（根据粒子最大径向距离调整）
max_r = np.max(r)
ax.set_xlim(-max_r, max_r)
ax.set_ylim(-max_r, max_r)
ax.set_zlim(-max_r, max_r)
ax.set_xlabel('X (m)')
ax.set_ylabel('Y (m)')
ax.set_zlabel('Z (m)')
ax.set_title('Charged Particle Motion in Planetary Magnetic Field')

# 动画更新函数
def update(frame):
    # 更新已走过的轨迹
    trajectory.set_data_3d(x[:frame+1], y[:frame+1], z[:frame+1])
    # 更新当前粒子位置
    particle.set_data_3d(x[frame], y[frame], z[frame])
    return trajectory, particle

# 创建动画：frames为时间点数量，interval为帧间隔（毫秒）
ani = animation.FuncAnimation(fig, update, frames=len(t), interval=50, blit=True)
plt.show()

4. 优化建议 #

• 若动画卡顿，可对时间序列采样（如每隔2个时间点取一次数据）：t_sampled = t[::2]，对应的数据也同步采样
• 调整interval参数控制播放速度：值越小，动画越快
• 可添加坐标轴标签、标题增强可读性，或调整视角ax.view_init(elev=30, azim=45)

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