用MATLAB写八大行星的运动程序

以下是一个简单的MATLAB程序，用于模拟八大行星的运动。这个程序假设行星之间没有相互作用，只受到太阳的引力影响。其中，轨道数据和质量数据来自NASA的公开数据（https://nssdc.gsfc.nasa.gov/planetary/factsheet/）

% 八大行星的质量（单位：kg）
mass = [3.30e23, 4.87e24, 5.97e24, 6.42e23, 1.90e27, 5.68e26, 8.68e25, 1.02e26];

% 八大行星的轨道半径（单位：m）
semi_major_axis = [57910000, 108200000, 149600000, 227940000, 778330000, 1429400000, 2870990000, 4495060000];

% 八大行星的轨道离心率
eccentricity = [0.2056, 0.0068, 0.0167, 0.0934, 0.0484, 0.0557, 0.0444, 0.2488];

% 八大行星的轨道倾角（单位：弧度）
inclination = [0.00005, 0.00005, 0.00005, 0.00005, 0.00005, 0.00005, 0.00005, 0.00005];

% 八大行星的轨道周期（单位：秒）
period = [87.969, 224.701, 365.256, 686.980, 4332.589, 10759.22, 30685.4, 60190];

% 八大行星的初始位置和速度（单位：m，m/s）
initial_pos = [semi_major_axis(1), 0, 0];
initial_vel = [0, sqrt(6.6743e-11 * mass(1) / semi_major_axis(1)), 0];

% 模拟时间
t_start = 0;
t_end = 5 * period(8);
dt = period(8) / 1000;
tspan = t_start:dt:t_end;

% 定义ODE函数，计算加速度
function dYdt = ode(t, Y, mu)
    r = norm(Y(1:3));
    dYdt(1:3, 1) = Y(4:6);
    dYdt(4:6, 1) = -mu / r^3 * Y(1:3);
end

% 使用ODE函数计算轨道
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-8);
[t, Y] = ode45(@(t, Y) ode(t, Y, 6.6743e-11 * mass(1)), tspan, [initial_pos, initial_vel], options);

% 绘制轨道图
figure;
hold on;
for i = 1:8
    plot3(Y(:, 3 * i - 2), Y(:, 3 * i - 1), Y(:, 3 * i), 'LineWidth', 1.5);
end
plot3(0, 0, 0, 'r.', 'MarkerSize', 20);
axis equal;
xlabel('x (m)');
ylabel('y (m)');
zlabel('z (m)');
title('Orbits of the eight planets (Sun at origin)');

这个程序使用ODE函数来计算每个行星的轨道，然后将结果绘制在三维坐标系中。在绘制过程中，我们还将太阳的位置用红点表示出来，以便更好地理解行星的运动轨迹。