Initial

field_audit.py

@@ -0,0 +1,235 @@
+"""Iron-clad field-format audit: firmware math vs datasheets, full wire round-trip,
+and dashboard contract check. Every assert is a proof; any failure = discrepancy."""
+import re
+import struct
+import sys
+
+sys.path.insert(0, "../stratoflights-tracker/gnuradio")
+import stratoflight_protocol as sp
+
+PASS = []
+def ok(name, detail=""):
+    PASS.append(name)
+    print(f"  ok  {name}" + (f"  ({detail})" if detail else ""))
+
+# ============================================================================
+# 1. MS5611: firmware integer math (sensors.cpp ms5611_compute, exact shifts)
+#    vs datasheet worked example (MS5611-01BA03.pdf p.7):
+#    C1..C6 = 40127,36924,23317,23282,33464,28312; D1=9085466, D2=8569150
+#    Expected: dT=2366, TEMP=2007 (=20.07 C), OFF=2420281617, SENS=1315097036,
+#    P=100009 (=1000.09 mbar = 100009 Pa)
+# ============================================================================
+def ms5611_compute_firmware(c, d1, d2):
+    # c[0]=C1 ... c[5]=C6, mirrors sensors.cpp lines 123-145 with C shift semantics
+    dt = d2 - (c[4] << 8)
+    temp = 2000 + ((dt * c[5]) >> 23)
+    off = (c[1] << 16) + ((c[3] * dt) >> 7)
+    sens = (c[0] << 15) + ((c[2] * dt) >> 8)
+    if temp < 2000:
+        t2 = (dt * dt) >> 31
+        off2 = 5 * (temp - 2000) * (temp - 2000) // 2
+        sens2 = 5 * (temp - 2000) * (temp - 2000) // 4
+        if temp < -1500:
+            off2 += 7 * (temp + 1500) * (temp + 1500)
+            sens2 += 11 * (temp + 1500) * (temp + 1500) // 2
+        temp -= t2
+        off -= off2
+        sens -= sens2
+    p = (((d1 * sens) >> 21) - off) >> 15
+    return dt, temp, off, sens, p
+
+C = [40127, 36924, 23317, 23282, 33464, 28312]
+dt, temp, off, sens, p = ms5611_compute_firmware(C, 9085466, 8569150)
+assert dt == 2366, dt
+assert temp == 2007, temp
+assert off == 2420281617, off
+assert sens == 1315097036, sens
+assert p == 100009, p
+ok("MS5611 firmware math == datasheet example", "P=100009 => 1000.09 mbar == 100009 Pa, unit is Pa exactly")
+
+# ============================================================================
+# 2. Barometric altitude formula (sensors.cpp line 154-156) sanity:
+#    ISA: 101325 Pa -> 0 m; 26436 Pa -> ~10000 m
+# ============================================================================
+alt0 = 44330.0 * (1.0 - (101325.0 / 101325.0) ** 0.1902949)
+alt10k = 44330.0 * (1.0 - (26436.0 / 101325.0) ** 0.1902949)
+assert abs(alt0) < 0.01
+assert abs(alt10k - 10000) < 30, alt10k
+ok("baro altitude formula ISA check", f"26436 Pa -> {alt10k:.0f} m (true 10000)")
+
+# ============================================================================
+# 3. SHT45 conversions (sensors.cpp lines 216-228) vs datasheet 4.6 + CRC 4.4
+#    CRC example from datasheet: CRC(0xBEEF) = 0x92
+# ============================================================================
+def sht45_crc8(msb, lsb):
+    crc = 0xFF
+    for b in (msb, lsb):
+        crc ^= b
+        for _ in range(8):
+            crc = ((crc << 1) ^ 0x31) & 0xFF if crc & 0x80 else (crc << 1) & 0xFF
+    return crc
+
+assert sht45_crc8(0xBE, 0xEF) == 0x92
+ok("SHT45 CRC-8 == datasheet example", "CRC(0xBEEF)=0x92")
+
+t_ticks = round((20.0 + 45.0) / 175.0 * 65535)   # 20 C
+rh_ticks = round((55.0 + 6.0) / 125.0 * 65535)   # 55 %RH
+t_c = -45.0 + 175.0 * t_ticks / 65535.0
+rh = -6.0 + 125.0 * rh_ticks / 65535.0
+assert abs(t_c - 20.0) < 0.01 and abs(rh - 55.0) < 0.01
+ok("SHT45 tick conversion (datasheet eq. 1,2)", f"{t_c:.3f} C, {rh:.3f} %RH")
+
+# ============================================================================
+# 4. NMEA coordinate math (gnss.cpp nmea_coord) with float32 atof as on AVR
+#    (avr-gcc double == 32-bit float)
+# ============================================================================
+import numpy as np
+def nmea_coord_f32(field, hemi):
+    v = np.float64(np.float32(field))   # atof returns 32-bit on AVR
+    deg = int(v / 100.0)
+    minutes = v - deg * 100.0
+    result = deg * 10000000 + int((minutes / 60.0) * 1e7)
+    return -result if hemi in "SW" else result
+
+lat = nmea_coord_f32("4807.038", "N")
+assert abs(lat / 1e7 - (48 + 7.038 / 60)) < 2e-6, lat   # < ~0.2 m
+lon_hi = nmea_coord_f32("12959.95695", "E")
+err_deg = abs(lon_hi / 1e7 - (129 + 59.95695 / 60))
+assert err_deg < 5e-5, err_deg  # float32 precision ceiling
+ok("NMEA ddmm.mmmm -> deg*1e7", f"lat exact; lon float32 err {err_deg*111320*0.48:.1f} m at 129.999 E")
+
+# ============================================================================
+# 5. Full wire round-trip: physical values -> beacon_t bytes (as firmware packs)
+#    -> on-air frame (PN9 whiten + CRC16) -> bit-by-bit Deframer -> parse_beacon
+#    -> webclient beacon_to_packet -> dashboard payload keys
+# ============================================================================
+phys = dict(lat=61.6820716, lon=129.8599520, gps_alt_m=15015.00, pressure_pa=11710,
+            baro_alt_m=14490.11, baro_temp_c=-9.34, humidity_pct=12.3,
+            sht_temp_c=-11.57, vbat_mv=3521, tmp_ext_c=-46.82, sats=9,
+            hh=12, mm=34, ss=56, cs=78, uptime=4321,
+            flags=sp.FLAGS[0][0:0])  # placeholder
+flags_val = 0x001B  # GPS_FIX|GPS_TIME_VALID|MS5611_OK|SHT45_OK
+
+beacon = struct.pack(
+    "<B6sIBBBBiiiIihHhHhBH",
+    sp.MSG_BEACON, b"YK0001",
+    phys["uptime"], phys["hh"], phys["mm"], phys["ss"], phys["cs"],
+    round(phys["lat"] * 1e7), round(phys["lon"] * 1e7),
+    round(phys["gps_alt_m"] * 100), phys["pressure_pa"],
+    round(phys["baro_alt_m"] * 100), round(phys["baro_temp_c"] * 100),
+    round(phys["humidity_pct"] * 10), round(phys["sht_temp_c"] * 100),
+    phys["vbat_mv"], round(phys["tmp_ext_c"] * 100), phys["sats"], flags_val,
+)
+assert len(beacon) == 48
+
+frame = sp.build_frame(beacon)          # preamble|sync|len|whitened payload|CRC16
+got = {}
+defr = sp.Deframer(on_packet=lambda pl, crc_ok: got.update(pl=pl, ok=crc_ok))
+body = frame[10:]                       # strip preamble(8)+sync(2); len byte onwards
+defr.reset()
+for byte in body:
+    for bit in range(7, -1, -1):        # MSB-first, as on air
+        defr.push_bit((byte >> bit) & 1)
+assert got["ok"], "CRC failed in deframer"
+name, fields = sp.parse_message(got["pl"])
+assert name == "BEACON"
+assert fields["callsign"] == "YK0001"
+assert abs(fields["lat"] - phys["lat"]) < 5e-8
+assert abs(fields["lon"] - phys["lon"]) < 5e-8
+assert fields["alt"] == phys["gps_alt_m"]
+assert fields["baro_alt"] == phys["baro_alt_m"]
+assert fields["pressure_pa"] == phys["pressure_pa"]
+assert fields["baro_temp_c"] == phys["baro_temp_c"]
+assert fields["humidity_pct"] == phys["humidity_pct"]
+assert fields["sht_temp_c"] == phys["sht_temp_c"]
+assert fields["vbat_mv"] == phys["vbat_mv"]
+assert fields["tmp_ext_c"] == phys["tmp_ext_c"]
+assert fields["sats"] == phys["sats"]
+assert fields["utc"] == "12:34:56.78"
+assert fields["uptime_s"] == phys["uptime"]
+assert set(fields["flag_names"]) == {"GPS_FIX", "GPS_TIME_VALID", "MS5611_OK", "SHT45_OK"}
+ok("wire round-trip beacon -> frame -> deframe -> parse", "all 15 fields bit-exact")
+
+# webclient mapping
+sys.path.insert(0, "../stratoflights-tracker/gnuradio")
+from stratoflight_webclient import beacon_to_packet
+pkt = beacon_to_packet(fields, timestamp=1783000000)
+assert pkt["lat"] == fields["lat"] and pkt["lon"] == fields["lon"] and pkt["alt"] == fields["alt"]
+payload_keys = set(pkt["payload"].keys())
+
+# dashboard contract: keys referenced by parsePacket in telemetry.ts
+ts = open("./src/lib/telemetry.ts").read()
+dash_keys = set(re.findall(r"p\.([a-z_0-9]+)", ts))
+missing = dash_keys - payload_keys
+assert not missing, f"dashboard expects keys the station never sends: {missing}"
+unused = payload_keys - dash_keys
+ok("webclient payload keys cover dashboard parser", f"unused extras: {sorted(unused) or 'none'}")
+
+# flags name contract
+ts_flags = set(re.findall(r"name: '([A-Z0-9_]+)'", ts))
+proto_flags = {n for _, n in sp.FLAGS}
+assert ts_flags == proto_flags, (ts_flags ^ proto_flags)
+ok("13 flag names identical firmware<->gnuradio<->dashboard")
+
+# ============================================================================
+# 6. TMP20: firmware linear inversion vs datasheet parabolic truth
+#    parabolic: V = -3.88e-6*T^2 - 1.15e-2*T + 1.8639  [V]
+#    firmware:  T = (1863.9 - mv) / 11.77
+# ============================================================================
+print("\nTMP20 systematic error (firmware linear vs datasheet parabolic):")
+worst = 0.0
+for t_true in range(-60, 41, 10):
+    v_mv = (-3.88e-6 * t_true**2 - 1.15e-2 * t_true + 1.8639) * 1000.0
+    t_fw = (1863.9 - v_mv) / 11.77
+    err = t_fw - t_true
+    worst = max(worst, abs(err))
+    print(f"  T={t_true:+4d} C -> V={v_mv:7.1f} mV -> firmware reads {t_fw:+7.2f} C  (err {err:+.2f})")
+print(f"  worst |err| in -60..+40: {worst:.2f} C  (fixable in-field via CMD_SET_CALIB quadratic)")
+
+# ============================================================================
+# 7. vbat scaling round-trip at the LOW_BATT threshold (board.h: 147/47, 10-bit)
+# ============================================================================
+for vbat_true in (3300, 3600, 5000):
+    v_pin = vbat_true * 47 / 147
+    counts = round(v_pin * 1023 / 2500)         # 2.5 V ref
+    v_pin_mv = counts * 2500 // 1023            # firmware integer math
+    vbat_fw = v_pin_mv * 147 / 47
+    assert abs(vbat_fw - vbat_true) < 12, (vbat_true, vbat_fw)
+print(f"\n  ok  vbat divider round-trip |err| < 12 mV (1 LSB ~ 7.6 mV after divider)")
+
+
+# ============================================================================
+# 8. JSON & encodings: the exact serialization chain
+#    station json.dumps -> Django json.loads -> Django json.dumps -> JSON.parse
+# ============================================================================
+import json
+
+# (a) full JSON round-trip of a real packet, twice (station hop + Django hop)
+wire1 = json.dumps(pkt)                          # station (webclient:106)
+assert wire1.isascii(), "ensure_ascii must keep WS frames pure ASCII"
+srv = json.loads(wire1)                          # Django consumers.py:27
+wire2 = json.dumps(srv)                          # Django broadcast :53
+browser = json.loads(wire2)                      # browser JSON.parse
+assert browser == pkt, "packet mutated across JSON hops"
+ok("JSON double round-trip station->Django->browser", "byte-level ASCII, deep-equal")
+
+# (b) strict-JSON check: browsers reject NaN/Infinity which Python allows.
+#     allow_nan=False raises if any leaf could emit them.
+json.dumps(pkt, allow_nan=False)
+ok("no NaN/Infinity possible in beacon payload", "all fields integer-derived")
+
+# (c) corrupted callsign bytes survive the chain (decode 'replace' + escaping)
+bad = bytearray(beacon)
+bad[1:7] = b"\xff\xfe\x80YK "                    # garbage callsign
+_, bad_fields = sp.parse_message(bytes(bad))
+bad_pkt = beacon_to_packet(bad_fields, timestamp=0)
+s = json.dumps(bad_pkt)
+assert s.isascii() and json.loads(s)["payload"]["callsign"] == bad_fields["callsign"]
+ok("garbage callsign bytes -> U+FFFD, still valid ASCII JSON")
+
+# (d) numeric exactness through IEEE754 doubles (JS numbers)
+assert round(json.loads(json.dumps(616820716 / 1e7)) * 1e7) == 616820716
+assert json.loads(json.dumps(4294967295)) == 4294967295   # u32 max < 2^53
+ok("lat*1e7 and u32 survive double round-trip exactly")
+
+print(f"\n=== {len(PASS)+1} proof groups passed ===")