mrtk_sbas.c

Go to the documentation of this file.

/*------------------------------------------------------------------------------
 * mrtk_sbas.c : SBAS message decoding and correction functions
 *
 * Copyright (C) 2026 H.SHIONO (MRTKLIB Project)
 * Copyright (C) 2023-2025 Cabinet Office, Japan
 * Copyright (C) 2024-2025 Lighthouse Technology & Consulting Co. Ltd.
 * Copyright (C) 2023-2025 Japan Aerospace Exploration Agency
 * Copyright (C) 2023-2025 TOSHIBA ELECTRONIC TECHNOLOGIES CORPORATION
 * Copyright (C) 2015- Mitsubishi Electric Corp.
 * Copyright (C) 2014 Geospatial Information Authority of Japan
 * Copyright (C) 2014 T.SUZUKI
 * Copyright (C) 2007-2023 T.TAKASU
 *
 * SPDX-License-Identifier: BSD-2-Clause
 *----------------------------------------------------------------------------*/
/**
 * @file mrtk_sbas.c
 * @brief MRTKLIB SBAS Module — SBAS message decoding and corrections.
 *
 * option : -DRRCENA  enable rrc correction
 *
 * references :
 *     [1] RTCA/DO-229C, Minimum operational performanc standards for global
 *         positioning system/wide area augmentation system airborne equipment,
 *         RTCA inc, November 28, 2001
 *     [2] IS-QZSS v.1.1, Quasi-Zenith Satellite System Navigation Service
 *         Interface Specification for QZSS, Japan Aerospace Exploration Agency,
 *         July 31, 2009
 */
#include "mrtklib/mrtk_sbas.h"
 
#include <math.h>
#include <stdlib.h>
#include <string.h>
 
#include "mrtklib/mrtk_atmos.h"
#include "mrtklib/mrtk_bits.h"
#include "mrtklib/mrtk_sys.h"
#include "mrtklib/mrtk_trace.h"
 
/*--- local constants (duplicated to avoid rtklib.h dependency) -------------*/
#define SYS_NONE 0x00
#define SYS_GPS 0x01
#define SYS_SBS 0x02
#define SYS_GLO 0x04
#define SYS_QZS 0x10
 
static const double PI = 3.1415926535897932;
static const double D2R = 3.1415926535897932 / 180.0;
static const double R2D = 180.0 / 3.1415926535897932;
static const double CLIGHT = 299792458.0;
static const double RE_WGS84 = 6378137.0;
 
#define P2_11 4.882812500000000E-04 /* 2^-11 */
#define P2_31 4.656612873077393E-10 /* 2^-31 */
#define P2_39 1.818989403545856E-12 /* 2^-39 */
#define MAXSBSAGEF 30.0             /* max age of SBAS fast correction (s) */
#define MAXSBSAGEL 1800.0           /* max age of SBAS long term corr (s) */
 
/*--- forward declarations for legacy functions resolved at link time -------*/
 
/* constants -----------------------------------------------------------------*/
 
#define WEEKOFFSET 1024 /* gps week offset for NovAtel OEM-3 */
 
/* sbas igp definition -------------------------------------------------------*/
static const int16_t x1[] = {-75, -65, -55, -50, -45, -40, -35, -30, -25, -20, -15, -10, -5, 0,
                             5,   10,  15,  20,  25,  30,  35,  40,  45,  50,  55,  65,  75, 85},
                     x2[] = {-55, -50, -45, -40, -35, -30, -25, -20, -15, -10, -5, 0,
                             5,   10,  15,  20,  25,  30,  35,  40,  45,  50,  55},
                     x3[] = {-75, -65, -55, -50, -45, -40, -35, -30, -25, -20, -15, -10, -5, 0,
                             5,   10,  15,  20,  25,  30,  35,  40,  45,  50,  55,  65,  75},
                     x4[] = {-85, -75, -65, -55, -50, -45, -40, -35, -30, -25, -20, -15, -10, -5,
                             0,   5,   10,  15,  20,  25,  30,  35,  40,  45,  50,  55,  65,  75},
                     x5[] = {-180, -175, -170, -165, -160, -155, -150, -145, -140, -135, -130, -125, -120, -115, -110,
                             -105, -100, -95,  -90,  -85,  -80,  -75,  -70,  -65,  -60,  -55,  -50,  -45,  -40,  -35,
                             -30,  -25,  -20,  -15,  -10,  -5,   0,    5,    10,   15,   20,   25,   30,   35,   40,
                             45,   50,   55,   60,   65,   70,   75,   80,   85,   90,   95,   100,  105,  110,  115,
                             120,  125,  130,  135,  140,  145,  150,  155,  160,  165,  170,  175},
                     x6[] = {-180, -170, -160, -150, -140, -130, -120, -110, -100, -90, -80, -70,
                             -60,  -50,  -40,  -30,  -20,  -10,  0,    10,   20,   30,  40,  50,
                             60,   70,   80,   90,   100,  110,  120,  130,  140,  150, 160, 170},
                     x7[] = {-180, -150, -120, -90, -60, -30, 0, 30, 60, 90, 120, 150},
                     x8[] = {-170, -140, -110, -80, -50, -20, 10, 40, 70, 100, 130, 160};
 
const sbsigpband_t igpband1[9][8] = {/* band 0-8 */
                                     {{-180, x1, 1, 28},
                                      {-175, x2, 29, 51},
                                      {-170, x3, 52, 78},
                                      {-165, x2, 79, 101},
                                      {-160, x3, 102, 128},
                                      {-155, x2, 129, 151},
                                      {-150, x3, 152, 178},
                                      {-145, x2, 179, 201}},
                                     {{-140, x4, 1, 28},
                                      {-135, x2, 29, 51},
                                      {-130, x3, 52, 78},
                                      {-125, x2, 79, 101},
                                      {-120, x3, 102, 128},
                                      {-115, x2, 129, 151},
                                      {-110, x3, 152, 178},
                                      {-105, x2, 179, 201}},
                                     {{-100, x3, 1, 27},
                                      {-95, x2, 28, 50},
                                      {-90, x1, 51, 78},
                                      {-85, x2, 79, 101},
                                      {-80, x3, 102, 128},
                                      {-75, x2, 129, 151},
                                      {-70, x3, 152, 178},
                                      {-65, x2, 179, 201}},
                                     {{-60, x3, 1, 27},
                                      {-55, x2, 28, 50},
                                      {-50, x4, 51, 78},
                                      {-45, x2, 79, 101},
                                      {-40, x3, 102, 128},
                                      {-35, x2, 129, 151},
                                      {-30, x3, 152, 178},
                                      {-25, x2, 179, 201}},
                                     {{-20, x3, 1, 27},
                                      {-15, x2, 28, 50},
                                      {-10, x3, 51, 77},
                                      {-5, x2, 78, 100},
                                      {0, x1, 101, 128},
                                      {5, x2, 129, 151},
                                      {10, x3, 152, 178},
                                      {15, x2, 179, 201}},
                                     {{20, x3, 1, 27},
                                      {25, x2, 28, 50},
                                      {30, x3, 51, 77},
                                      {35, x2, 78, 100},
                                      {40, x4, 101, 128},
                                      {45, x2, 129, 151},
                                      {50, x3, 152, 178},
                                      {55, x2, 179, 201}},
                                     {{60, x3, 1, 27},
                                      {65, x2, 28, 50},
                                      {70, x3, 51, 77},
                                      {75, x2, 78, 100},
                                      {80, x3, 101, 127},
                                      {85, x2, 128, 150},
                                      {90, x1, 151, 178},
                                      {95, x2, 179, 201}},
                                     {{100, x3, 1, 27},
                                      {105, x2, 28, 50},
                                      {110, x3, 51, 77},
                                      {115, x2, 78, 100},
                                      {120, x3, 101, 127},
                                      {125, x2, 128, 150},
                                      {130, x4, 151, 178},
                                      {135, x2, 179, 201}},
                                     {{140, x3, 1, 27},
                                      {145, x2, 28, 50},
                                      {150, x3, 51, 77},
                                      {155, x2, 78, 100},
                                      {160, x3, 101, 127},
                                      {165, x2, 128, 150},
                                      {170, x3, 151, 177},
                                      {175, x2, 178, 200}}};
const sbsigpband_t igpband2[2][5] =
    {/* band 9-10 */
     {{60, x5, 1, 72}, {65, x6, 73, 108}, {70, x6, 109, 144}, {75, x6, 145, 180}, {85, x7, 181, 192}},
     {{-60, x5, 1, 72}, {-65, x6, 73, 108}, {-70, x6, 109, 144}, {-75, x6, 145, 180}, {-85, x8, 181, 192}}};
/* extract field from line ---------------------------------------------------*/
static char* getfield(char* p, int pos) {
    for (pos--; pos > 0; pos--, p++) {
        if (!(p = strchr(p, ','))) {
            return NULL;
        }
    }
    return p;
}
/* variance of fast correction (udre=UDRE+1) ---------------------------------*/
static double varfcorr(int udre) {
    const double var[14] = {0.052,  0.0924, 0.1444, 0.283,  0.4678,  0.8315,   1.2992,
                            1.8709, 2.5465, 3.326,  5.1968, 20.7870, 230.9661, 2078.695};
    return 0 < udre && udre <= 14 ? var[udre - 1] : 0.0;
}
/* variance of ionosphere correction (give=GIVEI+1) --------------------------*/
static double varicorr(int give) {
    const double var[15] = {0.0084, 0.0333, 0.0749, 0.1331, 0.2079, 0.2994, 0.4075,  0.5322,
                            0.6735, 0.8315, 1.1974, 1.8709, 3.326,  20.787, 187.0826};
    return 0 < give && give <= 15 ? var[give - 1] : 0.0;
}
/* fast correction degradation -----------------------------------------------*/
static double degfcorr(int ai) {
    const double degf[16] = {0.00000, 0.00005, 0.00009, 0.00012, 0.00015, 0.00020, 0.00030, 0.00045,
                             0.00060, 0.00090, 0.00150, 0.00210, 0.00270, 0.00330, 0.00460, 0.00580};
    return 0 < ai && ai <= 15 ? degf[ai] : 0.0058;
}
/* decode type 1: prn masks --------------------------------------------------*/
static int decode_sbstype1(const sbsmsg_t* msg, sbssat_t* sbssat) {
    int i, n, sat;
 
    trace(NULL, 4, "decode_sbstype1:\n");
 
    for (i = 1, n = 0; i <= 210 && n < MAXSAT; i++) {
        if (getbitu(msg->msg, 13 + i, 1)) {
            if (i <= 37) {
                sat = satno(SYS_GPS, i); /*   0- 37: gps */
            } else if (i <= 61) {
                sat = satno(SYS_GLO, i - 37); /*  38- 61: glonass */
            } else if (i <= 119) {
                sat = 0; /*  62-119: future gnss */
            } else if (i <= 138) {
                sat = satno(SYS_SBS, i); /* 120-138: geo/waas */
            } else if (i <= 182) {
                sat = 0; /* 139-182: reserved */
            } else if (i <= 192) {
                sat = satno(SYS_SBS, i + 10); /* 183-192: qzss ref [2] */
            } else if (i <= 202) {
                sat = satno(SYS_QZS, i); /* 193-202: qzss ref [2] */
            } else {
                sat = 0; /* 203-   : reserved */
            }
            sbssat->sat[n++].sat = sat;
        }
    }
    sbssat->iodp = getbitu(msg->msg, 224, 2);
    sbssat->nsat = n;
 
    trace(NULL, 5, "decode_sbstype1: nprn=%d iodp=%d\n", n, sbssat->iodp);
    return 1;
}
/* decode type 2-5,0: fast corrections ---------------------------------------*/
static int decode_sbstype2(const sbsmsg_t* msg, sbssat_t* sbssat) {
    int i, j, iodf, type, udre;
    double prc, dt;
    gtime_t t0;
 
    trace(NULL, 4, "decode_sbstype2:\n");
 
    if (sbssat->iodp != (int)getbitu(msg->msg, 16, 2)) {
        return 0;
    }
 
    type = getbitu(msg->msg, 8, 6);
    iodf = getbitu(msg->msg, 14, 2);
 
    for (i = 0; i < 13; i++) {
        if ((j = 13 * ((type == 0 ? 2 : type) - 2) + i) >= sbssat->nsat) {
            break;
        }
        udre = getbitu(msg->msg, 174 + 4 * i, 4);
        t0 = sbssat->sat[j].fcorr.t0;
        prc = sbssat->sat[j].fcorr.prc;
        sbssat->sat[j].fcorr.t0 = gpst2time(msg->week, msg->tow);
        sbssat->sat[j].fcorr.prc = getbits(msg->msg, 18 + i * 12, 12) * 0.125f;
        sbssat->sat[j].fcorr.udre = udre + 1;
        dt = timediff(sbssat->sat[j].fcorr.t0, t0);
        if (t0.time == 0 || dt <= 0.0 || 18.0 < dt || sbssat->sat[j].fcorr.ai == 0) {
            sbssat->sat[j].fcorr.rrc = 0.0;
            sbssat->sat[j].fcorr.dt = 0.0;
        } else {
            sbssat->sat[j].fcorr.rrc = (sbssat->sat[j].fcorr.prc - prc) / dt;
            sbssat->sat[j].fcorr.dt = dt;
        }
        sbssat->sat[j].fcorr.iodf = iodf;
    }
    trace(NULL, 5, "decode_sbstype2: type=%d iodf=%d\n", type, iodf);
    return 1;
}
/* decode type 6: integrity info ---------------------------------------------*/
static int decode_sbstype6(const sbsmsg_t* msg, sbssat_t* sbssat) {
    int i, iodf[4], udre;
 
    trace(NULL, 4, "decode_sbstype6:\n");
 
    for (i = 0; i < 4; i++) {
        iodf[i] = getbitu(msg->msg, 14 + i * 2, 2);
    }
    for (i = 0; i < sbssat->nsat && i < MAXSAT; i++) {
        if (sbssat->sat[i].fcorr.iodf != iodf[i / 13]) {
            continue;
        }
        udre = getbitu(msg->msg, 22 + i * 4, 4);
        sbssat->sat[i].fcorr.udre = udre + 1;
    }
    trace(NULL, 5, "decode_sbstype6: iodf=%d %d %d %d\n", iodf[0], iodf[1], iodf[2], iodf[3]);
    return 1;
}
/* decode type 7: fast correction degradation factor -------------------------*/
static int decode_sbstype7(const sbsmsg_t* msg, sbssat_t* sbssat) {
    int i;
 
    trace(NULL, 4, "decode_sbstype7\n");
 
    if (sbssat->iodp != (int)getbitu(msg->msg, 18, 2)) {
        return 0;
    }
 
    sbssat->tlat = getbitu(msg->msg, 14, 4);
 
    for (i = 0; i < sbssat->nsat && i < MAXSAT; i++) {
        sbssat->sat[i].fcorr.ai = getbitu(msg->msg, 22 + i * 4, 4);
    }
    return 1;
}
/* decode type 9: geo navigation message -------------------------------------*/
static int decode_sbstype9(const sbsmsg_t* msg, nav_t* nav) {
    seph_t seph = {0};
    int i, sat, t;
 
    trace(NULL, 4, "decode_sbstype9:\n");
 
    if (!(sat = satno(SYS_SBS, msg->prn))) {
        trace(NULL, 2, "invalid prn in sbas type 9: prn=%3d\n", msg->prn);
        return 0;
    }
    t = (int)getbitu(msg->msg, 22, 13) * 16 - (int)msg->tow % 86400;
    if (t <= -43200) {
        t += 86400;
    } else if (t > 43200) {
        t -= 86400;
    }
    seph.sat = sat;
    seph.t0 = gpst2time(msg->week, msg->tow + t);
    seph.tof = gpst2time(msg->week, msg->tow);
    seph.sva = getbitu(msg->msg, 35, 4);
    seph.svh = seph.sva == 15 ? 1 : 0; /* unhealthy if ura==15 */
 
    seph.pos[0] = getbits(msg->msg, 39, 30) * 0.08;
    seph.pos[1] = getbits(msg->msg, 69, 30) * 0.08;
    seph.pos[2] = getbits(msg->msg, 99, 25) * 0.4;
    seph.vel[0] = getbits(msg->msg, 124, 17) * 0.000625;
    seph.vel[1] = getbits(msg->msg, 141, 17) * 0.000625;
    seph.vel[2] = getbits(msg->msg, 158, 18) * 0.004;
    seph.acc[0] = getbits(msg->msg, 176, 10) * 0.0000125;
    seph.acc[1] = getbits(msg->msg, 186, 10) * 0.0000125;
    seph.acc[2] = getbits(msg->msg, 196, 10) * 0.0000625;
 
    seph.af0 = getbits(msg->msg, 206, 12) * P2_31;
    seph.af1 = getbits(msg->msg, 218, 8) * P2_39 / 2.0;
 
    i = msg->prn - MINPRNSBS;
    if (!nav->seph || fabs(timediff(nav->seph[i].t0, seph.t0)) < 1E-3) { /* not change */
        return 0;
    }
    nav->seph[NSATSBS + i] = nav->seph[i]; /* previous */
    nav->seph[i] = seph;                   /* current */
 
    trace(NULL, 5, "decode_sbstype9: prn=%d\n", msg->prn);
    return 1;
}
/* decode type 18: ionospheric grid point masks ------------------------------*/
static int decode_sbstype18(const sbsmsg_t* msg, sbsion_t* sbsion) {
    const sbsigpband_t* p;
    int i, j, n, m, band = getbitu(msg->msg, 18, 4);
 
    trace(NULL, 4, "decode_sbstype18:\n");
 
    if (0 <= band && band <= 8) {
        p = igpband1[band];
        m = 8;
    } else if (9 <= band && band <= 10) {
        p = igpband2[band - 9];
        m = 5;
    } else {
        return 0;
    }
 
    sbsion[band].iodi = (int16_t)getbitu(msg->msg, 22, 2);
 
    for (i = 1, n = 0; i <= 201; i++) {
        if (!getbitu(msg->msg, 23 + i, 1)) {
            continue;
        }
        for (j = 0; j < m; j++) {
            if (i < p[j].bits || p[j].bite < i) {
                continue;
            }
            sbsion[band].igp[n].lat = band <= 8 ? p[j].y[i - p[j].bits] : p[j].x;
            sbsion[band].igp[n++].lon = band <= 8 ? p[j].x : p[j].y[i - p[j].bits];
            break;
        }
    }
    sbsion[band].nigp = n;
 
    trace(NULL, 5, "decode_sbstype18: band=%d nigp=%d\n", band, n);
    return 1;
}
/* decode half long term correction (vel code=0) -----------------------------*/
static int decode_longcorr0(const sbsmsg_t* msg, int p, sbssat_t* sbssat) {
    int i, n = getbitu(msg->msg, p, 6);
 
    trace(NULL, 4, "decode_longcorr0:\n");
 
    if (n == 0 || n > MAXSAT) {
        return 0;
    }
 
    sbssat->sat[n - 1].lcorr.iode = getbitu(msg->msg, p + 6, 8);
 
    for (i = 0; i < 3; i++) {
        sbssat->sat[n - 1].lcorr.dpos[i] = getbits(msg->msg, p + 14 + 9 * i, 9) * 0.125;
        sbssat->sat[n - 1].lcorr.dvel[i] = 0.0;
    }
    sbssat->sat[n - 1].lcorr.daf0 = getbits(msg->msg, p + 41, 10) * P2_31;
    sbssat->sat[n - 1].lcorr.daf1 = 0.0;
    sbssat->sat[n - 1].lcorr.t0 = gpst2time(msg->week, msg->tow);
 
    trace(NULL, 5, "decode_longcorr0:sat=%2d\n", sbssat->sat[n - 1].sat);
    return 1;
}
/* decode half long term correction (vel code=1) -----------------------------*/
static int decode_longcorr1(const sbsmsg_t* msg, int p, sbssat_t* sbssat) {
    int i, n = getbitu(msg->msg, p, 6), t;
 
    trace(NULL, 4, "decode_longcorr1:\n");
 
    if (n == 0 || n > MAXSAT) {
        return 0;
    }
 
    sbssat->sat[n - 1].lcorr.iode = getbitu(msg->msg, p + 6, 8);
 
    for (i = 0; i < 3; i++) {
        sbssat->sat[n - 1].lcorr.dpos[i] = getbits(msg->msg, p + 14 + i * 11, 11) * 0.125;
        sbssat->sat[n - 1].lcorr.dvel[i] = getbits(msg->msg, p + 58 + i * 8, 8) * P2_11;
    }
    sbssat->sat[n - 1].lcorr.daf0 = getbits(msg->msg, p + 47, 11) * P2_31;
    sbssat->sat[n - 1].lcorr.daf1 = getbits(msg->msg, p + 82, 8) * P2_39;
    t = (int)getbitu(msg->msg, p + 90, 13) * 16 - (int)msg->tow % 86400;
    if (t <= -43200) {
        t += 86400;
    } else if (t > 43200) {
        t -= 86400;
    }
    sbssat->sat[n - 1].lcorr.t0 = gpst2time(msg->week, msg->tow + t);
 
    trace(NULL, 5, "decode_longcorr1: sat=%2d\n", sbssat->sat[n - 1].sat);
    return 1;
}
/* decode half long term correction ------------------------------------------*/
static int decode_longcorrh(const sbsmsg_t* msg, int p, sbssat_t* sbssat) {
    trace(NULL, 4, "decode_longcorrh:\n");
 
    if (getbitu(msg->msg, p, 1) == 0) { /* vel code=0 */
        if (sbssat->iodp == (int)getbitu(msg->msg, p + 103, 2)) {
            return decode_longcorr0(msg, p + 1, sbssat) && decode_longcorr0(msg, p + 52, sbssat);
        }
    } else if (sbssat->iodp == (int)getbitu(msg->msg, p + 104, 2)) {
        return decode_longcorr1(msg, p + 1, sbssat);
    }
    return 0;
}
/* decode type 24: mixed fast/long term correction ---------------------------*/
static int decode_sbstype24(const sbsmsg_t* msg, sbssat_t* sbssat) {
    int i, j, iodf, blk, udre;
 
    trace(NULL, 4, "decode_sbstype24:\n");
 
    if (sbssat->iodp != (int)getbitu(msg->msg, 110, 2)) {
        return 0; /* check IODP */
    }
 
    blk = getbitu(msg->msg, 112, 2);
    iodf = getbitu(msg->msg, 114, 2);
 
    for (i = 0; i < 6; i++) {
        if ((j = 13 * blk + i) >= sbssat->nsat) {
            break;
        }
        udre = getbitu(msg->msg, 86 + 4 * i, 4);
 
        sbssat->sat[j].fcorr.t0 = gpst2time(msg->week, msg->tow);
        sbssat->sat[j].fcorr.prc = getbits(msg->msg, 14 + i * 12, 12) * 0.125f;
        sbssat->sat[j].fcorr.udre = udre + 1;
        sbssat->sat[j].fcorr.iodf = iodf;
    }
    return decode_longcorrh(msg, 120, sbssat);
}
/* decode type 25: long term satellite error correction ----------------------*/
static int decode_sbstype25(const sbsmsg_t* msg, sbssat_t* sbssat) {
    trace(NULL, 4, "decode_sbstype25:\n");
 
    return decode_longcorrh(msg, 14, sbssat) && decode_longcorrh(msg, 120, sbssat);
}
/* decode type 26: ionospheric delay corrections -----------------------------*/
static int decode_sbstype26(const sbsmsg_t* msg, sbsion_t* sbsion) {
    int i, j, block, delay, give, band = getbitu(msg->msg, 14, 4);
 
    trace(NULL, 4, "decode_sbstype26:\n");
 
    if (band > MAXBAND || sbsion[band].iodi != (int)getbitu(msg->msg, 217, 2)) {
        return 0;
    }
 
    block = getbitu(msg->msg, 18, 4);
 
    for (i = 0; i < 15; i++) {
        if ((j = block * 15 + i) >= sbsion[band].nigp) {
            continue;
        }
        give = getbitu(msg->msg, 22 + i * 13 + 9, 4);
 
        delay = getbitu(msg->msg, 22 + i * 13, 9);
        sbsion[band].igp[j].t0 = gpst2time(msg->week, msg->tow);
        sbsion[band].igp[j].delay = delay == 0x1FF ? 0.0f : delay * 0.125f;
        sbsion[band].igp[j].give = give + 1;
 
        if (sbsion[band].igp[j].give >= 16) {
            sbsion[band].igp[j].give = 0;
        }
    }
    trace(NULL, 5, "decode_sbstype26: band=%d block=%d\n", band, block);
    return 1;
}
/* update sbas corrections -----------------------------------------------------
 * update sbas correction parameters in navigation data with a sbas message
 * args   : sbsmg_t  *msg    I   sbas message
 *          nav_t    *nav    IO  navigation data
 * return : message type (-1: error or not supported type)
 * notes  : nav->seph must point to seph[NSATSBS*2] (array of seph_t)
 *               seph[prn-MINPRNSBS+1]          : sat prn current epehmeris
 *               seph[prn-MINPRNSBS+1+MAXPRNSBS]: sat prn previous epehmeris
 *-----------------------------------------------------------------------------*/
int sbsupdatecorr(const sbsmsg_t* msg, nav_t* nav) {
    int type = getbitu(msg->msg, 8, 6), stat = -1;
 
    trace(NULL, 3, "sbsupdatecorr: type=%d\n", type);
 
    if (msg->week == 0) {
        return -1;
    }
 
    switch (type) {
        case 0:
            stat = decode_sbstype2(msg, &nav->sbssat);
            break;
        case 1:
            stat = decode_sbstype1(msg, &nav->sbssat);
            break;
        case 2:
        case 3:
        case 4:
        case 5:
            stat = decode_sbstype2(msg, &nav->sbssat);
            break;
        case 6:
            stat = decode_sbstype6(msg, &nav->sbssat);
            break;
        case 7:
            stat = decode_sbstype7(msg, &nav->sbssat);
            break;
        case 9:
            stat = decode_sbstype9(msg, nav);
            break;
        case 18:
            stat = decode_sbstype18(msg, nav->sbsion);
            break;
        case 24:
            stat = decode_sbstype24(msg, &nav->sbssat);
            break;
        case 25:
            stat = decode_sbstype25(msg, &nav->sbssat);
            break;
        case 26:
            stat = decode_sbstype26(msg, nav->sbsion);
            break;
        case 63:
            break; /* null message */
 
            /*default: trace(NULL,2,"unsupported sbas message: type=%d\n",type); break;*/
    }
    return stat ? type : -1;
}
/* read sbas log file --------------------------------------------------------*/
static void readmsgs(const char* file, int sel, gtime_t ts, gtime_t te, sbs_t* sbs) {
    sbsmsg_t* sbs_msgs;
    int i, week, prn, ch, msg;
    uint32_t b;
    double tow, ep[6] = {0};
    char buff[256], *p;
    gtime_t time;
    FILE* fp;
 
    trace(NULL, 3, "readmsgs: file=%s sel=%d\n", file, sel);
 
    if (!(fp = fopen(file, "r"))) {
        trace(NULL, 2, "sbas message file open error: %s\n", file);
        return;
    }
    while (fgets(buff, sizeof(buff), fp)) {
        if (sscanf(buff, "%d %lf %d", &week, &tow, &prn) == 3 && (p = strstr(buff, ": "))) {
            p += 2; /* rtklib form */
        } else if (sscanf(buff, "%d %lf %lf %lf %lf %lf %lf %d", &prn, ep, ep + 1, ep + 2, ep + 3, ep + 4, ep + 5,
                          &msg) == 8) {
            /* ems (EGNOS Message Service) form */
            ep[0] += ep[0] < 70.0 ? 2000.0 : 1900.0;
            tow = time2gpst(epoch2time(ep), &week);
            p = buff + (msg >= 10 ? 25 : 24);
        } else if (!strncmp(buff, "#RAWWAASFRAMEA", 14)) { /* NovAtel OEM4/V */
            if (!(p = getfield(buff, 6))) {
                continue;
            }
            if (sscanf(p, "%d,%lf", &week, &tow) < 2) {
                continue;
            }
            if (!(p = strchr(++p, ';'))) {
                continue;
            }
            if (sscanf(++p, "%d,%d", &ch, &prn) < 2) {
                continue;
            }
            if (!(p = getfield(p, 4))) {
                continue;
            }
        } else if (!strncmp(buff, "$FRMA", 5)) { /* NovAtel OEM3 */
            if (!(p = getfield(buff, 2))) {
                continue;
            }
            if (sscanf(p, "%d,%lf,%d", &week, &tow, &prn) < 3) {
                continue;
            }
            if (!(p = getfield(p, 6))) {
                continue;
            }
            if (week < WEEKOFFSET) {
                week += WEEKOFFSET;
            }
        } else {
            continue;
        }
 
        if (sel != 0 && sel != prn) {
            continue;
        }
 
        time = gpst2time(week, tow);
 
        if (!screent(time, ts, te, 0.0)) {
            continue;
        }
 
        if (sbs->n >= sbs->nmax) {
            sbs->nmax = sbs->nmax == 0 ? 1024 : sbs->nmax * 2;
            if (!(sbs_msgs = (sbsmsg_t*)realloc(sbs->msgs, sbs->nmax * sizeof(sbsmsg_t)))) {
                trace(NULL, 1, "readsbsmsg malloc error: nmax=%d\n", sbs->nmax);
                free(sbs->msgs);
                sbs->msgs = NULL;
                sbs->n = sbs->nmax = 0;
                return;
            }
            sbs->msgs = sbs_msgs;
        }
        sbs->msgs[sbs->n].week = week;
        sbs->msgs[sbs->n].tow = (int)(tow + 0.5);
        sbs->msgs[sbs->n].prn = prn;
        for (i = 0; i < 29; i++) {
            sbs->msgs[sbs->n].msg[i] = 0;
        }
        for (i = 0; *(p - 1) && *p && i < 29; p += 2, i++) {
            if (sscanf(p, "%2X", &b) == 1) {
                sbs->msgs[sbs->n].msg[i] = (uint8_t)b;
            }
        }
        sbs->msgs[sbs->n++].msg[28] &= 0xC0;
    }
    fclose(fp);
}
/* compare sbas messages -----------------------------------------------------*/
static int cmpmsgs(const void* p1, const void* p2) {
    sbsmsg_t *q1 = (sbsmsg_t*)p1, *q2 = (sbsmsg_t*)p2;
    return q1->week != q2->week ? q1->week - q2->week
                                : (q1->tow < q2->tow ? -1 : (q1->tow > q2->tow ? 1 : q1->prn - q2->prn));
}
/* read sbas message file ------------------------------------------------------
 * read sbas message file
 * args   : char     *file   I   sbas message file (wind-card * is expanded)
 *          int      sel     I   sbas satellite prn number selection (0:all)
 *         (gtime_t  ts      I   start time)
 *         (gtime_t  te      I   end time  )
 *          sbs_t    *sbs    IO  sbas messages
 * return : number of sbas messages
 * notes  : sbas message are appended and sorted. before calling the funciton,
 *          sbs->n, sbs->nmax and sbs->msgs must be set properly. (initially
 *          sbs->n=sbs->nmax=0, sbs->msgs=NULL)
 *          only the following file extentions after wild card expanded are valid
 *          to read. others are skipped
 *          .sbs, .SBS, .ems, .EMS
 *-----------------------------------------------------------------------------*/
int sbsreadmsgt(const char* file, int sel, gtime_t ts, gtime_t te, sbs_t* sbs) {
    char *efiles[MAXEXFILE] = {0}, *ext;
    int i, n;
 
    trace(NULL, 3, "sbsreadmsgt: file=%s sel=%d\n", file, sel);
 
    for (i = 0; i < MAXEXFILE; i++) {
        if (!(efiles[i] = (char*)malloc(1024))) {
            for (i--; i >= 0; i--) {
                free(efiles[i]);
            }
            return 0;
        }
    }
    /* expand wild card in file path */
    n = expath(file, efiles, MAXEXFILE);
 
    for (i = 0; i < n; i++) {
        if (!(ext = strrchr(efiles[i], '.'))) {
            continue;
        }
        if (strcmp(ext, ".sbs") && strcmp(ext, ".SBS") && strcmp(ext, ".ems") && strcmp(ext, ".EMS")) {
            continue;
        }
 
        readmsgs(efiles[i], sel, ts, te, sbs);
    }
    for (i = 0; i < MAXEXFILE; i++) {
        free(efiles[i]);
    }
 
    /* sort messages */
    if (sbs->n > 0) {
        qsort(sbs->msgs, sbs->n, sizeof(sbsmsg_t), cmpmsgs);
    }
    return sbs->n;
}
int sbsreadmsg(const char* file, int sel, sbs_t* sbs) {
    gtime_t ts = {0}, te = {0};
 
    trace(NULL, 3, "sbsreadmsg: file=%s sel=%d\n", file, sel);
 
    return sbsreadmsgt(file, sel, ts, te, sbs);
}
/* output sbas messages --------------------------------------------------------
 * output sbas message record to output file in rtklib sbas log format
 * args   : FILE   *fp       I   output file pointer
 *          sbsmsg_t *sbsmsg I   sbas messages
 * return : none
 *-----------------------------------------------------------------------------*/
void sbsoutmsg(FILE* fp, sbsmsg_t* sbsmsg) {
    int i, prn = sbsmsg->prn, type = sbsmsg->msg[1] >> 2;
 
    trace(NULL, 4, "sbsoutmsg:\n");
 
    fprintf(fp, "%4d %6d %3d %2d : ", sbsmsg->week, sbsmsg->tow, prn, type);
    for (i = 0; i < 29; i++) {
        fprintf(fp, "%02X", sbsmsg->msg[i]);
    }
    fprintf(fp, "\n");
}
/* search igps ---------------------------------------------------------------*/
static void searchigp(gtime_t time, const double* pos, const sbsion_t* ion, const sbsigp_t** igp, double* x,
                      double* y) {
    int i, latp[2], lonp[4];
    double lat = pos[0] * R2D, lon = pos[1] * R2D;
    const sbsigp_t* p;
 
    trace(NULL, 4, "searchigp: pos=%.3f %.3f\n", pos[0] * R2D, pos[1] * R2D);
 
    if (lon >= 180.0) {
        lon -= 360.0;
    }
    if (-55.0 <= lat && lat < 55.0) {
        latp[0] = (int)floor(lat / 5.0) * 5;
        latp[1] = latp[0] + 5;
        lonp[0] = lonp[1] = (int)floor(lon / 5.0) * 5;
        lonp[2] = lonp[3] = lonp[0] + 5;
        *x = (lon - lonp[0]) / 5.0;
        *y = (lat - latp[0]) / 5.0;
    } else {
        latp[0] = (int)floor((lat - 5.0) / 10.0) * 10 + 5;
        latp[1] = latp[0] + 10;
        lonp[0] = lonp[1] = (int)floor(lon / 10.0) * 10;
        lonp[2] = lonp[3] = lonp[0] + 10;
        *x = (lon - lonp[0]) / 10.0;
        *y = (lat - latp[0]) / 10.0;
        if (75.0 <= lat && lat < 85.0) {
            lonp[1] = (int)floor(lon / 90.0) * 90;
            lonp[3] = lonp[1] + 90;
        } else if (-85.0 <= lat && lat < -75.0) {
            lonp[0] = (int)floor((lon - 50.0) / 90.0) * 90 + 40;
            lonp[2] = lonp[0] + 90;
        } else if (lat >= 85.0) {
            for (i = 0; i < 4; i++) {
                lonp[i] = (int)floor(lon / 90.0) * 90;
            }
        } else if (lat < -85.0) {
            for (i = 0; i < 4; i++) {
                lonp[i] = (int)floor((lon - 50.0) / 90.0) * 90 + 40;
            }
        }
    }
    for (i = 0; i < 4; i++) {
        if (lonp[i] == 180) {
            lonp[i] = -180;
        }
    }
    for (i = 0; i <= MAXBAND; i++) {
        for (p = ion[i].igp; p < ion[i].igp + ion[i].nigp; p++) {
            if (p->t0.time == 0) {
                continue;
            }
            if (p->lat == latp[0] && p->lon == lonp[0] && p->give > 0) {
                igp[0] = p;
            } else if (p->lat == latp[1] && p->lon == lonp[1] && p->give > 0) {
                igp[1] = p;
            } else if (p->lat == latp[0] && p->lon == lonp[2] && p->give > 0) {
                igp[2] = p;
            } else if (p->lat == latp[1] && p->lon == lonp[3] && p->give > 0) {
                igp[3] = p;
            }
            if (igp[0] && igp[1] && igp[2] && igp[3]) {
                return;
            }
        }
    }
}
/* sbas ionospheric delay correction -------------------------------------------
 * compute sbas ionosphric delay correction
 * args   : gtime_t  time    I   time
 *          nav_t    *nav    I   navigation data
 *          double   *pos    I   receiver position {lat,lon,height} (rad/m)
 *          double   *azel   I   satellite azimuth/elavation angle (rad)
 *          double   *delay  O   slant ionospheric delay (L1) (m)
 *          double   *var    O   variance of ionospheric delay (m^2)
 * return : status (1:ok, 0:no correction)
 * notes  : before calling the function, sbas ionosphere correction parameters
 *          in navigation data (nav->sbsion) must be set by callig
 *          sbsupdatecorr()
 *-----------------------------------------------------------------------------*/
int sbsioncorr(gtime_t time, const nav_t* nav, const double* pos, const double* azel, double* delay, double* var) {
    const double re = 6378.1363, hion = 350.0;
    int i, err = 0;
    double fp, posp[2], x = 0.0, y = 0.0, t, w[4] = {0};
    const sbsigp_t* igp[4] = {0}; /* {ws,wn,es,en} */
 
    trace(NULL, 4, "sbsioncorr: pos=%.3f %.3f azel=%.3f %.3f\n", pos[0] * R2D, pos[1] * R2D, azel[0] * R2D,
          azel[1] * R2D);
 
    *delay = *var = 0.0;
    if (pos[2] < -100.0 || azel[1] <= 0) {
        return 1;
    }
 
    /* ipp (ionospheric pierce point) position */
    fp = ionppp(pos, azel, re, hion, posp);
 
    /* search igps around ipp */
    searchigp(time, posp, nav->sbsion, igp, &x, &y);
 
    /* weight of igps */
    if (igp[0] && igp[1] && igp[2] && igp[3]) {
        w[0] = (1.0 - x) * (1.0 - y);
        w[1] = (1.0 - x) * y;
        w[2] = x * (1.0 - y);
        w[3] = x * y;
    } else if (igp[0] && igp[1] && igp[2]) {
        w[1] = y;
        w[2] = x;
        if ((w[0] = 1.0 - w[1] - w[2]) < 0.0) {
            err = 1;
        }
    } else if (igp[0] && igp[2] && igp[3]) {
        w[0] = 1.0 - x;
        w[3] = y;
        if ((w[2] = 1.0 - w[0] - w[3]) < 0.0) {
            err = 1;
        }
    } else if (igp[0] && igp[1] && igp[3]) {
        w[0] = 1.0 - y;
        w[3] = x;
        if ((w[1] = 1.0 - w[0] - w[3]) < 0.0) {
            err = 1;
        }
    } else if (igp[1] && igp[2] && igp[3]) {
        w[1] = 1.0 - x;
        w[2] = 1.0 - y;
        if ((w[3] = 1.0 - w[1] - w[2]) < 0.0) {
            err = 1;
        }
    } else {
        err = 1;
    }
 
    if (err) {
        trace(NULL, 2, "no sbas iono correction: lat=%3.0f lon=%4.0f\n", posp[0] * R2D, posp[1] * R2D);
        return 0;
    }
    for (i = 0; i < 4; i++) {
        if (!igp[i]) {
            continue;
        }
        t = timediff(time, igp[i]->t0);
        *delay += w[i] * igp[i]->delay;
        *var += w[i] * varicorr(igp[i]->give) * 9E-8 * fabs(t);
    }
    *delay *= fp;
    *var *= fp * fp;
 
    trace(NULL, 5, "sbsioncorr: dion=%7.2f sig=%7.2f\n", *delay, sqrt(*var));
    return 1;
}
/* get meterological parameters ----------------------------------------------*/
static void getmet(double lat, double* met) {
    static const double metprm[][10] = {/* lat=15,30,45,60,75 */
                                        {1013.25, 299.65, 26.31, 6.30E-3, 2.77, 0.00, 0.00, 0.00, 0.00E-3, 0.00},
                                        {1017.25, 294.15, 21.79, 6.05E-3, 3.15, -3.75, 7.00, 8.85, 0.25E-3, 0.33},
                                        {1015.75, 283.15, 11.66, 5.58E-3, 2.57, -2.25, 11.00, 7.24, 0.32E-3, 0.46},
                                        {1011.75, 272.15, 6.78, 5.39E-3, 1.81, -1.75, 15.00, 5.36, 0.81E-3, 0.74},
                                        {1013.00, 263.65, 4.11, 4.53E-3, 1.55, -0.50, 14.50, 3.39, 0.62E-3, 0.30}};
    int i, j;
    double a;
    lat = fabs(lat);
    if (lat <= 15.0) {
        for (i = 0; i < 10; i++) {
            met[i] = metprm[0][i];
        }
    } else if (lat >= 75.0) {
        for (i = 0; i < 10; i++) {
            met[i] = metprm[4][i];
        }
    } else {
        j = (int)(lat / 15.0);
        a = (lat - j * 15.0) / 15.0;
        for (i = 0; i < 10; i++) {
            met[i] = (1.0 - a) * metprm[j - 1][i] + a * metprm[j][i];
        }
    }
}
/* tropospheric delay correction -----------------------------------------------
 * compute sbas tropospheric delay correction (mops model)
 * args   : gtime_t time     I   time
 *          double   *pos    I   receiver position {lat,lon,height} (rad/m)
 *          double   *azel   I   satellite azimuth/elavation (rad)
 *          double   *var    O   variance of troposphric error (m^2)
 * return : slant tropospheric delay (m)
 *-----------------------------------------------------------------------------*/
double sbstropcorr(gtime_t time, const double* pos, const double* azel, double* var) {
    const double k1 = 77.604, k2 = 382000.0, rd = 287.054, gm = 9.784, g = 9.80665;
    static double pos_[3] = {0}, zh = 0.0, zw = 0.0;
    int i;
    double c, met[10], sinel = sin(azel[1]), h = pos[2], m;
 
    trace(NULL, 4, "sbstropcorr: pos=%.3f %.3f azel=%.3f %.3f\n", pos[0] * R2D, pos[1] * R2D, azel[0] * R2D,
          azel[1] * R2D);
 
    if (pos[2] < -100.0 || 10000.0 < pos[2] || azel[1] <= 0) {
        *var = 0.0;
        return 0.0;
    }
    if (zh == 0.0 || fabs(pos[0] - pos_[0]) > 1E-7 || fabs(pos[1] - pos_[1]) > 1E-7 || fabs(pos[2] - pos_[2]) > 1.0) {
        getmet(pos[0] * R2D, met);
        c = cos(2.0 * PI * (time2doy(time) - (pos[0] >= 0.0 ? 28.0 : 211.0)) / 365.25);
        for (i = 0; i < 5; i++) {
            met[i] -= met[i + 5] * c;
        }
        zh = 1E-6 * k1 * rd * met[0] / gm;
        zw = 1E-6 * k2 * rd / (gm * (met[4] + 1.0) - met[3] * rd) * met[2] / met[1];
        zh *= pow(1.0 - met[3] * h / met[1], g / (rd * met[3]));
        zw *= pow(1.0 - met[3] * h / met[1], (met[4] + 1.0) * g / (rd * met[3]) - 1.0);
        for (i = 0; i < 3; i++) {
            pos_[i] = pos[i];
        }
    }
    m = 1.001 / sqrt(0.002001 + sinel * sinel);
    *var = 0.12 * 0.12 * m * m;
    return (zh + zw) * m;
}
/* long term correction ------------------------------------------------------*/
static int sbslongcorr(gtime_t time, int sat, const sbssat_t* sbssat, double* drs, double* ddts) {
    const sbssatp_t* p;
    double t;
    int i;
 
    trace(NULL, 3, "sbslongcorr: sat=%2d\n", sat);
 
    for (p = sbssat->sat; p < sbssat->sat + sbssat->nsat; p++) {
        if (p->sat != sat || p->lcorr.t0.time == 0) {
            continue;
        }
        t = timediff(time, p->lcorr.t0);
        if (fabs(t) > MAXSBSAGEL) {
            trace(NULL, 2, "sbas long-term correction expired: %s sat=%2d t=%5.0f\n", time_str(time, 0), sat, t);
            return 0;
        }
        for (i = 0; i < 3; i++) {
            drs[i] = p->lcorr.dpos[i] + p->lcorr.dvel[i] * t;
        }
        *ddts = p->lcorr.daf0 + p->lcorr.daf1 * t;
 
        trace(NULL, 5, "sbslongcorr: sat=%2d drs=%7.2f%7.2f%7.2f ddts=%7.2f\n", sat, drs[0], drs[1], drs[2],
              *ddts * CLIGHT);
 
        return 1;
    }
    /* if sbas satellite without correction, no correction applied */
    if (satsys(sat, NULL) == SYS_SBS) {
        return 1;
    }
 
    trace(NULL, 2, "no sbas long-term correction: %s sat=%2d\n", time_str(time, 0), sat);
    return 0;
}
/* fast correction -----------------------------------------------------------*/
static int sbsfastcorr(gtime_t time, int sat, const sbssat_t* sbssat, double* prc, double* var) {
    const sbssatp_t* p;
    double t;
 
    trace(NULL, 3, "sbsfastcorr: sat=%2d\n", sat);
 
    for (p = sbssat->sat; p < sbssat->sat + sbssat->nsat; p++) {
        if (p->sat != sat) {
            continue;
        }
        if (p->fcorr.t0.time == 0) {
            break;
        }
        t = timediff(time, p->fcorr.t0) + sbssat->tlat;
 
        /* expire age of correction or UDRE==14 (not monitored) */
        if (fabs(t) > MAXSBSAGEF || p->fcorr.udre >= 15) {
            continue;
        }
        *prc = p->fcorr.prc;
#ifdef RRCENA
        if (p->fcorr.ai > 0 && fabs(t) <= 8.0 * p->fcorr.dt) {
            *prc += p->fcorr.rrc * t;
        }
#endif
        *var = varfcorr(p->fcorr.udre) + degfcorr(p->fcorr.ai) * t * t / 2.0;
 
        trace(NULL, 5, "sbsfastcorr: sat=%3d prc=%7.2f sig=%7.2f t=%5.0f\n", sat, *prc, sqrt(*var), t);
        return 1;
    }
    trace(NULL, 2, "no sbas fast correction: %s sat=%2d\n", time_str(time, 0), sat);
    return 0;
}
/* sbas satellite ephemeris and clock correction -------------------------------
 * correct satellite position and clock bias with sbas satellite corrections
 * args   : gtime_t time     I   reception time
 *          int    sat       I   satellite
 *          nav_t  *nav      I   navigation data
 *          double *rs       IO  sat position and corrected {x,y,z} (ecef) (m)
 *          double *dts      IO  sat clock bias and corrected (s)
 *          double *var      O   sat position and clock variance (m^2)
 * return : status (1:ok,0:no correction)
 * notes  : before calling the function, sbas satellite correction parameters
 *          in navigation data (nav->sbssat) must be set by callig
 *          sbsupdatecorr().
 *          satellite clock correction include long-term correction and fast
 *          correction.
 *          sbas clock correction is usually based on L1C/A code. TGD or DCB has
 *          to be considered for other codes
 *-----------------------------------------------------------------------------*/
int sbssatcorr(gtime_t time, int sat, const nav_t* nav, double* rs, double* dts, double* var) {
    double drs[3] = {0}, dclk = 0.0, prc = 0.0;
    int i;
 
    trace(NULL, 3, "sbssatcorr : sat=%2d\n", sat);
 
    /* sbas long term corrections */
    if (!sbslongcorr(time, sat, &nav->sbssat, drs, &dclk)) {
        return 0;
    }
    /* sbas fast corrections */
    if (!sbsfastcorr(time, sat, &nav->sbssat, &prc, var)) {
        return 0;
    }
    for (i = 0; i < 3; i++) {
        rs[i] += drs[i];
    }
 
    dts[0] += dclk + prc / CLIGHT;
 
    trace(NULL, 5, "sbssatcorr: sat=%2d drs=%6.3f %6.3f %6.3f dclk=%.3f %.3f var=%.3f\n", sat, drs[0], drs[1], drs[2],
          dclk, prc / CLIGHT, *var);
 
    return 1;
}
/* decode sbas message ---------------------------------------------------------
 * decode sbas message frame words and check crc
 * args   : gtime_t time     I   reception time
 *          int    prn       I   sbas satellite prn number
 *          uint32_t *word   I   message frame words (24bit x 10)
 *          sbsmsg_t *sbsmsg O   sbas message
 * return : status (1:ok,0:crc error)
 *-----------------------------------------------------------------------------*/
int sbsdecodemsg(gtime_t time, int prn, const uint32_t* words, sbsmsg_t* sbsmsg) {
    int i, j;
    uint8_t f[29];
    double tow;
 
    trace(NULL, 5, "sbsdecodemsg: prn=%d\n", prn);
 
    if (time.time == 0) {
        return 0;
    }
    tow = time2gpst(time, &sbsmsg->week);
    sbsmsg->tow = (int)(tow + DTTOL);
    sbsmsg->prn = prn;
    for (i = 0; i < 7; i++) {
        for (j = 0; j < 4; j++) {
            sbsmsg->msg[i * 4 + j] = (uint8_t)(words[i] >> ((3 - j) * 8));
        }
    }
    sbsmsg->msg[28] = (uint8_t)(words[7] >> 18) & 0xC0;
    for (i = 28; i > 0; i--) {
        f[i] = (sbsmsg->msg[i] >> 6) + (sbsmsg->msg[i - 1] << 2);
    }
    f[0] = sbsmsg->msg[0] >> 6;
 
    return rtk_crc24q(f, 29) == (words[7] & 0xFFFFFF); /* check crc */
}