feat: time-differenced ionosphere (ROT) and per-arc TEC detrending

computeIonoRate: bias-free sequential or standardized-baseline
differences (TECU/min) and second undivided differences, never
crossing arc boundaries; detrendIonoArcs anchors arcs at their first
observation. Ground-truth tested against a known dTEC/dt. Suggested by
Hans van der Marel.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: MiguelPuntoEs

CHANGELOG.md

@@ -1,5 +1,11 @@
 # Changelog
 
+## 1.6.0
+
+### New features
+
+- **Time-differenced ionosphere (rate of TEC)** — `computeIonoRate(result, intervalSec?, order?)`: sequential ΔSTEC/Δt in TECU/min with all biases cancelled; a standardized baseline (e.g. 60 s, all pairs) gives a sample-rate-independent picture; `order: 2` returns the second undivided difference (gradients removed, scintillation/noise amplified). Differences never cross arc boundaries, so cycle slips do not pollute the series. `detrendIonoArcs(result)` removes the per-arc bias by anchoring each arc at its first observation, showing pure TEC variation. `IonoSeries` gains `arcStarts`. Suggested by Hans van der Marel.
+
 ## 1.5.1
 
 ### Fixed

package.json

@@ -1,6 +1,6 @@
 {
   "name": "gnss-js",
-  "version": "1.5.1",
+  "version": "1.6.0",
   "description": "Comprehensive GNSS library for JavaScript — time scales, coordinates, RINEX parsing, RTCM3 decoding, orbit computation, and signal analysis",
   "type": "module",
   "sideEffects": false,

src/analysis/ionosphere.ts

@@ -58,6 +58,12 @@ export interface IonoSeries {
   tecuPerNs: number;
   /** Time series of slant TEC values (TECU, DCB-biased). */
   points: IonoPoint[];
+  /**
+   * Indices into `points` where a new continuous arc begins (gap,
+   * cycle slip, or geometry-free jump). Differencing and per-arc
+   * detrending must not cross these boundaries.
+   */
+  arcStarts: number[];
 }
 
 export interface IonoResult {
@@ -107,6 +113,7 @@ interface PairState {
 export class IonoAccumulator {
   private state = new Map<string, Map<string, PairState>>();
   private closed = new Map<string, Map<string, IonoPoint[]>>();
+  private closedArcStarts = new Map<string, Map<string, number[]>>();
   private interval: number;
   private obsIndices: Map<string, Map<string, { L: number; C: number | null }>>;
   private gloChannels: Record<string, number>;
@@ -254,6 +261,17 @@ export class IonoAccumulator {
         points = [];
         satArcs.set(pairKey, points);
       }
+      let satStarts = this.closedArcStarts.get(prn);
+      if (!satStarts) {
+        satStarts = new Map();
+        this.closedArcStarts.set(prn, satStarts);
+      }
+      let starts = satStarts.get(pairKey);
+      if (!starts) {
+        starts = [];
+        satStarts.set(pairKey, starts);
+      }
+      starts.push(points.length);
       for (let k = 0; k < arc.times.length; k++) {
         points.push({
           time: arc.times[k]!,
@@ -289,7 +307,7 @@ export class IonoAccumulator {
       }
       if (!bestKey) continue;
       const points = satArcs.get(bestKey)!;
-      points.sort((a, b) => a.time - b.time);
+      const arcStarts = this.closedArcStarts.get(prn)?.get(bestKey) ?? [0];
       const sys = prn[0]!;
       for (const p of points) {
         sum += p.stec;
@@ -308,6 +326,7 @@ export class IonoAccumulator {
         codes: this.pairCodes.get(`${sys}:${bestKey}`) ?? ['', ''],
         tecuPerNs,
         points,
+        arcStarts,
       });
     }
     series.sort((a, b) => a.prn.localeCompare(b.prn));
@@ -319,3 +338,132 @@ export class IonoAccumulator {
     };
   }
 }
+
+/* ================================================================== */
+/*  Time-differenced ionosphere (rate of TEC)                          */
+/* ================================================================== */
+
+export interface IonoRatePoint {
+  /** Epoch time (ms) — the later epoch of the differenced pair. */
+  time: number;
+  /** ROT in TECU/min (order 1) or undivided ΔΔSTEC in TECU (order 2). */
+  value: number;
+}
+
+export interface IonoRateSeries {
+  prn: string;
+  system: string;
+  points: IonoRatePoint[];
+}
+
+/**
+ * Sequential time differences of the slant TEC series — the biases
+ * (ambiguities, DCBs) cancel, leaving ionospheric rate of TEC, phase
+ * noise, and scintillation. Differences never cross arc boundaries,
+ * so cycle slips do not appear as outliers (arcs split there).
+ *
+ * @param intervalSec Standardized differencing baseline in seconds.
+ *   Omit for native-rate sequential differences; set e.g. 60 for a
+ *   sample-rate-independent picture (all pairs ~1 min apart are used).
+ * @param order 1 (default) = first difference in TECU/min;
+ *   2 = second undivided difference in TECU (native rate only,
+ *   gradients removed, noise and scintillation amplified).
+ */
+export function computeIonoRate(
+  result: IonoResult,
+  intervalSec?: number,
+  order: 1 | 2 = 1
+): IonoRateSeries[] {
+  const out: IonoRateSeries[] = [];
+  for (const s of result.series) {
+    const pts: IonoRatePoint[] = [];
+    const bounds = [...s.arcStarts, s.points.length];
+    for (let a = 0; a < bounds.length - 1; a++) {
+      const lo = bounds[a]!;
+      const hi = bounds[a + 1]!;
+      if (order === 2) {
+        // Second undivided difference on (approximately) uniform spacing
+        for (let i = lo + 2; i < hi; i++) {
+          const dt1 = s.points[i - 1]!.time - s.points[i - 2]!.time;
+          const dt2 = s.points[i]!.time - s.points[i - 1]!.time;
+          if (Math.abs(dt1 - dt2) > 0.1 * Math.max(dt1, dt2)) continue;
+          pts.push({
+            time: s.points[i]!.time,
+            value:
+              s.points[i]!.stec -
+              2 * s.points[i - 1]!.stec +
+              s.points[i - 2]!.stec,
+          });
+        }
+        continue;
+      }
+      if (intervalSec === undefined) {
+        for (let i = lo + 1; i < hi; i++) {
+          const dtMin = (s.points[i]!.time - s.points[i - 1]!.time) / 60_000;
+          if (dtMin <= 0) continue;
+          pts.push({
+            time: s.points[i]!.time,
+            value: (s.points[i]!.stec - s.points[i - 1]!.stec) / dtMin,
+          });
+        }
+      } else {
+        // All pairs ~intervalSec apart (two-pointer over the arc)
+        const targetMs = intervalSec * 1000;
+        let j = lo;
+        for (let i = lo; i < hi; i++) {
+          if (j <= i) j = i + 1;
+          while (j < hi && s.points[j]!.time - s.points[i]!.time < targetMs)
+            j++;
+          if (j >= hi) break;
+          const dtMs = s.points[j]!.time - s.points[i]!.time;
+          if (dtMs > 1.5 * targetMs) continue; // gap-ish, not a clean pair
+          pts.push({
+            time: s.points[j]!.time,
+            value: ((s.points[j]!.stec - s.points[i]!.stec) / dtMs) * 60_000,
+          });
+        }
+      }
+    }
+    if (pts.length > 0) out.push({ prn: s.prn, system: s.system, points: pts });
+  }
+  return out;
+}
+
+/**
+ * Remove the per-arc bias from a slant TEC result — each arc is
+ * shifted so its first observation reads zero, leaving only the TEC
+ * *variation* along the arc (Hans van der Marel's suggestion; using
+ * the first observation keeps the shape untouched at the small risk
+ * of anchoring on an outlier).
+ */
+export function detrendIonoArcs(result: IonoResult): IonoResult {
+  let sum = 0;
+  let count = 0;
+  let maxStec = -Infinity;
+  const series = result.series.map((s) => {
+    const bounds = [...s.arcStarts, s.points.length];
+    const points: IonoPoint[] = new Array<IonoPoint>(s.points.length);
+    for (let a = 0; a < bounds.length - 1; a++) {
+      const lo = bounds[a]!;
+      const hi = bounds[a + 1]!;
+      const anchor = s.points[lo]?.stec ?? 0;
+      for (let i = lo; i < hi; i++) {
+        points[i] = {
+          time: s.points[i]!.time,
+          stec: s.points[i]!.stec - anchor,
+        };
+      }
+    }
+    for (const p of points) {
+      sum += p.stec;
+      count++;
+      if (p.stec > maxStec) maxStec = p.stec;
+    }
+    return { ...s, points };
+  });
+  return {
+    series,
+    maxStec: count > 0 ? maxStec : 0,
+    meanStec: count > 0 ? sum / count : 0,
+  };
+}

test/iono-rate.test.ts

@@ -0,0 +1,92 @@
+/**
+ * Ground-truth tests for the time-differenced ionosphere helpers.
+ * A synthetic TEC(t) with known derivative pins the ROT computation;
+ * differencing across arc boundaries and sample-rate dependence are
+ * the two failure modes checked explicitly.
+ */
+import { describe, expect, it } from 'vitest';
+import type { IonoResult, IonoSeries } from '../src/analysis/ionosphere';
+import { computeIonoRate, detrendIonoArcs } from '../src/analysis/ionosphere';
+
+const T0 = Date.UTC(2026, 0, 1);
+
+/** TEC(t) = 20 + 5·sin(2π·t/7200): dTEC/dt peaks at ~0.26 TECU/min. */
+const tec = (tMs: number) =>
+  20 + 5 * Math.sin((2 * Math.PI * (tMs - T0)) / 7_200_000);
+const RATE_PEAK = ((5 * 2 * Math.PI) / 7200) * 60; // TECU/min
+
+function series(stepSec: number, n: number, arcStarts = [0]): IonoSeries {
+  return {
+    prn: 'G01',
+    system: 'G',
+    label: 'L1-L2',
+    codes: ['C1C', 'C2W'],
+    tecuPerNs: 2.85,
+    points: Array.from({ length: n }, (_, i) => {
+      const t = T0 + i * stepSec * 1000;
+      return { time: t, stec: tec(t) };
+    }),
+    arcStarts,
+  };
+}
+
+const asResult = (s: IonoSeries): IonoResult => ({
+  series: [s],
+  maxStec: 25,
+  meanStec: 20,
+});
+
+describe('computeIonoRate', () => {
+  it('native sequential differences track the true dTEC/dt', () => {
+    const rates = computeIonoRate(asResult(series(30, 240)));
+    const values = rates[0]!.points.map((p) => p.value);
+    const peak = Math.max(...values.map(Math.abs));
+    expect(peak).toBeGreaterThan(RATE_PEAK * 0.95);
+    expect(peak).toBeLessThan(RATE_PEAK * 1.05);
+  });
+
+  it('standardized 60 s baseline is sample-rate independent', () => {
+    const at30s = computeIonoRate(asResult(series(30, 240)), 60);
+    const at1s = computeIonoRate(asResult(series(1, 7200)), 60);
+    const peak30 = Math.max(...at30s[0]!.points.map((p) => Math.abs(p.value)));
+    const peak1 = Math.max(...at1s[0]!.points.map((p) => Math.abs(p.value)));
+    expect(peak30).toBeCloseTo(peak1, 1);
+    // ~1-min pairs from 1 s data: many more of them, same magnitude
+    expect(at1s[0]!.points.length).toBeGreaterThan(
+      at30s[0]!.points.length * 10
+    );
+  });
+
+  it('never differences across an arc boundary', () => {
+    // Two arcs with a 100 TECU step between them (a slip's signature):
+    // the step must NOT appear in the rate series.
+    const s = series(30, 100, [0, 50]);
+    for (let i = 50; i < 100; i++) s.points[i]!.stec += 100;
+    const rates = computeIonoRate(asResult(s));
+    expect(rates[0]!.points).toHaveLength(98); // 49 + 49, not 99
+    const peak = Math.max(...rates[0]!.points.map((p) => Math.abs(p.value)));
+    expect(peak).toBeLessThan(1); // the 100 TECU step would read ~200/min
+  });
+
+  it('second undivided difference is near zero for smooth TEC', () => {
+    const d2 = computeIonoRate(asResult(series(30, 240)), undefined, 2);
+    const peak = Math.max(...d2[0]!.points.map((p) => Math.abs(p.value)));
+    expect(peak).toBeLessThan(0.01); // curvature of the slow sinusoid only
+  });
+});
+
+describe('detrendIonoArcs', () => {
+  it('anchors every arc at its first observation', () => {
+    const s = series(30, 100, [0, 50]);
+    for (let i = 50; i < 100; i++) s.points[i]!.stec += 7; // arc bias
+    const out = detrendIonoArcs(asResult(s));
+    const pts = out.series[0]!.points;
+    expect(pts[0]!.stec).toBe(0);
+    expect(pts[50]!.stec).toBe(0);
+    // Variation within each arc is preserved
+    expect(pts[10]!.stec).toBeCloseTo(
+      s.points[10]!.stec - s.points[0]!.stec,
+      9
+    );
+  });
+});