source: ntrip/trunk/BNC/tests/tools/mkcrafted.py

#!/usr/bin/env python3
"""Write a BNC raw file with crafted, CRC-valid RTCM3 frames that use the
upper end of the field ranges allowed by the format (decoder robustness test).

Streams (one station ID each):
  CRAFT0IGS0  IGS-SSR VTEC (4076/201), decoded by SsrCorrIgs, followed by
              RTCM 1300/1301/1302 with 31-character names (BNC's RTCM scan
              lists them only for streams that also decode successfully);
              in the first epoch also two RTCM 1033 (antenna and receiver
              descriptors): a complete one and one whose receiver descriptor
              length (200) runs beyond the end of the message
  CRAFT0RTC0  RTCM-SSR VTEC (1264), decoded by SsrCorrRtcm
  CRAFT0NEW0  new RTCM-SSR metadata (60), not dispatched by BNC at present;
              decoded by ssr_harness.cpp with SsrCorrRtcmNew:
                - 7 entries (maximum), additional parameter of 40 and 0 bits
                - 31 entries (5-bit count, more than META_MAXNUMENTRIES=7)
              and satellite antenna messages (DF+010-DF+021):
                - BDS (83): C01, C17, C33, C49, C64 (mask bits in all four
                  16-bit quarters of DF394), individual blocks
                - GPS (80): G01, G02, G32 sharing one block and one frequency
                  block; nadir angle corrections incl. the invalid value

VTEC epochs (every 5 s):
  t0      2 layers: layer 1 degree 16/order 2, layer 2 degree 2/order 2
  t0+5    4 layers: layers 1-3 degree 1/order 1, layer 4 degree 16/order 16
  t0+10   1 layer degree 2/order 2 (decoding continues normally)
Degree and order are transmitted as N-1 in 4 bits, so 16 is the maximum
(DF474/DF475, IDF037/IDF038). Coefficient C(d,o) is (100*d + o)/200 TECU
and S(d,o) its negative, so every value in the .ssr output can be checked
by eye; the heights are 100 km * layer number.

With --rawfile, a raw file with hostile chunk headers is written instead
(test of the raw file reader, data/crafted_rawfile.raw):
  ../ESCAPE0  station ID with a path: RTCM-SSR VTEC whose .ssr file must be
              written inside the output directory
  CRAFT0RTC0  RTCM-SSR VTEC, 5 s later
  then a chunk header announcing 2000000000 bytes, which must be rejected

Usage: mkcrafted.py <out.raw>
       mkcrafted.py --rawfile <out.raw>
"""

import datetime
import os
import sys

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from bncraw import crc24q, write_chunks  # noqa: E402

T0 = datetime.datetime(2026, 10, 23, 12, 0, 0)  # GPS time, a Friday
NAME31 = "ABCDEFGHIJKLMNOPQRSTUVWXYZ01234"


class BitWriter:
    def __init__(self):
        self.bits = []

    def u(self, value, n):
        assert 0 <= value < (1 << n), (value, n)
        self.bits += [(value >> (n - 1 - i)) & 1 for i in range(n)]

    def s(self, value, n):
        assert -(1 << (n - 1)) <= value < (1 << (n - 1)), (value, n)
        self.u(value & ((1 << n) - 1), n)

    def text(self, string, nlen):
        self.u(len(string), nlen)
        for c in string.encode("ascii"):
            self.u(c, 8)

    def payload(self):
        bits = self.bits + [0] * (-len(self.bits) % 8)
        return bytes(int("".join(map(str, bits[i:i + 8])), 2)
                     for i in range(0, len(bits), 8))


def frame(payload):
    assert len(payload) <= 1023, len(payload)
    head = bytes([0xD3, len(payload) >> 8, len(payload) & 0xFF]) + payload
    return head + crc24q(head).to_bytes(3, "big")


def gps_sow(t):
    return (t.weekday() + 1) % 7 * 86400 + t.hour * 3600 + t.minute * 60 + t.second


def vtec_body(w, t, layers):
    """VTEC content after the message number(s): IDF003-IDF041 / DF385-DF478."""
    w.u(gps_sow(t), 20)      # epoch time 1 s
    w.u(2, 4)                # update interval 5 s
    w.u(0, 1)                # multiple message indicator
    w.u(1, 4)                # SSR IOD
    w.u(4711, 16)            # provider ID
    w.u(1, 4)                # solution ID
    w.u(20, 9)               # VTEC quality 1.0 TECU
    w.u(len(layers) - 1, 2)
    for n, (degree, order) in enumerate(layers, 1):
        w.u(10 * n, 8)       # height, 10 km units
        w.u(degree - 1, 4)
        w.u(order - 1, 4)
        for o in range(order + 1):
            for d in range(o, degree + 1):
                w.s(100 * d + o, 16)
        for o in range(1, order + 1):
            for d in range(o, degree + 1):
                w.s(-(100 * d + o), 16)


def vtec_igs(t, layers):
    w = BitWriter()
    w.u(4076, 12)
    w.u(1, 3)                # IGS SSR version
    w.u(201, 8)              # IGS message subtype: VTEC
    vtec_body(w, t, layers)
    return frame(w.payload())


def vtec_rtcm(t, layers):
    w = BitWriter()
    w.u(1264, 12)
    vtec_body(w, t, layers)
    return frame(w.payload())


def service_crs():
    w = BitWriter()
    w.u(1300, 12)
    w.text(NAME31, 5)
    w.u(2000, 16)            # coordinate epoch 1900 + 20.00
    return frame(w.payload())


def helmert():
    w = BitWriter()
    w.u(1301, 12)
    w.text(NAME31, 5)        # source name
    w.text(NAME31[::-1], 5)  # target name
    w.u(1, 8)                # system identification number
    w.u(0, 10)               # utilized transformation message indicator
    w.u(61000 - 44244, 16)   # MJD of t0
    for n in (23, 23, 23, 32, 32, 32, 25, 17, 17, 17, 17, 17, 17, 14):
        w.s(0, n)
    return frame(w.payload())


def rtcm_crs():
    w = BitWriter()
    w.u(1302, 12)
    w.text(NAME31, 5)
    w.u(1, 1)                # anchor: plate-fixed
    w.u(7, 5)                # plate number: Eurasia
    w.u(2, 3)                # number of database links
    w.text(NAME31, 5)
    w.text(NAME31[::-1], 5)
    return frame(w.payload())


def descriptors(strings, truncate=False):
    """RTCM 1033: antenna, setup ID, serial, receiver, firmware, serial."""
    w = BitWriter()
    w.u(1033, 12)
    w.u(0, 12)               # reference station ID
    ant, antser, rec, firm, recser = strings
    w.text(ant, 8)
    w.u(0, 8)                # antenna setup ID
    w.text(antser, 8)
    if truncate:
        w.u(200, 8)          # receiver descriptor length beyond the message
        for c in rec.encode("ascii"):
            w.u(c, 8)
    else:
        w.text(rec, 8)
        w.text(firm, 8)
        w.text(recser, 8)
    return frame(w.payload())


def metadata(entries):
    """entries: list of None (default model) or number of additional bits."""
    w = BitWriter()
    w.u(60, 12)
    w.u(1, 4)                # SSR IOD
    w.u(4711, 16)            # provider ID
    w.u(1, 4)                # solution ID
    w.u(len(entries), 5)
    for n, bits in enumerate(entries):
        w.u(n % 32, 5)       # model correction type
        w.u(1, 1)            # applied
        w.u(bits is not None, 1)
        if bits is not None:
            w.u(1, 3)        # non-default model identifier
        w.u(0, 1)            # no data IOD
        if bits is not None:
            w.u(1, 1)        # additional parameter follows
            w.u(bits, 8)
            for i in range(bits):
                w.u(i & 1, 1)
    return frame(w.payload())


def antenna(msg, sats, shared, nadir):
    """Satellite antenna message; sats: PRNs (1..64) of the mask."""
    w = BitWriter()
    w.u(msg, 12)
    w.u(4711, 16)            # provider ID
    w.u(5, 6)                # satellite antenna IOD
    w.u(1, 1)                # phase center information
    w.u(0, 1)                # group delay information
    w.u(1, 1)                # nadir angle dependent corrections follow
    w.u(3, 5)                # maximum off-nadir angle: corrections 0..3 deg
    w.u(1, 4)                # range extension: 3 + 1 = 4 bits
    w.u(shared, 1)           # satellite set indicator
    mask = 0
    for prn in sats:
        mask |= 1 << (64 - prn)
    w.u(mask, 64)
    for n in range(1 if shared else len(sats)):
        w.u(shared, 1)       # frequency set indicator
        w.u(0b101000, 6)     # frequencies 0 and 2
        for f in range(1 if shared else 2):
            w.u(1, 1)        # nadir correction follows
            w.s(100 * n + 10 * f + 1, 12)    # mm
            for deg in range(4):
                w.s(nadir[(n + f + deg) % len(nadir)], 4)   # mm, -8: invalid
    return frame(w.payload())


def rawfile(path):
    chunks = []
    for i, sta in enumerate(["../ESCAPE0", "CRAFT0RTC0"]):
        t = T0 + datetime.timedelta(seconds=5 * i)
        chunks.append((t.isoformat(), sta, "RTCM_3", vtec_rtcm(t, [(2, 2)])))
    write_chunks(path, chunks)
    t = T0 + datetime.timedelta(seconds=10)
    with open(path, "ab") as o:
        o.write(("%s CRAFT0RTC0 RTCM_3 2000000000\n" % t.isoformat()).encode())
        o.write(vtec_rtcm(t, [(2, 2)]) + b"\n")


def main():
    if len(sys.argv) == 3 and sys.argv[1] == "--rawfile":
        rawfile(sys.argv[2])
        return
    if len(sys.argv) != 2:
        sys.exit(__doc__)
    epochs = [
        [(16, 2), (2, 2)],
        [(1, 1), (1, 1), (1, 1), (16, 16)],
        [(2, 2)],
    ]
    chunks = []
    for i, layers in enumerate(epochs):
        t = T0 + datetime.timedelta(seconds=5 * i)
        stamp = t.isoformat()
        data = vtec_igs(t, layers) + service_crs() + helmert() + rtcm_crs()
        if i == 0:
            data += descriptors(["LEIAR25.R4      LEIT", "12345", "SEPT POLARX5",
                                 "5.5.0", "3001234"])
            data += descriptors(["TRUNCATED ANTENNA", "9", "RECEIVER", "", ""],
                                truncate=True)
        chunks.append((stamp, "CRAFT0IGS0", "RTCM_3", data))
        chunks.append((stamp, "CRAFT0RTC0", "RTCM_3", vtec_rtcm(t, layers)))
    stamp = (T0 + datetime.timedelta(seconds=15)).isoformat()
    chunks.append((stamp, "CRAFT0NEW0", "RTCM_3",
                   metadata([40, 0] + [None] * 5) + metadata([None] * 31)
                   + antenna(83, [1, 17, 33, 49, 64], 0, [-7, 0, 7, 3])
                   + antenna(80, [1, 2, 32], 1, [5, -8, -1, 2])))
    write_chunks(sys.argv[1], chunks)


if __name__ == "__main__":
    main()