fix(python): mixed SGP4/SDP4 propagate

.github/workflows/python.yaml

@@ -29,6 +29,22 @@ jobs:
           - os: ubuntu-latest
             python: '3.10'
             target: linux-x86_64
+          # Linux aarch64 (Graviton, Apple Silicon containers)
+          - os: ubuntu-24.04-arm
+            python: '3.14'
+            target: linux-aarch64
+          - os: ubuntu-24.04-arm
+            python: '3.13'
+            target: linux-aarch64
+          - os: ubuntu-24.04-arm
+            python: '3.12'
+            target: linux-aarch64
+          - os: ubuntu-24.04-arm
+            python: '3.11'
+            target: linux-aarch64
+          - os: ubuntu-24.04-arm
+            python: '3.10'
+            target: linux-aarch64
           # macOS ARM64
           - os: macos-latest
             python: '3.14'
@@ -76,13 +92,20 @@ jobs:
           echo "Python libdir: $PYTHON_LIBDIR"
           echo "Python lib: $PYTHON_LIB"
 
-          # On macOS, don't link against Python library - symbols resolve at load time
+          # On macOS, don't link against Python library - symbols resolve at load time.
+          # On Linux, pin the target explicitly: without -Dtarget Zig compiles for
+          # the runner's own CPU and glibc (AVX-512 baseline, glibc 2.39 on
+          # ubuntu-latest), and the wheel then SIGILLs on non-AVX-512 hosts and
+          # refuses to load on glibc < 2.36 (Amazon Linux 2023 has 2.34).
+          # glibc 2.28 = manylinux_2_28; a baseline CPU keeps the wheel portable
+          # while oma still dispatches to AVX2/AVX-512/SVE at runtime.
           if [[ "$RUNNER_OS" == "macOS" ]]; then
             zig build python-bindings \
               -Dpython-include="$PYTHON_INCLUDE" \
               -Doptimize=ReleaseFast
           else
             zig build python-bindings \
+              -Dtarget="$(uname -m)-linux-gnu.2.28" \
               -Dpython-include="$PYTHON_INCLUDE" \
               -Dpython-lib="$PYTHON_LIB" \
               -Dpython-lib-path="$PYTHON_LIBDIR" \

bindings/python/astroz/__init__.py

@@ -287,6 +287,55 @@ def _start_jd(start_time):
     return 2440587.5 + (start_time.timestamp() / 86400.0)
 
 
+def _gmst_rad(jd):
+    """Mirror of WorldCoordinateSystem.julianToGmst for an array of Julian dates."""
+    d = np.asarray(jd, dtype=np.float64) - 2451545.0
+    t = d / 36525.0
+    gmst = 280.46061837 + 360.98564736629 * d + 0.000387933 * t * t - t * t * t / 38710000.0
+    return np.deg2rad(np.mod(gmst, 360.0))
+
+
+_EARTH_ROTATION_RATE = 7.2921150e-5  # rad/s, constants.earth.rotationRate
+_WGS84_A = 6378.137
+_WGS84_E2 = 6.69437999014e-3
+
+
+def _teme_to_ecef_into(pos, vel, jd, sats):
+    """Rotate pos[:, sats] (and vel) from TEME to ECEF in place. Velocity is
+    the derivative of the ECEF track (rotation minus omega x r), matching the
+    Zig writeOutput path."""
+    theta = _gmst_rad(jd)[:, None]
+    c, s = np.cos(theta), np.sin(theta)
+    p = pos[:, sats]
+    x, y = p[..., 0].copy(), p[..., 1].copy()
+    p[..., 0] = c * x + s * y
+    p[..., 1] = c * y - s * x
+    if vel is not None:
+        v = vel[:, sats]
+        vx, vy = v[..., 0].copy(), v[..., 1].copy()
+        v[..., 0] = c * vx + s * vy + _EARTH_ROTATION_RATE * p[..., 1]
+        v[..., 1] = c * vy - s * vx - _EARTH_ROTATION_RATE * p[..., 0]
+
+
+def _ecef_to_geodetic_into(pos, sats):
+    """WGS84 ECEF km -> (lat deg, lon deg, alt km) in place; mirrors
+    WorldCoordinateSystem.ecefToGeodeticDeg (iterative, 10 rounds)."""
+    p = pos[:, sats]
+    x, y, z = p[..., 0].copy(), p[..., 1].copy(), p[..., 2].copy()
+    lon = np.arctan2(y, x)
+    r = np.sqrt(x * x + y * y)
+    lat = np.arctan2(z, r * (1.0 - _WGS84_E2))
+    for _ in range(10):
+        sin_lat = np.sin(lat)
+        n = _WGS84_A / np.sqrt(1.0 - _WGS84_E2 * sin_lat * sin_lat)
+        lat = np.arctan2(z + _WGS84_E2 * n * sin_lat, r)
+    sin_lat = np.sin(lat)
+    n = _WGS84_A / np.sqrt(1.0 - _WGS84_E2 * sin_lat * sin_lat)
+    p[..., 0] = np.degrees(lat)
+    p[..., 1] = np.degrees(lon)
+    p[..., 2] = r / np.cos(lat) - n
+
+
 def _pad_epochs_for_simd(epochs_arr):
     """Pad epoch array to multiple of 8 for SIMD alignment."""
     n = len(epochs_arr)
@@ -528,12 +577,12 @@ def propagate(
             output=output,
             reference_jd=start,
             time_major=True,
+            output_stride=n_sats,
         )
 
-    # SDP4 batch: positions [n_sgp4, n_sats), time-major with full-width stride
+    # SDP4 batch: positions [n_sgp4, n_sats), time-major with full-width stride.
+    # The SDP4 kernel emits TEME; ECEF / geodetic are applied here in numpy.
     if const._sdp4_satrecs:
-        if output != "teme":
-            raise NotImplementedError("deep-space (SDP4) satellites currently support output='teme' only")
         vel_out = vel if velocities else np.empty(shape, dtype=np.float64)
         _sdp4_batch_propagate_into(
             const._sdp4_satrecs,
@@ -544,6 +593,11 @@ def propagate(
             output_stride=n_sats,
             sat_offset=n_sgp4,
         )
+        if output != "teme":
+            sdp4 = slice(n_sgp4, n_sats)
+            _teme_to_ecef_into(pos, vel, start + times / 1440.0, sdp4)
+            if output == "geodetic":
+                _ecef_to_geodetic_into(pos, sdp4)
 
     return (pos, vel) if velocities else pos
 

src/Constellation.zig

@@ -494,15 +494,20 @@ inline fn writeOutput(
             const origIdx = origIndices[base + lane];
             const ob = outBase(layout, @as(usize, origIdx), timeIdx, ctx.numTimes, ctx.numSatellites);
             ctx.resultsPos[ob..][0..3].* = switch (mode) {
-                .geodetic => WCS.ecefToGeodetic(eciToEcef(results[lane][0], sinG, cosG)),
+                .geodetic => WCS.ecefToGeodeticDeg(eciToEcef(results[lane][0], sinG, cosG)),
                 .ecef => eciToEcef(results[lane][0], sinG, cosG),
                 .teme => results[lane][0],
             };
             if (hasVel) {
-                ctx.resultsVel.?[ob..][0..3].* = if (mode != .teme)
-                    eciToEcef(results[lane][1], sinG, cosG)
-                else
-                    results[lane][1];
+                ctx.resultsVel.?[ob..][0..3].* = if (mode != .teme) blk: {
+                    // Velocity in the rotating frame: rotate, then remove the
+                    // frame's own motion (omega x r). Without this a GEO bird
+                    // shows ~3 km/s in ECEF instead of ~0.
+                    const p = eciToEcef(results[lane][0], sinG, cosG);
+                    const v = eciToEcef(results[lane][1], sinG, cosG);
+                    const w = constants.earth.rotationRate;
+                    break :blk .{ v[0] + w * p[1], v[1] - w * p[0], v[2] };
+                } else results[lane][1];
             }
         }
     }
@@ -961,4 +966,26 @@ test "ECEF output mode" {
         try testing.expectApproxEqAbs(ey, posEcef[ob + 1], 1e-6);
         try testing.expectApproxEqAbs(posTeme[ob + 2], posEcef[ob + 2], 1e-6);
     }
+
+    // ECEF velocity is the time derivative of the ECEF track, so a 1 s
+    // finite difference must match it (LEO acceleration bounds the error at
+    // ~0.005 km/s over half a second) and the GEO bird must sit nearly still.
+    const dtDays = 1.0 / 86400.0;
+    const jd2 = [_]f64{ @floor(epoch), @floor(epoch) };
+    const fr2 = [_]f64{ epoch - @floor(epoch), epoch - @floor(epoch) + dtDays };
+    var pos2: [2 * 3 * 3]f64 = undefined;
+    var vel2: [2 * 3 * 3]f64 = undefined;
+    c.resetCarry();
+    try c.propagate(&jd2, &fr2, &pos2, &vel2, .ecef, .timeMajor);
+    for (0..3) |sat| {
+        const a = sat * 3;
+        const b = 9 + sat * 3;
+        for (0..3) |k| {
+            const fd = pos2[b + k] - pos2[a + k];
+            try testing.expectApproxEqAbs(fd, vel2[a + k], 0.01);
+        }
+    }
+    const geo = 3;
+    const geoSpeed = @sqrt(vel2[geo] * vel2[geo] + vel2[geo + 1] * vel2[geo + 1] + vel2[geo + 2] * vel2[geo + 2]);
+    try testing.expect(geoSpeed < 0.1);
 }