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

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

function [t,x,y,z,vx,vy,vz] = kepler_orbit(t0, tf, dt, x0, y0, z0, vx0, vy0, vz0, mu)
% t0: 初值时间
% tf: 终值时间
% dt: 步长
% x0, y0, z0: 初始位置
% vx0, vy0, vz0: 初始速度
% mu: 引力常数

% 初始化位置和速度
x = x0;
y = y0;
z = z0;
vx = vx0;
vy = vy0;
vz = vz0;

% 初始化时间
t = t0;

% 计算加速度
r = norm([x y z]);
ax = -mu*x/r^3;
ay = -mu*y/r^3;
az = -mu*z/r^3;

% 进行RK5积分
while t <= tf
    kx1 = vx;
    ky1 = vy;
    kz1 = vz;
    kvx1 = ax;
    kvy1 = ay;
    kvz1 = az;
    
    kx2 = vx + 0.5*dt*kvx1;
    ky2 = vy + 0.5*dt*kvy1;
    kz2 = vz + 0.5*dt*kvz1;
    [ax2, ay2, az2] = calc_acceleration(mu, x + 0.5*dt*kx1, y + 0.5*dt*ky1, z + 0.5*dt*kz1);
    kvx2 = ax2;
    kvy2 = ay2;
    kvz2 = az2;
    
    kx3 = vx + 0.5*dt*kvx2;
    ky3 = vy + 0.5*dt*kvy2;
    kz3 = vz + 0.5*dt*kvz2;
    [ax3, ay3, az3] = calc_acceleration(mu, x + 0.5*dt*kx2, y + 0.5*dt*ky2, z + 0.5*dt*kz2);
    kvx3 = ax3;
    kvy3 = ay3;
    kvz3 = az3;
    
    kx4 = vx + dt*kvx3;
    ky4 = vy + dt*kvy3;
    kz4 = vz + dt*kvz3;
    [ax4, ay4, az4] = calc_acceleration(mu, x + dt*kx3, y + dt*ky3, z + dt*kz3);
    kvx4 = ax4;
    kvy4 = ay4;
    kvz4 = az4;
    
    kx5 = vx + (1/6)*dt*(kvx1 + 2*kvx2 + 2*kvx3 + kvx4);
    ky5 = vy + (1/6)*dt*(kvy1 + 2*kvy2 + 2*kvy3 + kvy4);
    kz5 = vz + (1/6)*dt*(kvz1 + 2*kvz2 + 2*kvz3 + kvz4);
    [ax5, ay5, az5] = calc_acceleration(mu, x + (1/6)*dt*(kx1 + 2*kx2 + 2*kx3 + kx4), y + (1/6)*dt*(ky1 + 2*ky2 + 2*ky3 + ky4), z + (1/6)*dt*(kz1 + 2*kz2 + 2*kz3 + kz4));
    kvx5 = ax5;
    kvy5 = ay5;
    kvz5 = az5;
    
    x = x + (1/6)*dt*(kx1 + 2*kx2 + 2*kx3 + kx4);
    y = y + (1/6)*dt*(ky1 + 2*ky2 + 2*ky3 + ky4);
    z = z + (1/6)*dt*(kz1 + 2*kz2 + 2*kz3 + kz4);
    vx = vx + (1/6)*dt*(kvx1 + 2*kvx2 + 2*kvx3 + kvx4);
    vy = vy + (1/6)*dt*(kvy1 + 2*kvy2 + 2*kvy3 + kvy4);
    vz = vz + (1/6)*dt*(kvz1 + 2*kvz2 + 2*kvz3 + kvz4);
    
    t = t + dt;
end

end

function [ax, ay, az] = calc_acceleration(mu, x, y, z)
% 计算加速度
r = norm([x y z]);
ax = -mu*x/r^3;
ay = -mu*y/r^3;
az = -mu*z/r^3;
end

使用方法：

t0 = 0;
tf = 86400*365; % 一年的时间
dt = 86400; % 一天的时间
x0 = 1.4710e+11;
y0 = 0;
z0 = 0;
vx0 = 0;
vy0 = 3.0287e+04;
vz0 = 0;
mu = 1.3271e+20;

[t,x,y,z,vx,vy,vz] = kepler_orbit(t0, tf, dt, x0, y0, z0, vx0, vy0, vz0, mu);

% 绘制轨道图
plot3(x,y,z);
xlabel('x');
ylabel('y');
zlabel('z');
title('Kepler Orbit');