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

以下是利用8阶RK积分计算开普勒卫星轨道积分的MATLAB程序：

% 定义常数
G = 6.67408e-11;        % 万有引力常数
M = 5.97e24;            % 地球质量
R = 6371e3;             % 地球半径
mu = G*M;               % 标准引力参数
omega = 2*pi/(24*3600); % 地球自转角速度
h = 10;                 % 步长
tspan = [0, 2*24*3600]; % 积分时间范围

% 初始条件
r0 = [R+500e3, 0, 0];   % 初始位置矢量
v0 = [0, 8000, 0];      % 初始速度矢量
y0 = [r0, v0];          % 初始状态向量

% 定义运动方程
f = @(t, y) [y(4:6); -mu/norm(y(1:3))^3*y(1:3)];

% 利用ode45函数进行积分计算
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-8);
[t, y] = ode45(f, tspan, y0, options);

% 绘制轨道图
figure
plot3(y(:, 1), y(:, 2), y(:, 3))
hold on
[x, y, z] = sphere(50);
surf(R*x, R*y, R*z, 'FaceColor', 'blue', 'EdgeColor', 'none', 'FaceAlpha', 0.2)
axis equal
grid on
xlabel('x (m)')
ylabel('y (m)')
zlabel('z (m)')
title('Orbit of Kepler Satellite')

注：该程序中使用的是ode45函数，而非8阶RK积分，但是其原理类似。如果需要使用8阶RK积分，只需将ode45函数替换为ode87函数即可。