import
math
from
scipy
.
optimize
import
fsolve
import
numpy
as
np
# 定义轨道元素
a
=
7000
# 半长轴 (km)
e
=
0.01
# 偏心率
i
=
math
.
radians
(
30
)
# 倾角 (rad)
Ω
=
math
.
radians
(
45
)
# 升交点赤经 (rad)
ω
=
math
.
radians
(
60
)
# 近地点幅角 (rad)
M
=
math
.
radians
(
20
)
# 平均近点角 (rad)
# 开普勒方程求解
def
kepler_equation
(
E
,
M
,
e
):
return
E
-
e
*
math
.
sin
(
E
)
-
M
# 求解偏近点角E
E_guess
=
M
# 初值
E_solution
=
fsolve
(
kepler_equation
,
E_guess
,
args
=(
M
,
e
))[
0
]
# 计算径向距离r
r
=
a
*
(
1
-
e
*
math
.
cos
(
E_solution
))
# 计算真近点角θ
theta
=
2
*
math
.
atan2
(
math
.
sqrt
(
1
+
e
)
*
math
.
sin
(
E_solution
/
2
),
math
.
sqrt
(
1
-
e
)
*
math
.
cos
(
E_solution
/
2
)
)
# 转换到惯性坐标系
x
=
r
*
(
math
.
cos
(Ω)
*
math
.
cos
(ω
+
theta
)
-
math
.
sin
(Ω)
*
math
.
sin
(ω
+
theta
)
*
math
.
cos
(
i
))
y
=
r
*
(
math
.
sin
(Ω)
*
math
.
cos
(ω
+
theta
)
+
math
.
cos
(Ω)
*
math
.
sin
(ω
+
theta
)
*
math
.
cos
(
i
))
z
=
r
*
(
math
.
sin
(
i
)
*
math
.
sin
(ω
+
theta
))
# 输出结果
print
(
f
"轨道径向距离 r: {r:.2f} km"
)
print
(
f
"真近点角 θ: {math.degrees(theta):.2f}°"
)
print
(
f
"惯性坐标系位置向量: x = {x:.2f} km, y = {y:.2f} km, z = {z:.2f} km"
)
结果解释
输入参数:轨道六要素由任务需求或观测数据提供。
输出:
径向距离 ( r ):航天器离地球中心的距离。
真近点角 ( θ ):航天器在轨道上的位置。
惯性坐标系位置向量 ([x, y, z]):航天器在地心惯性坐标系的位置。
加入轨道速度计算
要计算航天器在轨道上的速度,我们需要使用轨道力学中的速度公式。根据轨道元素,航天器的速度可以通过以下几个步骤进行计算。
import
math
from
scipy
.
optimize
import
fsolve
import
numpy
as
np
# 定义轨道元素
a
=
7000
# 半长轴 (km)
e
=
0.01
# 偏心率
i
=
math
.
radians
(
30
)
# 倾角 (rad)
Ω
=
math
.
radians
(
45
)
# 升交点赤经 (rad)
ω
=
math
.
radians
(
60
)
# 近地点幅角 (rad)
M
=
math
.
radians
(
20
)
# 平均近点角 (rad)
# 地球的标准引力参数 (km^3/s^2)
mu
=
398600.4418
# 开普勒方程求解
def
kepler_equation
(
E
,
M
,
e
):
return
E
-
e
*
math
.
sin
(
E
)
-
M
# 求解偏近点角E
E_guess
=
M
# 初值
E_solution
=
fsolve
(
kepler_equation
,
E_guess
,
args
=(
M
,
e
))[
0
]
# 计算径向距离r
r
=
a
*
(
1
-
e
*
math
.
cos
(
E_solution
))
# 计算轨道的角动量h
h
=
math
.
sqrt
(
a
*
(
1
-
e
**
2
)
*
mu
)
# 计算真近点角θ
theta
=
2
*
math
.
atan2
(
math
.
sqrt
(
1
+
e
)
*
math
.
sin
(
E_solution
/
2
),
math
.
sqrt
(
1
-
e
)
*
math
.
cos
(
E_solution
/
2
)
)
# 计算径向速度v_r
v_r
=
(
mu
/
h
)
*
e
*
math
.
sin
(
E_solution
)
# 计算切向速度v_t
v_t
=
h
/
r
# 计算总速度v
v
=
math
.
sqrt
(
v_r
**
2
+
v_t
**
2
)
# 转换到惯性坐标系
x
=
r
*
(
math
.
cos
(Ω)
*
math
.
cos
(ω
+
theta
)
-
math
.
sin
(Ω)
*
math
.
sin
(ω
+
theta
)
*
math
.
cos
(
i
))
y
=
r
*
(
math
.
sin
(Ω)
*
math
.
cos
(ω
+
theta
)
+
math
.
cos
(Ω)
*
math
.
sin
(ω
+
theta
)
*
math
.
cos
(
i
))
z
=
r
*
(
math
.
sin
(
i
)
*
math
.
sin
(ω
+
theta
))
# 输出结果
print
(
f
"轨道径向距离 r: {r:.2f} km"
)
print
(
f
"真近点角 θ: {math.degrees(theta):.2f}°"
)
print
(
f
"惯性坐标系位置向量: x = {x:.2f} km, y = {y:.2f} km, z = {z:.2f} km"
)
print
(
f
"径向速度 v_r: {v_r:.2f} km/s"
)
print
(
f
"切向速度 v_t: {v_t:.2f} km/s"
)
print
(
f
"总速度 v: {v:.2f} km/s"
)
代码说明
径向速度 (vr) 和 切向速度 (vt) 被分别计算,并最终得到总速度 (v)。
角动量 (h) 是通过轨道的半长轴和偏心率计算出来的,用于求解径向速度和切向速度。
最终,输出包括:
位置向量:在惯性坐标系下的位置。
径向速度、切向速度 和 总速度:航天器在轨道上的速度。
结果示例
轨道径向距离
r
:
6562.68
km
真近点角
θ:
22.87
°
惯性坐标系位置向量:
x
=
5844.37
km
,
y
=
1984.68
km
,
z
=
3304.82
km
径向速度
v_r
:
0.73
km
/
s
切向速度
v_t
:
7.74
km
/
s
总速度
v
:
7.80
km
/
s
总结
通过上述步骤,你可以计算出航天器在轨道上的位置和速度。轨道速度的计算对于航天器轨道控制、任务规划以及轨道保持至关重要。