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
// JEOD_INV: TS.01 — `<SelfRef>` / `<SelfPlanet>` are runtime-resolved storage-boundary wildcards; see `docs/JEOD_invariants.md` row TS.01 and the lint at `tests/self_ref_self_planet_discipline.rs`.
//! Stage 6 of [`super::super::Simulation::step_internal`]: drag, SRP,
//! and gravity-gradient-torque interactions. Reads atmospheric state
//! (stage 5), gravity acceleration (stage 4), and source positions
//! (stages 2-2b); writes per-body force/torque outputs and returns
//! `sun_pos` / `moon_pos` so stages 7 and 9 don't have to re-resolve
//! them.
use glam::{DMat3, DVec3};
use astrodyn::{IntegOrigin, Position, RadiationForce, RootInertial, Vec3Ext};
use super::super::Simulation;
impl Simulation {
/// Stage 6 — drag, SRP (flat-plate or cannonball), gravity gradient
/// torque. Returns the resolved Sun and Moon inertial positions in
/// root-inertial coordinates (computed once at the top so subsequent
/// stages can reuse them). `sun_pos` / `moon_pos` are typed
/// `Position<RootInertial>` so SRP / solar-beta / earth-lighting
/// callers cannot silently mix them with integration-frame body
/// state — RF.10.
pub(super) fn compute_interactions(
&mut self,
dt: f64,
body_integ_origins: &[IntegOrigin],
) -> (
Option<Position<RootInertial>>,
Option<Position<RootInertial>>,
) {
// sun_pos is also used in stage 9 (solar beta, earth lighting); compute once here.
let sun_pos = self.sun_source.map(|idx| self.source_position_typed(idx));
let moon_pos = self.moon_source.map(|idx| self.source_position_typed(idx));
let source_frame_ids = &self.source_frame_ids;
let frame_tree = &self.frame_tree;
let root_fid = self.root_frame_id;
for (body_idx, body) in self.bodies.iter_mut().enumerate() {
// JEOD_INV: RF.10 — SRP / shadow consume root-inertial position
// (because `sun_position` is root-inertial). Drag, by contrast,
// is *not* in the shift class: `compute_ballistic_drag` subtracts
// `atmos.wind` (computed via `omega × atmosphere_position` in
// the atmosphere planet's frame) from the vehicle velocity,
// so the velocity must match the same planet-centered frame —
// which is the body's integration frame in realistic configs.
//
// The body position is held typed (`Position<RootInertial>`) so
// every subtraction `body_pos - sun_position` typechecks only
// when both sides agree on `RootInertial`. Mixing
// `body.trans.position` (`Position<IntegrationFrame>`) with
// `sun_position` (`Position<RootInertial>`) is a compile error.
let body_pos_typed = body
.trans
.to_inertial(&body_integ_origins[body_idx])
.position;
let inertial_pos = body_pos_typed.raw_si();
// Compute structural transform once (shared by drag and flat-plate SRP)
let t_inertial_body = body.rot.as_ref().map_or(DMat3::IDENTITY, |r| {
r.q_inertial_body.as_witness().left_quat_to_transformation()
});
let t_inertial_struct =
astrodyn::compute_t_inertial_struct(&body.t_struct_body, &t_inertial_body);
// Aerodynamic drag — uses planet-centered velocity (NOT a
// shift site; see comment above).
body.aero_force = None;
if let Some(ref drag_config) = body.drag {
body.aero_force = Some(astrodyn::compute_ballistic_drag(
drag_config,
&body.atmospheric_state,
body.trans.velocity.raw_si(),
&t_inertial_struct,
));
}
// Solar radiation pressure (flat-plate)
body.radiation_force = None;
if let Some(ref mut fps) = body.flat_plate_state {
fps.stage_inputs = None;
}
if let Some(sun_position_typed) = sun_pos {
// Local raw alias for the legacy DVec3-arg consumers
// (`compute_shadow_fraction`, `compute_cannonball_srp`,
// `compute_flat_plate_srp_thermal`). The typed value
// remains in scope for the structural guards below.
let sun_position = sun_position_typed.raw_si();
if let Some(ref mut fps) = body.flat_plate_state {
// Flat-plate SRP with thermal emission. Typed
// subtraction enforces matching frames (RF.10):
// `body_pos_typed - sun_position_typed` only typechecks
// when both are `Position<RootInertial>`, structurally
// catching the integration-frame mixing bug.
let sun_to_vehicle: Position<RootInertial> =
body_pos_typed - sun_position_typed;
let distance = sun_to_vehicle.raw_si().length();
let sun_to_vehicle = sun_to_vehicle.raw_si();
// Skip SRP (not the whole body) if too close to Sun
if distance >= 1.0 {
let flux_inertial_hat = sun_to_vehicle / distance;
let flux_mag = astrodyn::solar_flux_at_distance(distance);
// Shadow fraction (step-constant — matches JEOD's
// scheduled-class shadow evaluation in SIM_3_ORBIT).
let illum_factor = body
.shadow_body
.map(|(idx, radius)| {
astrodyn::compute_shadow_fraction(
inertial_pos,
sun_position,
{
let fid = source_frame_ids[idx].inertial;
if fid == root_fid {
DVec3::ZERO
} else {
frame_tree.get(fid).state.trans.position
}
},
radius,
astrodyn::SOLAR_RADIUS,
)
})
.unwrap_or(1.0);
// `mass.position` is the CoM in the vehicle's
// structural frame (raw `DVec3`). Re-tag at
// the boundary into the typed structural-
// frame phantom for the SRP kernel inputs;
// the kernel re-extracts the raw DVec3 via
// `.raw_si()` for in-frame arithmetic.
let center_grav_struct = body
.mass
.as_ref()
.map_or(DVec3::ZERO, |m| m.center_of_mass.raw_si())
.m_at::<astrodyn::StructuralFrame<astrodyn::SelfRef>>();
let center_grav = center_grav_struct.raw_si();
// Kinematic children skip the standard
// integrator (`integrate.rs` continues past
// them on `kinematic_only`), and the
// per-RK4-stage derivative-class SRP/thermal
// recompute lives inside that gated branch
// (`integrate_body_coupled`). Falling through
// to the derivative arm below would stash
// `stage_inputs` that nobody consumes — plate
// temperatures freeze and `radiation_force`
// stays `None`. Force the Scheduled
// single-shot path so the appendage's thermal
// state advances and a representative
// `RadiationForce` is published, mirroring the
// analogous Bevy fix in PR #287
// (`flat_plate_srp_system` for
// `KinematicChildC` entities). Only the
// orbital state is gated by composite-rigid-
// body integration; Sun illumination, plate
// flux, and thermal evolution are
// state-independent of which body in the
// chain ends up at the root.
let order = if body.kinematic_only {
astrodyn::ThermalIntegrationOrder::Scheduled
} else {
fps.integration_order
};
match order {
astrodyn::ThermalIntegrationOrder::Scheduled => {
// Scheduled-class: compute SRP force + Euler T
// update once per step (JEOD SIM_3_ORBIT).
let flux_struct_hat = t_inertial_struct * flux_inertial_hat;
let srp_result = astrodyn::compute_flat_plate_srp_thermal(
&fps.plates,
&fps.t_pow4_cached,
flux_struct_hat,
flux_mag,
center_grav,
illum_factor,
);
// Force: structural → inertial. Torque: stays structural.
let force_inertial =
t_inertial_struct.transpose() * srp_result.force;
body.radiation_force = Some(RadiationForce {
force: force_inertial,
torque: srp_result.torque,
});
fps.integrate_temperatures(&srp_result.temp_dots, dt);
}
astrodyn::ThermalIntegrationOrder::DerivativeFirstOrder
| astrodyn::ThermalIntegrationOrder::DerivativeRk4 => {
// Derivative-class: SRP force (and optionally T)
// recomputed per RK4 stage. Stash the step-start
// inputs on the plate state; Stage 8 consumes
// them via `integrate_body_coupled` below.
fps.stage_inputs = Some(astrodyn::FlatPlateStageInputs {
// Typed `sun_position` carries the
// RootInertial frame phantom into the
// RK4 derivative closure, where it
// can only be subtracted from a
// typed `Position<RootInertial>` —
// structurally refusing the bug
// shape `stage_trans.position
// - sun_position` as DVec3-DVec3.
sun_position: sun_position_typed,
illum_factor,
center_grav: center_grav_struct,
});
// `radiation_force` is left None here; Stage 8
// writes a representative final-stage value so
// `VehicleOutput` still reports SRP force.
}
}
}
} else if let Some((cx_area, albedo, diffuse)) = body.cannonball_srp {
let illum_factor = body
.shadow_body
.map(|(idx, radius)| {
astrodyn::compute_shadow_fraction(
inertial_pos,
sun_position,
{
let fid = source_frame_ids[idx].inertial;
if fid == root_fid {
DVec3::ZERO
} else {
frame_tree.get(fid).state.trans.position
}
},
radius,
astrodyn::SOLAR_RADIUS,
)
})
.unwrap_or(1.0);
let force = astrodyn::compute_cannonball_srp(
inertial_pos,
sun_position,
cx_area,
albedo,
diffuse,
illum_factor,
);
if force != DVec3::ZERO {
body.radiation_force = Some(RadiationForce {
force,
torque: DVec3::ZERO,
});
}
}
}
// Gravity gradient torque
body.gravity_torque = None;
if body.compute_gravity_torque {
if let (Some(ref rot), Some(ref mass)) = (&body.rot, &body.mass) {
let t_parent_this = rot
.q_inertial_body
.as_witness()
.left_quat_to_transformation();
let inertia_dmat = mass.inertia.as_dmat3();
body.gravity_torque = Some(astrodyn::compute_gravity_torque(
&body.gravity_accel.grav_grad,
&t_parent_this,
&inertia_dmat,
));
}
}
}
(sun_pos, moon_pos)
}
}