使用 Python 进行低轨卫星星座轨道力学计算的完整模拟方法。内容涵盖基础理论（牛顿万有引力、开普勒轨道）、轨道参数与坐标系统、二体问题与摄动模型（J2、大气阻力、太阳辐射压）。文章提供了完整的 Python 代码实现，包括轨道计算类、数值积分器、Walker 星座设计、可见性分析与链路预算工具。此外，还展示了 3D 可视化、并行计算优化及 GPU 加速方案，并通过 Starlink 星座模拟案例验证了框架的有效性。该框架适用于星座设计、任务规划、教育培训及科研分析。
本文详细介绍了使用 Python 进行低轨卫星星座轨道力学计算的完整模拟方法。我们将从基础理论出发，逐步构建一个完整的轨道模拟系统，涵盖卫星轨道参数计算、星座构型设计、轨道摄动模型以及可见性分析等关键方面。
低地球轨道（LEO）卫星星座已成为全球通信、地球观测和导航系统的关键基础设施。星座由数十至数千颗卫星组成，如 Starlink、OneWeb 等系统。精确的轨道力学计算对于星座设计、卫星管理和通信链路维护至关重要。
Python 凭借其强大的科学计算库（如 NumPy、SciPy）和可视化工具（如 Matplotlib），成为轨道力学模拟的理想选择。本文将构建一个完整的轨道模拟框架，涵盖从基础理论到实际应用的各个方面。
轨道力学的核心是牛顿万有引力定律：
F = G * m1 * m2 / r^2
对于地球与卫星系统：
a = -mu * r / r^3
其中 mu = GM 是地球引力常数，约为 3.986×10^14 m^3/s^2。
卫星在地球引力场中的运动方程为：
d^2r/dt^2 = -mu * r / r^3 + apert
其中 apert 表示各种摄动加速度。
卫星轨道通常用六个轨道要素描述：
轨道计算涉及多个坐标系统：
对于椭圆轨道，开普勒方程为：
M = E - e * sin(E)
其中 M 是平近点角，E 是偏近点角。
根据轨道要素计算卫星位置和速度：
import numpy as np
from math import sin, cos, sqrt, radians, degrees
class KeplerOrbit:
"""开普勒轨道计算类"""
def __init__(self, a, e, i, raan, arg_peri, nu, mu=3.986004418e14):
"""
初始化轨道要素
参数:
a: 半长轴 (m)
e: 偏心率
i: 轨道倾角 (度)
raan: 升交点赤经 (度)
arg_peri: 近地点幅角 (度)
nu: 真近点角 (度)
mu: 地球引力常数 (m^3/s^2)
"""
self.a = a
self.e = e
self.i = radians(i)
self.raan = radians(raan)
self.arg_peri = radians(arg_peri)
self.nu = radians(nu)
self.mu = mu
def calculate_position_velocity(self):
"""计算卫星在地心惯性坐标系中的位置和速度"""
# 1. 计算轨道平面内的位置和速度
# 半正焦弦 p
p = self.a * (1 - self.e**2)
# 轨道平面内的极坐标位置
r = p / (1 + self.e * cos(self.nu))
# 位置在轨道平面内的分量
r_perifocal = np.array([r * cos(self.nu), r * sin(self.nu), 0.0])
# 速度在轨道平面内的分量
v_perifocal = np.array([
-sin(self.nu) * sqrt(self.mu / p),
(self.e + cos(self.nu)) * sqrt(self.mu / p),
0.0
])
# 2. 转换矩阵：从轨道平面到地心惯性坐标系
# 绕 Z 轴旋转 (-raan)
R3_raan = np.array([
[cos(-self.raan), -sin(-self.raan), 0],
[sin(-self.raan), cos(-self.raan), 0],
[0, 0, 1]
])
# 绕 X 轴旋转 (-i)
R1_i = np.array([
[1, 0, 0],
[0, cos(-self.i), -sin(-self.i)],
[0, sin(-self.i), cos(-self.i)]
])
# 绕 Z 轴旋转 (-arg_peri)
R3_arg_peri = np.array([
[cos(-self.arg_peri), -sin(-self.arg_peri), 0],
[sin(-self.arg_peri), cos(-self.arg_peri), 0],
[0, 0, 1]
])
# 完整的转换矩阵 Q
Q = R3_raan @ R1_i @ R3_arg_peri
# 3. 转换到地心惯性坐标系
r_eci = Q @ r_perifocal
v_eci = Q @ v_perifocal
return r_eci, v_eci
def calculate_orbital_period(self):
"""计算轨道周期"""
return 2 * np.pi * sqrt(self.a**3 / self.mu)
def propagate_orbit(self, delta_t):
"""传播轨道到新时间点
参数:
delta_t: 时间增量 (秒)
"""
# 计算平近点角变化
n = sqrt(self.mu / self.a**3)
# 平均运动
M0 = self.eccentric_to_mean(self.true_to_eccentric(self.nu))
M = M0 + n * delta_t
# 解开普勒方程求偏近点角
E = self.solve_kepler(M, self.e)
# 计算新的真近点角
nu_new = self.eccentric_to_true(E)
# 更新轨道要素
return KeplerOrbit(
self.a, self.e, degrees(self.i), degrees(self.raan),
degrees(self.arg_peri), degrees(nu_new), self.mu
)
@staticmethod
def true_to_eccentric(nu, e):
"""真近点角转偏近点角"""
return 2 * np.arctan(np.sqrt((1 - e) / (1 + e)) * np.tan(nu / 2))
@staticmethod
def eccentric_to_true(E, e):
"""偏近点角转真近点角"""
return 2 * np.arctan(np.sqrt((1 + e) / (1 - e)) * np.tan(E / 2))
@staticmethod
def eccentric_to_mean(E, e):
"""偏近点角转平近点角"""
return E - e * np.sin(E)
@staticmethod
def solve_kepler(M, e, tol=1e-12, max_iter=100):
"""解开普勒方程 E - e*sin(E) = M
参数:
M: 平近点角
e: 偏心率
tol: 容差
max_iter: 最大迭代次数
返回:
E: 偏近点角
"""
# 初始估计
if e < 0.8:
E = M
else:
E = np.pi
# 牛顿 - 拉弗森迭代
for i in range(max_iter):
f = E - e * np.sin(E) - M
f_prime = 1 - e * np.cos(E)
delta = f / f_prime
E -= delta
if abs(delta) < tol:
break
return E
实际轨道受多种摄动影响，主要考虑：
地球非球形导致的主要摄动项：
U_J2 = (mu * J2 * Re^2) / (2 * r^3) * (3 * sin^2(phi) - 1)
其中 J2 ≈ 1.08263×10^-3，Re 为地球赤道半径。
class PerturbationModel:
"""轨道摄动模型"""
def __init__(self, J2=1.08263e-3, R_e=6378137.0):
self.J2 = J2
self.R_e = R_e
def J2_acceleration(self, r_eci):
"""计算 J2 摄动加速度
参数:
r_eci: 卫星在地心惯性坐标系中的位置向量 (m)
返回:
a_J2: J2 摄动加速度向量 (m/s^2)
"""
x, y, z = r_eci
r = np.linalg.norm(r_eci)
# 地球引力常数
mu = 3.986004418e14
# J2 摄动加速度分量
factor = 1.5 * self.J2 * mu * self.R_e**2 / r**5
a_x = factor * x * (5 * z**2 / r**2 - 1)
a_y = factor * y * (5 * z**2 / r**2 - 1)
a_z = factor * z * (5 * z**2 / r**2 - 3)
return np.array([a_x, a_y, a_z])
def atmospheric_drag(self, r_eci, v_eci, satellite_params):
"""大气阻力摄动
参数:
r_eci: 位置向量 (m)
v_eci: 速度向量 (m/s)
satellite_params: 卫星参数字典
返回:
a_drag: 大气阻力加速度 (m/s^2)
"""
# 提取卫星参数
Cd = satellite_params.get('drag_coefficient', 2.2)
A = satellite_params.get('cross_sectional_area', 1.0) # m^2
m = satellite_params.get(, )
r = np.linalg.norm(r_eci)
h = r - .R_e
rho0 =
H =
rho = rho0 * np.exp(-h / H)
v_rel = np.linalg.norm(v_eci)
a_drag_magnitude = - * Cd * (A/m) * rho * v_rel**
v_rel > :
a_drag = a_drag_magnitude * (-v_eci / v_rel)
:
a_drag = np.zeros()
a_drag
():
P_sun =
Cr = satellite_params.get(, )
A = satellite_params.get(, )
m = satellite_params.get(, )
sun_direction = np.array([, , ])
sat_to_sun = sun_direction - r_eci / np.linalg.norm(r_eci)
sat_to_sun = sat_to_sun / np.linalg.norm(sat_to_sun)
a_srp_magnitude = -Cr * (A/m) * P_sun
a_srp_magnitude * sat_to_sun
对于包含摄动的精确轨道计算，需要使用数值积分：
class NumericalIntegrator:
"""轨道数值积分器"""
def __init__(self, perturbation_model=None):
self.perturbation_model = perturbation_model or PerturbationModel()
self.mu = 3.986004418e14
def ode_equations(self, t, y, satellite_params):
"""轨道运动微分方程
参数:
t: 时间 (s)
y: 状态向量 [x, y, z, vx, vy, vz]
satellite_params: 卫星参数
返回:
dydt: 状态向量导数
"""
# 位置和速度
r = y[:3]
v = y[3:]
# 二体问题加速度
r_norm = np.linalg.norm(r)
a_two_body = -self.mu * r / r_norm**3
# 总加速度（初始化为二体加速度）
a_total = a_two_body
# 添加摄动加速度
if self.perturbation_model:
# J2 摄动
a_J2 = self.perturbation_model.J2_acceleration(r)
a_total += a_J2
# 大气阻力（仅适用于低轨）
r_norm_km = r_norm / 1000
if r_norm_km < 2000: # 高度小于 2000km
a_drag = self.perturbation_model.atmospheric_drag(r, v, satellite_params)
a_total += a_drag
# 太阳辐射压力
a_srp = self.perturbation_model.solar_radiation_pressure(r, satellite_params)
a_total += a_srp
# 状态导数
dydt = np.zeros(6)
dydt[:3] = v # 位置导数 = 速度
dydt[3:] = a_total
dydt
():
t0, tf = t_span
n_steps = ((tf - t0) / dt) +
t_array = np.linspace(t0, tf, n_steps)
state_history = np.zeros((n_steps, ))
state_history[] = initial_state
method == :
i (n_steps - ):
t = t_array[i]
y = state_history[i]
k1 = .ode_equations(t, y, satellite_params)
k2 = .ode_equations(t + dt/, y + dt/ * k1, satellite_params)
k3 = .ode_equations(t + dt/, y + dt/ * k2, satellite_params)
k4 = .ode_equations(t + dt, y + dt * k3, satellite_params)
state_history[i+] = y + dt/ * (k1 + *k2 + *k3 + k4)
method == :
i ():
t = t_array[i]
y = state_history[i]
k1 = .ode_equations(t, y, satellite_params)
k2 = .ode_equations(t + dt/, y + dt/ * k1, satellite_params)
k3 = .ode_equations(t + dt/, y + dt/ * k2, satellite_params)
k4 = .ode_equations(t + dt, y + dt * k3, satellite_params)
state_history[i+] = y + dt/ * (k1 + *k2 + *k3 + k4)
i (, n_steps - ):
t = t_array[i]
y = state_history[i]
f_0 = .ode_equations(t_array[i], state_history[i], satellite_params)
f_1 = .ode_equations(t_array[i-], state_history[i-], satellite_params)
f_2 = .ode_equations(t_array[i-], state_history[i-], satellite_params)
f_3 = .ode_equations(t_array[i-], state_history[i-], satellite_params)
state_history[i+] = y + dt/ * (*f_0 - *f_1 + *f_2 - *f_3)
t_array, state_history
Walker 星座是最常用的星座构型，由三个参数描述：
class WalkerConstellation:
"""Walker 星座设计类"""
def __init__(self, T, P, F, altitude, inclination, raan0=0):
"""
初始化 Walker 星座
参数:
T: 卫星总数
P: 轨道面数
F: 相位因子 (0 到 P-1 之间的整数)
altitude: 轨道高度 (km)
inclination: 轨道倾角 (度)
raan0: 第一个轨道面的升交点赤经 (度)
"""
self.T = T
self.P = P
self.F = F
self.altitude = altitude # km
self.inclination = inclination # 度
self.raan0 = raan0 # 度
# 计算每个轨道面的卫星数
self.S = T // P # 每个轨道面的卫星数
# 计算轨道半长轴
R_e = 6371.0 # 地球半径 (km)
self.a = (R_e + altitude) * 1000 # 转换为米
# 生成星座中的所有卫星
self.satellites = self.generate_constellation()
def generate_constellation(self):
"""生成 Walker 星座中的所有卫星"""
satellites = []
# 每个轨道面
for p in range(self.P):
# 轨道面升交点赤经
raan = self.raan0 + p * (360.0 / self.P)
# 每个轨道面内的卫星
for s in (.S):
phase = ( / .T) * (s * .P + p * .F)
nu = phase
orbit = KeplerOrbit(
a=.a, e=, i=.inclination,
raan=raan, arg_peri=,
nu=nu
)
satellites.append({
: p,
: s,
: phase,
: orbit
})
satellites
():
satellite = .satellites[satellite_idx]
time != :
orbit = satellite[].propagate_orbit(time)
:
orbit = satellite[]
orbit.calculate_position_velocity()
():
coverage_stats = {
: (ground_points),
: min_elevation,
: [],
: [],
: []
}
lat, lon ground_points:
time_step =
total_time = *
n_steps = total_time // time_step
coverage_mask = np.zeros(n_steps, dtype=)
step (n_steps):
t = step * time_step
visible =
sat .satellites:
orbit = sat[].propagate_orbit(t)
r_eci, _ = orbit.calculate_position_velocity()
elevation = .calculate_elevation(r_eci, lat, lon, t)
elevation >= min_elevation:
visible =
coverage_mask[step] = visible
coverage_percentage = np.(coverage_mask) / n_steps *
gaps = .find_coverage_gaps(coverage_mask, time_step)
gaps:
avg_gap = np.mean(gaps)
max_gap = np.(gaps)
:
avg_gap =
max_gap =
coverage_stats[].append(coverage_percentage)
coverage_stats[].append(avg_gap)
coverage_stats[].append(max_gap)
coverage_stats
():
omega_e =
lat_rad = np.radians(lat)
lon_rad = np.radians(lon)
R_e =
lon_rad_rotated = lon_rad + omega_e * t
r_ground_eci = np.array([
R_e * np.cos(lat_rad) * np.cos(lon_rad_rotated),
R_e * np.cos(lat_rad) * np.sin(lon_rad_rotated),
R_e * np.sin(lat_rad)
])
r_rel = r_eci - r_ground_eci
n_ground = r_ground_eci / np.linalg.norm(r_ground_eci)
r_rel_norm = np.linalg.norm(r_rel)
sin_elevation = np.dot(r_rel, n_ground) / r_rel_norm
elevation = np.degrees(np.arcsin(sin_elevation))
elevation
():
gaps = []
current_gap =
covered coverage_mask:
covered:
current_gap += time_step
:
current_gap > :
gaps.append(current_gap)
current_gap =
current_gap > :
gaps.append(current_gap)
gaps
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd
from datetime import datetime, timedelta
class SatelliteOrbitSimulator:
"""卫星轨道模拟器"""
def __init__(self, constellation_config=None):
"""
初始化轨道模拟器
参数:
constellation_config: 星座配置字典
"""
self.constellation_config = constellation_config
self.constellation = None
self.integrator = NumericalIntegrator()
# 地球参数
self.R_e = 6378137.0 # 赤道半径 (m)
self.mu = 3.986004418e14 # 地球引力常数
# 时间系统
self.jd2000 = 2451545.0 # 2000 年 1 月 1 日 12:00 UT 的儒略日
def create_walker_constellation(self, T, P, F, altitude, inclination):
"""创建 Walker 星座"""
self.constellation = WalkerConstellation(T, P, F, altitude, inclination)
return self.constellation
def simulate_constellation(self, duration_hours=24, time_step=):
.constellation :
ValueError()
total_time = duration_hours *
n_steps = (total_time / time_step) +
simulation_data = {
: np.linspace(, total_time, n_steps),
: [],
: [],
: []
}
n_satellites = (.constellation.satellites)
i (n_satellites):
simulation_data[].append(np.zeros((n_steps, )))
simulation_data[].append(np.zeros((n_steps, )))
simulation_data[].append(np.zeros((n_steps, )))
step_idx, t (simulation_data[]):
sat_idx (n_satellites):
r_eci, v_eci = .constellation.get_satellite_position_velocity(sat_idx, t)
simulation_data[][sat_idx][step_idx] = r_eci
simulation_data[][sat_idx][step_idx] = v_eci
lat, lon = .eci_to_geodetic(r_eci, t)
simulation_data[][sat_idx][step_idx] = [lat, lon]
simulation_data
():
x, y, z = r_eci
omega_e =
theta = omega_e * t
x_rot = x * np.cos(theta) + y * np.sin(theta)
y_rot = -x * np.sin(theta) + y * np.cos(theta)
z_rot = z
lon = np.degrees(np.arctan2(y_rot, x_rot))
r = np.sqrt(x_rot** + y_rot** + z_rot**)
lat_geocentric = np.degrees(np.arcsin(z_rot / r))
lat = lat_geocentric
lat, lon
():
.constellation :
ValueError()
total_time = duration_hours *
time_step =
n_steps = (total_time / time_step)
coverage_data = {
: ground_stations,
: np.linspace(, total_time, n_steps),
: np.zeros(((ground_stations), n_steps)),
: np.zeros(((ground_stations), (.constellation.satellites), n_steps))
}
t_idx, t (coverage_data[]):
gs_idx, (lat, lon, _) (ground_stations):
best_elevation = -
sat_idx, sat (.constellation.satellites):
r_eci, _ = .constellation.get_satellite_position_velocity(sat_idx, t)
elevation = .constellation.calculate_elevation(r_eci, lat, lon, t)
coverage_data[][gs_idx, sat_idx, t_idx] = elevation
elevation > best_elevation:
best_elevation = elevation
best_elevation >= min_elevation:
coverage_data[][gs_idx, t_idx] =
coverage_data
class OrbitAnalyzer:
"""轨道分析工具"""
@staticmethod
def calculate_orbit_elements(r_eci, v_eci, mu=3.986004418e14):
"""从位置和速度计算轨道要素
参数:
r_eci: 位置向量 (m)
v_eci: 速度向量 (m/s)
mu: 引力常数
返回:
elements: 轨道要素字典
"""
# 位置和速度的模
r = np.linalg.norm(r_eci)
v = np.linalg.norm(v_eci)
# 角动量向量
h = np.cross(r_eci, v_eci)
h_norm = np.linalg.norm(h)
# 节点向量
n = np.cross([0, 0, 1], h)
n_norm = np.linalg.norm(n)
# 偏心率向量
e_vec = ((v**2 - mu/r) * r_eci - np.dot(r_eci, v_eci) * v_eci) / mu
e = np.linalg.norm(e_vec)
# 轨道能量
energy = v**2/2 - mu/r
# 半长轴
if e != 1:
a = -mu/(2*energy)
else:
a = float('inf')
# 轨道倾角
i = np.arccos(h[2]/h_norm)
# 升交点赤经
if n_norm != 0:
raan = np.arccos(n[0]/n_norm)
if n[1] < 0:
raan = 2*np.pi - raan
else:
raan = 0
# 近地点幅角
if n_norm != 0 and e > 0:
arg_peri = np.arccos(np.dot(n, e_vec)/(n_norm*e))
e_vec[] < :
arg_peri = *np.pi - arg_peri
:
arg_peri =
e > :
nu = np.arccos(np.dot(e_vec, r_eci)/(e*r))
np.dot(r_eci, v_eci) < :
nu = *np.pi - nu
:
nu =
e < :
E = * np.arctan(np.sqrt((-e)/(+e)) * np.tan(nu/))
M = E - e*np.sin(E)
:
M =
{
: a,
: e,
: np.degrees(i),
: np.degrees(raan),
: np.degrees(arg_peri),
: np.degrees(nu),
: M,
: *np.pi*np.sqrt(a**/mu) a> ()
}
():
n_stations = (coverage_data[])
n_time_steps = (coverage_data[])
analysis = {
: {},
: {}
}
gs_idx, (lat, lon, name) (coverage_data[]):
visibility = coverage_data[][gs_idx]
coverage_percentage = np.(visibility) / n_time_steps *
gaps = []
current_gap =
v visibility:
v == :
current_gap +=
:
current_gap > :
gaps.append(current_gap)
current_gap =
current_gap > :
gaps.append(current_gap)
avg_gap = np.mean(gaps) gaps
max_gap = np.(gaps) gaps
analysis[][name] = {
: lat,
: lon,
: coverage_percentage,
: avg_gap,
: max_gap,
: gaps
}
all_visibility = coverage_data[]
global_coverage = np.mean(all_visibility) *
simultaneous_coverage = np.(all_visibility, axis=)
analysis[] = {
: global_coverage,
: np.(simultaneous_coverage),
: np.(simultaneous_coverage),
: np.mean(simultaneous_coverage)
}
analysis
class OrbitVisualizer:
"""轨道可视化工具"""
def __init__(self):
self.fig = None
self.ax = None
def plot_3d_orbit(self, simulation_data, title="卫星轨道", show_earth=True):
"""绘制 3D 轨道图"""
self.fig = plt.figure(figsize=(12, 10))
self.ax = self.fig.add_subplot(111, projection='3d')
# 绘制地球
if show_earth:
self._plot_earth()
# 绘制每颗卫星的轨道
colors = plt.cm.rainbow(np.linspace(0, 1, len(simulation_data['satellite_positions'])))
for idx, positions in enumerate(simulation_data['satellite_positions']):
# 提取坐标
x = positions[:, 0] / 1000 # 转换为 km
y = positions[:, 1] / 1000
z = positions[:, 2] / 1000
# 绘制轨道
self.ax.plot(x, y, z, color=colors[idx], linewidth=1.0, alpha=0.7, label=f'Sat {idx+}')
.ax.set_xlabel()
.ax.set_ylabel()
.ax.set_zlabel()
.ax.set_title(title)
max_range = ([np.ptp(x) x [.ax.get_xlim(), .ax.get_ylim(), .ax.get_zlim()]] ) /
mid_x = (.ax.get_xlim()[] + .ax.get_xlim()[]) *
mid_y = (.ax.get_ylim()[] + .ax.get_ylim()[]) *
mid_z = (.ax.get_zlim()[] + .ax.get_zlim()[]) *
.ax.set_xlim(mid_x - max_range, mid_x + max_range)
.ax.set_ylim(mid_y - max_range, mid_y + max_range)
.ax.set_zlim(mid_z - max_range, mid_z + max_range)
plt.legend(loc=, bbox_to_anchor=(, ))
plt.tight_layout()
.fig, .ax
():
R_e =
u = np.linspace(, * np.pi, )
v = np.linspace(, np.pi, )
x = R_e * np.outer(np.cos(u), np.sin(v))
y = R_e * np.outer(np.sin(u), np.sin(v))
z = R_e * np.outer(np.ones(np.size(u)), np.cos(v))
.ax.plot_surface(x, y, z, color=, alpha=, edgecolor=)
theta = np.linspace(, *np.pi, )
x_eq = R_e * np.cos(theta)
y_eq = R_e * np.sin(theta)
.ax.plot(x_eq, y_eq, , , alpha=, linewidth=)
():
fig, axes = plt.subplots(, , figsize=(, ))
ax axes:
ax.set_xlim(-, )
ax.set_ylim(-, )
ax.set_xlabel()
ax.set_ylabel()
ax.grid(, alpha=)
ax.set_xticks(np.arange(-, , ))
ax.set_yticks(np.arange(-, , ))
colors = plt.cm.rainbow(np.linspace(, , (simulation_data[])))
idx, ground_track (simulation_data[]):
lons = ground_track[:, ]
lats = ground_track[:, ]
lons_wrapped = np.where(lons > , lons - , lons)
axes[].plot(lons_wrapped, lats, color=colors[idx], linewidth=, alpha=)
axes[].scatter(lons_wrapped[::], lats[::],
color=colors[idx], s=, alpha=)
axes[].set_title()
axes[].set_title()
plt.suptitle(title)
plt.tight_layout()
fig, axes
():
station_name:
name, data analysis_results[].items():
name == station_name:
station_data = data
:
ValueError()
fig, axes = plt.subplots(, , figsize=(, ))
ax = axes[, ]
gs_idx, (lat, lon, name_gs) (coverage_data[]):
name_gs == station_name:
visibility = coverage_data[][gs_idx]
time_hours = coverage_data[] /
ax.plot(time_hours, visibility, , linewidth=)
ax.fill_between(time_hours, , visibility, alpha=)
ax.set_xlabel()
ax.set_ylabel()
ax.set_title()
ax.set_ylim(-, )
ax.grid(, alpha=)
ax = axes[, ]
gs_idx, (lat, lon, name_gs) (coverage_data[]):
name_gs == station_name:
n_satellites = coverage_data[].shape[]
time_hours = coverage_data[] /
sat_idx (n_satellites):
elevation = coverage_data[][gs_idx, sat_idx, :]
ax.plot(time_hours, elevation, linewidth=, alpha=)
ax.set_xlabel()
ax.set_ylabel()
ax.set_title()
ax.grid(, alpha=)
ax = axes[, ]
station_data[]:
gaps_minutes = np.array(station_data[])
ax.hist(gaps_minutes, bins=, edgecolor=, alpha=)
ax.set_xlabel()
ax.set_ylabel()
ax.set_title()
ax.grid(, alpha=)
ax = axes[, ]
gs_idx, (lat, lon, name_gs) (coverage_data[]):
name_gs == station_name:
all_elevations = coverage_data[][gs_idx, :, :].flatten()
all_elevations = all_elevations[all_elevations > -]
ax.hist(all_elevations, bins=, edgecolor=, alpha=)
ax.set_xlabel()
ax.set_ylabel()
ax.set_title()
ax.grid(, alpha=)
plt.suptitle(, fontsize=)
plt.tight_layout()
:
fig, axes = plt.subplots(, , figsize=(, ))
ax = axes[, ]
station_names = []
coverage_percentages = []
name, data analysis_results[].items():
station_names.append(name)
coverage_percentages.append(data[])
y_pos = np.arange((station_names))
bars = ax.barh(y_pos, coverage_percentages, alpha=)
ax.set_yticks(y_pos)
ax.set_yticklabels(station_names)
ax.set_xlabel()
ax.set_title()
ax.grid(, alpha=, axis=)
ax = axes[, ]
time_hours = coverage_data[] /
all_visibility = coverage_data[]
simultaneous_coverage = np.(all_visibility, axis=)
ax.plot(time_hours, simultaneous_coverage, , linewidth=)
ax.fill_between(time_hours, , simultaneous_coverage, alpha=)
ax.set_xlabel()
ax.set_ylabel()
ax.set_title()
ax.grid(, alpha=)
ax = axes[, ]
metrics = analysis_results[]
metric_names = [, , , ]
metric_values = [
metrics[],
metrics[],
metrics[],
metrics[]
]
x_pos = np.arange((metric_names))
bars = ax.bar(x_pos, metric_values, alpha=)
ax.set_xticks(x_pos)
ax.set_xticklabels(metric_names, rotation=, ha=)
ax.set_ylabel()
ax.set_title()
ax.grid(, alpha=, axis=)
bar, value (bars, metric_values):
height = bar.get_height()
ax.text(bar.get_x() + bar.get_width()/, height, , ha=, va=)
ax = axes[, ]
lats = []
lons = []
coverage_values = []
name, data analysis_results[].items():
lats.append(data[])
lons.append(data[])
coverage_values.append(data[])
scatter = ax.scatter(lons, lats, c=coverage_values, cmap=, s=, edgecolor=)
ax.set_xlabel()
ax.set_ylabel()
ax.set_title()
ax.grid(, alpha=)
plt.colorbar(scatter, ax=ax, label=)
plt.suptitle(, fontsize=)
plt.tight_layout()
fig, axes
class OrbitAnimator:
"""轨道动画生成器"""
def __init__(self, simulation_data):
self.simulation_data = simulation_data
self.fig = None
self.ax = None
def create_3d_animation(self, output_file='orbit_animation.mp4', fps=30, dpi=100):
"""创建 3D 轨道动画"""
from matplotlib.animation import FuncAnimation
import matplotlib.pyplot as plt
# 设置图形
self.fig = plt.figure(figsize=(12, 10))
self.ax = self.fig.add_subplot(111, projection='3d')
# 绘制地球
self._plot_earth()
# 初始化卫星点
n_satellites = len(self.simulation_data['satellite_positions'])
satellites_points = []
satellites_lines = []
colors = plt.cm.rainbow(np.linspace(0, 1, n_satellites))
for i in range(n_satellites):
# 卫星当前位置点
point, = self.ax.plot([], [], [], 'o', color=colors[i], markersize=6, label=f'Sat {i+1}')
satellites_points.append(point)
line, = .ax.plot([], [], [], color=colors[i], linewidth=, alpha=)
satellites_lines.append(line)
.ax.set_xlabel()
.ax.set_ylabel()
.ax.set_zlabel()
.ax.set_title()
all_positions = np.vstack(.simulation_data[])
max_range = np.(np.(all_positions)) / *
.ax.set_xlim(-max_range, max_range)
.ax.set_ylim(-max_range, max_range)
.ax.set_zlim(-max_range, max_range)
n_frames = (.simulation_data[])
frame_step = (, n_frames // (fps * ))
():
idx = frame * frame_step
idx >= n_frames:
satellites_points + satellites_lines
i (n_satellites):
positions = .simulation_data[][i]
x = positions[idx, ] /
y = positions[idx, ] /
z = positions[idx, ] /
satellites_points[i].set_data([x], [y])
satellites_points[i].set_3d_properties([z])
start_idx = (, idx - )
x_line = positions[start_idx:idx+, ] /
y_line = positions[start_idx:idx+, ] /
z_line = positions[start_idx:idx+, ] /
satellites_lines[i].set_data(x_line, y_line)
satellites_lines[i].set_3d_properties(z_line)
time_seconds = .simulation_data[][idx]
hours = (time_seconds // )
minutes = ((time_seconds % ) // )
.ax.set_title()
satellites_points + satellites_lines
anim = FuncAnimation(.fig, update, frames=n_frames//frame_step, interval=/fps, blit=)
()
anim.save(output_file, writer=, fps=fps, dpi=dpi)
()
plt.close(.fig)
anim
():
R_e =
u = np.linspace(, * np.pi, )
v = np.linspace(, np.pi, )
x = R_e * np.outer(np.cos(u), np.sin(v))
y = R_e * np.outer(np.sin(u), np.sin(v))
z = R_e * np.outer(np.ones(np.size(u)), np.cos(v))
.ax.plot_surface(x, y, z, color=, alpha=, edgecolor=)
class VisibilityPredictor:
"""卫星可见性预测器"""
def __init__(self, constellation, min_elevation=10):
self.constellation = constellation
self.min_elevation = min_elevation
self.earth_radius = 6378137.0 # 地球半径 (m)
self.mu = 3.986004418e14 # 地球引力常数
def predict_visibility(self, ground_point, start_time, duration_hours=24, time_step=60):
"""预测卫星可见性
参数:
ground_point: 地面点 (lat, lon) (度)
start_time: 开始时间 (datetime 对象)
duration_hours: 预测时长 (小时)
time_step: 时间步长 (秒)
返回:
visibility_schedule: 可见性时间表
"""
lat, lon = ground_point
lat_rad = np.radians(lat)
lon_rad = np.radians(lon)
# 时间数组
total_time = duration_hours * 3600
time_array = np.arange(0, total_time, time_step)
# 初始化可见性数组
n_satellites = len(self.constellation.satellites)
visibility = np.zeros((n_satellites, len(time_array)))
elevations = np.zeros((n_satellites, len(time_array)))
# 计算每个时间步
for t_idx, t in enumerate(time_array):
# 地球自转角度
omega_e = 7.2921150e-5 # rad/s
theta = omega_e * t
# 地面点 ECI 坐标
x_ground = self.earth_radius * np.cos(lat_rad) * np.cos(lon_rad + theta)
y_ground = self.earth_radius * np.cos(lat_rad) * np.sin(lon_rad + theta)
z_ground = .earth_radius * np.sin(lat_rad)
ground_pos = np.array([x_ground, y_ground, z_ground])
sat_idx, sat (.constellation.satellites):
r_eci, _ = .constellation.get_satellite_position_velocity(sat_idx, t)
r_rel = r_eci - ground_pos
r_rel_norm = np.linalg.norm(r_rel)
ground_norm = ground_pos / np.linalg.norm(ground_pos)
sin_elevation = np.dot(r_rel, ground_norm) / r_rel_norm
elevation = np.degrees(np.arcsin(sin_elevation))
elevations[sat_idx, t_idx] = elevation
elevation >= .min_elevation:
visibility[sat_idx, t_idx] =
visibility_schedule = {
: time_array,
: [start_time + timedelta(seconds=(t)) t time_array],
: visibility,
: elevations,
: ground_point,
: .min_elevation
}
visibility_schedule
():
visibility = visibility_schedule[]
time_array = visibility_schedule[]
n_satellites, n_times = visibility.shape
analysis = {
: n_satellites,
: time_array[-] / ,
: time_array[] - time_array[] (time_array) > ,
: {},
: {}
}
sat_idx (n_satellites):
vis_times = visibility[sat_idx]
n_visible = np.(vis_times)
percentage = n_visible / n_times *
visible_periods = []
in_period =
start_idx =
t_idx, visible (vis_times):
visible in_period:
in_period =
start_idx = t_idx
visible in_period:
in_period =
visible_periods.append({
: time_array[start_idx],
: time_array[t_idx-],
: time_array[t_idx-] - time_array[start_idx]
})
in_period:
visible_periods.append({
: time_array[start_idx],
: time_array[-],
: time_array[-] - time_array[start_idx]
})
analysis[][sat_idx] = {
: (n_visible),
: percentage,
: visible_periods,
: np.mean([p[] p visible_periods]) visible_periods ,
: np.([p[] p visible_periods]) visible_periods
}
any_sat_visible = np.(visibility, axis=)
any_visible_percentage = np.(any_sat_visible) / n_times *
simultaneous = np.(visibility, axis=)
analysis[] = {
: any_visible_percentage,
: np.mean(simultaneous),
: np.(simultaneous),
: np.(simultaneous),
: ._find_coverage_gaps(any_sat_visible, time_array)
}
analysis
():
gaps = []
current_gap_start =
i, visible (any_sat_visible):
visible current_gap_start :
current_gap_start = time_array[i]
visible current_gap_start :
gap_duration = time_array[i] - current_gap_start
gaps.append({
: current_gap_start,
: time_array[i],
: gap_duration
})
current_gap_start =
current_gap_start :
gap_duration = time_array[-] - current_gap_start
gaps.append({
: current_gap_start,
: time_array[-],
: gap_duration
})
gaps
():
visibility = visibility_schedule[][sat_idx]
time_array = visibility_schedule[]
i (current_time_idx, (visibility)):
visibility[i]:
start_idx = i
j (start_idx, (visibility)):
visibility[j]:
end_idx = j -
:
end_idx = (visibility) -
elevations = visibility_schedule[][sat_idx]
max_elevation_idx = np.argmax(elevations[start_idx:end_idx+]) + start_idx
{
: sat_idx,
: time_array[start_idx],
: time_array[end_idx],
: time_array[end_idx] - time_array[start_idx],
: time_array[max_elevation_idx],
: elevations[max_elevation_idx],
: elevations[start_idx],
: elevations[end_idx]
}
class LinkBudgetAnalyzer:
"""链路预算分析器"""
def __init__(self, frequency=10e9): # 默认 10GHz
self.frequency = frequency # Hz
self.c = 299792458 # 光速 (m/s)
self.k = 1.380649e-23 # 玻尔兹曼常数
def calculate_free_space_loss(self, distance):
"""计算自由空间损耗
参数:
distance: 距离 (m)
返回:
fs_loss: 自由空间损耗 (dB)
"""
wavelength = self.c / self.frequency
fs_loss = 20 * np.log10(distance) + 20 * np.log10(self.frequency) + 20 * np.log10(4 * np.pi / self.c)
return fs_loss
def calculate_link_margin(self, tx_power, tx_gain, rx_gain, distance, system_temperature=500, data_rate=100e6, additional_losses=3.0):
"""计算链路余量
参数:
tx_power: 发射功率 (dBW)
tx_gain: 发射天线增益 (dBi)
rx_gain: 接收天线增益 (dBi)
distance: 距离 (m)
system_temperature: 系统噪声温度 (K)
data_rate: 数据速率 (bps)
additional_losses: 附加损耗 (dB)
返回:
link_margin: 链路余量 (dB)
"""
# 自由空间损耗
fs_loss = self.calculate_free_space_loss(distance)
# 接收功率
rx_power = tx_power + tx_gain + rx_gain - fs_loss - additional_losses
# 噪声功率
noise_power = 10 * np.log10(.k * system_temperature * data_rate)
cnr = rx_power - noise_power
required_ebn0 =
link_margin = cnr - required_ebn0 - * np.log10(data_rate)
{
: rx_power,
: noise_power,
: cnr,
: link_margin,
: fs_loss
}
():
results = {}
sat_idx, sat (constellation.satellites):
time_step =
total_time = *
time_array = np.arange(, total_time, time_step)
link_margins = []
distances = []
t time_array:
r_eci, _ = constellation.get_satellite_position_velocity(sat_idx, t)
lat, lon = ground_station
lat_rad = np.radians(lat)
lon_rad = np.radians(lon)
omega_e =
theta = omega_e * t
R_e =
x_ground = R_e * np.cos(lat_rad) * np.cos(lon_rad + theta)
y_ground = R_e * np.cos(lat_rad) * np.sin(lon_rad + theta)
z_ground = R_e * np.sin(lat_rad)
ground_pos = np.array([x_ground, y_ground, z_ground])
distance = np.linalg.norm(r_eci - ground_pos)
distances.append(distance)
link_result = .calculate_link_margin(
tx_power=satellite_params[],
tx_gain=satellite_params[],
rx_gain=ground_station_params[],
distance=distance,
system_temperature=ground_station_params[],
data_rate=satellite_params[]
)
link_margins.append(link_result[])
results[sat_idx] = {
: np.(link_margins),
: np.(link_margins),
: np.mean(link_margins),
: np.(distances),
: np.(distances),
: np.mean(distances),
: np.(np.array(link_margins) > ) / (link_margins) * ,
: {
: time_array,
: link_margins,
: distances
}
}
results
from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor
import multiprocessing as mp
class ParallelOrbitSimulator:
"""并行轨道模拟器"""
def __init__(self, n_workers=None):
self.n_workers = n_workers or mp.cpu_count()
def simulate_constellation_parallel(self, constellation, duration_hours=24, time_step=60):
"""并行模拟星座"""
total_time = duration_hours * 3600
n_steps = int(total_time / time_step) + 1
n_satellites = len(constellation.satellites)
# 初始化结果数组
positions = np.zeros((n_satellites, n_steps, 3))
velocities = np.zeros((n_satellites, n_steps, 3))
# 并行计算每颗卫星的轨道
with ProcessPoolExecutor(max_workers=self.n_workers) as executor:
# 准备参数
args = [(sat_idx, constellation, total_time, time_step) for sat_idx in range(n_satellites)]
# 并行执行
results = list(executor.map(self._simulate_single_satellite, args))
# 整理结果
for sat_idx, (pos, vel) in enumerate(results):
positions[sat_idx] = pos
velocities[sat_idx] = vel
return {
: np.linspace(, total_time, n_steps),
: positions,
: velocities
}
():
sat_idx, constellation, total_time, time_step = args
n_steps = (total_time / time_step) +
positions = np.zeros((n_steps, ))
velocities = np.zeros((n_steps, ))
step (n_steps):
t = step * time_step
r_eci, v_eci = constellation.get_satellite_position_velocity(sat_idx, t)
positions[step] = r_eci
velocities[step] = v_eci
positions, velocities
():
total_time = duration_hours *
time_step =
n_steps = (total_time / time_step)
n_points = (ground_points)
visibility = np.zeros((n_points, n_steps))
ProcessPoolExecutor(max_workers=.n_workers) executor:
args = [(point_idx, constellation, ground_points, total_time, time_step, min_elevation) point_idx (n_points)]
results = (executor.(._calculate_single_point_coverage, args))
point_idx, vis (results):
visibility[point_idx] = vis
{
: np.linspace(, total_time, n_steps),
: visibility
}
():
point_idx, constellation, ground_points, total_time, time_step, min_elevation = args
lat, lon = ground_points[point_idx]
n_steps = (total_time / time_step)
visibility = np.zeros(n_steps)
step (n_steps):
t = step * time_step
visible =
sat_idx ((constellation.satellites)):
r_eci, _ = constellation.get_satellite_position_velocity(sat_idx, t)
elevation = constellation.calculate_elevation(r_eci, lat, lon, t)
elevation >= min_elevation:
visible =
visibility[step] = visible
visibility
try:
import cupy as cp
class GPUOrbitSimulator:
"""GPU 加速轨道模拟器"""
def __init__(self):
self.device = cp.cuda.Device(0)
self.mu = cp.float64(3.986004418e14)
def propagate_orbits_gpu(self, initial_states, delta_t):
"""GPU 加速轨道传播
参数:
initial_states: 初始状态数组 (n_satellites, 6)
delta_t: 时间增量 (标量或数组)
返回:
final_states: 最终状态数组
"""
# 将数据转移到 GPU
states_gpu = cp.array(initial_states, dtype=cp.float64)
if cp.isscalar(delta_t):
delta_t_gpu = cp.float64(delta_t)
else:
delta_t_gpu = cp.array(delta_t, dtype=cp.float64)
# 定义 GPU 核函数
kernel = cp.ElementwiseKernel(
in_params='T x, T y, T z, T vx, T vy, T vz, T dt',
out_params='T fx, T fy, T fz, T fvx, T fvy, T fvz',
''' // 位置和速度
T r = sqrt(x*x + y*y + z*z); T r3 = r * r * r;
// 二体问题加速度
T ax = -mu * x / r3; T ay = -mu * y / r3; T az = -mu * z / r3;
// RK4 积分
// k1
T k1_vx = ax; T k1_vy = ay; T k1_vz = az; T k1_x = vx; T k1_y = vy; T k1_z = vz;
// k2
T x2 = x + 0.5 * dt * k1_x; T y2 = y + 0.5 * dt * k1_y; T z2 = z + 0.5 * dt * k1_z;
T vx2 = vx + 0.5 * dt * k1_vx; T vy2 = vy + 0.5 * dt * k1_vy; T vz2 = vz + 0.5 * dt * k1_vz;
T r2 = sqrt(x2*x2 + y2*y2 + z2*z2); T r3_2 = r2 * r2 * r2;
T k2_vx = -mu * x2 / r3_2; T k2_vy = -mu * y2 / r3_2; T k2_vz = -mu * z2 / r3_2;
T k2_x = vx2; T k2_y = vy2; T k2_z = vz2;
// k3
T x3 = x + 0.5 * dt * k2_x; T y3 = y + 0.5 * dt * k2_y; T z3 = z + 0.5 * dt * k2_z;
T vx3 = vx + 0.5 * dt * k2_vx; T vy3 = vy + 0.5 * dt * k2_vy; T vz3 = vz + 0.5 * dt * k2_vz;
T r3 = sqrt(x3*x3 + y3*y3 + z3*z3); T r3_3 = r3 * r3 * r3;
T k3_vx = -mu * x3 / r3_3; T k3_vy = -mu * y3 / r3_3; T k3_vz = -mu * z3 / r3_3;
T k3_x = vx3; T k3_y = vy3; T k3_z = vz3;
// k4
T x4 = x + dt * k3_x; T y4 = y + dt * k3_y; T z4 = z + dt * k3_z;
T vx4 = vx + dt * k3_vx; T vy4 = vy + dt * k3_vy; T vz4 = vz + dt * k3_vz;
T r4 = sqrt(x4*x4 + y4*y4 + z4*z4); T r3_4 = r4 * r4 * r4;
T k4_vx = -mu * x4 / r3_4; T k4_vy = -mu * y4 / r3_4; T k4_vz = -mu * z4 / r3_4;
T k4_x = vx4; T k4_y = vy4; T k4_z = vz4;
// 最终状态
fx = x + dt/6.0 * (k1_x + 2.0*k2_x + 2.0*k3_x + k4_x);
fy = y + dt/6.0 * (k1_y + 2.0*k2_y + 2.0*k3_y + k4_y);
fz = z + dt/6.0 * (k1_z + 2.0*k2_z + 2.0*k3_z + k4_z);
fvx = vx + dt/6.0 * (k1_vx + 2.0*k2_vx + 2.0*k3_vx + k4_vx);
fvy = vy + dt/6.0 * (k1_vy + 2.0*k2_vy + 2.0*k3_vy + k4_vy);
fvz = vz + dt/6.0 * (k1_vz + 2.0*k2_vz + 2.0*k3_vz + k4_vz);
''',
name=,
preamble=
)
x = states_gpu[:, ]
y = states_gpu[:, ]
z = states_gpu[:, ]
vx = states_gpu[:, ]
vy = states_gpu[:, ]
vz = states_gpu[:, ]
fx, fy, fz, fvx, fvy, fvz = kernel(x, y, z, vx, vy, vz, delta_t_gpu)
final_states = cp.stack([fx, fy, fz, fvx, fvy, fvz], axis=)
cp.asnumpy(final_states)
ImportError:
()
def simulate_starlink_constellation():
"""模拟 Starlink 星座"""
print("开始模拟 Starlink 星座...")
# 创建模拟器
simulator = SatelliteOrbitSimulator()
# Starlink Phase 1 参数（简化）
# 轨道壳 1: 550km, 53 度倾角，72 个轨道面，每面 22 颗卫星
constellation_1 = simulator.create_walker_constellation(
T=1584, # 72 * 22
P=72, F=1, altitude=550, inclination=53.0
)
# 轨道壳 2: 540km, 53.2 度倾角，72 个轨道面，每面 22 颗卫星
constellation_2 = simulator.create_walker_constellation(
T=1584, P=72, F=1, altitude=540, inclination=53.2
)
# 模拟运行（简化版，使用较小时长和步长）
print("模拟星座运行...")
simulation_data = simulator.simulate_constellation(
duration_hours=2, # 2 小时模拟
time_step=300 # 5 分钟步长
)
# 定义地面站
ground_stations = [
(40.7128, -74.0060, 'New York'),
(51.5074, -0.1278, 'London'),
(35.6762, 139.6503, 'Tokyo'),
(-33.8688, 151.2093, 'Sydney'),
(-23.5505, -46.6333, 'Sao Paulo'),
(28.6139, 77.2090, 'Delhi'),
(, , )
]
()
coverage_data = simulator.calculate_coverage(
ground_stations=ground_stations,
duration_hours=,
min_elevation=
)
analyzer = OrbitAnalyzer()
analysis_results = analyzer.analyze_constellation_coverage(coverage_data)
visualizer = OrbitVisualizer()
()
fig_3d, ax_3d = visualizer.plot_3d_orbit(
simulation_data, title=
)
plt.savefig(, dpi=, bbox_inches=)
fig_cov, axes_cov = visualizer.plot_coverage_analysis(coverage_data, analysis_results)
plt.savefig(, dpi=, bbox_inches=)
()
()
()
()
()
name, data analysis_results[].items():
()
()
global_metrics = analysis_results[]
()
()
plt.show()
{
: constellation_1,
: constellation_2,
: simulation_data,
: coverage_data,
: analysis_results
}
():
()
simulator = SatelliteOrbitSimulator()
constellation = simulator.create_walker_constellation(
T=, P=, F=, altitude=, inclination=
)
ground_station = (, -)
satellite_params = {
: ,
: ,
:
}
ground_station_params = {
: ,
:
}
link_analyzer = LinkBudgetAnalyzer(frequency=)
link_results = link_analyzer.analyze_constellation_link(
constellation, ground_station, satellite_params, ground_station_params
)
()
sat_idx, result link_results.items():
()
()
()
()
()
()
()
link_results
():
()
time
constellation_sizes = [, , , ]
simulation_times = []
size constellation_sizes:
()
simulator = SatelliteOrbitSimulator()
constellation = simulator.create_walker_constellation(
T=size, P=(np.sqrt(size)), F=, altitude=, inclination=
)
start_time = time.time()
simulation_data = simulator.simulate_constellation(
duration_hours=, time_step=
)
elapsed_time = time.time() - start_time
simulation_times.append(elapsed_time)
()
()
plt.figure(figsize=(, ))
plt.plot(constellation_sizes, simulation_times, , linewidth=, markersize=)
plt.xlabel()
plt.ylabel()
plt.title()
plt.grid(, alpha=)
plt.xscale()
plt.yscale()
x_log = np.log10(constellation_sizes)
y_log = np.log10(simulation_times)
coeffs = np.polyfit(x_log, y_log, )
trend_line = **coeffs[] * np.array(constellation_sizes)**coeffs[]
plt.plot(constellation_sizes, trend_line, , alpha=, label=)
plt.legend()
plt.savefig(, dpi=, bbox_inches=)
plt.show()
constellation_sizes, simulation_times
def main():
"""主函数：运行完整的星座模拟"""
print("=" * 60)
print("Python 卫星通信模拟：低轨星座的轨道力学计算")
print("=" * 60)
# 1. 模拟 Starlink 星座
starlink_results = simulate_starlink_constellation()
# 2. 链路预算分析
link_results = analyze_link_budget_example()
# 3. 性能基准测试
perf_results = performance_benchmark()
# 4. 创建动画（可选）
print("\n=== 生成轨道动画 ===")
animator = OrbitAnimator(starlink_results['simulation_data'])
try:
animator.create_3d_animation(
output_file='constellation_animation.mp4',
fps=30, dpi=100
)
print("动画生成完成!")
except Exception as e:
print(f"动画生成失败：{e}")
print("请确保已安装 ffmpeg")
print("\n" + "=" * 60)
print("模拟完成!")
print("=" * 60)
return {
'starlink': starlink_results,
'link_budget': link_results,
'performance': perf_results
}
if __name__ == "__main__":
# 运行主程序
results = main()
本文构建了一个完整的 Python 低轨卫星星座轨道力学计算框架，具有以下特点：
在开发过程中遇到的主要挑战及解决方案：
本框架可应用于：
完整代码可通过开源社区获取，安装和使用方法：
依赖库包括：
详细的开普勒轨道要素与笛卡尔坐标转换公式...
典型低轨卫星星座参数参考...
通过本文介绍的 Python 卫星轨道模拟框架，研究人员和工程师可以快速构建、分析和优化低轨卫星星座系统。该框架结合了理论深度和工程实用性，为卫星通信系统的设计和分析提供了强大的工具支持。随着商业航天和卫星互联网的快速发展，这样的模拟工具将变得越来越重要。
微信公众号「极客日志」，在微信中扫描左侧二维码关注。展示文案：极客日志 zeeklog
使用加密算法（如AES、TripleDES、Rabbit或RC4）加密和解密文本明文。 在线工具，加密/解密文本在线工具，online
解析常见 curl 参数并生成 fetch、axios、PHP curl 或 Python requests 示例代码。 在线工具，curl 转代码在线工具，online
将字符串编码和解码为其 Base64 格式表示形式即可。 在线工具，Base64 字符串编码/解码在线工具，online
将字符串、文件或图像转换为其 Base64 表示形式。 在线工具，Base64 文件转换器在线工具，online
将 Markdown（GFM）转为 HTML 片段，浏览器内 marked 解析；与 HTML转Markdown 互为补充。 在线工具，Markdown转HTML在线工具，online
将 HTML 片段转为 GitHub Flavored Markdown，支持标题、列表、链接、代码块与表格等；浏览器内处理，可链接预填。 在线工具，HTML转Markdown在线工具，online