使用matlab代码实现轨道六根数向位置速度矢量的转换

以下是使用Matlab代码实现轨道六根数向位置速度矢量的转换的示例：

% 定义常数
mu = 398600.4418; % 地球引力常数

% 定义轨道六根数
a = 7000; % 半长轴（km）
e = 0.1; % 偏心率
i = 30; % 轨道倾角（deg）
Omega = 45; % 升交点赤经（deg）
omega = 60; % 近地点幅角（deg）
M = 0; % 平近点角（deg）

% 计算轨道周期
T = 2 * pi * sqrt(a^3 / mu);

% 计算平近点角速度
n = 2 * pi / T;

% 计算真近点角
E = M;
while abs(E - e * sin(E) - M) > 1e-8
    E = E - (E - e * sin(E) - M) / (1 - e * cos(E));
end
theta = 2 * atan(sqrt((1 + e) / (1 - e)) * tan(E / 2));

% 计算升交点经度
Omega_rad = deg2rad(Omega);
omega_rad = deg2rad(omega);
i_rad = deg2rad(i);
theta_rad = deg2rad(theta);
x = a * (cos(theta) - e);
y = a * sqrt(1 - e^2) * sin(theta);
r_pqw = [x; y; 0];
v_pqw = sqrt(mu / a) / (1 - e * cos(theta)) * [-sin(theta); sqrt(1 - e^2) * cos(theta); 0];
R_3_Omega = [cos(Omega_rad) sin(Omega_rad) 0; -sin(Omega_rad) cos(Omega_rad) 0; 0 0 1];
R_1_i = [1 0 0; 0 cos(i_rad) sin(i_rad); 0 -sin(i_rad) cos(i_rad)];
R_3_omega = [cos(omega_rad) sin(omega_rad) 0; -sin(omega_rad) cos(omega_rad) 0; 0 0 1];
Q_pX = R_3_Omega * R_1_i * R_3_omega;
r_XYZ = Q_pX * r_pqw;
v_XYZ = Q_pX * v_pqw;

% 显示结果
disp("位置矢量（km）：");
disp(r_XYZ);
disp("速度矢量（km/s）：");
disp(v_XYZ);

其中，输入的六根数为半长轴 `a`（单位：千米）、偏心率 `e`、轨道倾角 `i`（单位：度）、升交点赤经 `Omega`（单位：度）、近地点幅角 `omega`（单位：度）和平近点角 `M`（单位：度）。输出的位置矢量和速度矢量均为三维向量，单位分别为千米和千米每秒。