请帮我写一个3维模拟三体问题的代码

当然，我可以帮你写一个简单的3维模拟三体问题的代码。以下是一个使用Python编写的示例代码：

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

# 定义常量
G = 6.67430e-11  # 万有引力常量

# 定义初始条件
m1 = 1.989e30  # 太阳质量
m2 = 5.972e24  # 地球质量
m3 = 7.348e22  # 月球质量

r1 = np.array([0, 0, 0])  # 太阳位置
r2 = np.array([147e9, 0, 0])  # 地球位置（距离太阳平均距离）
r3 = np.array([147e9 + 3.844e8, 0, 0])  # 月球位置（距离地球平均距离）

v1 = np.array([0, 0, 0])  # 太阳速度
v2 = np.array([0, 29.29e3, 0])  # 地球速度（绕太阳的速度）
v3 = np.array([0, 29.29e3 + 1.022e3, 0])  # 月球速度（绕地球的速度）

# 定义运动方程
def acceleration(r1, r2, r3):
    r12 = np.linalg.norm(r2 - r1)
    r13 = np.linalg.norm(r3 - r1)
    r23 = np.linalg.norm(r3 - r2)

    a1 = G * m2 * (r2 - r1) / r12**3 + G * m3 * (r3 - r1) / r13**3
    a2 = G * m1 * (r1 - r2) / r12**3 + G * m3 * (r3 - r2) / r23**3
    a3 = G * m1 * (r1 - r3) / r13**3 + G * m2 * (r2 - r3) / r23**3

    return a1, a2, a3

# 定义时间步长和模拟时间
dt = 60  # 时间步长（秒）
t_sim = 60 * 60 * 24 * 30  # 模拟时间（秒）

# 初始化数组
num_steps = int(t_sim / dt)
r1_arr = np.zeros((num_steps, 3))
r2_arr = np.zeros((num_steps, 3))
r3_arr = np.zeros((num_steps, 3))

# 进行模拟
for i in range(num_steps):
    a1, a2, a3 = acceleration(r1, r2, r3)

    # 更新位置和速度
    r1 += v1 * dt + 0.5 * a1 * dt**2
    r2 += v2 * dt + 0.5 * a2 * dt**2
    r3 += v3 * dt + 0.5 * a3 * dt**2

    v1 += a1 * dt
    v2 += a2 * dt
    v3 += a3 * dt

    r1_arr[i] = r1
    r2_arr[i] = r2
    r3_arr[i] = r3

# 绘制结果
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(r1_arr[:, 0], r1_arr[:, 1], r1_arr[:, 2], label='Sun')
ax.plot(r2_arr[:, 0], r2_arr[:, 1], r2_arr[:, 2], label='Earth')
ax.plot(r3_arr[:, 0], r3_arr[:, 1], r3_arr[:, 2], label='Moon')
ax.set_xlabel('X [m]')
ax.set_ylabel('Y [m]')
ax.set_zlabel('Z [m]')
ax.legend()
plt.show()

这段代码使用欧拉方法对三个天体（太阳、地球、月球）的运动进行了模拟，并绘制了它们在三维空间中的轨迹。你可以根据需要调整初始条件、时间步长和模拟时间来进行更详细的模拟。希望对你有所帮助！