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

以下是一个简单的 MATLAB 程序，用于使用 6 阶 RK 积分公式计算开普勒卫星的轨道积分：

function [t, r, v] = kepler_orbit(r0, v0, t0, tf, h)
% Inputs:
% r0: 初始位置向量
% v0: 初始速度向量
% t0: 初始时间
% tf: 结束时间
% h: 步长
% Outputs:
% t: 时间向量
% r: 位置矩阵
% v: 速度矩阵

% 常数
G = 6.67428e-11;
M = 5.97e24;
mu = G*M;

% 初始条件
r1 = r0;
v1 = v0;
t1 = t0;

% 初始化输出变量
t = t0:h:tf;
n = length(t);
r = zeros(3, n);
v = zeros(3, n);
r(:, 1) = r1;
v(:, 1) = v1;

% 迭代
for i = 1:n-1
    k1_r = v1;
    k1_v = -mu/norm(r1)^3*r1;
    k2_r = v1 + h/2*k1_v;
    k2_v = -mu/norm(r1+h/2*k1_r)^3*(r1+h/2*k1_r);
    k3_r = v1 + h/2*k2_v;
    k3_v = -mu/norm(r1+h/2*k2_r)^3*(r1+h/2*k2_r);
    k4_r = v1 + h*k3_v;
    k4_v = -mu/norm(r1+h*k3_r)^3*(r1+h*k3_r);
    r1 = r1 + h/6*(k1_r+2*k2_r+2*k3_r+k4_r);
    v1 = v1 + h/6*(k1_v+2*k2_v+2*k3_v+k4_v);
    r(:, i+1) = r1;
    v(:, i+1) = v1;
    t1 = t1 + h;
end

end

这个函数需要初始位置和速度向量，初始和结束时间以及步长作为输入，并返回轨道的时间、位置和速度向量作为输出。你可以像这样调用这个函数：

r0 = [10000; 0; 0];
v0 = [0; 8000; 0];
t0 = 0;
tf = 10000;
h = 10;
[t, r, v] = kepler_orbit(r0, v0, t0, tf, h);

请注意，这只是一个简单的例子，不包括任何误差控制或其他特殊功能。