快速入门

从零开始生成第一条周期轨道。

创建 CR3BP 系统

from e2m2e.algorithm.dynamics import CR3BP_System
from e2m2e.data.constants import Datum

# 创建地月 CR3BP 系统(μ 取 DE421 基准,ADR 0022)
system = CR3BP_System(
    mu=Datum.DE421.mu,
    primary="Earth",
    secondary="Moon",
)._with_default_scales()

# 计算平动点
system.compute_libration_points()
print(f"L1 = {system.L1}")
print(f"L2 = {system.L2}")

# 查看系统信息
system.info()

生成 DRO 轨道族

以下示例从种子轨道出发,通过微分修正 + 延拓生成 DRO 轨道族:

from e2m2e.algorithm.dynamics import CR3BP_System, CR3BP_Dynamics
from e2m2e.data.constants import Datum
from e2m2e.data.types.orbit import Orbit
from e2m2e.algorithm.solver import DifferentialCorrection, Continuation
import numpy as np

# 1. 创建系统与动力学(μ 取 DE421 基准,ADR 0022)
system = CR3BP_System(
    mu=Datum.DE421.mu, primary="Earth", secondary="Moon"
)._with_default_scales()
system.compute_libration_points()
dynamics = CR3BP_Dynamics(system)

# 2. 种子轨道(DRO 初始猜测)
initial_state = [0.79188556619742, 0.0, 0.0, 0.0, 0.53682, 0.0]
seed_orbit = Orbit(states=[initial_state], times=[0], system=system)
seed_orbit.period = 6.307  # 周期初猜(TU),修正中迭代更新

# 3. 微分修正:固定 x0 的 2D 对称策略
corrector = DifferentialCorrection(dynamics)
corrector.setup_2D_symmetric_x_fixed_x0(x0=initial_state[0])
result = corrector.iterate_correction(initial_guess=seed_orbit)
seed_dro = result.orbit  # 修正后的轨道(None 表示失败)

if seed_dro is not None:
    print(f"修正成功,周期 = {seed_dro.period:.6f}")

# 4. 延拓生成轨道族
continuation = Continuation(corrector=corrector)
cont_result = continuation.natural_continuation(
    seed_orbit=seed_dro,
    param_range=(0.14, 0.9),
    step_size=0.005,
)
print(f"轨道族包含 {len(cont_result.family.orbits)} 条轨道")

生成 Halo 轨道

from e2m2e.algorithm.family.halo_initial_guess import compute_halo_initial_guess

# Richardson 三阶解析近似生成初始猜测
z0 = 0.001  # z 方向振幅(小振幅种子,Richardson 近似精度高)
guess = compute_halo_initial_guess(system.mu, z0, L=1, halo_class=0)

initial_state = np.array([
    guess["x0"], 0.0, z0,
    0.0, guess["vy0"], 0.0,
])

# Halo 微分修正策略
corrector = DifferentialCorrection(dynamics)
corrector.setup_halo_orbit_fixed_z0(z0=z0, libration_point=1)

initial_guess = Orbit(
    states=initial_state.reshape(1, -1),
    times=np.array([0.0]),
    system=system,
)
initial_guess.period = guess["T_half"] * 2

halo_result = corrector.iterate_correction(initial_guess=initial_guess)
halo = halo_result.orbit
if halo is not None:
    print(f"Halo 周期: {halo.period:.6f}")

多重打靶

沿用上一节微分修正得到的 seed_dro 作为打靶初值,因此以下示例 以 seed_dro 修正成功(不为 None)为前提:

from e2m2e.algorithm.solver import MultipleShooting, sample_patch_points
from e2m2e.data.templates import ConvergenceState

ms = MultipleShooting(dynamics=dynamics)
t_patch, state_patch = sample_patch_points(seed_dro, n_points=5)

result = ms.correct(
    t_patch=t_patch,
    state_patch=state_patch,
    max_iter=50,
    tolerance=1e-10,
    var_time=True,
)

if result.status == ConvergenceState.CONVERGED:
    print(f"收敛,最大残差 {result.max_residual:.2e}")

转移轨道设计

from e2m2e.algorithm.transfer import Transfer

transfer = Transfer(dynamics)
result = transfer.set_orbit(start=dro_orbit, end=ro_orbit).optimize(
    initial_guess={"alpha": 1.0, "transfer_time": 15.0, "t_ins": 5.0},
    alpha_range=(0.5, 2.5),
)

力模型组合(星历传播)

from e2m2e.algorithm.coordinate import (
    CelestialBodyOrigin,
    CoordinateSystem,
    ICRSAxes,
)
from e2m2e.algorithm.dynamics import EphemerisSystem
from e2m2e.data.kernels.manager import SPICEManager
from e2m2e.algorithm.forces import ForceModel, GravityField, DragModel
from e2m2e.algorithm.forces.atmosphere import ExponentialAtmosphere

# 加载 SPICE 内核
spice = SPICEManager()
spice.load_kernel("kernels/de440s.bsp")

# 创建星历系统(frame 默认为 J2000)
eph_system = EphemerisSystem(
    bodies=["EARTH", "MOON", "SUN"],
    spice=spice, origin="EARTH",
)

# ForceModel 要求系统已设置坐标系
eph_system.coordinate_system = CoordinateSystem(
    axes=ICRSAxes(), origin=CelestialBodyOrigin(body="EARTH", spice=spice)
)

# 组合力模型
fm = ForceModel(eph_system)
fm.add_force(GravityField("EARTH", degree=2, order=0), name="j2")
fm.add_force(
    DragModel(atmosphere=ExponentialAtmosphere(), area=10.0, mass=1000.0),
    name="drag",
)

# 传播 LEO 轨道
r = 6378.137 + 400.0
v = np.sqrt(398600.4415 / r)
state0 = np.array([r, 0.0, 0.0, 0.0, v, 0.0])

et0 = spice.utc_to_et("2025-06-21T11:00:06")
result = fm.propagate(state0, (et0, et0 + 86400.0))

下一步