利用4阶RK积分计算开普勒卫星轨道积分的MATLAB程序

以下是一个简单的 MATLAB 程序，使用4阶RK积分方法来计算开普勒卫星轨道积分。程序包括输入轨道参数、计算轨道和绘制结果的过程。

%% 输入轨道参数
a = 7000;   % 长半轴（km）
e = 0.1;    % 离心率
i = 30;     % 轨道倾角（度）
O = 0;      % 升交点赤经（度）
w = 0;      % 近地点角（度）
M = 0;      % 平近点角（度）
tspan = [0 1]*3600; % 计算时间范围（秒）

%% 计算轨道
mu = 3.986e5; % 地球引力常数（km^3/s^2）
n = sqrt(mu/a^3); % 平均运动（rad/s）
E0 = M2E(M, e); % 初始偏近点角（弧度）
y0 = [a*(cos(E0)-e); a*sqrt(1-e^2)*sin(E0); 0; n]; % 初始状态向量
options = odeset('RelTol',1e-8,'AbsTol',1e-8); % ODE求解选项
[t,y] = ode45(@(t,y) keplerODE(t,y,mu), tspan, y0, options); % 用 RK4 求解 ODE

%% 绘制结果
figure
plot3(y(:,1),y(:,2),y(:,3))
axis equal
xlabel('X (km)')
ylabel('Y (km)')
zlabel('Z (km)')
title('Orbit of Kepler Satellite')

%% 帮助函数
function E = M2E(M, e)
    E0 = M;
    E = M + e*sin(E0);
    while abs(E - E0) > 1e-8
        E0 = E;
        E = M + e*sin(E0);
    end
end

function dydt = keplerODE(t,y,mu)
    r = norm(y(1:3)); % 卫星距离地心的距离
    dydt = [y(4); y(5); y(6); -mu/r^3*y(1); -mu/r^3*y(2); -mu/r^3*y(3)];
end

这个程序定义了输入轨道参数，例如长半轴、离心率、轨道倾角等。它使用 `M2E` 函数将平近点角转换为偏近点角，并使用 `ode45` 函数计算轨道。最后，它通过绘制卫星的轨道来可视化结果。