matlab求密切轨道根数

Williamlliw

已于 2023-11-19 17:09:59 修改

阅读量413

点赞数 9

文章标签： matlab 学习

于 2023-06-29 17:00:55 首次发布

本文链接：https://blog.csdn.net/m0_62150241/article/details/131455860

版权

文章介绍了如何在不考虑奇异情况下，基于航天器的位置和速度矢量计算轨道根数，包括升交点赤经、轨道倾角、近地点幅角、真近点角、轨道半长轴和离心率。通过角动量和几何关系简化了升交点赤经的计算，并详细说明了近地点幅角的求解过程，涉及MATLAB代码实现。

摘要由CSDN通过智能技术生成

航天器在轨道上运行时，根据一组确定的位置矢量r和速度矢量v可以确定一组确定的轨道根数，包括升交点赤经 $\Omega$ ，轨道倾角i，近地点幅角 w，真近点角 $\theta$ ，轨道半长轴 a，离心率 p。其中 $\Omega$ 、w、 $\theta$ 的范围均为0 ~ 360°，i 的范围是 0 ~ 180°。

根据教材中的公式把代码敲出来并不难，但教材中往往不会提及如何将反三角函数求得的0 ~ 180°角化到定义的 0 ~ 360° 中；此外，当出现奇异(如赤道面轨道，圆轨道)时的处理方法也常常不会提及，因此在此做一整理。

建立春分点惯性坐标系中。设中心天体的质量与引力常数之积为μ，该时刻航天器位置速度向量如下。

$\{\overrightarrow{r} \}_i = \{x_0\ y_0\ z_0\}^T$

$\{\overrightarrow{v} \}_i = \{v_0\ v_0\ v_0\}^T$

可计算得角动量

$\{\overrightarrow{H}\}_i = [\overrightarrow{r}]_i^X\{\overrightarrow{v}\}_i = \{H_x\ H_y\ H_z\}^T$

并求出半长轴和离心率

$E = \dfrac{1}{2}v^2 - \dfrac{\mu }{r} =-\dfrac{\mu }{2a}\ \therefore a = \dfrac{r}{2 - rv^2/\mu}$

$p = \dfrac{H^2}{\mu } \ \therefore e = \sqrt{1-\dfrac{H^2}{a\mu }}$

不考虑奇异情况下，计算其他参数的步骤如下。

$\cos i = \dfrac{H_z}{H}\ \therefore i = \arccos \dfrac{H_z}{H}$

下面是一个很精彩的几何关系，可以大大简化升交点赤经的计算，并能很方便地分析出 $\Omega$ 在 0 ~ 360° 间的范围。记角动量 H 在 XOY 面的投影为向量 E，向量 E 绕 z 轴逆时针旋转 90° 即为升交点N的方向，记为向量 N'，原因请读者自行推导。

$\{\overrightarrow{E} \}_i = \{H_x\ H_y\}^T$

$\{\overrightarrow{N'} \}_i = \{-H_y\ H_x\}^T$

从x轴顺时针转到向量 N' 的极角即为升交点赤经 $\Omega$ ，用 matlab 的 atan2 函数可以很快求出 $\Omega$ 。

下面计算真近点角f的余弦

$r = \dfrac{p}{1+e\cos f}\ \therefore \cos f = \dfrac{p/r-1}{e}$

然后鉴别 f 在 0 ~ 360° 中的区间，在不考虑奇异，航天器总是在一半的轨道上加速，在另一半的轨道上减速，对应引力与速度的夹角一半是锐角，另一半是钝角。记 $\theta$ 为航天器位置矢量r与速度矢量v间的夹角。若 $\theta$ ≤ 90°，航天器减速，在远离近地点段， $f\subset [0, 180^\circ]$ 。若 $\theta$ > 90°，航天器加速，在接近近地点段， $f\subset (180^\circ, 360^\circ)$ 。

$\theta \leqslant 90^\circ,\ f = \arccos\dfrac{p/r-1}{e}$

$\theta > 90^\circ,\ f = -\arccos\dfrac{p/r-1}{e}+360^\circ$

最后计算近地点幅角 w。记轨道节点单位节点矢量为N，

$\{\overrightarrow{N} \}_i = \dfrac{1}{\sqrt{H_x^2+H_y^2}}\{-H_y\ H_x\ 0\}^T$

则轨道幅角 u 的余弦值为

$\cos u = \dfrac{\overrightarrow{r}\cdot \overrightarrow{N}}{\overrightarrow{r}}$

然后甄别u的范围，当航天器在赤道面上方时，u 在 0 ~ 180° 间；反之则在 180° ~ 360°间，有

$z_0\geqslant 0,u = \arccos\dfrac{\overrightarrow{r}\cdot \overrightarrow{N}}{\overrightarrow{r}}$

$z_0\geqslant 0,u = -\arccos\dfrac{\overrightarrow{r}\cdot \overrightarrow{N}}{\overrightarrow{r}}+360^\circ$

近地点幅角w = u - f，当其值小于0时加360°，可得6个密切轨道要素。

在不考虑奇异情况下的密切轨道计算代码如下：

% 定义参数
clear;
miu = 3.986004415e5; % 引力常数
% 初始位置和速度，修改该处初值值以得到不同结果
r0 = [-4109.2; 4789.7; 2257.6];
v0 = [6.8; 4.4; 0.1];
r0m = norm(r0);
v0m = norm(v0);
% 半长轴
a = r0m/(2 - r0m*v0m^2/miu);
H = cross(r0, v0);
p = norm(H)^2/miu;
% 离心率
e = sqrt(1 - p/a);
i = acos(H(3)/norm(H));
% 升交点赤经
omega = atan2(H(1), -H(2));
omega = mod(omega, 2*pi);
N = [-H(2); H(1); 0]/sqrt(H(2)^2 + H(1)^2);
% 近地点幅角
f = acos((p-r0m)/(r0m*e)); % 暂时
u = acos(r0'*N/norm(r0)); % 暂时
if r0(3) < 0
    u = -u + 2*pi;
end
thetaRV = acos(r0'*v0/(r0m * v0m));
if thetaRV > pi/2
    f = -f + 2*pi;
end
% 近地点幅角
w = u - f;
w = mod(w, 2*pi);
E = 2*atan(sqrt((1-e)/(1+e)) * tan(f/2));
E = mod(E, 2*pi);
M = E - e*sin(E);
fprintf('omega = %.2f °, i = %.2f °\na = %.2f km, e = %.2f\nw = %.2f °, M = %.2f °\n', omega/pi*180, i/pi*180, a, e, w/pi*180, M/pi*180);

点赞超过10就更一下考虑奇异和非椭圆轨道时的轨道根数代码