Matlab定义开普勒轨道函数Kepler_orbit

以下是一个简单的Matlab函数，用于计算开普勒轨道的位置和速度：

function [r, v] = Kepler_orbit(mu, r0, v0, t)
% mu: 标准引力常数
% r0: 初始位置矢量
% v0: 初始速度矢量
% t: 时间向量

% 计算初始角动量和半径
h = cross(r0, v0);
r = norm(r0);

% 计算初始能量和半通径
E = 0.5*norm(v0)^2 - mu/r;
a = -mu/(2*E);

% 计算偏心率和真近点角
e = sqrt(1 - norm(h)^2/(mu*a));
theta0 = acos(dot(e, r0)/(e*r));

% 计算平近点角和平均角速度
M0 = theta0 - e*sin(theta0);
n = sqrt(mu/a^3);

% 计算解析解
M = M0 + n*t;
E = keplerEq(M, e);
theta = 2*atan(sqrt((1+e)/(1-e))*tan(E/2));
r = a*(1 - e^2)/(1 + e*cos(theta));
p = a*(1 - e^2);
h = sqrt(mu*p);
v = (1/h)*cross(h, r);

end

function E = keplerEq(M, e)
% 计算开普勒方程的解
tol = 1e-10;
E = M;
dE = 1;
while abs(dE) > tol
    dE = (M - E + e*sin(E))/(1 - e*cos(E));
    E = E + dE;
end
end

使用方法：

1. 定义标准引力常数mu、初始位置矢量r0、初始速度矢量v0和时间向量t。
2. 调用函数Kepler_orbit(mu, r0, v0, t)。
3. 返回值是位置向量r和速度向量v，分别与时间向量t对应。

注意：该函数仅适用于开普勒问题，即单个质点绕一个质量远远大于它的天体运动的情况。如果涉及到多个天体的相互作用，则需要使用更复杂的模型。