由卫星的轨道参数得到卫星的位置和速度矢量代码

以下是一个示例代码，可以从卫星的轨道参数中计算出卫星的位置和速度矢量：

import numpy as np

# 定义常数
G = 6.6743e-11  # 地球引力常数
M = 5.97e24  # 地球质量
R = 6371000  # 地球半径

# 卫星轨道参数
a = 42164000  # 长半轴（m）
e = 0.5  # 偏心率
i = 45  # 轨道倾角（deg）
Omega = 30  # 升交点赤经（deg）
omega = 60  # 近地点幅角（deg）
M0 = 0  # 初始平近点角（deg）
t = 3600  # 时间（s）

# 计算轨道周期和平均角速度
T = 2 * np.pi * np.sqrt(a**3 / (G * M))
n = 2 * np.pi / T

# 计算真近点角
E0 = M0
E1 = M0 + e * np.sin(E0)
while abs(E1 - E0) > 1e-8:
    E0 = E1
    E1 = M0 + e * np.sin(E0)

theta = 2 * np.arctan(np.sqrt((1 + e) / (1 - e)) * np.tan(E1 / 2))

# 计算卫星在轨道平面上的位置
r = a * (1 - e**2) / (1 + e * np.cos(theta))
x = r * np.cos(theta)
y = r * np.sin(theta)

# 计算升交点处的卫星速度矢量
vx = np.sqrt(G * M / a) * (-np.sin(theta))
vy = np.sqrt(G * M / a) * (e + np.cos(theta))

# 将卫星位置和速度转换到地心惯性系
R3_Omega = np.array([[np.cos(np.deg2rad(Omega)), np.sin(np.deg2rad(Omega)), 0],
                     [-np.sin(np.deg2rad(Omega)), np.cos(np.deg2rad(Omega)), 0],
                     [0, 0, 1]])
R1_i = np.array([[1, 0, 0],
                 [0, np.cos(np.deg2rad(i)), np.sin(np.deg2rad(i))],
                 [0, -np.sin(np.deg2rad(i)), np.cos(np.deg2rad(i))]])
R3_omega = np.array([[np.cos(np.deg2rad(omega)), np.sin(np.deg2rad(omega)), 0],
                     [-np.sin(np.deg2rad(omega)), np.cos(np.deg2rad(omega)), 0],
                     [0, 0, 1]])

P = np.dot(R3_Omega, np.dot(R1_i, R3_omega))
r_ECI = np.dot(P, np.array([x, y, 0]))
v_ECI = np.dot(P, np.array([vx, vy, 0]))

# 计算卫星在地心惯性系中的位置和速度矢量
r_ECEF = np.dot(np.array([[np.cos(n * t), np.sin(n * t), 0],
                          [-np.sin(n * t), np.cos(n * t), 0],
                          [0, 0, 1]]), r_ECI)
v_ECEF = np.dot(np.array([[np.cos(n * t), np.sin(n * t), 0],
                          [-np.sin(n * t), np.cos(n * t), 0],
                          [0, 0, 1]]), v_ECI - np.cross(np.array([0, 0, n]), r_ECI))

# 计算卫星在地心固定系中的位置和速度矢量
GMST = 280.4606 + 360.9856473 * t / 86400
if GMST > 360:
    GMST -= 360
elif GMST < 0:
    GMST += 360

R3_GMST = np.array([[np.cos(np.deg2rad(GMST)), np.sin(np.deg2rad(GMST)), 0],
                    [-np.sin(np.deg2rad(GMST)), np.cos(np.deg2rad(GMST)), 0],
                    [0, 0, 1]])

r_ECEF_fixed = np.dot(R3_GMST, r_ECEF)
v_ECEF_fixed = np.dot(R3_GMST, v_ECEF)

# 将卫星位置和速度矢量转换为经纬度和海拔高度
lon = np.arctan2(r_ECEF_fixed[1], r_ECEF_fixed[0])
lat = np.arctan2(r_ECEF_fixed[2], np.sqrt(r_ECEF_fixed[0]**2 + r_ECEF_fixed[1]**2))
alt = np.sqrt(r_ECEF_fixed[0]**2 + r_ECEF_fixed[1]**2 + r_ECEF_fixed[2]**2) - R

print('经度（deg）：', np.rad2deg(lon))
print('纬度（deg）：', np.rad2deg(lat))
print('高度（m）：', alt)
print('速度（m/s）：', np.linalg.norm(v_ECEF_fixed))

这个代码使用了 Kepler 六个轨道参数来计算卫星的位置和速度矢量，然后将其转换到地心惯性系、地心固定系和经纬度高度坐标系。注意，这个代码中使用了一些简化，例如假设地球是一个完美的球体，没有考虑大气阻力等因素。实际应用中，可能需要更加复杂的模型来计算卫星的位置和速度。