大气密度模型
e2m2e 提供可插拔的大气密度模型，供阻力力模型 DragModel 依赖注入使用。
当前实现基于 US Standard Atmosphere 1976 分段指数模型，覆盖 0-1000 km 高度范围。
核心概念
ExponentialAtmosphere：USSA76 分段指数大气密度模型，提供
density(altitude)接口，支持 F10.7 / Ap 一阶修正。
ExponentialAtmosphere
US Standard Atmosphere 1976 分段指数大气密度模型。
在每个高度层内使用 ρ(h) = ρ₀ · exp(-(h - h₀) / H) 计算密度，
层间标高由相邻断点密度比推导，确保密度连续且单调递减。
F10.7 太阳射电通量和 Ap 地磁指数通过线性乘法因子对基准密度做一阶修正。
from e2m2e.algorithm.forces.atmosphere import ExponentialAtmosphere
# 默认参数：F10.7=150 sfu（中等太阳活动），Ap=15（中等地磁活动）
atm = ExponentialAtmosphere()
# 查询不同高度的密度
rho_surface = atm.density(0.0) # 1.225 kg/m³
rho_100km = atm.density(100.0) # 5.604e-7 kg/m³
rho_400km = atm.density(400.0) # 2.803e-12 kg/m³
rho_1000km = atm.density(1000.0) # 0.0（超出模型上限）
# 调整太阳活动参数
atm_high = ExponentialAtmosphere(f107=200, ap=50)
rho_high = atm_high.density(400.0) # 密度高于默认值
模型范围与边界行为
高度范围：0 km 至 1000 km。
高于 1000 km：返回 0.0（阻力可忽略）。
低于 0 km：钳到 0 km，返回地表密度（避免负高度导致指数爆炸）。
参数说明
参数 |
说明 |
默认值 |
|---|---|---|
|
F10.7 太阳射电通量（sfu），反映太阳极紫外辐射强度 |
150.0 |
|
Ap 地磁指数，反映地磁活动强度 |
15.0 |
F10.7 和 Ap 越高，大气热膨胀越显著，同高度密度越大。 模型对两者的修正为线性乘法因子：
\[\text{factor} = \left(1 + 0.5 \cdot \frac{f107 - 150}{150}\right)
\cdot \left(1 + 0.1 \cdot \frac{ap - 15}{15}\right)\]
与 DragModel 配合：LEO 轨道衰减分析
以下示例展示 LEO 轨道衰减分析工作流：
选择指数大气模型、配置 DragModel、设置初始状态、运行轨道衰减传播。
import numpy as np
from e2m2e.data.kernels.manager import SPICEManager
from e2m2e.algorithm.dynamics.ephemeris_system import EphemerisSystem
from e2m2e.algorithm.coordinate.coordinate_system import CoordinateSystem
from e2m2e.algorithm.coordinate.standard_axes import ICRSAxes
from e2m2e.algorithm.coordinate.standard_origins import CelestialBodyOrigin
from e2m2e.algorithm.forces.atmosphere import ExponentialAtmosphere
from e2m2e.algorithm.forces import ForceModel, GravityField, DragModel
# 1. 准备星历系统（ICRF + 地球中心）
spice = SPICEManager()
spice.load_kernel("path/to/de440.bsp")
system = EphemerisSystem(bodies=["EARTH"], spice=spice, origin="EARTH")
system.coordinate_system = CoordinateSystem(
axes=ICRSAxes(),
origin=CelestialBodyOrigin(body="EARTH", spice=spice),
)
# 2. 选择大气模型并配置阻力
atmosphere = ExponentialAtmosphere(f107=150.0, ap=4.0)
drag = DragModel(
atmosphere=atmosphere,
body="EARTH",
cd=2.2, # 阻力系数
area=10.0, # 迎风面积 m²
mass=1000.0, # 航天器质量 kg
)
# 3. 构建力模型：J2 重力 + 大气阻力
fm = ForceModel(system)
fm.add_force(GravityField("EARTH", degree=2, order=0), name="j2")
fm.add_force(drag, name="drag")
# 4. 设置 LEO 初始状态（400 km 圆轨道，km / km/s）
r = 6378.137 + 400.0
v = np.sqrt(398600.4415 / r)
state0 = np.array([r, 0.0, 0.0, 0.0, v, 0.0])
# 5. 传播 1 天，观察轨道衰减
et0 = spice.utc_to_et("2025-06-21T11:00:06")
t_span = (et0, et0 + 86400.0)
t_eval = np.linspace(et0, et0 + 86400.0, 200)
result = fm.propagate(state0, t_span, t_eval=t_eval, max_steps=200_000)
# 6. 分析结果：半长轴随时间衰减
states = result["states"]
times = result["time"]
def semi_major_axis(state, mu=398600.4415):
r = np.linalg.norm(state[:3])
v = np.linalg.norm(state[3:6])
energy = v**2 / 2.0 - mu / r
return -mu / (2.0 * energy)
a_history = np.array([semi_major_axis(s) for s in states])
delta_a = a_history[-1] - a_history[0]
print(f"初始半长轴: {a_history[0]:.2f} km")
print(f"1 天后半长轴: {a_history[-1]:.2f} km")
print(f"半长轴衰减: {delta_a:.4f} km")
配置驱动方式（推荐）
上述工作流也可通过配置字典一键构建，便于序列化与复用：
config = {
"version": 1,
"forces": [
{
"name": "j2",
"type": "GravityField",
"enabled": True,
"params": {"body": "EARTH", "degree": 2, "order": 0},
},
{
"name": "drag",
"type": "DragModel",
"enabled": True,
"params": {
"body": "EARTH",
"cd": 2.2,
"area": 10.0,
"mass": 1000.0,
"atmosphere": {
"type": "ExponentialAtmosphere",
"params": {"f107": 150.0, "ap": 4.0},
},
},
},
],
}
fm = ForceModel.from_config(config, system)
# round-trip 契约：序列化后的配置可原样重建，且重建结果再序列化不变
fm2 = ForceModel.from_config(fm.to_config(), system)
assert fm2.to_config() == fm.to_config()
注意 to_config() 输出的是规范化配置，会补出手写配置省略的默认值键
（如 GravityField 的 input_frame、gravity_file），
因此不能直接断言 fm.to_config() == config。