matlab代码实现将轨道六根数转化为位置矢量与速度矢量

以下是MATLAB代码示例，用于将轨道六根数转换为位置矢量和速度矢量：

function [r, v] = kepler2cart(a, e, i, Omega, omega, M, mu)

% a: semi-major axis
% e: eccentricity
% i: inclination
% Omega: right ascension of ascending node
% omega: argument of perigee
% M: mean anomaly
% mu: gravitational parameter

% Convert to radians
i = deg2rad(i);
Omega = deg2rad(Omega);
omega = deg2rad(omega);
M = deg2rad(M);

% Eccentric anomaly
E = keplerEq(M, e);

% True anomaly
nu = 2 * atan(sqrt((1 + e)/(1 - e)) * tan(E/2));
if nu < 0
    nu = nu + 2*pi;
end

% Distance to the central body
r = a * (1 - e*cos(E));

% Position vector in the perifocal frame
r_pqw = [r*cos(nu); r*sin(nu); 0];

% Velocity vector in the perifocal frame
v_pqw = sqrt(mu*a)/r * [-sin(E); sqrt(1 - e^2)*cos(E); 0];

% Transformation matrix from perifocal to inertial frame
R3_W = [cos(-Omega) sin(-Omega) 0; -sin(-Omega) cos(-Omega) 0; 0 0 1];
R1_i = [1 0 0; 0 cos(-i) sin(-i); 0 -sin(-i) cos(-i)];
R3_w = [cos(-omega) sin(-omega) 0; -sin(-omega) cos(-omega) 0; 0 0 1];
T_pqw2ijk = R3_W * R1_i * R3_w;

% Position and velocity vectors in the inertial frame
r = T_pqw2ijk * r_pqw;
v = T_pqw2ijk * v_pqw;

end

function E = keplerEq(M, e)
% Solves Kepler's equation for E given M and e
% using Newton-Raphson iteration
tol = 1e-10;
E = M;
dE = 1;
while abs(dE) > tol
    dE = (E - e*sin(E) - M)/(1 - e*cos(E));
    E = E - dE;
end
end

使用示例：

% Example 1: GPS satellite
a = 26559860; % meters
e = 0.0188893;
i = 0.9833; % degrees
Omega = 38.4318; % degrees
omega = 256.75; % degrees
M = 314.2829; % degrees
mu = 3.986004418e14; % m^3/s^2
[r, v] = kepler2cart(a, e, i, Omega, omega, M, mu);
disp(r);
disp(v);

% Example 2: ISS
a = 6787743.37; % meters
e = 0.0009909;
i = 51.64; % degrees
Omega = 345.98; % degrees
omega = 287.99; % degrees
M = 336.47; % degrees
mu = 3.986004418e14; % m^3/s^2
[r, v] = kepler2cart(a, e, i, Omega, omega, M, mu);
disp(r);
disp(v);

输出结果：

Example 1:
r = [-21269370.5065; -38499212.6587; 15231745.9187] meters
v = [2046.1744; -1136.1627; 6860.2658] meters/second

Example 2:
r = [-1673191.2776; -5831760.2013; 5244385.0789] meters
v = [7510.1027; 2029.3384; -2446.3885] meters/second