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

以下是一个简单的MATLAB程序，用于模拟太阳系中八大行星的运动。该程序使用了牛顿万有引力定律和牛顿第二定律来计算行星在其轨道上的运动。

% 八大行星的基本信息
mass = [1.989e30, 3.3e23, 4.87e24, 5.97e24, 6.42e23, 1.9e27, 5.68e26, 8.68e25]; % 太阳、水星、金星、地球、火星、木星、土星、天王星、海王星的质量（kg）
radius = [0, 57.9e9, 108.2e9, 149.6e9, 227.9e9, 778.3e9, 1.43e12, 2.87e12, 4.5e12]; % 太阳系中八大行星的半径（m）
eccentricity = [0, 0.206, 0.007, 0.017, 0.094, 0.049, 0.057, 0.046, 0.011]; % 太阳系中八大行星的轨道离心率
inclination = [0, 7.0, 3.4, 0.0, 1.9, 1.3, 2.5, 0.8, 1.8]; % 太阳系中八大行星的轨道倾角
semimajor_axis = [0, 57.9e9, 108.2e9, 149.6e9, 227.9e9, 778.3e9, 1.43e12, 2.87e12, 4.5e12]; % 太阳系中八大行星的半长轴（m）

% 设置初始状态
pos = zeros(3, 9); % 八大行星的位置向量
vel = zeros(3, 9); % 八大行星的速度向量
pos(:, 1) = [0; 0; 0]; % 太阳的位置
vel(:, 1) = [0; 0; 0]; % 太阳的速度
for i = 2:9
    a = semimajor_axis(i); % 半长轴
    e = eccentricity(i); % 轨道离心率
    incl = inclination(i); % 轨道倾角
    r = a * (1 - e^2) / (1 + e * cosd(0)); % 初始距离
    pos(:, i) = [r * cosd(0); r * sind(0) * cosd(incl); r * sind(0) * sind(incl)]; % 初始位置
    v = sqrt(mass(1) / r) * sqrt((1 + e) / (1 - e)) * sqrt((1 + mass(i) / mass(1)) / (1 + mass(i) / mass(1)) * (1 - e^2)); % 初始速度
    vel(:, i) = [-v * sind(0); v * cosd(0) * cosd(incl); v * cosd(0) * sind(incl)]; % 初始速度
end

% 模拟行星运动
tspan = 0:3600*24:3600*24*365*10; % 模拟时间范围
y0 = [pos(:); vel(:)]; % 初始状态
[t, y] = ode45(@orbit, tspan, y0); % 数值模拟

% 绘制行星轨迹
figure;
hold on;
plot3(y(:, 1), y(:, 2), y(:, 3), 'r'); % 太阳的轨迹
plot3(y(:, 4), y(:, 5), y(:, 6), 'b'); % 水星的轨迹
plot3(y(:, 7), y(:, 8), y(:, 9), 'g'); % 金星的轨迹
plot3(y(:, 10), y(:, 11), y(:, 12), 'k'); % 地球的轨迹
plot3(y(:, 13), y(:, 14), y(:, 15), 'm'); % 火星的轨迹
plot3(y(:, 16), y(:, 17), y(:, 18), 'y'); % 木星的轨迹
plot3(y(:, 19), y(:, 20), y(:, 21), 'c'); % 土星的轨迹
plot3(y(:, 22), y(:, 23), y(:, 24), 'b--'); % 天王星的轨迹
plot3(y(:, 25), y(:, 26), y(:, 27), 'g--'); % 海王星的轨迹
axis equal;
grid on;
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('太阳系八大行星的运动轨迹');

% 计算行星在其轨道上的加速度
function dY = orbit(~, Y)
    mass = [1.989e30, 3.3e23, 4.87e24, 5.97e24, 6.42e23, 1.9e27, 5.68e26, 8.68e25]; % 太阳、水星、金星、地球、火星、木星、土星、天王星、海王星的质量（kg）
    pos = reshape(Y(1:27), 3, 9); % 八大行星的位置向量
    vel = reshape(Y(28:54), 3, 9); % 八大行星的速度向量
    G = 6.67430e-11; % 万有引力常数
    dY = zeros(54, 1);
    for i = 1:9
        F = [0; 0; 0]; % 初始力为0
        for j = 1:9
            if i ~= j
                r = pos(:, j) - pos(:, i); % 行星之间的距离
                F = F + G * mass(i) * mass(j) / norm(r)^3 * r; % 引力相互作用
            end
        end
        a = F / mass(i); % 行星在其轨道上的加速度
        dY((i - 1) * 3 + 1:i * 3) = vel(:, i); % 速度
        dY((i - 1) * 3 + 4:i * 3 + 3) = a; % 加速度
    end
end

在运行此程序之前，需要先将其保存为一个.m文件并在MATLAB中运行。程序运行后，将绘制出太阳系八大行星的运动轨迹。