开普勒 计算卫星轨道 代码

以下是开普勒计算卫星轨道的 Python 代码：

import math

# Define constants
G = 6.67259e-11  # Gravitational constant (m^3 kg^-1 s^-2)
M = 5.97219e24   # Mass of the Earth (kg)
R = 6371.0e3     # Radius of the Earth (m)

# Define initial conditions
a = float(input("Enter semi-major axis (m): "))
e = float(input("Enter eccentricity: "))
i = float(input("Enter inclination (deg): "))
Omega = float(input("Enter right ascension of ascending node (deg): "))
omega = float(input("Enter argument of perigee (deg): "))
M0 = float(input("Enter mean anomaly at epoch (deg): "))
t0 = float(input("Enter epoch time (s): "))

# Convert angles to radians
i = math.radians(i)
Omega = math.radians(Omega)
omega = math.radians(omega)
M0 = math.radians(M0)

# Calculate derived parameters
n = math.sqrt(G * M / a**3)  # Mean motion (rad/s)
T = 2 * math.pi / n         # Period (s)
t = float(input("Enter time since epoch (s): "))
M = M0 + n * t             # Mean anomaly (rad)

# Use Newton-Raphson iteration to solve Kepler's equation for E
E0 = M
E1 = M + e * math.sin(E0)
while abs(E1 - E0) > 1e-8:
    E0 = E1
    E1 = M + e * math.sin(E0)
    
E = E1

# Calculate true anomaly (nu) and distance from Earth (r)
nu = 2 * math.atan(math.sqrt((1 + e) / (1 - e)) * math.tan(E / 2))
r = a * (1 - e**2) / (1 + e * math.cos(nu))

# Calculate position and velocity vectors in the orbital plane
x_orb = r * math.cos(nu)
y_orb = r * math.sin(nu)
vx_orb = math.sqrt(G * M / a) * (-math.sin(E))
vy_orb = math.sqrt(G * M / a) * (math.sqrt(1 - e**2) * math.cos(E))

# Transform position and velocity vectors to ECI coordinates
x_ECI = x_orb * (math.cos(omega) * math.cos(Omega) - math.sin(omega) * math.sin(Omega) * math.cos(i)) - \
        y_orb * (math.sin(omega) * math.cos(Omega) + math.cos(omega) * math.sin(Omega) * math.cos(i))
y_ECI = x_orb * (math.cos(omega) * math.sin(Omega) + math.sin(omega) * math.cos(Omega) * math.cos(i)) + \
        y_orb * (math.cos(omega) * math.cos(Omega) * math.cos(i) - math.sin(omega) * math.sin(Omega))
z_ECI = x_orb * (math.sin(omega) * math.sin(i)) + \
        y_orb * (math.cos(omega) * math.sin(i))
vx_ECI = vx_orb * (math.cos(omega) * math.cos(Omega) - math.sin(omega) * math.sin(Omega) * math.cos(i)) - \
         vy_orb * (math.sin(omega) * math.cos(Omega) + math.cos(omega) * math.sin(Omega) * math.cos(i))
vy_ECI = vx_orb * (math.cos(omega) * math.sin(Omega) + math.sin(omega) * math.cos(Omega) * math.cos(i)) + \
         vy_orb * (math.cos(omega) * math.cos(Omega) * math.cos(i) - math.sin(omega) * math.sin(Omega))
vz_ECI = vx_orb * (math.sin(omega) * math.sin(i)) + \
         vy_orb * (math.cos(omega) * math.sin(i))

# Print results
print("Position vector (m):")
print("[{:.2f}, {:.2f}, {:.2f}]".format(x_ECI, y_ECI, z_ECI))
print("Velocity vector (m/s):")
print("[{:.2f}, {:.2f}, {:.2f}]".format(vx_ECI, vy_ECI, vz_ECI))

这个代码使用了开普勒定律和牛顿-拉弗森迭代法来计算卫星的位置和速度，并将结果从轨道坐标系转换到 ECI 坐标系。你需要输入卫星的初始条件和时间，然后程序会输出卫星在 ECI 坐标系下的位置和速度。