astrodyn_runner 0.1.1

Standalone arena-state simulation harness driving the astrodyn pipeline without Bevy ECS
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
//! Supporting types for [`super::Simulation`].
//!
//! - Public surface declared here: [`VehicleOutput`],
//!   [`ContactPairConfig`] (re-exported through `simulation::mod` and
//!   `crate::lib` for API stability).
//! - Crate-internal: [`SimBody`], [`GravityData`].
//!
//! `SourceFrameIds` was lifted to `astrodyn::source_frames` (issue #71)
//! so the Bevy adapter can build source frames against the same
//! structure `astrodyn_runner` uses. `DetachedSubtreeState` lives in
//! `astrodyn::subtree` (issue #253 Task C) — pure rigid-body
//! kinematics, no `Simulation` dependency. It is re-exported from
//! `simulation::mod` so consumers reach it via
//! `astrodyn_runner::DetachedSubtreeState` regardless.
//!
//! JEOD_INV: TS.01 — runner storage boundary. The `<SelfRef>` /
//! `<SelfPlanet>` wildcards on `SimBody`, `VehicleConfig`, and
//! `VehicleOutput` fields below are the runner-internal analogue of the
//! Bevy `Component` boundary in `src/components.rs`: per-body storage
//! whose vehicle/planet identity is decided at runtime by the body's
//! position in the runner's `Vec<SimBody>` and its `gravity_controls`
//! source index, not by a compile-time mission-crate phantom. System
//! code (`step::*`, `mass_tree`, `frame_attach`) reads these fields and
//! threads `<V: Vehicle>` / `<P: Planet>` parameters through the
//! kernel calls that consume them. See
//! `tests/self_ref_self_planet_discipline.rs` for the workspace lint.

use glam::{DMat3, DVec3};

use astrodyn::{
    AerodynamicForce, AtmosphereState, DragConfig, DynamicsConfig, EulerSequence, FrameDerivatives,
    FrameSwitchConfig, GeodeticState, GravityAccelerationTyped, GravityControls, GravitySource,
    LvlhFrame, MassPropertiesTyped, OrbitalElements, RadiationForce, RootInertial, RotationModel,
    RotationalStateTyped, SelfPlanet, SelfRef, SrpModel, TotalForce, TranslationalStateTyped,
    VehicleConfig,
};
use astrodyn::{ContactFacet, FrameId, GroundFacet, IntegrationFrame, MassBodyId, MassPointState};

/// Registration of a contact interaction between two bodies.
///
/// Port of JEOD's `ContactPair` — the two facets are registered with the
/// `Contact` manager, which runs `check_contact()` at each derivative
/// evaluation (see `contact.sm`). In this runner, registered pairs are
/// evaluated at every RK4 stage inside [`super::Simulation::step`] when
/// any pairs are present.
#[derive(Debug, Clone)]
pub struct ContactPairConfig {
    /// Index of body A (the "subject" in JEOD terminology).
    pub body_a: usize,
    /// Facet on body A (shape positions in A's structural frame).
    pub facet_a: ContactFacet,
    /// Index of body B (the "target" in JEOD terminology).
    pub body_b: usize,
    /// Facet on body B (shape positions in B's structural frame).
    pub facet_b: ContactFacet,
}

/// Registration of a ground-contact interaction between a vehicle and a
/// planetary surface.
///
/// Symmetric to [`ContactPairConfig`] but bodyless on the ground side —
/// the ground doesn't integrate, has no Newton's-third-law reaction
/// applied to it, and is queried per-step from the [`GroundFacet`]'s
/// terrain model. The vehicle facet/material must match the ground
/// facet's material exactly (single-pair `SpringPairInteraction`
/// semantics).
#[derive(Debug, Clone)]
pub struct GroundContactPairConfig {
    /// Index of the vehicle body.
    pub body_a: usize,
    /// Vehicle facet (shape positions in the body's structural frame).
    pub vehicle_facet: ContactFacet,
    /// Ground facet (terrain, alt_offset, material).
    pub ground_facet: GroundFacet,
    /// JEOD initialization-time impulse, computed at registration via
    /// the `Phase::Initialization` evaluator (`facet_pos_body == 0`).
    /// Consumed at stage 1 of the first integration step (RK4 weight
    /// 1/6) and cleared to `None` thereafter — mirrors
    /// `ContactSurface::collect_forces_torques` zeroing `facet.force`
    /// after stage 1 in JEOD.
    ///
    /// Note that the init-phase evaluator essentially always reports
    /// contact for any realistic planet radius (because `|rel_state|`
    /// is O(1–2 m) — the facet's body-frame surface extent — while
    /// `|ground|` is O(R)). So this field is `Some(...)` for every
    /// successfully registered pair on initialization, regardless of
    /// the vehicle's actual altitude; the impulsive launch JEOD's CSV
    /// records is precisely this initialization-state effect. The
    /// field becomes `None` only after the first integration step
    /// consumes it.
    pub pending_initial_impulse: Option<GroundContactImpulse>,
}

/// Impulsive contact contribution from JEOD's pre-propagation
/// `GroundInteraction::initialize` call (force on the vehicle in
/// inertial coords + body-frame torque about CoM).
#[derive(Debug, Clone, Copy)]
pub struct GroundContactImpulse {
    pub force_inertial: DVec3,
    pub torque_body: DVec3,
}

/// Attachment of a body to a non-body reference frame (port of JEOD's
/// `DynBody::frame_attach` member, populated by `attach_to_frame`).
///
/// Captures the parent ref-frame ID and the rigid-body offset between
/// the parent frame and the attached body's composite-body frame at the
/// instant of attach. The runner's per-step pass derives the body's
/// state by composing the parent frame's current state with this fixed
/// offset (see `Simulation::propagate_frame_attached_state`); the body
/// stays glued to the parent frame as the frame moves under ephemeris,
/// planet rotation, or kinematic-joint drives.
///
/// JEOD source:
/// `models/dynamics/dyn_body/src/dyn_body_attach.cc:271-379` (the
/// three `attach_to_frame` overloads); the captured offset corresponds
/// to JEOD's `frame_attach.attach_offset` (`X_pframe_to_struct`).
#[derive(Debug, Clone, Copy)]
pub struct FrameAttachState {
    /// `FrameId` of the parent reference frame. Resolved from a
    /// frame-tree lookup at attach time and never reparented while the
    /// attachment holds.
    pub parent_frame_id: FrameId,
    /// Rigid-body offset from the parent frame to this body's
    /// composite-body frame, expressed in parent-frame coordinates.
    /// Frozen at attach time; changes only when the attachment is
    /// released and re-established.
    pub attach_offset: MassPointState,
}

/// Gravity-specific data associated with a source (decoupled from frame tree).
///
/// The frame tree stores position/velocity/rotation state; this struct stores
/// the physical gravity model data that lives alongside it. The `velocity`
/// field stores source velocity for relativistic corrections — for central
/// bodies at the root frame, the tree node has zero velocity but the source
/// may still have physical velocity (e.g., Sun orbiting the barycenter).
pub(crate) struct GravityData {
    /// Physical gravity source (mu, model: PointMass or SphericalHarmonics).
    pub source: GravitySource,
    /// Source velocity in the inertial frame (m/s). Used for relativistic
    /// corrections. Stored here rather than in the tree because the root
    /// frame's velocity must be zero (it's the reference origin).
    pub velocity: DVec3,
    /// Tidal ΔC20 to add to the base C20 coefficient. Updated each step.
    pub delta_c20: f64,
    /// Tidal configuration. When `Some`, the simulation computes ΔC20 each step.
    pub tidal_config: Option<astrodyn::TidalConfig>,
    /// Rotation model for updating planet-fixed frame each step.
    pub rotation_model: RotationModel,
    /// Sidereal angular velocity (rad/s) for the planet-fixed frame's
    /// `ang_vel_this`. Sourced from `PlanetConfig::omega` at setup time.
    /// JEOD sets this as `[0, 0, planet_omega]` in `planet_rnp.cc`.
    pub planet_omega: f64,
}

/// Read-only view of vehicle state after stepping.
///
/// Returned by [`super::Simulation::body`]. Contains the current integrated state
/// plus any derived states that were configured.
#[derive(Debug, Clone)]
pub struct VehicleOutput {
    /// Current translational state (position, velocity) in the integration frame.
    ///
    /// Carries the `IntegrationFrame` phantom (mirrors `SimBody.trans`'s
    /// integration-frame storage). Consumers that need root-inertial
    /// coordinates apply `body.trans.to_inertial(&integ_origin)` first;
    /// see `JEOD_invariants.md` row RF.10 for the shift discipline.
    pub trans: TranslationalStateTyped<IntegrationFrame>,
    /// Frame ID of the current integration frame in the simulation's frame tree.
    pub integ_frame_id: FrameId,
    /// Current rotational state (attitude, body-frame angular velocity).
    /// `None` for 3-DOF.
    pub rot: Option<RotationalStateTyped<SelfRef>>,
    /// Total translational acceleration in the integration frame (m/s²) at the
    /// end of the last `step()`. Sum of gravity and non-gravity contributions —
    /// mirrors JEOD's `derivs.trans_accel`. Zero before the first `step()`.
    pub trans_accel: DVec3,
    /// Total rotational acceleration in the body frame (rad/s²) at the end of
    /// the last `step()` — mirrors JEOD's `derivs.rot_accel`. `None` for 3-DOF
    /// bodies; zero before the first `step()`.
    pub rot_accel: Option<DVec3>,
    /// Orbital elements from the latest step.
    pub orbital_elements: Option<OrbitalElements<SelfPlanet>>,
    /// Euler angles `[phi, theta, psi]` from the latest step.
    pub euler_angles: Option<[f64; 3]>,
    /// LVLH frame from the latest step.
    pub lvlh_frame: Option<LvlhFrame>,
    /// Geodetic state (latitude, longitude, altitude).
    pub geodetic_state: Option<GeodeticState>,
    /// Solar beta angle (radians).
    pub solar_beta: Option<f64>,
    /// Earth lighting state (sun/moon occlusion, albedo).
    pub earth_lighting: Option<astrodyn::EarthLightingState>,
}

/// Internal per-body simulation state. Combines user config with bookkeeping
/// and output fields. Not exposed publicly — users interact through
/// [`VehicleConfig`] (input) and [`VehicleOutput`] (output).
pub(crate) struct SimBody {
    // ── Config (from VehicleConfig) ──
    /// Translational state in this body's integration frame.
    ///
    /// `IntegrationFrame` is *kind-distinct* from `RootInertial` so that
    /// root-inertial consumers (gravity, relativistic, SRP, solar beta,
    /// earth lighting — the *shift sites*, which mix body state with
    /// root-inertial source positions for Sun, Moon, or gravity sources)
    /// cannot silently take the integration-frame value — they must call
    /// `body.trans.to_inertial(&integ_origin)` first. Planet-inertial
    /// consumers (atmosphere, drag velocity, LVLH, geodetic, orbital
    /// elements — *non-shift sites*, which operate within a single
    /// planet's inertial frame) take `body.trans.position.raw_si()`
    /// directly: the body's integration frame is that planet's inertial
    /// frame in realistic configs, so applying the shift would change the
    /// planet-relative coordinates and produce wrong physics. See
    /// issue #255 and `JEOD_invariants.md` RF.10 for the split.
    pub trans: TranslationalStateTyped<IntegrationFrame>,
    pub rot: Option<RotationalStateTyped<SelfRef>>,
    pub mass: Option<MassPropertiesTyped<SelfRef>>,
    /// If this body participates in a mass tree, its node ID.
    pub mass_body_id: Option<MassBodyId>,
    /// When `true`, this body's `trans`/`rot` are derived each step from
    /// its mass-tree parent via `propagate_state_via_storage` rather
    /// than integrated. Mirrors the Bevy adapter's `KinematicChildC`
    /// marker. The body must be a non-root node in the mass tree;
    /// integration is skipped, force accumulation still runs (matching
    /// JEOD's `compute_point_derivative` flag, which lets a child body
    /// accumulate forces without integrating them — useful for
    /// post-step acceleration introspection).
    ///
    /// JEOD precedent: `DynBody::propagate_state_from_structure`
    /// (`models/dynamics/dyn_body/src/dyn_body_propagate_state.cc`)
    /// derives child states from the parent's structure each step;
    /// only the root integrates (`DB.17`).
    pub kinematic_only: bool,
    /// When `Some`, this body is attached to a non-body reference frame
    /// (port of JEOD `DynBody::attach_to_frame`,
    /// `models/dynamics/dyn_body/src/dyn_body_attach.cc:271-379`). The
    /// body's `trans` / `rot` are derived each step from the parent
    /// reference frame's state composed with the captured offset, and
    /// translational + rotational integration is suppressed (mirrors
    /// `dyn_body_integration.cc:309-333`'s `frame_attach.isAttached()`
    /// branch). Distinct from [`kinematic_only`](Self::kinematic_only),
    /// which targets a parent **body** in the mass tree;
    /// `frame_attach` targets a parent **reference frame**. JEOD_INV:
    /// DB.21 — only unattached bodies integrate.
    pub frame_attach: Option<FrameAttachState>,
    pub config: DynamicsConfig,
    pub gravity_controls: GravityControls<usize>,
    pub integrator: astrodyn::IntegratorType,
    pub drag: Option<DragConfig>,
    pub flat_plate_state: Option<astrodyn::FlatPlateState<astrodyn::SelfRef>>,
    pub cannonball_srp: Option<(f64, f64, f64)>,
    pub shadow_body: Option<(usize, f64)>,
    pub t_struct_body: DMat3,
    pub compute_gravity_torque: bool,
    /// Per-body atmospheric state slot.
    ///
    /// Always present; the runner's atmosphere stage gates execution on
    /// `body.drag.is_some()` (the same predicate that decides whether
    /// the slot was previously primed). This mirrors the Bevy adapter's
    /// `AtmosphericStateC<P>`, which `DragConfigC` auto-inserts via
    /// `#[require(...)]` — both consumers now carry the typed slot
    /// unconditionally and use a separate config-presence predicate to
    /// decide whether to fill it.
    pub atmospheric_state: AtmosphereState<SelfPlanet>,
    pub external_force: DVec3,
    pub external_torque: DVec3,
    /// Externally applied force in the body's structural frame (N).
    ///
    /// JEOD's `Force` is collected in the body's structural frame and
    /// rotated to inertial at force-collection time
    /// (`models/dynamics/dyn_body/src/dyn_body_collect.cc:219-221`).
    /// Tier 3 sims that schedule struct-frame force events
    /// (`SIM_verif_attach_detach`'s `RUN_compute_child_derivative`) need
    /// this entry point so the inertial-frame contribution tracks the
    /// body's current attitude across each integration step.
    pub external_force_struct: DVec3,
    /// Externally applied torque in the body's structural frame (N·m).
    ///
    /// Mirrors [`external_force_struct`](Self::external_force_struct);
    /// rotated to body frame at force-collection time via the body's
    /// structural-to-body transform.
    pub external_torque_struct: DVec3,

    // ── Frame switching ──
    pub integ_frame_id: FrameId,
    pub body_frame_id: FrameId,
    pub frame_switches: Vec<FrameSwitchConfig>,

    // ── Bookkeeping (written each step, not user-visible) ──
    pub gravity_accel: GravityAccelerationTyped<RootInertial>,
    pub total_force: TotalForce,
    pub frame_derivs: FrameDerivatives,
    pub aero_force: Option<AerodynamicForce>,
    pub radiation_force: Option<RadiationForce>,
    pub gravity_torque: Option<DVec3>,

    // ── Derived state config ──
    pub orbital_elements_source: Option<usize>,
    pub euler_sequence: Option<EulerSequence>,
    pub compute_lvlh: bool,
    pub geodetic_planet: Option<(usize, f64, f64)>,
    pub compute_solar_beta: bool,
    pub earth_lighting_config: Option<(f64, f64, f64)>,

    // ── Derived state outputs ──
    pub orbital_elements: Option<OrbitalElements<SelfPlanet>>,
    pub euler_angles: Option<[f64; 3]>,
    pub lvlh_frame: Option<LvlhFrame>,
    pub geodetic_state: Option<GeodeticState>,
    pub solar_beta: Option<f64>,
    pub earth_lighting: Option<astrodyn::EarthLightingState>,

    // ── Integrator state ──
    pub gj_state: Option<astrodyn::GaussJacksonState>,
    pub abm4_state: Option<astrodyn::Abm4State>,
}

impl SimBody {
    /// Convert a user-facing VehicleConfig into an internal SimBody.
    pub(crate) fn from_config(
        config: VehicleConfig,
        integ_frame_id: FrameId,
        body_frame_id: FrameId,
    ) -> Self {
        let has_rot = config.rot.is_some();
        let dynamics_config = DynamicsConfig {
            translational_dynamics: true,
            rotational_dynamics: has_rot,
            three_dof: !has_rot,
        };

        let (flat_plate_state, cannonball_srp) = match config.srp {
            Some(SrpModel::FlatPlate(fps)) => (Some(fps), None),
            Some(SrpModel::Cannonball {
                cx_area,
                albedo,
                diffuse,
            }) => (None, Some((cx_area, albedo, diffuse))),
            None => (None, None),
        };

        let shadow_body = config.shadow_body.map(|sb| (sb.source_idx, sb.radius));

        Self {
            // VehicleConfig::trans is typed `<RootInertial>`; the
            // runner re-tags as `<IntegrationFrame>` because integration
            // runs in the body's integration frame (numerically
            // coincides with root for root-integrated bodies). See #255.
            trans: config.trans.relabel_to::<IntegrationFrame>(),
            // VehicleConfig stores typed `<SelfRef>` rot/mass at the
            // mission boundary; SimBody now carries the same typed
            // siblings end-to-end (issue #397), so this is a direct
            // hand-off — no untyped demotion.
            rot: config.rot,
            mass: config.mass,
            mass_body_id: None,
            kinematic_only: false,
            frame_attach: None,
            config: dynamics_config,
            gravity_controls: config.gravity_controls,
            integrator: config.integrator,
            drag: config.drag,
            flat_plate_state,
            cannonball_srp,
            shadow_body,
            t_struct_body: config.t_struct_body,
            compute_gravity_torque: config.compute_gravity_gradient,
            atmospheric_state: AtmosphereState::<SelfPlanet>::default(),
            external_force: config.external_force.raw_si(),
            external_torque: config.external_torque.raw_si(),
            external_force_struct: DVec3::ZERO,
            external_torque_struct: DVec3::ZERO,

            integ_frame_id,
            body_frame_id,
            frame_switches: config.frame_switches,

            gravity_accel: GravityAccelerationTyped::<RootInertial>::default(),
            total_force: TotalForce::default(),
            frame_derivs: FrameDerivatives::default(),
            aero_force: None,
            radiation_force: None,
            gravity_torque: None,

            orbital_elements_source: config.derived.orbital_elements_source,
            euler_sequence: config.derived.euler_sequence,
            compute_lvlh: config.derived.lvlh,
            geodetic_planet: config
                .derived
                .geodetic
                .map(|g| (g.source_idx, g.r_eq, g.r_pol)),
            compute_solar_beta: config.derived.solar_beta,
            earth_lighting_config: config
                .derived
                .earth_lighting
                .map(|e| (e.earth_radius, e.moon_radius, e.sun_radius)),

            orbital_elements: None,
            euler_angles: None,
            lvlh_frame: None,
            geodetic_state: None,
            solar_beta: None,
            earth_lighting: None,

            gj_state: None,
            abm4_state: None,
        }
    }

    /// Create a VehicleOutput view of the current state.
    pub(crate) fn output(&self) -> VehicleOutput {
        VehicleOutput {
            // VehicleOutput::trans is typed against the integration
            // frame, matching `SimBody.trans`. Direct copy — no phantom
            // change, no demotion.
            trans: self.trans,
            integ_frame_id: self.integ_frame_id,
            rot: self.rot,
            trans_accel: self.frame_derivs.trans_accel,
            rot_accel: self.rot.map(|_| self.frame_derivs.rot_accel),
            orbital_elements: self.orbital_elements.clone(),
            euler_angles: self.euler_angles,
            lvlh_frame: self.lvlh_frame,
            geodetic_state: self.geodetic_state,
            solar_beta: self.solar_beta,
            earth_lighting: self.earth_lighting.clone(),
        }
    }
}