fix: place coned rotor at its true CM before the parallel-axis shift

The parallel-axis theorem I_O = I_C + m(|r|^2 I - r outer r) only holds
when the tensor is taken about the body centre of mass. The rotor tensor
was formed about the hub apex, but a coned rotor's mass sits off the hub
plane along the shaft by span*sin(cone), so its CM is offset from the
apex. Using the apex both as the tensor reference and as the shift origin
double-counts that offset and drops the cross terms for a rotor with real
overhang and shaft tilt.

Compute the axial first moment (exact per-segment Simpson), take the
shaft-axis CM offset as first-moment / mass, reduce the axial second
moment to that CM (a mass-weighted variance, always non-negative), and
place the rotor body at apex + offset along the shaft axis. A flat rotor
(sin(cone) = 0) leaves the offset and axial term zero, so the planar
I_polar/2 split about the apex is unchanged.

For a coned rotor this raises the tower-top vertical CM by the precone
term real ElastoDyn also carries (~0.06 m on IEA-22); the deck adapter
approximates the blades as a point mass at the apex and omits it, so the
cross-check test now documents that ~1 % divergence.

Author: SMI-Lab-Inha

src/pybmodes/io/windio.py

@@ -93,6 +93,24 @@ def _second_moment(w: np.ndarray, c: np.ndarray, x: np.ndarray) -> float:
     return float(np.sum(dx / 6.0 * (f0 + 4.0 * fm + f1)))
 
 
+def _first_moment(w: np.ndarray, c: np.ndarray, x: np.ndarray) -> float:
+    """Exact ``integral w(x)*c(x) dx`` for piecewise-linear ``w`` and ``c``.
+
+    The integrand is quadratic on each segment (a product of two linear
+    factors), which Simpson's rule integrates exactly. Used for the rotor's
+    axial first moment, i.e. its shaft-axis CM offset, so a coned rotor body
+    is placed at its centre of mass before the parallel-axis shift.
+    """
+    w = np.asarray(w, dtype=float)
+    c = np.asarray(c, dtype=float)
+    x = np.asarray(x, dtype=float)
+    dx = np.diff(x)
+    w0, w1 = w[:-1], w[1:]
+    c0, c1 = c[:-1], c[1:]
+    fm = 0.5 * (w0 + w1) * 0.5 * (c0 + c1)
+    return float(np.sum(dx / 6.0 * (w0 * c0 + 4.0 * fm + w1 * c1)))
+
+
 def _require_yaml():
     """Import PyYAML or raise the documented friendly error.
 
@@ -812,7 +830,7 @@ def _blade_reference_axis(six: dict, comp: dict, blade_component: str) -> dict:
 def _blade_span_mass_inertia(
     six: dict, ref: dict, blade_component: str,
     hub_r: float, cone_cos: float, cone_sin: float,
-) -> tuple[float, float, float]:
+) -> tuple[float, float, float, float]:
     """Integrate a blade's span into ``(mass, polar, axial)`` second moments.
 
     ``six`` is
@@ -825,11 +843,14 @@ def _blade_span_mass_inertia(
     (prebend / sweep) are optional and default to zero (a straight span).
 
     Returns the single-blade mass, the polar second moment about the rotor
-    axis ``∫ (dm/ds) · r(s)² ds`` (radial ``r = hub_r + z·cone_cos``), and
-    the axial second moment ``∫ (dm/ds) · a(s)² ds`` (axial offset
-    ``a = z·cone_sin``). The caller builds the rotor's polar inertia from
-    the former and the coned transverse (diametral) inertia
-    ``polar/2 + N_bl·axial`` from both (issue #130).
+    axis ``∫ (dm/ds) · r(s)² ds`` (radial ``r = hub_r + z·cone_cos``), the
+    axial second moment ``∫ (dm/ds) · a(s)² ds`` (axial offset
+    ``a = z·cone_sin``), and the axial first moment ``∫ (dm/ds) · a(s) ds``.
+    The caller builds the rotor's polar inertia from the second radial
+    moment, the coned transverse (diametral) inertia from the radial and
+    axial second moments, and places the rotor body at its coned centre of
+    mass (shaft-axis offset ``axial_first / mass``) before the parallel-axis
+    shift (issue #130).
     """
     im = _require_mapping(six.get("inertia_matrix"), "blade six_x_six.inertia_matrix")
     grid = np.asarray(_require_key(im, "grid", "blade inertia_matrix"), dtype=float)
@@ -910,7 +931,8 @@ def _blade_span_mass_inertia(
     axial = span * cone_sin
     polar_second_moment = _second_moment(mpl, radial, s)
     axial_second_moment = _second_moment(mpl, axial, s)
-    return mass, polar_second_moment, axial_second_moment
+    axial_first_moment = _first_moment(mpl, axial, s)
+    return mass, polar_second_moment, axial_second_moment, axial_first_moment
 
 
 def read_windio_rna(
@@ -1071,6 +1093,7 @@ def _nac_geom(key: str) -> float:
     m_blade_each = 0.0
     i_polar_each = 0.0
     i_axial_each = 0.0
+    m_axial1_each = 0.0
     if n_blades > 0:
         blade = _require_mapping(
             comps.get(component_blade), f"components.{component_blade}"
@@ -1084,21 +1107,35 @@ def _nac_geom(key: str) -> float:
             f"components.{component_blade}.elastic_properties_mb.six_x_six",
         )
         ref = _blade_reference_axis(six, blade, component_blade)
-        m_blade_each, i_polar_each, i_axial_each = _blade_span_mass_inertia(
+        (
+            m_blade_each,
+            i_polar_each,
+            i_axial_each,
+            m_axial1_each,
+        ) = _blade_span_mass_inertia(
             six, ref, component_blade, hub_r, cone_cos, cone_sin,
         )
     m_blades = n_blades * m_blade_each
 
     # Rotor inertia from the spanwise blade mass (issue #130). In the shaft
-    # frame the rotor tensor about the hub is diag([I_polar, I_diam, I_diam]):
-    # the polar term ``I_polar = N_bl·∫dm·r²`` (radial), and the transverse
-    # (diametral) term ``I_diam = I_polar/2 + N_bl·∫dm·a²``. The ``I_polar/2``
-    # is the perpendicular-axis split for a planar rotor; the ``N_bl·∫dm·a²``
-    # adds the coned rotor's shaft-axis (``a = span·sin(cone)``) offset, which
-    # a flat-disc assumption would drop. Lumping the blades as a bare point
-    # mass at the apex would drop the whole tensor.
+    # frame the rotor tensor is diag([I_polar, I_diam, I_diam]): the polar
+    # term ``I_polar = N_bl·∫dm·r²`` (radial ``r = hub_r + span·cos(cone)``),
+    # and the transverse (diametral) term. For a coned rotor the blade mass
+    # sits off the hub plane by ``a = span·sin(cone)``, so the rotor CM shifts
+    # along the shaft axis by ``offset = ∫dm·a / m_blades`` and the tensor is
+    # formed about that CM: ``I_diam = I_polar/2 + N_bl·∫dm·a² − m_blades·
+    # offset²`` (the axial second moment reduced to the CM by the parallel-
+    # axis theorem, a mass-weighted variance that stays non-negative). The
+    # body is placed at that CM below so the shift to the tower top is exact.
+    # A flat rotor (``sin(cone) = 0``) leaves the offset and the axial term
+    # zero, recovering the planar ``I_polar/2`` split about the apex. Lumping
+    # the blades as a bare point mass at the apex would drop the whole tensor.
     i_polar_rotor = n_blades * i_polar_each
-    i_diam_rotor = 0.5 * i_polar_rotor + n_blades * i_axial_each
+    cm_axial_offset = (m_axial1_each / m_blade_each) if m_blade_each > 0.0 else 0.0
+    i_axial_cm = max(
+        0.0, n_blades * i_axial_each - m_blades * cm_axial_offset * cm_axial_offset
+    )
+    i_diam_rotor = 0.5 * i_polar_rotor + i_axial_cm
     i_blades = np.diag([i_polar_rotor, i_diam_rotor, i_diam_rotor])
 
     # Rotor apex relative to the tower top. An upwind rotor sits at negative
@@ -1124,10 +1161,17 @@ def _nac_geom(key: str) -> float:
     i_hub = r_tilt @ i_hub @ r_tilt.T
     i_blades = r_tilt @ i_blades @ r_tilt.T
 
+    # The coned rotor CM sits ``cm_axial_offset`` beyond the apex along the
+    # shaft axis (r_tilt's first column, the outboard unit vector the overhang
+    # already runs along), so a precone displaces the blade mass further from
+    # the tower in the same sense for both rotor orientations.
+    shaft_axis = r_tilt[:, 0]
+    rotor_cm = apex + cm_axial_offset * shaft_axis
+
     bodies = [
         (m_nac, r_nac, i_nac),
         (m_hub, apex, i_hub),
-        (m_blades, apex.copy(), i_blades),
+        (m_blades, rotor_cm, i_blades),
     ]
 
     m_total = float(sum(m for m, _, _ in bodies))

tests/test_windio_rna.py

@@ -200,6 +200,35 @@ def _mk(cone: float):
     assert coned.mass == pytest.approx(flat.mass)
 
 
+def test_rna_coned_cm_shift_with_uptilt(tmp_path: pathlib.Path) -> None:
+    """A coned rotor's CM sits off the hub apex along the (tilted) shaft, so
+    with uptilt the tower-top vertical CM shifts by m_blades·offset·sin(uptilt)
+    over the total mass. Pins that the rotor body is placed at its coned CM,
+    not the apex, before the parallel-axis shift (Codex review on #130)."""
+    pytest.importorskip("yaml")
+    import numpy as np
+
+    from pybmodes.io.windio import read_windio_rna
+
+    uptilt = 0.10
+
+    def _mk(cone: float):
+        o = _base_ontology()
+        o["components"]["hub"]["cone_angle"] = cone
+        o["components"]["nacelle"]["drivetrain"]["uptilt"] = uptilt
+        return read_windio_rna(_write(o, tmp_path, f"tc_{cone}.yaml"))
+
+    flat = _mk(0.0)
+    coned = _mk(0.15)
+    # base blade: uniform 100 kg/m over 0..50 m -> mass-weighted mean span 25 m,
+    # so the coned shaft-axis CM offset is sin(cone)·25.
+    offset = float(np.sin(0.15)) * 25.0
+    m_blades = 3 * 5000.0
+    expected = m_blades * offset * float(np.sin(uptilt)) / flat.mass
+    assert coned.mass == pytest.approx(flat.mass)
+    assert coned.cm_axial - flat.cm_axial == pytest.approx(expected, rel=1e-6)
+
+
 def test_rna_diagonal_inertia_accepted(tmp_path: pathlib.Path) -> None:
     """A diagonal 3-vector system_inertia [Ixx, Iyy, Izz] is accepted as a
     diagonal tensor, matching the equivalent zero-product 6-vector (Codex
@@ -558,9 +587,15 @@ def test_from_windio_with_monopile_lumped_rna_cal(tmp_path: pathlib.Path) -> Non
     reason="IEA-22 WindIO ontology + ElastoDyn deck not present",
 )
 def test_rna_iea22_matches_elastodyn_deck() -> None:
-    """The WindIO auto-RNA mass and CM match the IEA-22 ElastoDyn deck's
-    tower-top assembly to <0.5 % (issue #82). Inertia is WindIO-native and
-    is not expected to byte-match the deck's NacYIner (ecosystem drift)."""
+    """The WindIO auto-RNA mass matches the IEA-22 ElastoDyn deck's tower-top
+    assembly to <0.5 % (issue #82), and its CM agrees to ~1 %. The small CM
+    divergence is a documented, intentional one: the deck adapter lumps the
+    blades as a point mass at the rotor apex, whereas the WindIO rigid-rotor
+    path (issue #130) places the coned rotor at its true centre of mass,
+    which sits ~2.6 m outboard of the apex for the 4 deg precone and so
+    raises the tower-top vertical CM by ~0.06 m through the 6 deg shaft tilt
+    (real ElastoDyn carries the same precone term). Inertia is WindIO-native
+    and is not expected to byte-match the deck's NacYIner (ecosystem drift)."""
     pytest.importorskip("yaml")
     from pybmodes.io._elastodyn.adapter import _tower_top_assembly_mass
     from pybmodes.io.elastodyn_reader import (
@@ -576,7 +611,10 @@ def test_rna_iea22_matches_elastodyn_deck() -> None:
     )
 
     assert rna.mass == pytest.approx(deck.mass, rel=5e-3)       # <0.5 %
-    assert rna.cm_axial == pytest.approx(deck.cm_axial, rel=5e-3)
+    # The precone CM offset the deck adapter drops raises the WindIO CM
+    # slightly; it stays within ~1.5 % and lies above the apex-lumped value.
+    assert rna.cm_axial == pytest.approx(deck.cm_axial, rel=1.5e-2)
+    assert rna.cm_axial >= deck.cm_axial
     # WindIO-native inertia carries the fore-aft rotary term the same order
     # of magnitude as the deck, but not byte-identical (documented drift).
     assert rna.iyy > 0.0