diff --git a/.github/workflows/python.yaml b/.github/workflows/python.yaml index 7dcb527..b084959 100644 --- a/.github/workflows/python.yaml +++ b/.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" \ diff --git a/bindings/python/astroz/__init__.py b/bindings/python/astroz/__init__.py index 763c4b5..877d338 100644 --- a/bindings/python/astroz/__init__.py +++ b/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 diff --git a/src/Constellation.zig b/src/Constellation.zig index e4617ab..28b6cd0 100644 --- a/src/Constellation.zig +++ b/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); }