From ac56911d9999c5ff6b39b99f52b61780a5b4f42b Mon Sep 17 00:00:00 2001
From: Ryan Pavlick <218223+ryanpavlick@users.noreply.github.com>
Date: Sun, 3 May 2026 23:05:23 -0400
Subject: [PATCH] Add local and body frame transforms

---
 README.md                              |  63 ++++-
 src/pymap3d/__init__.py                | 108 +++++++-
 src/pymap3d/body.py                    | 331 +++++++++++++++++++++++++
 src/pymap3d/enu.py                     |  27 +-
 src/pymap3d/frames.py                  | 260 +++++++++++++++++++
 src/pymap3d/ned.py                     |  19 ++
 src/pymap3d/nwu.py                     | 141 +++++++++++
 src/pymap3d/sez.py                     | 141 +++++++++++
 src/pymap3d/tests/test_body.py         | 125 ++++++++++
 src/pymap3d/tests/test_local_frames.py |  94 +++++++
 10 files changed, 1306 insertions(+), 3 deletions(-)
 create mode 100644 src/pymap3d/body.py
 create mode 100644 src/pymap3d/frames.py
 create mode 100644 src/pymap3d/nwu.py
 create mode 100644 src/pymap3d/sez.py
 create mode 100644 src/pymap3d/tests/test_body.py
 create mode 100644 src/pymap3d/tests/test_local_frames.py

diff --git a/README.md b/README.md
index a4bcbb5..0520c40 100644
--- a/README.md
+++ b/README.md
@@ -88,11 +88,28 @@ converted to the desired coordinate system:
 
 ```
 aer2ecef  aer2enu  aer2geodetic  aer2ned
+         aer2nwu  aer2sez
+body2ecef body2ecefv body2enu body2geodetic body2ned
+         body2nwu body2sez
 ecef2aer  ecef2enu  ecef2enuv  ecef2geodetic  ecef2ned  ecef2nedv
+         ecef2body ecef2bodyv
+         ecef2nwu  ecef2nwuv  ecef2sez  ecef2sezv
 ecef2eci  eci2ecef eci2aer aer2eci geodetic2eci eci2geodetic
 enu2aer  enu2ecef   enu2geodetic
-geodetic2aer  geodetic2ecef  geodetic2enu  geodetic2ned
+         enu2body  enu2body_matrix
+         ecef2enu_matrix  enu2ecef_matrix
+geodetic2aer  geodetic2body  geodetic2ecef  geodetic2enu  geodetic2ned
+              geodetic2nwu  geodetic2sez
 ned2aer  ned2ecef   ned2geodetic
+         ned2body  ned2body_matrix
+nwu2body nwu2body_matrix
+sez2body sez2body_matrix
+         ecef2ned_matrix  ned2ecef_matrix
+nwu2aer  nwu2ecef   nwu2geodetic
+         ecef2nwu_matrix  nwu2ecef_matrix
+sez2aer  sez2ecef   sez2geodetic
+         ecef2sez_matrix  sez2ecef_matrix
+body_matrix body2enu_matrix body2ned_matrix body2nwu_matrix body2sez_matrix
 azel2radec radec2azel
 lookAtSpheroid
 track2 departure meanm
@@ -123,8 +140,52 @@ Abbreviations:
 * [ECI: Earth-centered Inertial using IERS](https://www.iers.org/IERS/EN/Home/home_node.html) via `astropy`
 * [ENU: East North Up](https://en.wikipedia.org/wiki/Axes_conventions#Ground_reference_frames:_ENU_and_NED)
 * [NED: North East Down](https://en.wikipedia.org/wiki/North_east_down)
+* NWU: North West Up
+* SEZ: South East Zenith
 * [radec: right ascension, declination](https://en.wikipedia.org/wiki/Right_ascension)
 
+### Local Frames
+
+PyMap3D includes several local tangent-plane frames:
+
+* `ENU` East, North, Up
+* `NED` North, East, Down
+* `NWU` North, West, Up
+* `SEZ` South, East, Zenith
+
+For vector-only workflows, explicit 3x3 rotation matrices are also available:
+
+```python
+import pymap3d as pm
+
+r_ecef_to_enu = pm.ecef2enu_matrix(lat0, lon0)
+r_enu_to_ecef = pm.enu2ecef_matrix(lat0, lon0)
+```
+
+### Body Frame
+
+PyMap3D also includes an aerospace body frame:
+
+* `x` forward
+* `y` right
+* `z` down
+
+The default body attitude uses yaw / pitch / roll with `sequence="zyx"`:
+
+* positive yaw turns right
+* positive pitch is nose up
+* positive roll is right wing down
+
+```python
+import pymap3d as pm
+
+xb, yb, zb = pm.ned2body(n, e, d, yaw, pitch, roll)
+n, e, d = pm.body2ned(xb, yb, zb, yaw, pitch, roll)
+```
+
+The composed body transforms can use any of the supported local parent frames:
+`frame="ned"`, `frame="enu"`, `frame="nwu"`, or `frame="sez"`.
+
 ### Ellipsoid
 
 Numerous functions in pymap3d use an ellipsoid model.
diff --git a/src/pymap3d/__init__.py b/src/pymap3d/__init__.py
index 3108384..8614f68 100644
--- a/src/pymap3d/__init__.py
+++ b/src/pymap3d/__init__.py
@@ -32,6 +32,31 @@
 __version__ = "3.2.0"
 
 from .aer import aer2ecef, aer2geodetic, ecef2aer, geodetic2aer
+from .body import (
+    body2ecef,
+    body2ecefv,
+    body2enu,
+    body2enu_matrix,
+    body2geodetic,
+    body2ned,
+    body2ned_matrix,
+    body2nwu,
+    body2nwu_matrix,
+    body2sez,
+    body2sez_matrix,
+    body_matrix,
+    ecef2body,
+    ecef2bodyv,
+    enu2body,
+    enu2body_matrix,
+    geodetic2body,
+    ned2body,
+    ned2body_matrix,
+    nwu2body,
+    nwu2body_matrix,
+    sez2body,
+    sez2body_matrix,
+)
 from .dca import (
     enu2dca,
     dca2enu,
@@ -57,16 +82,50 @@
     uvw2enu,
 )
 from .ellipsoid import Ellipsoid
-from .enu import aer2enu, enu2aer, enu2geodetic, geodetic2enu, enu2ecefv
+from .enu import (
+    aer2enu,
+    enu2aer,
+    enu2geodetic,
+    geodetic2enu,
+    enu2ecefv,
+    ecef2enu_matrix,
+    enu2ecef_matrix,
+)
 from .ned import (
     aer2ned,
     ecef2ned,
     ecef2nedv,
+    ecef2ned_matrix,
     geodetic2ned,
+    ned2ecef_matrix,
     ned2aer,
     ned2ecef,
     ned2geodetic,
 )
+from .nwu import (
+    aer2nwu,
+    ecef2nwu,
+    ecef2nwuv,
+    ecef2nwu_matrix,
+    geodetic2nwu,
+    nwu2aer,
+    nwu2ecef,
+    nwu2ecefv,
+    nwu2ecef_matrix,
+    nwu2geodetic,
+)
+from .sez import (
+    aer2sez,
+    ecef2sez,
+    ecef2sezv,
+    ecef2sez_matrix,
+    geodetic2sez,
+    sez2aer,
+    sez2ecef,
+    sez2ecefv,
+    sez2ecef_matrix,
+    sez2geodetic,
+)
 from .sidereal import datetime2sidereal, greenwichsrt
 from .spherical import geodetic2spherical, spherical2geodetic
 from .timeconv import str2dt
@@ -97,6 +156,29 @@
     "aer2geodetic",
     "ecef2aer",
     "geodetic2aer",
+    "body_matrix",
+    "ned2body_matrix",
+    "body2ned_matrix",
+    "enu2body_matrix",
+    "body2enu_matrix",
+    "ned2body",
+    "body2ned",
+    "enu2body",
+    "body2enu",
+    "sez2body",
+    "body2sez",
+    "nwu2body",
+    "body2nwu",
+    "sez2body_matrix",
+    "body2sez_matrix",
+    "nwu2body_matrix",
+    "body2nwu_matrix",
+    "ecef2body",
+    "body2ecef",
+    "ecef2bodyv",
+    "body2ecefv",
+    "geodetic2body",
+    "body2geodetic",
     "ecef2enu",
     "ecef2enuv",
     "ecef2geodetic",
@@ -112,13 +194,37 @@
     "enu2geodetic",
     "enu2ecefv",
     "geodetic2enu",
+    "ecef2enu_matrix",
+    "enu2ecef_matrix",
     "aer2ned",
     "ecef2ned",
     "ecef2nedv",
+    "ecef2ned_matrix",
     "geodetic2ned",
     "ned2aer",
     "ned2ecef",
+    "ned2ecef_matrix",
     "ned2geodetic",
+    "aer2sez",
+    "ecef2sez",
+    "ecef2sezv",
+    "ecef2sez_matrix",
+    "geodetic2sez",
+    "sez2aer",
+    "sez2ecef",
+    "sez2ecefv",
+    "sez2ecef_matrix",
+    "sez2geodetic",
+    "aer2nwu",
+    "ecef2nwu",
+    "ecef2nwuv",
+    "ecef2nwu_matrix",
+    "geodetic2nwu",
+    "nwu2aer",
+    "nwu2ecef",
+    "nwu2ecefv",
+    "nwu2ecef_matrix",
+    "nwu2geodetic",
     "datetime2sidereal",
     "greenwichsrt",
     "geodetic2spherical",
diff --git a/src/pymap3d/body.py b/src/pymap3d/body.py
new file mode 100644
index 0000000..d1057ff
--- /dev/null
+++ b/src/pymap3d/body.py
@@ -0,0 +1,331 @@
+"""Transforms involving an aerospace body frame."""
+
+from __future__ import annotations
+
+from .ecef import ecef2enu, enu2ecef
+from .ellipsoid import Ellipsoid
+from .enu import enu2geodetic, geodetic2enu
+from .frames import (
+    _body2ecef_rotation,
+    _body2enu_rotation,
+    _body2ned_rotation,
+    _body2nwu_rotation,
+    _body2sez_rotation,
+    _ecef2body_rotation,
+    _enu2body_rotation,
+    _matvec3,
+    _ned2body_rotation,
+    _nwu2body_rotation,
+    _sez2body_rotation,
+)
+from .ned import ecef2ned, ned2ecef, ned2geodetic, geodetic2ned
+from .nwu import ecef2nwu, geodetic2nwu, nwu2ecef, nwu2geodetic
+from .sez import ecef2sez, geodetic2sez, sez2ecef, sez2geodetic
+
+__all__ = [
+    "body_matrix",
+    "ned2body_matrix",
+    "body2ned_matrix",
+    "enu2body_matrix",
+    "body2enu_matrix",
+    "sez2body_matrix",
+    "body2sez_matrix",
+    "nwu2body_matrix",
+    "body2nwu_matrix",
+    "ned2body",
+    "body2ned",
+    "enu2body",
+    "body2enu",
+    "sez2body",
+    "body2sez",
+    "nwu2body",
+    "body2nwu",
+    "ecef2body",
+    "body2ecef",
+    "ecef2bodyv",
+    "body2ecefv",
+    "geodetic2body",
+    "body2geodetic",
+]
+
+
+def body_matrix(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """
+    Rotation matrix that maps NED coordinates into the body frame.
+
+    The body frame uses aerospace axes: x forward, y right, z down.
+    Positive yaw turns right, positive pitch is nose up, and positive
+    roll is right wing down.
+    """
+
+    return _ned2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+
+
+def ned2body_matrix(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps NED coordinates into the body frame."""
+
+    return _ned2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+
+
+def body2ned_matrix(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps body coordinates into NED."""
+
+    return _body2ned_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+
+
+def enu2body_matrix(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps ENU coordinates into the body frame."""
+
+    return _enu2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+
+
+def body2enu_matrix(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps body coordinates into ENU."""
+
+    return _body2enu_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+
+
+def sez2body_matrix(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps SEZ coordinates into the body frame."""
+
+    return _sez2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+
+
+def body2sez_matrix(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps body coordinates into SEZ."""
+
+    return _body2sez_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+
+
+def nwu2body_matrix(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps NWU coordinates into the body frame."""
+
+    return _nwu2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+
+
+def body2nwu_matrix(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps body coordinates into NWU."""
+
+    return _body2nwu_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+
+
+def ned2body(n, e, d, yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Convert NED coordinates into body coordinates."""
+
+    return _matvec3(ned2body_matrix(yaw, pitch, roll, sequence=sequence, deg=deg), (n, e, d))
+
+
+def body2ned(xb, yb, zb, yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Convert body coordinates into NED."""
+
+    return _matvec3(body2ned_matrix(yaw, pitch, roll, sequence=sequence, deg=deg), (xb, yb, zb))
+
+
+def enu2body(e, n, u, yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Convert ENU coordinates into body coordinates."""
+
+    return _matvec3(enu2body_matrix(yaw, pitch, roll, sequence=sequence, deg=deg), (e, n, u))
+
+
+def body2enu(xb, yb, zb, yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Convert body coordinates into ENU."""
+
+    return _matvec3(body2enu_matrix(yaw, pitch, roll, sequence=sequence, deg=deg), (xb, yb, zb))
+
+
+def sez2body(s, e, z, yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Convert SEZ coordinates into body coordinates."""
+
+    return _matvec3(sez2body_matrix(yaw, pitch, roll, sequence=sequence, deg=deg), (s, e, z))
+
+
+def body2sez(xb, yb, zb, yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Convert body coordinates into SEZ."""
+
+    return _matvec3(body2sez_matrix(yaw, pitch, roll, sequence=sequence, deg=deg), (xb, yb, zb))
+
+
+def nwu2body(n, w, u, yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Convert NWU coordinates into body coordinates."""
+
+    return _matvec3(nwu2body_matrix(yaw, pitch, roll, sequence=sequence, deg=deg), (n, w, u))
+
+
+def body2nwu(xb, yb, zb, yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Convert body coordinates into NWU."""
+
+    return _matvec3(body2nwu_matrix(yaw, pitch, roll, sequence=sequence, deg=deg), (xb, yb, zb))
+
+
+def ecef2body(
+    x,
+    y,
+    z,
+    lat0,
+    lon0,
+    h0,
+    yaw,
+    pitch,
+    roll,
+    ell: Ellipsoid | None = None,
+    frame: str = "ned",
+    sequence: str = "zyx",
+    deg: bool = True,
+):
+    """Convert a target ECEF point into body coordinates."""
+
+    if frame == "ned":
+        local = ecef2ned(x, y, z, lat0, lon0, h0, ell, deg=deg)
+        return ned2body(*local, yaw, pitch, roll, sequence=sequence, deg=deg)
+    if frame == "enu":
+        local = ecef2enu(x, y, z, lat0, lon0, h0, ell, deg=deg)
+        return enu2body(*local, yaw, pitch, roll, sequence=sequence, deg=deg)
+    if frame == "sez":
+        local = ecef2sez(x, y, z, lat0, lon0, h0, ell, deg=deg)
+        return sez2body(*local, yaw, pitch, roll, sequence=sequence, deg=deg)
+    if frame == "nwu":
+        local = ecef2nwu(x, y, z, lat0, lon0, h0, ell, deg=deg)
+        return nwu2body(*local, yaw, pitch, roll, sequence=sequence, deg=deg)
+
+    raise ValueError("frame must be one of 'ned', 'enu', 'sez', or 'nwu'")
+
+
+def body2ecef(
+    xb,
+    yb,
+    zb,
+    lat0,
+    lon0,
+    h0,
+    yaw,
+    pitch,
+    roll,
+    ell: Ellipsoid | None = None,
+    frame: str = "ned",
+    sequence: str = "zyx",
+    deg: bool = True,
+):
+    """Convert body coordinates into a target ECEF point."""
+
+    if frame == "ned":
+        local = body2ned(xb, yb, zb, yaw, pitch, roll, sequence=sequence, deg=deg)
+        return ned2ecef(*local, lat0, lon0, h0, ell, deg=deg)
+    if frame == "enu":
+        local = body2enu(xb, yb, zb, yaw, pitch, roll, sequence=sequence, deg=deg)
+        return enu2ecef(*local, lat0, lon0, h0, ell, deg=deg)
+    if frame == "sez":
+        local = body2sez(xb, yb, zb, yaw, pitch, roll, sequence=sequence, deg=deg)
+        return sez2ecef(*local, lat0, lon0, h0, ell, deg=deg)
+    if frame == "nwu":
+        local = body2nwu(xb, yb, zb, yaw, pitch, roll, sequence=sequence, deg=deg)
+        return nwu2ecef(*local, lat0, lon0, h0, ell, deg=deg)
+
+    raise ValueError("frame must be one of 'ned', 'enu', 'sez', or 'nwu'")
+
+
+def ecef2bodyv(
+    x,
+    y,
+    z,
+    lat0,
+    lon0,
+    yaw,
+    pitch,
+    roll,
+    frame: str = "ned",
+    sequence: str = "zyx",
+    deg: bool = True,
+):
+    """Convert an ECEF vector into body coordinates."""
+
+    rotation = _ecef2body_rotation(
+        lat0, lon0, yaw, pitch, roll, sequence=sequence, frame=frame, deg=deg
+    )
+    return _matvec3(rotation, (x, y, z))
+
+
+def body2ecefv(
+    xb,
+    yb,
+    zb,
+    lat0,
+    lon0,
+    yaw,
+    pitch,
+    roll,
+    frame: str = "ned",
+    sequence: str = "zyx",
+    deg: bool = True,
+):
+    """Convert a body-frame vector into ECEF."""
+
+    rotation = _body2ecef_rotation(
+        lat0, lon0, yaw, pitch, roll, sequence=sequence, frame=frame, deg=deg
+    )
+    return _matvec3(rotation, (xb, yb, zb))
+
+
+def geodetic2body(
+    lat,
+    lon,
+    h,
+    lat0,
+    lon0,
+    h0,
+    yaw,
+    pitch,
+    roll,
+    ell: Ellipsoid | None = None,
+    frame: str = "ned",
+    sequence: str = "zyx",
+    deg: bool = True,
+):
+    """Convert a target geodetic point into body coordinates."""
+
+    if frame == "ned":
+        local = geodetic2ned(lat, lon, h, lat0, lon0, h0, ell, deg=deg)
+        return ned2body(*local, yaw, pitch, roll, sequence=sequence, deg=deg)
+    if frame == "enu":
+        local = geodetic2enu(lat, lon, h, lat0, lon0, h0, ell, deg=deg)
+        return enu2body(*local, yaw, pitch, roll, sequence=sequence, deg=deg)
+    if frame == "sez":
+        local = geodetic2sez(lat, lon, h, lat0, lon0, h0, ell, deg=deg)
+        return sez2body(*local, yaw, pitch, roll, sequence=sequence, deg=deg)
+    if frame == "nwu":
+        local = geodetic2nwu(lat, lon, h, lat0, lon0, h0, ell, deg=deg)
+        return nwu2body(*local, yaw, pitch, roll, sequence=sequence, deg=deg)
+
+    raise ValueError("frame must be one of 'ned', 'enu', 'sez', or 'nwu'")
+
+
+def body2geodetic(
+    xb,
+    yb,
+    zb,
+    lat0,
+    lon0,
+    h0,
+    yaw,
+    pitch,
+    roll,
+    ell: Ellipsoid | None = None,
+    frame: str = "ned",
+    sequence: str = "zyx",
+    deg: bool = True,
+):
+    """Convert body coordinates into a target geodetic point."""
+
+    if frame == "ned":
+        local = body2ned(xb, yb, zb, yaw, pitch, roll, sequence=sequence, deg=deg)
+        return ned2geodetic(*local, lat0, lon0, h0, ell, deg=deg)
+    if frame == "enu":
+        local = body2enu(xb, yb, zb, yaw, pitch, roll, sequence=sequence, deg=deg)
+        return enu2geodetic(*local, lat0, lon0, h0, ell, deg=deg)
+    if frame == "sez":
+        local = body2sez(xb, yb, zb, yaw, pitch, roll, sequence=sequence, deg=deg)
+        return sez2geodetic(*local, lat0, lon0, h0, ell, deg=deg)
+    if frame == "nwu":
+        local = body2nwu(xb, yb, zb, yaw, pitch, roll, sequence=sequence, deg=deg)
+        return nwu2geodetic(*local, lat0, lon0, h0, ell, deg=deg)
+
+    raise ValueError("frame must be one of 'ned', 'enu', 'sez', or 'nwu'")
diff --git a/src/pymap3d/enu.py b/src/pymap3d/enu.py
index 9fa0bbb..2cd1265 100644
--- a/src/pymap3d/enu.py
+++ b/src/pymap3d/enu.py
@@ -11,9 +11,18 @@
 
 from .ecef import ecef2geodetic, enu2ecef, geodetic2ecef, uvw2enu
 from .ellipsoid import Ellipsoid
+from .frames import _ecef2enu_rotation, _enu2ecef_rotation
 from .mathfun import atan2, cos, degrees, hypot, radians, sin
 
-__all__ = ["enu2aer", "aer2enu", "enu2geodetic", "geodetic2enu", "enu2ecefv"]
+__all__ = [
+    "enu2aer",
+    "aer2enu",
+    "enu2geodetic",
+    "geodetic2enu",
+    "enu2ecefv",
+    "ecef2enu_matrix",
+    "enu2ecef_matrix",
+]
 
 
 def enu2aer(e, n, u, deg: bool = True) -> tuple:
@@ -241,3 +250,19 @@ def enu2ecefv(e, n, u, lat0, lon0, deg: bool = True) -> tuple:
     z = cos(lat0) * n + sin(lat0) * u
 
     return x, y, z
+
+
+def ecef2enu_matrix(lat0, lon0, deg: bool = True):
+    """
+    Rotation matrix that maps ECEF vectors into ENU coordinates.
+    """
+
+    return _ecef2enu_rotation(lat0, lon0, deg=deg)
+
+
+def enu2ecef_matrix(lat0, lon0, deg: bool = True):
+    """
+    Rotation matrix that maps ENU vectors into ECEF coordinates.
+    """
+
+    return _enu2ecef_rotation(lat0, lon0, deg=deg)
diff --git a/src/pymap3d/frames.py b/src/pymap3d/frames.py
new file mode 100644
index 0000000..eb14a07
--- /dev/null
+++ b/src/pymap3d/frames.py
@@ -0,0 +1,260 @@
+"""Shared local-frame rotation helpers."""
+
+from __future__ import annotations
+
+from .mathfun import cos, radians, sin
+
+
+def _matmul3(a, b):
+    """Multiply two 3x3 matrices represented as tuples."""
+
+    return tuple(
+        tuple(sum(a[i][k] * b[k][j] for k in range(3)) for j in range(3)) for i in range(3)
+    )
+
+
+def _transpose3(matrix):
+    """Transpose a 3x3 matrix represented as tuples."""
+
+    return tuple(tuple(matrix[j][i] for j in range(3)) for i in range(3))
+
+
+def _matvec3(matrix, vector):
+    """Apply a 3x3 matrix to a 3-vector."""
+
+    return tuple(sum(matrix[i][j] * vector[j] for j in range(3)) for i in range(3))
+
+
+def _rotation_x(angle):
+    """Active rotation matrix about x."""
+
+    cos_angle = cos(angle)
+    sin_angle = sin(angle)
+    return (
+        (1.0, 0.0, 0.0),
+        (0.0, cos_angle, -sin_angle),
+        (0.0, sin_angle, cos_angle),
+    )
+
+
+def _rotation_y(angle):
+    """Active rotation matrix about y."""
+
+    cos_angle = cos(angle)
+    sin_angle = sin(angle)
+    return (
+        (cos_angle, 0.0, sin_angle),
+        (0.0, 1.0, 0.0),
+        (-sin_angle, 0.0, cos_angle),
+    )
+
+
+def _rotation_z(angle):
+    """Active rotation matrix about z."""
+
+    cos_angle = cos(angle)
+    sin_angle = sin(angle)
+    return (
+        (cos_angle, -sin_angle, 0.0),
+        (sin_angle, cos_angle, 0.0),
+        (0.0, 0.0, 1.0),
+    )
+
+
+def _ecef2enu_rotation(lat0, lon0, deg: bool = True):
+    """Rotation matrix that maps ECEF vectors into ENU coordinates."""
+
+    if deg:
+        lat0 = radians(lat0)
+        lon0 = radians(lon0)
+
+    sin_lat = sin(lat0)
+    cos_lat = cos(lat0)
+    sin_lon = sin(lon0)
+    cos_lon = cos(lon0)
+
+    return (
+        (-sin_lon, cos_lon, 0.0),
+        (-sin_lat * cos_lon, -sin_lat * sin_lon, cos_lat),
+        (cos_lat * cos_lon, cos_lat * sin_lon, sin_lat),
+    )
+
+
+def _enu2ecef_rotation(lat0, lon0, deg: bool = True):
+    """Rotation matrix that maps ENU vectors into ECEF coordinates."""
+
+    return _transpose3(_ecef2enu_rotation(lat0, lon0, deg=deg))
+
+
+def _ecef2ned_rotation(lat0, lon0, deg: bool = True):
+    """Rotation matrix that maps ECEF vectors into NED coordinates."""
+
+    return _matmul3(_ENU2NED, _ecef2enu_rotation(lat0, lon0, deg=deg))
+
+
+def _ned2ecef_rotation(lat0, lon0, deg: bool = True):
+    """Rotation matrix that maps NED vectors into ECEF coordinates."""
+
+    return _transpose3(_ecef2ned_rotation(lat0, lon0, deg=deg))
+
+
+def _ecef2sez_rotation(lat0, lon0, deg: bool = True):
+    """Rotation matrix that maps ECEF vectors into SEZ coordinates."""
+
+    return _matmul3(_ENU2SEZ, _ecef2enu_rotation(lat0, lon0, deg=deg))
+
+
+def _sez2ecef_rotation(lat0, lon0, deg: bool = True):
+    """Rotation matrix that maps SEZ vectors into ECEF coordinates."""
+
+    return _transpose3(_ecef2sez_rotation(lat0, lon0, deg=deg))
+
+
+def _ecef2nwu_rotation(lat0, lon0, deg: bool = True):
+    """Rotation matrix that maps ECEF vectors into NWU coordinates."""
+
+    return _matmul3(_ENU2NWU, _ecef2enu_rotation(lat0, lon0, deg=deg))
+
+
+def _nwu2ecef_rotation(lat0, lon0, deg: bool = True):
+    """Rotation matrix that maps NWU vectors into ECEF coordinates."""
+
+    return _transpose3(_ecef2nwu_rotation(lat0, lon0, deg=deg))
+
+
+def _ned2body_rotation(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """
+    Rotation matrix that maps NED vectors into body coordinates.
+
+    The body frame uses aerospace axes: x forward, y right, z down.
+    Positive yaw turns right, positive pitch is nose up, and positive
+    roll is right wing down.
+    """
+
+    if sequence != "zyx":
+        raise ValueError("only sequence='zyx' is currently supported")
+
+    if deg:
+        yaw = radians(yaw)
+        pitch = radians(pitch)
+        roll = radians(roll)
+
+    body2ned = _matmul3(_rotation_z(yaw), _matmul3(_rotation_y(pitch), _rotation_x(roll)))
+    return _transpose3(body2ned)
+
+
+def _body2ned_rotation(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps body vectors into NED coordinates."""
+
+    return _transpose3(_ned2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg))
+
+
+def _enu2body_rotation(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps ENU vectors into body coordinates."""
+
+    return _matmul3(_ned2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg), _ENU2NED)
+
+
+def _body2enu_rotation(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps body vectors into ENU coordinates."""
+
+    return _transpose3(_enu2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg))
+
+
+def _sez2body_rotation(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps SEZ vectors into body coordinates."""
+
+    return _matmul3(_enu2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg), _transpose3(_ENU2SEZ))
+
+
+def _body2sez_rotation(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps body vectors into SEZ coordinates."""
+
+    return _transpose3(_sez2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg))
+
+
+def _nwu2body_rotation(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps NWU vectors into body coordinates."""
+
+    return _matmul3(_enu2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg), _transpose3(_ENU2NWU))
+
+
+def _body2nwu_rotation(yaw, pitch, roll, sequence: str = "zyx", deg: bool = True):
+    """Rotation matrix that maps body vectors into NWU coordinates."""
+
+    return _transpose3(_nwu2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg))
+
+
+def _ecef2body_rotation(
+    lat0,
+    lon0,
+    yaw,
+    pitch,
+    roll,
+    sequence: str = "zyx",
+    frame: str = "ned",
+    deg: bool = True,
+):
+    """Rotation matrix that maps ECEF vectors into body coordinates."""
+
+    if frame == "ned":
+        local2body = _ned2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+        ecef2local = _ecef2ned_rotation(lat0, lon0, deg=deg)
+    elif frame == "enu":
+        local2body = _enu2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+        ecef2local = _ecef2enu_rotation(lat0, lon0, deg=deg)
+    elif frame == "sez":
+        local2body = _sez2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+        ecef2local = _ecef2sez_rotation(lat0, lon0, deg=deg)
+    elif frame == "nwu":
+        local2body = _nwu2body_rotation(yaw, pitch, roll, sequence=sequence, deg=deg)
+        ecef2local = _ecef2nwu_rotation(lat0, lon0, deg=deg)
+    else:
+        raise ValueError("frame must be one of 'ned', 'enu', 'sez', or 'nwu'")
+
+    return _matmul3(local2body, ecef2local)
+
+
+def _body2ecef_rotation(
+    lat0,
+    lon0,
+    yaw,
+    pitch,
+    roll,
+    sequence: str = "zyx",
+    frame: str = "ned",
+    deg: bool = True,
+):
+    """Rotation matrix that maps body vectors into ECEF coordinates."""
+
+    return _transpose3(
+        _ecef2body_rotation(
+            lat0,
+            lon0,
+            yaw,
+            pitch,
+            roll,
+            sequence=sequence,
+            frame=frame,
+            deg=deg,
+        )
+    )
+
+
+_ENU2NED = (
+    (0.0, 1.0, 0.0),
+    (1.0, 0.0, 0.0),
+    (0.0, 0.0, -1.0),
+)
+
+_ENU2SEZ = (
+    (0.0, -1.0, 0.0),
+    (1.0, 0.0, 0.0),
+    (0.0, 0.0, 1.0),
+)
+
+_ENU2NWU = (
+    (0.0, 1.0, 0.0),
+    (-1.0, 0.0, 0.0),
+    (0.0, 0.0, 1.0),
+)
diff --git a/src/pymap3d/ned.py b/src/pymap3d/ned.py
index 94cb98a..589a871 100644
--- a/src/pymap3d/ned.py
+++ b/src/pymap3d/ned.py
@@ -5,6 +5,7 @@
 from .ecef import ecef2enu, ecef2enuv, ecef2geodetic, enu2ecef
 from .ellipsoid import Ellipsoid
 from .enu import aer2enu, enu2aer, geodetic2enu
+from .frames import _ecef2ned_rotation, _ned2ecef_rotation
 
 __all__ = [
     "aer2ned",
@@ -14,6 +15,8 @@
     "ecef2ned",
     "geodetic2ned",
     "ecef2nedv",
+    "ecef2ned_matrix",
+    "ned2ecef_matrix",
 ]
 
 
@@ -303,3 +306,19 @@ def ecef2nedv(x, y, z, lat0, lon0, deg: bool = True) -> tuple[float, float, floa
     e, n, u = ecef2enuv(x, y, z, lat0, lon0, deg=deg)
 
     return n, e, -u
+
+
+def ecef2ned_matrix(lat0, lon0, deg: bool = True):
+    """
+    Rotation matrix that maps ECEF vectors into NED coordinates.
+    """
+
+    return _ecef2ned_rotation(lat0, lon0, deg=deg)
+
+
+def ned2ecef_matrix(lat0, lon0, deg: bool = True):
+    """
+    Rotation matrix that maps NED vectors into ECEF coordinates.
+    """
+
+    return _ned2ecef_rotation(lat0, lon0, deg=deg)
diff --git a/src/pymap3d/nwu.py b/src/pymap3d/nwu.py
new file mode 100644
index 0000000..44a58a2
--- /dev/null
+++ b/src/pymap3d/nwu.py
@@ -0,0 +1,141 @@
+"""Transforms involving NWU North West Up."""
+
+from __future__ import annotations
+
+from .ecef import ecef2enu, ecef2enuv, enu2ecef
+from .ellipsoid import Ellipsoid
+from .enu import aer2enu, enu2aer, enu2ecefv, geodetic2enu, enu2geodetic
+from .frames import _ecef2nwu_rotation, _nwu2ecef_rotation
+
+__all__ = [
+    "aer2nwu",
+    "nwu2aer",
+    "nwu2geodetic",
+    "nwu2ecef",
+    "nwu2ecefv",
+    "ecef2nwu",
+    "ecef2nwuv",
+    "geodetic2nwu",
+    "ecef2nwu_matrix",
+    "nwu2ecef_matrix",
+]
+
+
+def aer2nwu(az, elev, slantRange, deg: bool = True) -> tuple:
+    """
+    Convert azimuth, elevation, range to North, West, Up.
+    """
+
+    e, n, u = aer2enu(az, elev, slantRange, deg=deg)
+    return n, -e, u
+
+
+def nwu2aer(n, w, u, deg: bool = True) -> tuple:
+    """
+    Convert North, West, Up to azimuth, elevation, range.
+    """
+
+    return enu2aer(-w, n, u, deg=deg)
+
+
+def nwu2geodetic(
+    n,
+    w,
+    u,
+    lat0,
+    lon0,
+    h0,
+    ell: Ellipsoid | None = None,
+    deg: bool = True,
+) -> tuple:
+    """
+    Convert North, West, Up to target geodetic coordinates.
+    """
+
+    return enu2geodetic(-w, n, u, lat0, lon0, h0, ell, deg=deg)
+
+
+def nwu2ecef(
+    n,
+    w,
+    u,
+    lat0,
+    lon0,
+    h0,
+    ell: Ellipsoid | None = None,
+    deg: bool = True,
+) -> tuple:
+    """
+    Convert North, West, Up to target ECEF coordinates.
+    """
+
+    return enu2ecef(-w, n, u, lat0, lon0, h0, ell, deg=deg)
+
+
+def nwu2ecefv(n, w, u, lat0, lon0, deg: bool = True) -> tuple:
+    """
+    Convert an NWU vector into ECEF coordinates.
+    """
+
+    return enu2ecefv(-w, n, u, lat0, lon0, deg=deg)
+
+
+def ecef2nwu(
+    x,
+    y,
+    z,
+    lat0,
+    lon0,
+    h0,
+    ell: Ellipsoid | None = None,
+    deg: bool = True,
+) -> tuple:
+    """
+    Convert target ECEF coordinates to North, West, Up.
+    """
+
+    e, n, u = ecef2enu(x, y, z, lat0, lon0, h0, ell, deg=deg)
+    return n, -e, u
+
+
+def ecef2nwuv(x, y, z, lat0, lon0, deg: bool = True) -> tuple:
+    """
+    Convert an ECEF vector into North, West, Up.
+    """
+
+    e, n, u = ecef2enuv(x, y, z, lat0, lon0, deg=deg)
+    return n, -e, u
+
+
+def geodetic2nwu(
+    lat,
+    lon,
+    h,
+    lat0,
+    lon0,
+    h0,
+    ell: Ellipsoid | None = None,
+    deg: bool = True,
+) -> tuple:
+    """
+    Convert target geodetic coordinates to North, West, Up.
+    """
+
+    e, n, u = geodetic2enu(lat, lon, h, lat0, lon0, h0, ell, deg=deg)
+    return n, -e, u
+
+
+def ecef2nwu_matrix(lat0, lon0, deg: bool = True):
+    """
+    Rotation matrix that maps ECEF vectors into NWU coordinates.
+    """
+
+    return _ecef2nwu_rotation(lat0, lon0, deg=deg)
+
+
+def nwu2ecef_matrix(lat0, lon0, deg: bool = True):
+    """
+    Rotation matrix that maps NWU vectors into ECEF coordinates.
+    """
+
+    return _nwu2ecef_rotation(lat0, lon0, deg=deg)
diff --git a/src/pymap3d/sez.py b/src/pymap3d/sez.py
new file mode 100644
index 0000000..acc1f8a
--- /dev/null
+++ b/src/pymap3d/sez.py
@@ -0,0 +1,141 @@
+"""Transforms involving SEZ South East Zenith."""
+
+from __future__ import annotations
+
+from .ecef import ecef2enu, ecef2enuv, enu2ecef
+from .ellipsoid import Ellipsoid
+from .enu import aer2enu, enu2aer, enu2ecefv, geodetic2enu, enu2geodetic
+from .frames import _ecef2sez_rotation, _sez2ecef_rotation
+
+__all__ = [
+    "aer2sez",
+    "sez2aer",
+    "sez2geodetic",
+    "sez2ecef",
+    "sez2ecefv",
+    "ecef2sez",
+    "ecef2sezv",
+    "geodetic2sez",
+    "ecef2sez_matrix",
+    "sez2ecef_matrix",
+]
+
+
+def aer2sez(az, elev, slantRange, deg: bool = True) -> tuple:
+    """
+    Convert azimuth, elevation, range to South, East, Zenith.
+    """
+
+    e, n, u = aer2enu(az, elev, slantRange, deg=deg)
+    return -n, e, u
+
+
+def sez2aer(s, e, z, deg: bool = True) -> tuple:
+    """
+    Convert South, East, Zenith to azimuth, elevation, range.
+    """
+
+    return enu2aer(e, -s, z, deg=deg)
+
+
+def sez2geodetic(
+    s,
+    e,
+    z,
+    lat0,
+    lon0,
+    h0,
+    ell: Ellipsoid | None = None,
+    deg: bool = True,
+) -> tuple:
+    """
+    Convert South, East, Zenith to target geodetic coordinates.
+    """
+
+    return enu2geodetic(e, -s, z, lat0, lon0, h0, ell, deg=deg)
+
+
+def sez2ecef(
+    s,
+    e,
+    z,
+    lat0,
+    lon0,
+    h0,
+    ell: Ellipsoid | None = None,
+    deg: bool = True,
+) -> tuple:
+    """
+    Convert South, East, Zenith to target ECEF coordinates.
+    """
+
+    return enu2ecef(e, -s, z, lat0, lon0, h0, ell, deg=deg)
+
+
+def sez2ecefv(s, e, z, lat0, lon0, deg: bool = True) -> tuple:
+    """
+    Convert a SEZ vector into ECEF coordinates.
+    """
+
+    return enu2ecefv(e, -s, z, lat0, lon0, deg=deg)
+
+
+def ecef2sez(
+    x,
+    y,
+    z,
+    lat0,
+    lon0,
+    h0,
+    ell: Ellipsoid | None = None,
+    deg: bool = True,
+) -> tuple:
+    """
+    Convert target ECEF coordinates to South, East, Zenith.
+    """
+
+    e, n, u = ecef2enu(x, y, z, lat0, lon0, h0, ell, deg=deg)
+    return -n, e, u
+
+
+def ecef2sezv(x, y, z, lat0, lon0, deg: bool = True) -> tuple:
+    """
+    Convert an ECEF vector into South, East, Zenith.
+    """
+
+    e, n, u = ecef2enuv(x, y, z, lat0, lon0, deg=deg)
+    return -n, e, u
+
+
+def geodetic2sez(
+    lat,
+    lon,
+    h,
+    lat0,
+    lon0,
+    h0,
+    ell: Ellipsoid | None = None,
+    deg: bool = True,
+) -> tuple:
+    """
+    Convert target geodetic coordinates to South, East, Zenith.
+    """
+
+    e, n, u = geodetic2enu(lat, lon, h, lat0, lon0, h0, ell, deg=deg)
+    return -n, e, u
+
+
+def ecef2sez_matrix(lat0, lon0, deg: bool = True):
+    """
+    Rotation matrix that maps ECEF vectors into SEZ coordinates.
+    """
+
+    return _ecef2sez_rotation(lat0, lon0, deg=deg)
+
+
+def sez2ecef_matrix(lat0, lon0, deg: bool = True):
+    """
+    Rotation matrix that maps SEZ vectors into ECEF coordinates.
+    """
+
+    return _sez2ecef_rotation(lat0, lon0, deg=deg)
diff --git a/src/pymap3d/tests/test_body.py b/src/pymap3d/tests/test_body.py
new file mode 100644
index 0000000..8f61f33
--- /dev/null
+++ b/src/pymap3d/tests/test_body.py
@@ -0,0 +1,125 @@
+from math import radians
+
+import pymap3d as pm
+import pytest
+from pytest import approx
+
+lla0 = (42, -82, 200)
+ypr0 = (30, 12, -8)
+
+
+def _matvec3(matrix, vector):
+    return tuple(sum(matrix[i][j] * vector[j] for j in range(3)) for i in range(3))
+
+
+def _matmul3(a, b):
+    return tuple(
+        tuple(sum(a[i][k] * b[k][j] for k in range(3)) for j in range(3)) for i in range(3)
+    )
+
+
+def _assert_matrix_close(actual, expected):
+    for row_actual, row_expected in zip(actual, expected):
+        assert row_actual == approx(row_expected)
+
+
+def test_body_matrix_identity():
+    identity = ((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))
+
+    _assert_matrix_close(pm.body_matrix(0, 0, 0), identity)
+    _assert_matrix_close(pm.ned2body_matrix(0, 0, 0), identity)
+    _assert_matrix_close(pm.body2ned_matrix(0, 0, 0), identity)
+
+
+def test_body_rotation_conventions():
+    assert pm.ned2body(1, 0, 0, 0, 0, 0) == approx((1, 0, 0))
+    assert pm.ned2body(0, 1, 0, 90, 0, 0) == approx((1, 0, 0), abs=1e-12)
+    assert pm.ned2body(0, 0, -1, 0, 90, 0) == approx((1, 0, 0), abs=1e-12)
+    assert pm.ned2body(0, 0, 1, 0, 0, 90) == approx((0, 1, 0), abs=1e-12)
+
+
+def test_body_matrix_inverse_and_radians():
+    yaw, pitch, roll = ypr0
+    matrix = pm.ned2body_matrix(yaw, pitch, roll)
+    inverse = pm.body2ned_matrix(yaw, pitch, roll)
+    identity = ((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))
+
+    _assert_matrix_close(_matmul3(matrix, inverse), identity)
+    _assert_matrix_close(
+        matrix,
+        pm.ned2body_matrix(radians(yaw), radians(pitch), radians(roll), deg=False),
+    )
+
+
+def test_ned_enu_body_consistency():
+    enu = (5, 3, 2)
+    ned = (enu[1], enu[0], -enu[2])
+    sez = (-enu[1], enu[0], enu[2])
+    nwu = (enu[1], -enu[0], enu[2])
+    yaw, pitch, roll = ypr0
+
+    assert pm.enu2body(*enu, yaw, pitch, roll) == approx(pm.ned2body(*ned, yaw, pitch, roll))
+    assert pm.sez2body(*sez, yaw, pitch, roll) == approx(pm.enu2body(*enu, yaw, pitch, roll))
+    assert pm.nwu2body(*nwu, yaw, pitch, roll) == approx(pm.enu2body(*enu, yaw, pitch, roll))
+    assert pm.body2enu(*pm.enu2body(*enu, yaw, pitch, roll), yaw, pitch, roll) == approx(enu)
+    assert pm.body2ned(*pm.ned2body(*ned, yaw, pitch, roll), yaw, pitch, roll) == approx(ned)
+    assert pm.body2sez(*pm.sez2body(*sez, yaw, pitch, roll), yaw, pitch, roll) == approx(sez)
+    assert pm.body2nwu(*pm.nwu2body(*nwu, yaw, pitch, roll), yaw, pitch, roll) == approx(nwu)
+
+
+def test_body_matrices_match_vector_transforms():
+    local_ned = (4, -2, 1)
+    local_enu = (local_ned[1], local_ned[0], -local_ned[2])
+    local_sez = (-local_enu[1], local_enu[0], local_enu[2])
+    local_nwu = (local_enu[1], -local_enu[0], local_enu[2])
+    yaw, pitch, roll = ypr0
+
+    assert _matvec3(pm.ned2body_matrix(yaw, pitch, roll), local_ned) == approx(
+        pm.ned2body(*local_ned, yaw, pitch, roll)
+    )
+    assert _matvec3(pm.enu2body_matrix(yaw, pitch, roll), local_enu) == approx(
+        pm.enu2body(*local_enu, yaw, pitch, roll)
+    )
+    assert _matvec3(pm.sez2body_matrix(yaw, pitch, roll), local_sez) == approx(
+        pm.sez2body(*local_sez, yaw, pitch, roll)
+    )
+    assert _matvec3(pm.nwu2body_matrix(yaw, pitch, roll), local_nwu) == approx(
+        pm.nwu2body(*local_nwu, yaw, pitch, roll)
+    )
+
+
+@pytest.mark.parametrize("frame", ["ned", "enu", "sez", "nwu"])
+def test_ecef_body_roundtrip(frame):
+    yaw, pitch, roll = ypr0
+    target_ecef = pm.aer2ecef(33, 40, 1200, *lla0)
+
+    body = pm.ecef2body(*target_ecef, *lla0, yaw, pitch, roll, frame=frame)
+    assert pm.body2ecef(*body, *lla0, yaw, pitch, roll, frame=frame) == approx(target_ecef)
+
+
+@pytest.mark.parametrize("frame", ["ned", "enu", "sez", "nwu"])
+def test_ecef_body_vector_roundtrip(frame):
+    yaw, pitch, roll = ypr0
+    ecef_vector = (5, 3, 2)
+
+    body_vector = pm.ecef2bodyv(*ecef_vector, *lla0[:2], yaw, pitch, roll, frame=frame)
+    assert pm.body2ecefv(*body_vector, *lla0[:2], yaw, pitch, roll, frame=frame) == approx(
+        ecef_vector
+    )
+
+
+@pytest.mark.parametrize("frame", ["ned", "enu", "sez", "nwu"])
+def test_geodetic_body_roundtrip(frame):
+    yaw, pitch, roll = ypr0
+    target_lla = pm.aer2geodetic(33, 40, 1200, *lla0)
+
+    body = pm.geodetic2body(*target_lla, *lla0, yaw, pitch, roll, frame=frame)
+    assert pm.body2geodetic(*body, *lla0, yaw, pitch, roll, frame=frame) == approx(target_lla)
+
+
+def test_body_invalid_options():
+    with pytest.raises(ValueError):
+        pm.body_matrix(0, 0, 0, sequence="xyz")
+
+    with pytest.raises(ValueError):
+        pm.ecef2body(1, 2, 3, *lla0, 0, 0, 0, frame="xyz")
diff --git a/src/pymap3d/tests/test_local_frames.py b/src/pymap3d/tests/test_local_frames.py
new file mode 100644
index 0000000..0b933ff
--- /dev/null
+++ b/src/pymap3d/tests/test_local_frames.py
@@ -0,0 +1,94 @@
+from math import radians
+
+import pymap3d as pm
+from pytest import approx
+
+lla0 = (42, -82, 200)
+aer0 = (33, 70, 1000)
+
+
+def _matvec3(matrix, vector):
+    return tuple(sum(matrix[i][j] * vector[j] for j in range(3)) for i in range(3))
+
+
+def _matmul3(a, b):
+    return tuple(
+        tuple(sum(a[i][k] * b[k][j] for k in range(3)) for j in range(3)) for i in range(3)
+    )
+
+
+def _assert_matrix_close(actual, expected):
+    for row_actual, row_expected in zip(actual, expected):
+        assert row_actual == approx(row_expected)
+
+
+def test_sez_conversions():
+    enu = pm.aer2enu(*aer0)
+    sez = (-enu[1], enu[0], enu[2])
+    xyz = pm.aer2ecef(*aer0, *lla0)
+
+    assert pm.aer2sez(*aer0) == approx(sez)
+    assert pm.sez2aer(*sez) == approx(aer0)
+    assert pm.ecef2sez(*xyz, *lla0) == approx(sez)
+    assert pm.sez2ecef(*sez, *lla0) == approx(xyz)
+
+    lla1 = pm.aer2geodetic(*aer0, *lla0)
+    assert pm.geodetic2sez(*lla1, *lla0) == approx(sez)
+    assert pm.sez2geodetic(*sez, *lla0) == approx(lla1)
+
+
+def test_nwu_conversions():
+    enu = pm.aer2enu(*aer0)
+    nwu = (enu[1], -enu[0], enu[2])
+    xyz = pm.aer2ecef(*aer0, *lla0)
+
+    assert pm.aer2nwu(*aer0) == approx(nwu)
+    assert pm.nwu2aer(*nwu) == approx(aer0)
+    assert pm.ecef2nwu(*xyz, *lla0) == approx(nwu)
+    assert pm.nwu2ecef(*nwu, *lla0) == approx(xyz)
+
+    lla1 = pm.aer2geodetic(*aer0, *lla0)
+    assert pm.geodetic2nwu(*lla1, *lla0) == approx(nwu)
+    assert pm.nwu2geodetic(*nwu, *lla0) == approx(lla1)
+
+
+def test_vector_frame_conversions():
+    ecef_vector = (5, 3, 2)
+    enu_vector = pm.ecef2enuv(*ecef_vector, *lla0[:2])
+    sez_vector = (-enu_vector[1], enu_vector[0], enu_vector[2])
+    nwu_vector = (enu_vector[1], -enu_vector[0], enu_vector[2])
+
+    assert pm.ecef2sezv(*ecef_vector, *lla0[:2]) == approx(sez_vector)
+    assert pm.sez2ecefv(*sez_vector, *lla0[:2]) == approx(ecef_vector)
+    assert pm.ecef2nwuv(*ecef_vector, *lla0[:2]) == approx(nwu_vector)
+    assert pm.nwu2ecefv(*nwu_vector, *lla0[:2]) == approx(ecef_vector)
+
+
+def test_frame_matrices_match_vector_transforms():
+    ecef_vector = (5, 3, 2)
+    lat0, lon0 = lla0[:2]
+
+    frame_cases = [
+        (pm.ecef2enu_matrix(lat0, lon0), pm.enu2ecef_matrix(lat0, lon0), pm.ecef2enuv(*ecef_vector, lat0, lon0)),
+        (pm.ecef2ned_matrix(lat0, lon0), pm.ned2ecef_matrix(lat0, lon0), pm.ecef2nedv(*ecef_vector, lat0, lon0)),
+        (pm.ecef2sez_matrix(lat0, lon0), pm.sez2ecef_matrix(lat0, lon0), pm.ecef2sezv(*ecef_vector, lat0, lon0)),
+        (pm.ecef2nwu_matrix(lat0, lon0), pm.nwu2ecef_matrix(lat0, lon0), pm.ecef2nwuv(*ecef_vector, lat0, lon0)),
+    ]
+
+    identity = ((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))
+
+    for forward, inverse, expected in frame_cases:
+        assert _matvec3(forward, ecef_vector) == approx(expected)
+        assert _matvec3(inverse, expected) == approx(ecef_vector)
+        _assert_matrix_close(_matmul3(forward, inverse), identity)
+
+
+def test_frame_matrices_rad_inputs():
+    lat0, lon0 = lla0[:2]
+    rlat0 = radians(lat0)
+    rlon0 = radians(lon0)
+
+    _assert_matrix_close(pm.ecef2enu_matrix(lat0, lon0), pm.ecef2enu_matrix(rlat0, rlon0, deg=False))
+    _assert_matrix_close(pm.ecef2ned_matrix(lat0, lon0), pm.ecef2ned_matrix(rlat0, rlon0, deg=False))
+    _assert_matrix_close(pm.ecef2sez_matrix(lat0, lon0), pm.ecef2sez_matrix(rlat0, rlon0, deg=False))
+    _assert_matrix_close(pm.ecef2nwu_matrix(lat0, lon0), pm.ecef2nwu_matrix(rlat0, rlon0, deg=False))