Add SIM_verif_attach_mass RUN_09 (non-identity root orientation) (18→19/21) (#99)

crates/astrodyn_dynamics/src/mass_body.rs

@@ -962,34 +962,60 @@ impl MassTree {
         }
     }
 
-    /// Composite inertia — port of JEOD `mass_calc_composite_inertia.cc`.
+    /// Composite inertia — port of JEOD `mass_calc_composite_inertia.cc`
+    /// (with the body-frame offset transforms of `mass_update.cc:99-107`).
     ///
-    /// Starts with the core body's inertia shifted to the composite CoM via
-    /// the parallel axis theorem, then adds each child's composite inertia
-    /// (rotated to this body's structural frame) plus the child's parallel
-    /// axis contribution.
+    /// Computed in **this body's body frame**: JEOD's comment is "the core and
+    /// composite masses share a common body frame," and
+    /// `composite_properties.T_parent_this == core_properties.T_parent_this`
+    /// (`mass.cc:203`). The core inertia is already body-frame (the
+    /// `StructCG`/`Struct` init specs are rotated struct→body in
+    /// `mass_properties_init.cc:103/119`), so it enters unrotated; the
+    /// parallel-axis offsets are differences of struct-frame CoMs rotated into
+    /// the body frame by `T = t_parent_this` (`mass_update.cc:101/107`). Each
+    /// child's composite inertia, expressed in the *child's* body frame, is
+    /// rotated into this body frame by `composite_wrt_pbdy.T_parent_this`
+    /// (`mass_attach.cc:519`); for a direct child the body→body rotation is
+    /// `r = T · Sᵀ · T_childᵀ` (`S = structure_point.t_parent_this`,
+    /// parent-struct→child-struct; `T_child = child composite t_parent_this`),
+    /// and the rotated tensor is `r · I · rᵀ`.
+    ///
+    /// Reduces to the pure struct-frame computation (`Sᵀ · I · S`, struct
+    /// offsets) when every `t_parent_this` is identity — the common case and
+    /// every attach-mass RUN except RUN_09. For `yaw_180` + diagonal inertia
+    /// (Apollo) the body-frame result is bit-identical to the struct-frame one.
     fn calc_composite_inertia(&mut self, id: MassBodyId) {
         let cm = self.nodes[id].composite_properties.position;
-
-        // Core contribution: inertia + point-mass shift from core CoM to
-        // composite CoM (JEOD mass_calc_composite_inertia.cc lines 61-64).
+        // JEOD_INV: MA.25 — composite inertia in the body frame: body-frame core
+        // unrotated, offsets rotated struct→body by `t_parent_this`, children
+        // rotated child-body→parent-body.
+        // Struct→body for this composite body (shares the core's body frame).
+        let t_parent = self.nodes[id].composite_properties.t_parent_this;
+
+        // Core contribution: body-frame inertia (unrotated) plus the point-mass
+        // shift from the core CoM to the composite CoM, the offset rotated
+        // struct→body (JEOD mass_update.cc:107 + mass_calc_composite_inertia.cc:61-64).
         let core = &self.nodes[id].core_properties;
-        let core_offset = core.position - cm;
-        let mut composite_inertia = core.inertia + point_mass_inertia(core.mass, core_offset);
+        let core_offset_body = t_parent * (core.position - cm);
+        let mut composite_inertia = core.inertia + point_mass_inertia(core.mass, core_offset_body);
 
-        // Child contributions (lines 67-84).
+        // Child contributions (JEOD mass_calc_composite_inertia.cc:67-84).
         for &cid in &self.children[id] {
             let child = &self.nodes[cid];
-            let child_offset = child.composite_wrt_pstr.position - cm;
-
-            // Rotate child's composite inertia from child struct frame to
-            // parent struct frame: T^T * I_child * T
-            // This is JEOD's transpose_transform_matrix.
-            let t = child.structure_point.t_parent_this;
-            let rotated_inertia = t.transpose() * child.composite_properties.inertia * t;
-
-            composite_inertia +=
-                rotated_inertia + point_mass_inertia(child.composite_properties.mass, child_offset);
+            // Offset in this body's body frame (JEOD mass_update.cc:101).
+            let child_offset_body = t_parent * (child.composite_wrt_pstr.position - cm);
+
+            // Rotate the child's composite inertia from the child's body frame
+            // into this body's body frame: `r = child_body → parent_body
+            // = T · Sᵀ · T_childᵀ`, rotated tensor `r · I · rᵀ` (JEOD's
+            // `transpose_transform_matrix(composite_wrt_pbdy.T_parent_this, …)`).
+            let s = child.structure_point.t_parent_this;
+            let t_child = child.composite_properties.t_parent_this;
+            let r = t_parent * s.transpose() * t_child.transpose();
+            let rotated_inertia = r * child.composite_properties.inertia * r.transpose();
+
+            composite_inertia += rotated_inertia
+                + point_mass_inertia(child.composite_properties.mass, child_offset_body);
         }
 
         self.nodes[id].composite_properties.inertia = composite_inertia;

crates/astrodyn_runner/tests/runner_attach_detach_momentum.rs

@@ -1193,3 +1193,143 @@ fn from_builder_preserves_attached_bodies_initial_state() {
          got {parent_composite_mass}, expected {total_core_mass}"
     );
 }
+
+/// Inertia of a point mass `m` at offset `r` (parallel-axis term),
+/// computed independently of the kernel's `point_mass_inertia`.
+fn point_mass(m: f64, r: DVec3) -> DMat3 {
+    let outer = DMat3::from_cols(r * r.x, r * r.y, r * r.z);
+    DMat3::from_diagonal(DVec3::splat(r.length_squared())) * m - outer * m
+}
+
+/// Max per-column L2 distance between two matrices.
+fn mat3_max_col_diff(a: DMat3, b: DMat3) -> f64 {
+    let d = a - b;
+    [d.x_axis, d.y_axis, d.z_axis]
+        .into_iter()
+        .map(|c| c.length())
+        .fold(0.0_f64, f64::max)
+}
+
+/// A composite whose **parent has a non-identity, non-180° struct→body
+/// orientation** must carry its inertia in the body frame end-to-end:
+/// `recompute_composites` builds the composite in the parent body frame and
+/// `sync_body_mass_from_tree` hands it to the integrated body unchanged, so the
+/// rotational integrator (Euler's equation, body-frame `inertia · ω`) consumes
+/// a body-frame tensor.
+///
+/// `tier3_sim_attach_mass::RUN_09` cross-validates the body-frame composite
+/// *value* against JEOD's `mass.out`. This test guards the full Simulation
+/// `attach → sync_body_mass_from_tree → body.mass → integrate` pipeline for a
+/// **general** orientation — the case every existing trajectory scenario is
+/// blind to, because the only non-identity orientations in the suite are
+/// Apollo's `yaw_180`, which is inertia-invariant on its diagonal tensors.
+///
+/// The reference is derived independently via the frame-invariance identity:
+/// the same composite built in the **struct** frame (core rotated body→struct
+/// by `Sᵀ·I·S`, struct-frame parallel-axis offsets) and conjugated by `S`
+/// equals the body-frame composite. A sensitivity guard asserts the pipeline
+/// result is *not* the unrotated struct composite, so a regression that dropped
+/// the body-frame rotation fails loudly rather than silently.
+#[test]
+fn runner_attach_composite_inertia_is_body_frame() {
+    // General struct→body rotation (0.5 rad about a tilted axis): a yaw_180 or
+    // identity would hide the struct/body distinction this test exists to pin.
+    let s = DMat3::from_axis_angle(DVec3::new(1.0, 2.0, 3.0).normalize(), 0.5);
+
+    // Parent: asymmetric inertia in the BODY frame (a StructCG init would have
+    // already rotated it to body), with the non-identity struct→body transform.
+    let parent_body_inertia = DMat3::from_diagonal(DVec3::new(150.0, 200.0, 250.0));
+    let parent_mass = SimMassProperties::with_inertia(2.0, parent_body_inertia, DVec3::ZERO)
+        .with_t_parent_this(s);
+    // Child: identity orientation, body-frame inertia, attached at an off-axis
+    // struct offset so the composite is genuinely asymmetric (off-diagonal).
+    let child_body_inertia = DMat3::from_diagonal(DVec3::new(40.0, 50.0, 60.0));
+    let child_mass = SimMassProperties::with_inertia(1.0, child_body_inertia, DVec3::ZERO);
+    let offset = DVec3::new(2.0, 1.0, -0.5);
+
+    let omega0 = DVec3::new(0.05, 0.02, -0.01);
+    let parent_trans = TranslationalState {
+        position: DVec3::new(7e6, 0.0, 0.0),
+        velocity: DVec3::new(0.0, 7600.0, 0.0),
+    };
+    let parent_rot = Some(RotationalState {
+        quaternion: JeodQuat::identity(),
+        ang_vel_body: omega0,
+    });
+    let child_trans = TranslationalState {
+        position: DVec3::new(7e6, 0.0, 0.0) + offset,
+        velocity: DVec3::new(0.0, 7600.0, 0.0),
+    };
+    let child_rot = Some(RotationalState {
+        quaternion: JeodQuat::identity(),
+        ang_vel_body: omega0,
+    });
+
+    let (mut sim, parent_idx, child_idx, _pid, _cid) = build_pair(
+        0.5,
+        parent_mass,
+        parent_trans,
+        parent_rot,
+        child_mass,
+        child_trans,
+        child_rot,
+    );
+    sim.attach(child_idx, parent_idx, offset, DMat3::IDENTITY);
+
+    // Pipeline result: the integrated parent's composite inertia, body-frame.
+    let pipeline_body = sim
+        .body_mass(parent_idx)
+        .expect("parent mass after attach")
+        .inertia
+        .as_dmat3();
+
+    // Independent reference: build the composite in the struct frame and
+    // conjugate by S. Composite CoM (struct) = m_c/(m_p+m_c) along the offset.
+    let cm_struct = offset * (child_mass.mass / (parent_mass.mass + child_mass.mass));
+    let parent_core_struct = s.transpose() * parent_body_inertia * s;
+    let struct_composite = parent_core_struct
+        + point_mass(parent_mass.mass, -cm_struct)
+        + child_body_inertia
+        + point_mass(child_mass.mass, offset - cm_struct);
+    let expected_body = s * struct_composite * s.transpose();
+
+    assert!(
+        mat3_max_col_diff(pipeline_body, expected_body) < 1e-9,
+        "pipeline composite inertia must be the body-frame composite \
+         (S·I_struct·Sᵀ); got {pipeline_body:?}, expected {expected_body:?}"
+    );
+    // Sensitivity: for this general S the body-frame composite differs
+    // substantially from the unrotated struct composite — proves the
+    // struct→body rotation actually happened along the full pipeline.
+    assert!(
+        mat3_max_col_diff(pipeline_body, struct_composite) > 1.0,
+        "body-frame composite must differ from the struct-frame composite for a \
+         general orientation (else the body-frame rotation was dropped)"
+    );
+
+    // Smoke: the integrator propagates torque-free (mu = 0, no external torque)
+    // with the body-frame inertia. Body-frame angular momentum magnitude
+    // |I·ω| is conserved for a torque-free rigid body.
+    let omega_after_attach = sim
+        .body(parent_idx)
+        .rot
+        .expect("parent rot after attach")
+        .ang_vel_body
+        .raw_si();
+    let h0 = (pipeline_body * omega_after_attach).length();
+    for _ in 0..200 {
+        sim.step().expect("torque-free step");
+    }
+    let omega_final = sim
+        .body(parent_idx)
+        .rot
+        .expect("parent rot after stepping")
+        .ang_vel_body
+        .raw_si();
+    let hf = (pipeline_body * omega_final).length();
+    assert!(
+        omega_final.is_finite() && (hf - h0).abs() <= 1e-6 * h0.max(1.0),
+        "torque-free body-frame angular momentum |I·ω| must be conserved: \
+         |H0|={h0}, |Hf|={hf}"
+    );
+}

crates/astrodyn_verif_jeod/test_data/attach_mass_09_mass.out

@@ -0,0 +1,101 @@
+
+
+=============================================================
+
+Parent
+-------------------------------------------------------------
+
+Body Area
+Offset Vector [m]:
+             0.000000            0.000000            0.000000
+T_struct_struct [-]:
+             1.000000            0.000000            0.000000
+             0.000000            1.000000            0.000000
+             0.000000            0.000000            1.000000
+-------------------------------------------------------------
+Mass Properties
+M.P. CM vector [m]:
+             0.000000            0.000000            0.000000
+M.P. Mass [kg]:
+             1.000000
+M.P. Ib tensor [kgM2]:
+             0.354167           -0.025516            0.025516
+            -0.025516            0.239583           -0.156250
+             0.025516           -0.156250            0.239583
+M.P. T_struct_body [q]:
+             0.866025           -0.500000            0.000000
+             0.353553            0.612372            0.707107
+            -0.353553           -0.612372            0.707107
+-------------------------------------------------------------
+Composite Mass Properties
+C.M.P. CM vector [m]:
+            -0.500000            0.000000            0.000000
+C.M.P. Mass [kg]:
+             2.000000
+C.M.P. Ib tensor [kgM2]:
+             0.833333           -0.204124            0.204124
+            -0.204124            0.916667           -0.250000
+             0.204124           -0.250000            0.916667
+C.M.P. T_struct_body [q]:
+             0.866025           -0.500000            0.000000
+             0.353553            0.612372            0.707107
+            -0.353553           -0.612372            0.707107
+-------------------------------------------------------------
+Derived Items
+C.M.P. Inverse mass [1/kg]:
+             0.500000
+C.M.P. Inverse inertia tensor [1/(kgM2)]:
+             0.000000            0.000000            0.000000
+             0.000000            0.000000            0.000000
+             0.000000            0.000000            0.000000
+-------------------------------------------------------------
+
+=============================================================
+
+Child1
+-------------------------------------------------------------
+
+Body Area
+Offset Vector [m]:
+            -1.000000            0.000000            0.000000
+T_struct_struct [-]:
+             1.000000            0.000000            0.000000
+             0.000000            1.000000            0.000000
+             0.000000            0.000000            1.000000
+-------------------------------------------------------------
+Mass Properties
+M.P. CM vector [m]:
+             0.000000            0.000000            0.000000
+M.P. Mass [kg]:
+             1.000000
+M.P. Ib tensor [kgM2]:
+             0.333333            0.000000            0.000000
+             0.000000            0.416667            0.000000
+             0.000000            0.000000            0.083333
+M.P. T_struct_body [q]:
+             1.000000            0.000000            0.000000
+             0.000000            1.000000            0.000000
+             0.000000            0.000000            1.000000
+-------------------------------------------------------------
+Composite Mass Properties
+C.M.P. CM vector [m]:
+             0.000000            0.000000            0.000000
+C.M.P. Mass [kg]:
+             1.000000
+C.M.P. Ib tensor [kgM2]:
+             0.333333            0.000000            0.000000
+             0.000000            0.416667            0.000000
+             0.000000            0.000000            0.083333
+C.M.P. T_struct_body [q]:
+             1.000000            0.000000            0.000000
+             0.000000            1.000000            0.000000
+             0.000000            0.000000            1.000000
+-------------------------------------------------------------
+Derived Items
+C.M.P. Inverse mass [1/kg]:
+             1.000000
+C.M.P. Inverse inertia tensor [1/(kgM2)]:
+             0.000000            0.000000            0.000000
+             0.000000            0.000000            0.000000
+             0.000000            0.000000            0.000000
+-------------------------------------------------------------