MRTKLIB 0.7.14
Modernized GNSS Positioning Library
Loading...
Searching...
No Matches
mrtk_rtkpos.c
Go to the documentation of this file.
1/*------------------------------------------------------------------------------
2 * mrtk_rtkpos.c : RTK positioning functions
3 *
4 * Copyright (C) 2026 H.SHIONO (MRTKLIB Project)
5 * Copyright (C) 2023-2025 Cabinet Office, Japan
6 * Copyright (C) 2024-2025 Lighthouse Technology & Consulting Co. Ltd.
7 * Copyright (C) 2023-2025 Japan Aerospace Exploration Agency
8 * Copyright (C) 2023-2025 TOSHIBA ELECTRONIC TECHNOLOGIES CORPORATION
9 * Copyright (C) 2015- Mitsubishi Electric Corp.
10 * Copyright (C) 2014 Geospatial Information Authority of Japan
11 * Copyright (C) 2014 T.SUZUKI
12 * Copyright (C) 2013- Tim Everett (rtklibexplorer) [demo5 RTK enhancements]
13 * Copyright (C) 2007-2023 T.TAKASU
14 *
15 * SPDX-License-Identifier: BSD-2-Clause
16 *----------------------------------------------------------------------------*/
17#include "mrtklib/mrtk_rtkpos.h"
18
19#include <math.h>
20#include <stdarg.h>
21#include <stdio.h>
22#include <stdlib.h>
23#include <string.h>
24
25#include "mrtklib/mrtk_antenna.h"
26#include "mrtklib/mrtk_atmos.h"
27#include "mrtklib/mrtk_bias_sinex.h"
28#include "mrtklib/mrtk_bits.h"
29#include "mrtklib/mrtk_clas.h"
30#include "mrtklib/mrtk_const.h"
31#include "mrtklib/mrtk_coords.h"
32#include "mrtklib/mrtk_fcb.h"
33#include "mrtklib/mrtk_lambda.h"
34#include "mrtklib/mrtk_madoca.h"
35#include "mrtklib/mrtk_mat.h"
36#include "mrtklib/mrtk_peph.h"
37#include "mrtklib/mrtk_ppp.h"
38#include "mrtklib/mrtk_ppp_rtk.h"
39#include "mrtklib/mrtk_spp.h"
40#include "mrtklib/mrtk_sys.h"
41#include "mrtklib/mrtk_tides.h"
42#include "mrtklib/mrtk_trace.h"
43#include "mrtklib/mrtk_vrs.h"
44
45/* local constants (duplicated from rtklib.h to avoid circular dependency) */
46#define CLIGHT 299792458.0 /* speed of light (m/s) */
47#define D2R (3.1415926535897932384626433832795 / 180.0)
48#define R2D (180.0 / 3.1415926535897932384626433832795)
49#define FREQ1 1.57542E9 /* L1/E1/B1C frequency (Hz) */
50
51#define SYS_NONE 0x00 /* navigation system: none */
52#define SYS_GPS 0x01 /* navigation system: GPS */
53#define SYS_SBS 0x02 /* navigation system: SBAS */
54#define SYS_GLO 0x04 /* navigation system: GLONASS */
55#define SYS_GAL 0x08 /* navigation system: Galileo */
56#define SYS_QZS 0x10 /* navigation system: QZSS */
57#define SYS_CMP 0x20 /* navigation system: BeiDou */
58#define SYS_IRN 0x40 /* navigation system: IRNSS */
59
60#define EFACT_GPS 1.0 /* error factor: GPS */
61#define EFACT_GLO 1.5 /* error factor: GLONASS */
62#define EFACT_SBS 3.0 /* error factor: SBAS */
63
64#define INT_SWAP_STAT 86400.0 /* swap interval of solution status (s) */
65#define SSR_VENDOR_L6 0 /* SSR vendor: L6 */
66#define SSR_VENDOR_RTCM 1 /* SSR vendor: RTCM */
67#define MAXAGESSRL6 60.0 /* max age of SSR L6 (s) */
68
69/* forward declarations (implemented in rtkcmn.c/ppp_iono.c, resolved at link time) */
70extern int miono_get_corr(const double* rr, nav_t* nav);
71
72/* constants/macros ----------------------------------------------------------*/
73
74#define SQR(x) ((x) * (x))
75#define SQRT(x) ((x) <= 0.0 || (x) != (x) ? 0.0 : sqrt(x))
76#define MIN(x, y) ((x) <= (y) ? (x) : (y))
77#define ROUND(x) (int)floor((x) + 0.5)
78
79#define VAR_POS SQR(30.0) /* initial variance of receiver pos (m^2) */
80#define VAR_VEL SQR(10.0) /* initial variance of receiver vel ((m/s)^2) */
81#define VAR_ACC SQR(10.0) /* initial variance of receiver acc ((m/ss)^2) */
82#define VAR_HWBIAS SQR(1.0) /* initial variance of h/w bias ((m/MHz)^2) */
83#define VAR_GRA SQR(0.001) /* initial variance of gradient (m^2) */
84#define INIT_ZWD 0.15 /* initial zwd (m) */
85
86#define PRN_HWBIAS 1E-6 /* process noise of h/w bias (m/MHz/sqrt(s)) */
87#define GAP_RESION 120 /* gap to reset ionosphere parameters (epochs) */
88#define MAXACC 30.0 /* max accel for doppler slip detection (m/s^2) */
89
90#define VAR_HOLDAMB 0.001 /* constraint to hold ambiguity (cycle^2) */
91
92#define TTOL_MOVEB (1.0 + 2 * DTTOL)
93/* time sync tolerance for moving-baseline (s) */
94
95/* number of parameters (pos,ionos,tropos,hw-bias,phase-bias,real,estimated) */
96#define NF(opt) ((opt)->ionoopt == IONOOPT_IFLC ? 1 : (opt)->nf)
97#define NP(opt) ((opt)->dynamics == 0 ? 3 : 9)
98#define NI(opt) (((opt)->ionoopt != IONOOPT_EST && (opt)->ionoopt != IONOOPT_EST_ADPT) ? 0 : MAXSAT)
99#define NT(opt) ((opt)->tropopt < TROPOPT_EST ? 0 : ((opt)->tropopt < TROPOPT_ESTG ? 2 : 6))
100#define NL(opt) ((opt)->glomodear != 2 ? 0 : NFREQGLO)
101#define NB(opt) ((opt)->mode <= PMODE_DGPS ? 0 : MAXSAT * NF(opt))
102#define NR(opt) (NP(opt) + NI(opt) + NT(opt) + NL(opt))
103#define NX(opt) (NR(opt) + NB(opt))
104
105/* state variable index */
106#define II(s, opt) (NP(opt) + (s) - 1) /* ionos (s:satellite no) */
107#define IT(r, opt) (NP(opt) + NI(opt) + NT(opt) / 2 * (r)) /* tropos (r:0=rov,1:ref) */
108#define IL(f, opt) (NP(opt) + NI(opt) + NT(opt) + (f)) /* receiver h/w bias */
109#define IB(s, f, opt) (NR(opt) + MAXSAT * (f) + (s) - 1) /* phase bias (s:satno,f:freq) */
110
111/* polynomial coefficients for adaptive AR ratio threshold based on sat-pair count
112 (fitted to LAMBDA reliability curves from TU Delft; valid for 1–50 pairs) */
113static const double ar_poly_coeffs[3][5] = {
114 {-1.94058448e-01, -7.79023476e+00, 1.24231120e+02, -4.03126050e+02, 3.50413202e+02},
115 {6.42237302e-01, -8.39813962e+00, 2.92107285e+01, -2.37577308e+01, -1.14307128e+00},
116 {-2.22600390e-02, 3.23169103e-01, -1.39837429e+00, 2.19282996e+00, -5.34583971e-02}};
117
118/* global variables ----------------------------------------------------------*/
119static int statlevel = 0; /* rtk status output level (0:off) */
120static FILE* fp_stat = NULL; /* rtk status file pointer */
121static char file_stat[1024] = ""; /* rtk status file original path */
122static gtime_t time_stat = {0}; /* rtk status file time */
123
124/* open solution status file ---------------------------------------------------
125 * open solution status file and set output level
126 * args : char *file I rtk status file
127 * int level I rtk status level (0: off)
128 * return : status (1:ok,0:error)
129 * notes : file can constain time keywords (%Y,%y,%m...) defined in reppath().
130 * The time to replace keywords is based on UTC of CPU time.
131 * output : solution status file record format
132 *
133 * $POS,week,tow,stat,posx,posy,posz,posxf,posyf,poszf
134 * week/tow : gps week no/time of week (s)
135 * stat : solution status
136 * posx/posy/posz : position x/y/z ecef (m) float
137 * posxf/posyf/poszf : position x/y/z ecef (m) fixed
138 *
139 * $VELACC,week,tow,stat,vele,veln,velu,acce,accn,accu,velef,velnf,veluf,accef,accnf,accuf
140 * week/tow : gps week no/time of week (s)
141 * stat : solution status
142 * vele/veln/velu : velocity e/n/u (m/s) float
143 * acce/accn/accu : acceleration e/n/u (m/s^2) float
144 * velef/velnf/veluf : velocity e/n/u (m/s) fixed
145 * accef/accnf/accuf : acceleration e/n/u (m/s^2) fixed
146 *
147 * $CLK,week,tow,stat,clk1,clk2,clk3,clk4
148 * week/tow : gps week no/time of week (s)
149 * stat : solution status
150 * clk1 : receiver clock bias GPS (ns)
151 * clk2 : receiver clock bias GLO-GPS (ns)
152 * clk3 : receiver clock bias GAL-GPS (ns)
153 * clk4 : receiver clock bias BDS-GPS (ns)
154 *
155 * $ION,week,tow,stat,sat,az,el,ion,ion-fixed
156 * week/tow : gps week no/time of week (s)
157 * stat : solution status
158 * sat : satellite id
159 * az/el : azimuth/elevation angle(deg)
160 * ion : vertical ionospheric delay L1 (m) float
161 * ion-fixed: vertical ionospheric delay L1 (m) fixed
162 *
163 * $TROP,week,tow,stat,rcv,ztd,ztdf
164 * week/tow : gps week no/time of week (s)
165 * stat : solution status
166 * rcv : receiver (1:rover,2:base station)
167 * ztd : zenith total delay (m) float
168 * ztdf : zenith total delay (m) fixed
169 *
170 * $HWBIAS,week,tow,stat,frq,bias,biasf
171 * week/tow : gps week no/time of week (s)
172 * stat : solution status
173 * frq : frequency (1:L1,2:L2,...)
174 * bias : h/w bias coefficient (m/MHz) float
175 * biasf : h/w bias coefficient (m/MHz) fixed
176 *
177 * $SAT,week,tow,sat,frq,az,el,resp,resc,vsat,snr,fix,slip,lock,outc,slipc,rejc
178 * week/tow : gps week no/time of week (s)
179 * sat/frq : satellite id/frequency (1:L1,2:L2,...)
180 * az/el : azimuth/elevation angle (deg)
181 * resp : pseudorange residual (m)
182 * resc : carrier-phase residual (m)
183 * vsat : valid data flag (0:invalid,1:valid)
184 * snr : signal strength (dbHz)
185 * fix : ambiguity flag (0:no data,1:float,2:fixed,3:hold,4:ppp)
186 * slip : cycle-slip flag (bit1:slip,bit2:parity unknown)
187 * lock : carrier-lock count
188 * outc : data outage count
189 * slipc : cycle-slip count
190 * rejc : data reject (outlier) count
191 *
192 *-----------------------------------------------------------------------------*/
193extern int rtkopenstat(const char* file, int level) {
194 gtime_t time = utc2gpst(timeget());
195 char path[1024];
196
197 trace(NULL, 3, "rtkopenstat: file=%s level=%d\n", file, level);
198
199 if (level <= 0) {
200 return 0;
201 }
202
203 reppath(file, path, time, "", "");
204
205 if (!(fp_stat = fopen(path, "w"))) {
206 trace(NULL, 1, "rtkopenstat: file open error path=%s\n", path);
207 return 0;
208 }
209 strcpy(file_stat, file);
210 time_stat = time;
211 statlevel = level;
212 return 1;
213}
214/* close solution status file --------------------------------------------------
215 * close solution status file
216 * args : none
217 * return : none
218 *-----------------------------------------------------------------------------*/
219extern void rtkclosestat(void) {
220 trace(NULL, 3, "rtkclosestat:\n");
221
222 if (fp_stat) {
223 fclose(fp_stat);
224 }
225 fp_stat = NULL;
226 file_stat[0] = '\0';
227 statlevel = 0;
228}
229/* write solution status to buffer -------------------------------------------*/
230extern int rtkoutstat(rtk_t* rtk, char* buff) {
231 ssat_t* ssat;
232 double tow, pos[3], vel[3], acc[3], vela[3] = {0}, acca[3] = {0}, xa[3];
233 int i, j, week, est, nfreq, nf = NF(&rtk->opt);
234 char id[32], *p = buff;
235
236 if (rtk->sol.stat <= SOLQ_NONE) {
237 return 0;
238 }
239 /* write ppp solution status to buffer */
240 if (rtk->opt.mode >= PMODE_PPP_KINEMA) {
241 return pppoutstat(rtk, buff);
242 }
243 est = rtk->opt.mode >= PMODE_DGPS;
244 nfreq = est ? nf : 1;
245 tow = time2gpst(rtk->sol.time, &week);
246
247 /* receiver position */
248 if (est) {
249 for (i = 0; i < 3; i++) {
250 xa[i] = i < rtk->na ? rtk->xa[i] : 0.0;
251 }
252 p += sprintf(p, "$POS,%d,%.3f,%d,%.4f,%.4f,%.4f,%.4f,%.4f,%.4f\n", week, tow, rtk->sol.stat, rtk->x[0],
253 rtk->x[1], rtk->x[2], xa[0], xa[1], xa[2]);
254 } else {
255 p += sprintf(p, "$POS,%d,%.3f,%d,%.4f,%.4f,%.4f,%.4f,%.4f,%.4f\n", week, tow, rtk->sol.stat, rtk->sol.rr[0],
256 rtk->sol.rr[1], rtk->sol.rr[2], 0.0, 0.0, 0.0);
257 }
258 /* receiver velocity and acceleration */
259 if (est && rtk->opt.dynamics) {
260 ecef2pos(rtk->sol.rr, pos);
261 ecef2enu(pos, rtk->x + 3, vel);
262 ecef2enu(pos, rtk->x + 6, acc);
263 if (rtk->na >= 6) {
264 ecef2enu(pos, rtk->xa + 3, vela);
265 }
266 if (rtk->na >= 9) {
267 ecef2enu(pos, rtk->xa + 6, acca);
268 }
269 p += sprintf(p, "$VELACC,%d,%.3f,%d,%.4f,%.4f,%.4f,%.5f,%.5f,%.5f,%.4f,%.4f,%.4f,%.5f,%.5f,%.5f\n", week, tow,
270 rtk->sol.stat, vel[0], vel[1], vel[2], acc[0], acc[1], acc[2], vela[0], vela[1], vela[2], acca[0],
271 acca[1], acca[2]);
272 } else {
273 ecef2pos(rtk->sol.rr, pos);
274 ecef2enu(pos, rtk->sol.rr + 3, vel);
275 p += sprintf(p, "$VELACC,%d,%.3f,%d,%.4f,%.4f,%.4f,%.5f,%.5f,%.5f,%.4f,%.4f,%.4f,%.5f,%.5f,%.5f\n", week, tow,
276 rtk->sol.stat, vel[0], vel[1], vel[2], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
277 }
278 /* receiver clocks */
279 p += sprintf(p, "$CLK,%d,%.3f,%d,%d,%.3f,%.3f,%.3f,%.3f\n", week, tow, rtk->sol.stat, 1, rtk->sol.dtr[0] * 1E9,
280 rtk->sol.dtr[1] * 1E9, rtk->sol.dtr[2] * 1E9, rtk->sol.dtr[3] * 1E9);
281
282 /* ionospheric parameters */
283 if (est && rtk->opt.ionoopt == IONOOPT_EST) {
284 for (i = 0; i < MAXSAT; i++) {
285 ssat = rtk->ssat + i;
286 if (!ssat->vs) {
287 continue;
288 }
289 satno2id(i + 1, id);
290 j = II(i + 1, &rtk->opt);
291 xa[0] = j < rtk->na ? rtk->xa[j] : 0.0;
292 p += sprintf(p, "$ION,%d,%.3f,%d,%s,%.1f,%.1f,%.4f,%.4f\n", week, tow, rtk->sol.stat, id,
293 ssat->azel[0] * R2D, ssat->azel[1] * R2D, rtk->x[j], xa[0]);
294 }
295 }
296 /* tropospheric parameters */
297 if (est && (rtk->opt.tropopt == TROPOPT_EST || rtk->opt.tropopt == TROPOPT_ESTG)) {
298 for (i = 0; i < 2; i++) {
299 j = IT(i, &rtk->opt);
300 xa[0] = j < rtk->na ? rtk->xa[j] : 0.0;
301 p += sprintf(p, "$TROP,%d,%.3f,%d,%d,%.4f,%.4f\n", week, tow, rtk->sol.stat, i + 1, rtk->x[j], xa[0]);
302 }
303 }
304 /* receiver h/w bias */
305 if (est && rtk->opt.glomodear == 2) {
306 for (i = 0; i < nfreq; i++) {
307 j = IL(i, &rtk->opt);
308 xa[0] = j < rtk->na ? rtk->xa[j] : 0.0;
309 p += sprintf(p, "$HWBIAS,%d,%.3f,%d,%d,%.4f,%.4f\n", week, tow, rtk->sol.stat, i + 1, rtk->x[j], xa[0]);
310 }
311 }
312 return (int)(p - buff);
313}
314/* swap solution status file -------------------------------------------------*/
315static void swapsolstat(void) {
316 gtime_t time = utc2gpst(timeget());
317 char path[1024];
318
319 if ((int)(time2gpst(time, NULL) / INT_SWAP_STAT) == (int)(time2gpst(time_stat, NULL) / INT_SWAP_STAT)) {
320 return;
321 }
322 time_stat = time;
323
324 if (!reppath(file_stat, path, time, "", "")) {
325 return;
326 }
327 if (fp_stat) {
328 fclose(fp_stat);
329 }
330
331 if (!(fp_stat = fopen(path, "w"))) {
332 trace(NULL, 2, "swapsolstat: file open error path=%s\n", path);
333 return;
334 }
335 trace(NULL, 3, "swapsolstat: path=%s\n", path);
336}
337/* output solution status ----------------------------------------------------*/
338static void outsolstat(rtk_t* rtk) {
339 ssat_t* ssat;
340 double tow;
341 char buff[MAXSOLMSG + 1], id[32];
342 int i, j, n, week, nfreq, nf = NF(&rtk->opt);
343
344 if (statlevel <= 0 || !fp_stat || !rtk->sol.stat) {
345 return;
346 }
347
348 trace(NULL, 3, "outsolstat:\n");
349
350 /* swap solution status file */
351 swapsolstat();
352
353 /* write solution status */
354 n = rtkoutstat(rtk, buff);
355 buff[n] = '\0';
356
357 fputs(buff, fp_stat);
358
359 if (rtk->sol.stat == SOLQ_NONE || statlevel <= 1) {
360 return;
361 }
362
363 tow = time2gpst(rtk->sol.time, &week);
364 nfreq = rtk->opt.mode >= PMODE_DGPS ? nf : 1;
365
366 /* write residuals and status */
367 for (i = 0; i < MAXSAT; i++) {
368 ssat = rtk->ssat + i;
369 if (!ssat->vs) {
370 continue;
371 }
372 satno2id(i + 1, id);
373 for (j = 0; j < nfreq; j++) {
374 fprintf(fp_stat, "$SAT,%d,%.3f,%s,%d,%.1f,%.1f,%.4f,%.4f,%d,%.1f,%d,%d,%d,%d,%d,%d\n", week, tow, id, j + 1,
375 ssat->azel[0] * R2D, ssat->azel[1] * R2D, ssat->resp[j], ssat->resc[j], ssat->vsat[j],
376 ssat->snr[j] * SNR_UNIT, ssat->fix[j], ssat->slip[j] & 3, ssat->lock[j], ssat->outc[j],
377 ssat->slipc[j], ssat->rejc[j]);
378 }
379 }
380}
381/* save error message --------------------------------------------------------*/
382static void errmsg(rtk_t* rtk, const char* format, ...) {
383 char buff[256], tstr[32];
384 int n;
385 va_list ap;
386 time2str(rtk->sol.time, tstr, 2);
387 n = sprintf(buff, "%s: ", tstr + 11);
388 va_start(ap, format);
389 n += vsprintf(buff + n, format, ap);
390 va_end(ap);
391 n = n < MAXERRMSG - rtk->neb ? n : MAXERRMSG - rtk->neb;
392 memcpy(rtk->errbuf + rtk->neb, buff, n);
393 rtk->neb += n;
394 trace(NULL, 2, "%s", buff);
395}
396/* single-differenced observable ---------------------------------------------*/
397static double sdobs(const obsd_t* obs, int i, int j, int k) {
398 double pi = (k < NFREQ) ? obs[i].L[k] : obs[i].P[k - NFREQ];
399 double pj = (k < NFREQ) ? obs[j].L[k] : obs[j].P[k - NFREQ];
400 return pi == 0.0 || pj == 0.0 ? 0.0 : pi - pj;
401}
402/* single-differenced geometry-free linear combination of phase --------------*/
403static double gfobs(const obsd_t* obs, int i, int j, int k, const nav_t* nav) {
404 double freq1, freq2, L1, L2;
405
406 freq1 = sat2freq(obs[i].sat, obs[i].code[0], nav);
407 freq2 = sat2freq(obs[i].sat, obs[i].code[k], nav);
408 L1 = sdobs(obs, i, j, 0);
409 L2 = sdobs(obs, i, j, k);
410 if (freq1 == 0.0 || freq2 == 0.0 || L1 == 0.0 || L2 == 0.0) {
411 return 0.0;
412 }
413 return L1 * CLIGHT / freq1 - L2 * CLIGHT / freq2;
414}
415/* single-differenced measurement error variance -----------------------------*/
416/* select IFLC second frequency (L5 for Galileo if nf>=3, L2 otherwise) -----*/
417static int seliflc(int optnf, int sys) { return (optnf == 2 || sys != SYS_GAL) ? 1 : 2; }
418static double varerr(int sat, int sys, double el, double snr_rover, double snr_base, double bl, double dt, int f,
419 const prcopt_t* opt) {
420 double a, b, c = opt->err[3] * bl / 1E4, d = CLIGHT * opt->sclkstab * dt, fact;
421 double sinel = sin(el), var;
422 int nf = NF(opt), frq, code;
423
424 frq = f % nf;
425 code = (f < nf) ? 0 : 1;
426 /* code uses eratio[frq] directly; phase scaled relative to eratio[0] */
427 if (code) {
428 fact = opt->eratio[frq];
429 } else {
430 fact = opt->eratio[frq] / opt->eratio[0];
431 }
432 if (fact <= 0.0) {
433 fact = opt->eratio[0];
434 }
435
436 /* per-constellation error factor */
437 switch (sys) {
438 case SYS_GPS:
439 fact *= EFACT_GPS;
440 break;
441 case SYS_GLO:
442 fact *= EFACT_GLO;
443 break;
444 case SYS_GAL:
445 fact *= EFACT_GAL;
446 break;
447 case SYS_SBS:
448 fact *= EFACT_SBS;
449 break;
450 case SYS_QZS:
451 fact *= EFACT_QZS;
452 break;
453 case SYS_CMP:
454 fact *= EFACT_CMP;
455 break;
456 case SYS_IRN:
457 fact *= EFACT_IRN;
458 break;
459 default:
460 fact *= EFACT_GPS;
461 break;
462 }
463 a = fact * opt->err[1];
464 b = fact * opt->err[2];
465
466 var = 2.0 * (a * a + b * b / sinel / sinel + c * c) + d * d;
467
468 if (opt->err[6] > 0.0) { /* add SNR term */
469 double e = fact * opt->err[6];
470 double snr_max = opt->err[5];
471 double dr = snr_max - snr_rover, db = snr_max - snr_base;
472 var += e * e * (pow(10.0, 0.1 * (dr > 0.0 ? dr : 0.0)) + pow(10.0, 0.1 * (db > 0.0 ? db : 0.0)));
473 }
474 /* IFLC noise amplification applied to entire variance (SQR(3)=9×) */
475 var *= (opt->ionoopt == IONOOPT_IFLC) ? SQR(3.0) : 1.0;
476 return var;
477}
478/* baseline length -----------------------------------------------------------*/
479static double baseline(const double* ru, const double* rb, double* dr) {
480 int i;
481 for (i = 0; i < 3; i++) {
482 dr[i] = ru[i] - rb[i];
483 }
484 return norm(dr, 3);
485}
486/* initialize state and covariance -------------------------------------------*/
487static void initx(rtk_t* rtk, double xi, double var, int i) {
488 int j;
489 rtk->x[i] = xi;
490 for (j = 0; j < rtk->nx; j++) {
491 rtk->P[i + j * rtk->nx] = rtk->P[j + i * rtk->nx] = i == j ? var : 0.0;
492 }
493}
494/* select common satellites between rover and reference station --------------*/
495static int selsat(const obsd_t* obs, double* azel, int nu, int nr, const prcopt_t* opt, int* sat, int* iu, int* ir) {
496 int i, j, k = 0;
497
498 trace(NULL, 3, "selsat : nu=%d nr=%d\n", nu, nr);
499
500 for (i = 0, j = nu; i < nu && j < nu + nr; i++, j++) {
501 if (obs[i].sat < obs[j].sat) {
502 j--;
503 } else if (obs[i].sat > obs[j].sat) {
504 i--;
505 } else if (azel[1 + j * 2] >= opt->elmin) { /* elevation at base station */
506 sat[k] = obs[i].sat;
507 iu[k] = i;
508 ir[k++] = j;
509 trace(NULL, 4, "(%2d) sat=%3d iu=%2d ir=%2d\n", k - 1, obs[i].sat, i, j);
510 }
511 }
512 return k;
513}
514/* temporal update of position/velocity/acceleration -------------------------*/
515static void udpos(rtk_t* rtk, double tt) {
516 double *F, *P, *FP, *x, *xp, pos[3], Q[9] = {0}, Qv[9], var = 0.0;
517 int i, j, *ix, nx;
518
519 trace(NULL, 3, "udpos : tt=%.3f\n", tt);
520
521 /* fixed mode */
522 if (rtk->opt.mode == PMODE_FIXED) {
523 for (i = 0; i < 3; i++) {
524 initx(rtk, rtk->opt.ru[i], 1E-8, i);
525 }
526 return;
527 }
528 /* initialize position for first epoch */
529 if (norm(rtk->x, 3) <= 0.0) {
530 for (i = 0; i < 3; i++) {
531 initx(rtk, rtk->sol.rr[i], VAR_POS, i);
532 }
533 if (rtk->opt.dynamics) {
534 for (i = 3; i < 6; i++) {
535 initx(rtk, rtk->sol.rr[i], VAR_VEL, i);
536 }
537 for (i = 6; i < 9; i++) {
538 initx(rtk, 1E-6, VAR_ACC, i);
539 }
540 }
541 }
542 /* static mode */
543 if (rtk->opt.mode == PMODE_STATIC) {
544 return;
545 }
546
547 /* kinmatic mode without dynamics */
548 if (!rtk->opt.dynamics) {
549 for (i = 0; i < 3; i++) {
550 initx(rtk, rtk->sol.rr[i], VAR_POS, i);
551 }
552 return;
553 }
554 /* check variance of estimated position */
555 for (i = 0; i < 3; i++) {
556 var += rtk->P[i + i * rtk->nx];
557 }
558 var /= 3.0;
559
560 if (var > VAR_POS) {
561 /* reset position with large variance */
562 for (i = 0; i < 3; i++) {
563 initx(rtk, rtk->sol.rr[i], VAR_POS, i);
564 }
565 for (i = 3; i < 6; i++) {
566 initx(rtk, rtk->sol.rr[i], VAR_VEL, i);
567 }
568 for (i = 6; i < 9; i++) {
569 initx(rtk, 1E-6, VAR_ACC, i);
570 }
571 trace(NULL, 2, "reset rtk position due to large variance: var=%.3f\n", var);
572 return;
573 }
574 /* generate valid state index */
575 ix = imat(rtk->nx, 1);
576 for (i = nx = 0; i < rtk->nx; i++) {
577 if (rtk->x[i] != 0.0 && rtk->P[i + i * rtk->nx] > 0.0) {
578 ix[nx++] = i;
579 }
580 }
581 if (nx < 9) {
582 free(ix);
583 return;
584 }
585 /* state transition of position/velocity/acceleration */
586 F = eye(nx);
587 P = mat(nx, nx);
588 FP = mat(nx, nx);
589 x = mat(nx, 1);
590 xp = mat(nx, 1);
591
592 for (i = 0; i < 6; i++) {
593 F[i + (i + 3) * nx] = tt;
594 }
595 /* include accel coupling only when position variance is low (converged) */
596 if (var < rtk->opt.thresar[1]) {
597 for (i = 0; i < 3; i++) {
598 F[i + (i + 6) * nx] = (tt >= 0.0 ? 1.0 : -1.0) * SQR(tt) / 2.0;
599 }
600 } else {
601 trace(NULL, 3, "pos var too high for accel term: var=%.4f\n", var);
602 }
603 for (i = 0; i < nx; i++) {
604 x[i] = rtk->x[ix[i]];
605 for (j = 0; j < nx; j++) {
606 P[i + j * nx] = rtk->P[ix[i] + ix[j] * rtk->nx];
607 }
608 }
609 /* x=F*x, P=F*P*F+Q */
610 matmul("NN", nx, 1, nx, 1.0, F, x, 0.0, xp);
611 matmul("NN", nx, nx, nx, 1.0, F, P, 0.0, FP);
612 matmul("NT", nx, nx, nx, 1.0, FP, F, 0.0, P);
613
614 for (i = 0; i < nx; i++) {
615 rtk->x[ix[i]] = xp[i];
616 for (j = 0; j < nx; j++) {
617 rtk->P[ix[i] + ix[j] * rtk->nx] = P[i + j * nx];
618 }
619 }
620 /* process noise added to only acceleration */
621 Q[0] = Q[4] = SQR(rtk->opt.prn[3]) * fabs(tt);
622 Q[8] = SQR(rtk->opt.prn[4]) * fabs(tt);
623 ecef2pos(rtk->x, pos);
624 covecef(pos, Q, Qv);
625 for (i = 0; i < 3; i++) {
626 for (j = 0; j < 3; j++) {
627 rtk->P[i + 6 + (j + 6) * rtk->nx] += Qv[i + j * 3];
628 }
629 }
630 free(ix);
631 free(F);
632 free(P);
633 free(FP);
634 free(x);
635 free(xp);
636}
637/* temporal update of ionospheric parameters ---------------------------------*/
638static void udion(rtk_t* rtk, double tt, double bl, const int* sat, int ns) {
639 double el, fact;
640 int i, j;
641
642 trace(NULL, 3, "udion : tt=%.3f bl=%.0f ns=%d\n", tt, bl, ns);
643
644 for (i = 1; i <= MAXSAT; i++) {
645 j = II(i, &rtk->opt);
646 if (rtk->x[j] != 0.0 && rtk->ssat[i - 1].outc[0] > GAP_RESION && rtk->ssat[i - 1].outc[1] > GAP_RESION) {
647 rtk->x[j] = 0.0;
648 }
649 }
650 for (i = 0; i < ns; i++) {
651 j = II(sat[i], &rtk->opt);
652
653 if (rtk->x[j] == 0.0) {
654 initx(rtk, 1E-6, SQR(rtk->opt.std[1] * bl / 1E4), j);
655 } else {
656 /* elevation dependent factor of process noise */
657 el = rtk->ssat[sat[i] - 1].azel[1];
658 fact = cos(el);
659 rtk->P[j + j * rtk->nx] += SQR(rtk->opt.prn[1] * bl / 1E4 * fact) * fabs(tt);
660 }
661 }
662}
663/* temporal update of tropospheric parameters --------------------------------*/
664static void udtrop(rtk_t* rtk, double tt, double bl) {
665 int i, j, k;
666
667 trace(NULL, 3, "udtrop : tt=%.3f\n", tt);
668
669 for (i = 0; i < 2; i++) {
670 j = IT(i, &rtk->opt);
671
672 if (rtk->x[j] == 0.0) {
673 initx(rtk, INIT_ZWD, SQR(rtk->opt.std[2]), j); /* initial zwd */
674
675 if (rtk->opt.tropopt >= TROPOPT_ESTG) {
676 for (k = 0; k < 2; k++) {
677 initx(rtk, 1E-6, VAR_GRA, ++j);
678 }
679 }
680 } else {
681 rtk->P[j + j * rtk->nx] += SQR(rtk->opt.prn[2]) * fabs(tt);
682
683 if (rtk->opt.tropopt >= TROPOPT_ESTG) {
684 for (k = 0; k < 2; k++) {
685 rtk->P[++j * (1 + rtk->nx)] += SQR(rtk->opt.prn[2] * 0.3) * fabs(tt);
686 }
687 }
688 }
689 }
690}
691/* temporal update of receiver h/w biases ------------------------------------*/
692static void udrcvbias(rtk_t* rtk, double tt) {
693 int i, j;
694
695 trace(NULL, 3, "udrcvbias: tt=%.3f\n", tt);
696
697 for (i = 0; i < NFREQGLO; i++) {
698 j = IL(i, &rtk->opt);
699
700 if (rtk->x[j] == 0.0) {
701 initx(rtk, 1E-6, VAR_HWBIAS, j);
702 }
703 /* hold to fixed solution */
704 else if (rtk->nfix >= rtk->opt.minfix && rtk->sol.ratio > rtk->opt.thresar[0]) {
705 initx(rtk, rtk->xa[j], rtk->Pa[j + j * rtk->na], j);
706 } else {
707 rtk->P[j + j * rtk->nx] += SQR(PRN_HWBIAS) * fabs(tt);
708 }
709 }
710}
711/* detect cycle slip by LLI --------------------------------------------------*/
712/* detect cycle slip by observation code change (e.g. C1C→C1S switch) --------*/
713static void detslp_code(rtk_t* rtk, const obsd_t* obs, int i, int rcv) {
714 int f, sat = obs[i].sat, nf = rtk->opt.nf < NFREQ ? rtk->opt.nf : NFREQ;
715 for (f = 0; f < nf; f++) {
716 uint8_t code = obs[i].code[f];
717 if (code == CODE_NONE) {
718 continue;
719 }
720 uint8_t ccode = rtk->ssat[sat - 1].codeprev[f][rcv - 1];
721 if (code != ccode) {
722 rtk->ssat[sat - 1].codeprev[f][rcv - 1] = code;
723 /* if previously tracked, flag as cycle slip */
724 if (ccode != CODE_NONE) {
725 rtk->ssat[sat - 1].slip[f] |= 1;
726 errmsg(rtk,
727 "slip detected code change (sat=%2d rcv=%d F=%d "
728 "%s->%s)\n",
729 sat, rcv, f + 1, code2obs(ccode), code2obs(code));
730 }
731 }
732 }
733}
734/* detect cycle slip by LLI --------------------------------------------------*/
735static void detslp_ll(rtk_t* rtk, const obsd_t* obs, int i, int rcv) {
736 uint32_t slip, LLI;
737 int f, sat = obs[i].sat;
738
739 trace(NULL, 3, "detslp_ll: i=%d rcv=%d\n", i, rcv);
740
741 for (f = 0; f < rtk->opt.nf; f++) {
742 if (obs[i].L[f] == 0.0 || fabs(timediff(obs[i].time, rtk->ssat[sat - 1].pt[rcv - 1][f])) < DTTOL) {
743 continue;
744 }
745 /* restore previous LLI */
746 if (rcv == 1) {
747 LLI = getbitu(&rtk->ssat[sat - 1].slip[f], 0, 2); /* rover */
748 } else {
749 LLI = getbitu(&rtk->ssat[sat - 1].slip[f], 2, 2); /* base */
750 }
751
752 /* detect slip by cycle slip flag in LLI */
753 if (rtk->tt >= 0.0) { /* forward */
754 if (obs[i].LLI[f] & 1) {
755 errmsg(rtk, "slip detected forward (sat=%2d rcv=%d F=%d LLI=%x)\n", sat, rcv, f + 1, obs[i].LLI[f]);
756 }
757 slip = obs[i].LLI[f];
758 } else { /* backward */
759 if (LLI & 1) {
760 errmsg(rtk, "slip detected backward (sat=%2d rcv=%d F=%d LLI=%x)\n", sat, rcv, f + 1, LLI);
761 }
762 slip = LLI;
763 }
764 /* detect slip by parity unknown flag transition in LLI */
765 if (((LLI & 2) && !(obs[i].LLI[f] & 2)) || (!(LLI & 2) && (obs[i].LLI[f] & 2))) {
766 errmsg(rtk, "slip detected half-cyc (sat=%2d rcv=%d F=%d LLI=%x->%x)\n", sat, rcv, f + 1, LLI,
767 obs[i].LLI[f]);
768 slip |= 1;
769 }
770 /* save current LLI */
771 if (rcv == 1) {
772 setbitu(&rtk->ssat[sat - 1].slip[f], 0, 2, obs[i].LLI[f]);
773 } else {
774 setbitu(&rtk->ssat[sat - 1].slip[f], 2, 2, obs[i].LLI[f]);
775 }
776
777 /* save slip and half-cycle valid flag */
778 rtk->ssat[sat - 1].slip[f] |= (uint8_t)slip;
779 rtk->ssat[sat - 1].half[f] = (obs[i].LLI[f] & 2) ? 0 : 1;
780 }
781}
782/* detect cycle slip by geometry free phase jump -----------------------------*/
783static void detslp_gf(rtk_t* rtk, const obsd_t* obs, int i, int j, const nav_t* nav) {
784 int k, sat = obs[i].sat;
785 double g0, g1;
786
787 trace(NULL, 3, "detslp_gf: i=%d j=%d\n", i, j);
788
789 /* skip if disabled or slip already detected on any frequency */
790 if (rtk->opt.thresslip == 0.0) {
791 return;
792 }
793 for (k = 0; k < rtk->opt.nf; k++) {
794 if (rtk->ssat[sat - 1].slip[k] & 1) {
795 return;
796 }
797 }
798
799 for (k = 1; k < rtk->opt.nf; k++) {
800 if ((g1 = gfobs(obs, i, j, k, nav)) == 0.0) {
801 continue; /* try next freq, not abort */
802 }
803
804 g0 = rtk->ssat[sat - 1].gf[k - 1];
805 rtk->ssat[sat - 1].gf[k - 1] = g1;
806
807 if (g0 != 0.0 && fabs(g1 - g0) > rtk->opt.thresslip) {
808 rtk->ssat[sat - 1].slip[0] |= 1;
809 rtk->ssat[sat - 1].slip[k] |= 1;
810 errmsg(rtk, "slip detected GF jump (sat=%2d L1-L%d GF=%.3f %.3f)\n", sat, k + 1, g0, g1);
811 }
812 }
813}
814/* detect cycle slip by doppler and phase difference -------------------------*/
815/* detect cycle slip by doppler-phase difference with mean clock-jump removal --
816 * Mean of all-sat doppler differences removes common receiver clock jump,
817 * allowing per-satellite slip detection without false positives.
818 * Controlled by opt.thresdop (cycle/s); set to 0 to disable. */
819static void detslp_dop(rtk_t* rtk, const obsd_t* obs, const int* ix, int ns, int rcv, const nav_t* nav) {
820 (void)nav;
821 int i, ii, f, sat, ndop = 0, nf = MIN(rtk->opt.nf, NFREQ); /* clamp to array bounds */
822 double dph, dpt, mean_dop = 0.0;
823 double dopdif[MAXSAT][NFREQ], tt[MAXSAT][NFREQ];
824
825 trace(NULL, 4, "detslp_dop: rcv=%d\n", rcv);
826
827 if (rtk->opt.thresdop <= 0.0) {
828 return; /* disabled */
829 }
830
831 /* compute per-sat doppler-phase difference rate (cycle/s) */
832 for (i = 0; i < ns; i++) {
833 ii = ix[i];
834 sat = obs[ii].sat;
835 for (f = 0; f < nf; f++) {
836 dopdif[i][f] = 0.0;
837 tt[i][f] = 0.0;
838 if (obs[ii].L[f] == 0.0 || obs[ii].D[f] == 0.0 || rtk->ssat[sat - 1].ph[rcv - 1][f] == 0.0) {
839 continue;
840 }
841 if (fabs(tt[i][f] = timediff(obs[ii].time, rtk->ssat[sat - 1].pt[rcv - 1][f])) < DTTOL) {
842 continue;
843 }
844 dph = (obs[ii].L[f] - rtk->ssat[sat - 1].ph[rcv - 1][f]) / tt[i][f];
845 dpt = -obs[ii].D[f];
846 dopdif[i][f] = dph - dpt;
847 /* include in mean only if not an extreme outlier (3× threshold) */
848 if (fabs(dopdif[i][f]) < 3.0 * rtk->opt.thresdop) {
849 mean_dop += dopdif[i][f];
850 ndop++;
851 }
852 }
853 }
854 if (ndop == 0) {
855 return; /* no valid obs — likely large clock error */
856 }
857 mean_dop /= ndop; /* mean ≈ common clock-jump component */
858
859 /* flag slip if mean-removed difference exceeds threshold */
860 for (i = 0; i < ns; i++) {
861 sat = obs[ix[i]].sat;
862 for (f = 0; f < nf; f++) {
863 if (dopdif[i][f] == 0.0) {
864 continue;
865 }
866 if (fabs(dopdif[i][f] - mean_dop) > rtk->opt.thresdop) {
867 rtk->ssat[sat - 1].slip[f] |= 1;
868 errmsg(rtk,
869 "slip detected doppler (sat=%2d rcv=%d dL%d=%.3f "
870 "off=%.3f tt=%.2f)\n",
871 sat, rcv, f + 1, dopdif[i][f] - mean_dop, mean_dop, tt[i][f]);
872 }
873 }
874 }
875}
876/* temporal update of phase biases -------------------------------------------*/
877static void udbias(rtk_t* rtk, double tt, const obsd_t* obs, const int* sat, const int* iu, const int* ir, int ns,
878 const nav_t* nav) {
879 double cp, pr, cp1, cp2, pr1, pr2, *bias, offset, freqi, freq1, freq2, C1, C2;
880 int i, j, k, slip, rejc, reset, nf = NF(&rtk->opt);
881
882 trace(NULL, 3, "udbias : tt=%.3f ns=%d\n", tt, ns);
883
884 /* clear cycle-slip flags for all common sats */
885 for (i = 0; i < ns; i++) {
886 for (k = 0; k < rtk->opt.nf; k++) {
887 rtk->ssat[sat[i] - 1].slip[k] &= 0xFC;
888 }
889 }
890 /* detect cycle slip by doppler-phase difference (processes all sats at once
891 to remove common clock-jump via mean subtraction) */
892 detslp_dop(rtk, obs, iu, ns, 1, nav);
893 detslp_dop(rtk, obs, ir, ns, 2, nav);
894
895 for (i = 0; i < ns; i++) {
896 /* detect cycle slip by observation code change */
897 detslp_code(rtk, obs, iu[i], 1);
898 detslp_code(rtk, obs, ir[i], 2);
899
900 /* detect cycle slip by LLI */
901 detslp_ll(rtk, obs, iu[i], 1);
902 detslp_ll(rtk, obs, ir[i], 2);
903
904 /* detect cycle slip by geometry-free phase jump */
905 detslp_gf(rtk, obs, iu[i], ir[i], nav);
906
907 /* update half-cycle valid flag */
908 for (k = 0; k < nf; k++) {
909 rtk->ssat[sat[i] - 1].half[k] = !((obs[iu[i]].LLI[k] & 2) || (obs[ir[i]].LLI[k] & 2));
910 }
911 }
912 for (k = 0; k < nf; k++) {
913 /* reset phase-bias if instantaneous AR or expire obs outage counter */
914 for (i = 1; i <= MAXSAT; i++) {
915 reset = ++rtk->ssat[i - 1].outc[k] > (uint32_t)rtk->opt.maxout;
916
917 if (rtk->opt.modear == ARMODE_INST && rtk->x[IB(i, k, &rtk->opt)] != 0.0) {
918 initx(rtk, 0.0, 0.0, IB(i, k, &rtk->opt));
919 } else if (reset && rtk->x[IB(i, k, &rtk->opt)] != 0.0) {
920 initx(rtk, 0.0, 0.0, IB(i, k, &rtk->opt));
921 trace(NULL, 3, "udbias : obs outage counter overflow (sat=%3d L%d n=%d)\n", i, k + 1,
922 rtk->ssat[i - 1].outc[k]);
923 rtk->ssat[i - 1].outc[k] = 0;
924 }
925 if (rtk->opt.modear != ARMODE_INST && reset) {
926 rtk->ssat[i - 1].lock[k] = -rtk->opt.minlock;
927 }
928 }
929 /* reset phase-bias if detecting cycle slip */
930 for (i = 0; i < ns; i++) {
931 j = IB(sat[i], k, &rtk->opt);
932 rtk->P[j + j * rtk->nx] += rtk->opt.prn[0] * rtk->opt.prn[0] * fabs(tt);
933 slip = rtk->ssat[sat[i] - 1].slip[k];
934 if (rtk->opt.ionoopt == IONOOPT_IFLC) {
935 int f2i = seliflc(rtk->opt.nf, rtk->ssat[sat[i] - 1].sys);
936 slip |= rtk->ssat[sat[i] - 1].slip[f2i];
937 }
938 rejc = rtk->ssat[sat[i] - 1].rejc[k];
939 if (rtk->opt.modear == ARMODE_INST || (!(slip & 1) && rejc < 2)) {
940 continue;
941 }
942 rtk->x[j] = 0.0;
943 rtk->ssat[sat[i] - 1].rejc[k] = 0;
944 rtk->ssat[sat[i] - 1].lock[k] = -rtk->opt.minlock;
945 }
946 bias = zeros(ns, 1);
947
948 /* estimate approximate phase-bias by phase - code */
949 for (i = j = 0, offset = 0.0; i < ns; i++) {
950 if (rtk->opt.ionoopt != IONOOPT_IFLC) {
951 cp = sdobs(obs, iu[i], ir[i], k); /* cycle */
952 pr = sdobs(obs, iu[i], ir[i], k + NFREQ);
953 freqi = sat2freq(sat[i], obs[iu[i]].code[k], nav);
954 if (cp == 0.0 || pr == 0.0 || freqi == 0.0) {
955 continue;
956 }
957
958 bias[i] = cp - pr * freqi / CLIGHT;
959 } else {
960 int f2b = seliflc(rtk->opt.nf, rtk->ssat[sat[i] - 1].sys);
961 cp1 = sdobs(obs, iu[i], ir[i], 0);
962 cp2 = sdobs(obs, iu[i], ir[i], f2b);
963 pr1 = sdobs(obs, iu[i], ir[i], NFREQ);
964 pr2 = sdobs(obs, iu[i], ir[i], NFREQ + f2b);
965 freq1 = sat2freq(sat[i], obs[iu[i]].code[0], nav);
966 freq2 = sat2freq(sat[i], obs[iu[i]].code[f2b], nav);
967 if (cp1 == 0.0 || cp2 == 0.0 || pr1 == 0.0 || pr2 == 0.0 || freq1 == 0.0 || freq2 <= 0.0) {
968 continue;
969 }
970
971 C1 = SQR(freq1) / (SQR(freq1) - SQR(freq2));
972 C2 = -SQR(freq2) / (SQR(freq1) - SQR(freq2));
973 bias[i] = (C1 * cp1 * CLIGHT / freq1 + C2 * cp2 * CLIGHT / freq2) - (C1 * pr1 + C2 * pr2);
974 }
975 if (rtk->x[IB(sat[i], k, &rtk->opt)] != 0.0) {
976 offset += bias[i] - rtk->x[IB(sat[i], k, &rtk->opt)];
977 j++;
978 }
979 }
980 /* correct phase-bias offset to ensure phase-code coherency */
981 if (j > 0) {
982 for (i = 1; i <= MAXSAT; i++) {
983 if (rtk->x[IB(i, k, &rtk->opt)] != 0.0) {
984 rtk->x[IB(i, k, &rtk->opt)] += offset / j;
985 }
986 }
987 }
988 /* set initial states of phase-bias */
989 for (i = 0; i < ns; i++) {
990 if (bias[i] == 0.0 || rtk->x[IB(sat[i], k, &rtk->opt)] != 0.0) {
991 continue;
992 }
993 initx(rtk, bias[i], SQR(rtk->opt.std[0]), IB(sat[i], k, &rtk->opt));
994 }
995 free(bias);
996 }
997}
998/* temporal update of states --------------------------------------------------*/
999static void udstate(rtk_t* rtk, const obsd_t* obs, const int* sat, const int* iu, const int* ir, int ns,
1000 const nav_t* nav) {
1001 double tt = rtk->tt, bl, dr[3];
1002
1003 trace(NULL, 3, "udstate : ns=%d\n", ns);
1004
1005 /* temporal update of position/velocity/acceleration */
1006 udpos(rtk, tt);
1007
1008 /* temporal update of ionospheric parameters */
1009 if (rtk->opt.ionoopt >= IONOOPT_EST) {
1010 bl = baseline(rtk->x, rtk->rb, dr);
1011 udion(rtk, tt, bl, sat, ns);
1012 }
1013 /* temporal update of tropospheric parameters */
1014 if (rtk->opt.tropopt >= TROPOPT_EST) {
1015 udtrop(rtk, tt, bl);
1016 }
1017 /* temporal update of eceiver h/w bias */
1018 if (rtk->opt.glomodear == 2 && (rtk->opt.navsys & SYS_GLO)) {
1019 udrcvbias(rtk, tt);
1020 }
1021 /* temporal update of phase-bias */
1022 if (rtk->opt.mode > PMODE_DGPS) {
1023 udbias(rtk, tt, obs, sat, iu, ir, ns, nav);
1024 }
1025}
1026/* UD (undifferenced) phase/code residual for satellite ----------------------*/
1027static void zdres_sat(int base, double r, const obsd_t* obs, const nav_t* nav, const double* azel, const double* dant,
1028 const prcopt_t* opt, double* y, double* freq) {
1029 double freq1, freq2, C1, C2, dant_if;
1030 int i, nf = NF(opt), f2;
1031
1032 if (opt->ionoopt == IONOOPT_IFLC) { /* iono-free linear combination */
1033 f2 = seliflc(nf, satsys(obs->sat, NULL));
1034 freq1 = sat2freq(obs->sat, obs->code[0], nav);
1035 freq2 = sat2freq(obs->sat, obs->code[f2], nav);
1036 if (freq1 == 0.0 || freq2 == 0.0) {
1037 return;
1038 }
1039
1040 if (testsnr(base, 0, azel[1], obs->SNR[0] * SNR_UNIT, &opt->snrmask) ||
1041 testsnr(base, f2, azel[1], obs->SNR[f2] * SNR_UNIT, &opt->snrmask)) {
1042 return;
1043 }
1044
1045 C1 = SQR(freq1) / (SQR(freq1) - SQR(freq2));
1046 C2 = -SQR(freq2) / (SQR(freq1) - SQR(freq2));
1047 dant_if = C1 * dant[0] + C2 * dant[f2];
1048
1049 if (obs->L[0] != 0.0 && obs->L[f2] != 0.0) {
1050 y[0] = C1 * obs->L[0] * CLIGHT / freq1 + C2 * obs->L[f2] * CLIGHT / freq2 - r - dant_if;
1051 }
1052 if (obs->P[0] != 0.0 && obs->P[f2] != 0.0) {
1053 y[1] = C1 * obs->P[0] + C2 * obs->P[f2] - r - dant_if;
1054 }
1055 freq[0] = 1.0;
1056 } else {
1057 for (i = 0; i < nf; i++) {
1058 if ((freq[i] = sat2freq(obs->sat, obs->code[i], nav)) == 0.0) {
1059 continue;
1060 }
1061
1062 /* check SNR mask */
1063 if (testsnr(base, i, azel[1], obs->SNR[i] * SNR_UNIT, &opt->snrmask)) {
1064 continue;
1065 }
1066 /* residuals = observable - pseudorange */
1067 if (obs->L[i] != 0.0) {
1068 y[i] = obs->L[i] * CLIGHT / freq[i] - r - dant[i];
1069 }
1070 if (obs->P[i] != 0.0) {
1071 y[i + nf] = obs->P[i] - r - dant[i];
1072 }
1073 }
1074 }
1075}
1076/* UD (undifferenced) phase/code residuals -----------------------------------*/
1077static int zdres(int base, const obsd_t* obs, int n, const double* rs, const double* dts, const double* var,
1078 const int* svh, const nav_t* nav, const double* rr, const prcopt_t* opt, int index, double* y,
1079 double* e, double* azel, double* freq) {
1080 double r, rr_[3], pos[3], dant[NFREQPCV] = {0}, disp[3];
1081 double zhd, zazel[] = {0.0, 90.0 * D2R};
1082 int i, nf = NF(opt);
1083
1084 trace(NULL, 3, "zdres : n=%d\n", n);
1085
1086 for (i = 0; i < n * nf * 2; i++) {
1087 y[i] = 0.0;
1088 }
1089
1090 if (norm(rr, 3) <= 0.0) {
1091 return 0; /* no receiver position */
1092 }
1093
1094 for (i = 0; i < 3; i++) {
1095 rr_[i] = rr[i];
1096 }
1097
1098 /* earth tide correction */
1099 if (opt->tidecorr) {
1100 tidedisp(gpst2utc(obs[0].time), rr_, opt->tidecorr, &nav->erp, opt->odisp[base], disp);
1101 for (i = 0; i < 3; i++) {
1102 rr_[i] += disp[i];
1103 }
1104 }
1105 ecef2pos(rr_, pos);
1106
1107 for (i = 0; i < n; i++) {
1108 /* compute geometric-range and azimuth/elevation angle */
1109 if ((r = geodist(rs + i * 6, rr_, e + i * 3)) <= 0.0) {
1110 continue;
1111 }
1112 if (satazel(pos, e + i * 3, azel + i * 2) < opt->elmin) {
1113 continue;
1114 }
1115
1116 /* test elevation mask for rover */
1117 if (!base && testelmask(azel + i * 2, opt->elmaskopt)) {
1118 continue;
1119 }
1120
1121 /* excluded satellite? */
1122 if (satexclude(obs[i].sat, var[i], svh[i], opt)) {
1123 continue;
1124 }
1125
1126 /* satellite clock-bias */
1127 r += -CLIGHT * dts[i * 2];
1128
1129 /* troposphere delay model (hydrostatic) */
1130 zhd = tropmodel(obs[0].time, pos, zazel, 0.0);
1131 r += tropmapf(obs[i].time, pos, azel + i * 2, NULL) * zhd;
1132
1133 /* receiver antenna phase center correction */
1134 antmodel(opt->pcvr + index, opt->antdel[index], azel + i * 2, opt->posopt[1], dant);
1135
1136 /* UD phase/code residual for satellite */
1137 zdres_sat(base, r, obs + i, nav, azel + i * 2, dant, opt, y + i * nf * 2, freq + i * nf);
1138 }
1139 trace(NULL, 4, "rr_=%.3f %.3f %.3f\n", rr_[0], rr_[1], rr_[2]);
1140 trace(NULL, 4, "pos=%.9f %.9f %.3f\n", pos[0] * R2D, pos[1] * R2D, pos[2]);
1141 for (i = 0; i < n; i++) {
1142 trace(NULL, 4, "sat=%2d %13.3f %13.3f %13.3f %13.10f %6.1f %5.1f\n", obs[i].sat, rs[i * 6], rs[1 + i * 6],
1143 rs[2 + i * 6], dts[i * 2], azel[i * 2] * R2D, azel[1 + i * 2] * R2D);
1144 }
1145 trace(NULL, 4, "y=\n");
1146 tracemat(NULL, 4, y, nf * 2, n, 13, 3);
1147
1148 return 1;
1149}
1150/* test valid observation data -----------------------------------------------*/
1151static int validobs(int i, int j, int f, int nf, double* y) {
1152 /* if no phase observable, psudorange is also unusable */
1153 return y[f + i * nf * 2] != 0.0 && y[f + j * nf * 2] != 0.0 &&
1154 (f < nf || (y[f - nf + i * nf * 2] != 0.0 && y[f - nf + j * nf * 2] != 0.0));
1155}
1156/* DD (double-differenced) measurement error covariance ----------------------*/
1157static void ddcov(const int* nb, int n, const double* Ri, const double* Rj, int nv, double* R) {
1158 int i, j, k = 0, b;
1159
1160 trace(NULL, 3, "ddcov : n=%d\n", n);
1161
1162 for (i = 0; i < nv * nv; i++) {
1163 R[i] = 0.0;
1164 }
1165 for (b = 0; b < n; k += nb[b++]) {
1166 for (i = 0; i < nb[b]; i++) {
1167 for (j = 0; j < nb[b]; j++) {
1168 R[k + i + (k + j) * nv] = Ri[k + i] + (i == j ? Rj[k + i] : 0.0);
1169 }
1170 }
1171 }
1172 trace(NULL, 5, "R=\n");
1173 tracemat(NULL, 5, R, nv, nv, 8, 6);
1174}
1175/* baseline length constraint ------------------------------------------------*/
1176static int constbl(rtk_t* rtk, const double* x, const double* P, double* v, double* H, double* Ri, double* Rj,
1177 int index) {
1178 const double thres = 0.1; /* threshold for nonliearity (v.2.3.0) */
1179 double xb[3], b[3], bb, var = 0.0;
1180 int i;
1181
1182 trace(NULL, 3, "constbl : \n");
1183
1184 /* no constraint */
1185 if (rtk->opt.baseline[0] <= 0.0) {
1186 return 0;
1187 }
1188
1189 /* time-adjusted baseline vector and length */
1190 for (i = 0; i < 3; i++) {
1191 xb[i] = rtk->rb[i];
1192 b[i] = x[i] - xb[i];
1193 }
1194 bb = norm(b, 3);
1195
1196 /* approximate variance of solution */
1197 if (P) {
1198 for (i = 0; i < 3; i++) {
1199 var += P[i + i * rtk->nx];
1200 }
1201 var /= 3.0;
1202 }
1203 /* check nonlinearity */
1204 if (var > SQR(thres * bb)) {
1205 trace(NULL, 3, "constbl : equation nonlinear (bb=%.3f var=%.3f)\n", bb, var);
1206 return 0;
1207 }
1208 /* constraint to baseline length */
1209 v[index] = rtk->opt.baseline[0] - bb;
1210 if (H) {
1211 for (i = 0; i < 3; i++) {
1212 H[i + index * rtk->nx] = b[i] / bb;
1213 }
1214 }
1215 Ri[index] = 0.0;
1216 Rj[index] = SQR(rtk->opt.baseline[1]);
1217
1218 trace(NULL, 4, "baseline len v=%13.3f R=%8.6f %8.6f\n", v[index], Ri[index], Rj[index]);
1219
1220 return 1;
1221}
1222/* precise tropspheric model -------------------------------------------------*/
1223static double prectrop(gtime_t time, const double* pos, int r, const double* azel, const prcopt_t* opt, const double* x,
1224 double* dtdx) {
1225 double m_w = 0.0, cotz, grad_n, grad_e;
1226 int i = IT(r, opt);
1227
1228 /* wet mapping function */
1229 tropmapf(time, pos, azel, &m_w);
1230
1231 if (opt->tropopt >= TROPOPT_ESTG && azel[1] > 0.0) {
1232 /* m_w=m_0+m_0*cot(el)*(Gn*cos(az)+Ge*sin(az)): ref [6] */
1233 cotz = 1.0 / tan(azel[1]);
1234 grad_n = m_w * cotz * cos(azel[0]);
1235 grad_e = m_w * cotz * sin(azel[0]);
1236 m_w += grad_n * x[i + 1] + grad_e * x[i + 2];
1237 dtdx[1] = grad_n * x[i];
1238 dtdx[2] = grad_e * x[i];
1239 } else {
1240 dtdx[1] = dtdx[2] = 0.0;
1241 }
1242 dtdx[0] = m_w;
1243 return m_w * x[i];
1244}
1245/* test satellite system (m=0:GPS/SBS,1:GLO,2:GAL,3:BDS,4:QZS,5:IRN) ---------*/
1246static int test_sys(int sys, int m) {
1247 switch (sys) {
1248 case SYS_GPS:
1249 return m == 0;
1250 case SYS_SBS:
1251 return m == 0;
1252 case SYS_GLO:
1253 return m == 1;
1254 case SYS_GAL:
1255 return m == 2;
1256 case SYS_CMP:
1257 return m == 3;
1258 case SYS_QZS:
1259 return m == 4;
1260 case SYS_IRN:
1261 return m == 5;
1262 }
1263 return 0;
1264}
1265/* DD (double-differenced) phase/code residuals ------------------------------*/
1266static int ddres(rtk_t* rtk, const nav_t* nav, const obsd_t* obs, double dt, const double* x, const double* P,
1267 const int* sat, double* y, double* e, double* azel, double* freq, const int* iu, const int* ir, int ns,
1268 double* v, double* H, double* R, int* vflg) {
1269 prcopt_t* opt = &rtk->opt;
1270 double bl, dr[3], posu[3], posr[3], didxi = 0.0, didxj = 0.0, *im;
1271 double *tropr, *tropu, *dtdxr, *dtdxu, *Ri, *Rj, freqi, freqj, *Hi = NULL;
1272 int i, j, k, m, f, nv = 0, nb[NFREQ * 4 * 2 + 2] = {0}, b = 0, sysi, sysj, nf = NF(opt);
1273
1274 trace(NULL, 3, "ddres : dt=%.1f nx=%d ns=%d\n", dt, rtk->nx, ns);
1275
1276 bl = baseline(x, rtk->rb, dr);
1277 ecef2pos(x, posu);
1278 ecef2pos(rtk->rb, posr);
1279
1280 Ri = mat(ns * nf * 2 + 2, 1);
1281 Rj = mat(ns * nf * 2 + 2, 1);
1282 im = mat(ns, 1);
1283 tropu = mat(ns, 1);
1284 tropr = mat(ns, 1);
1285 dtdxu = mat(ns, 3);
1286 dtdxr = mat(ns, 3);
1287
1288 for (i = 0; i < MAXSAT; i++) {
1289 for (j = 0; j < NFREQ; j++) {
1290 rtk->ssat[i].resp[j] = rtk->ssat[i].resc[j] = 0.0;
1291 }
1292 }
1293 /* compute factors of ionospheric and tropospheric delay */
1294 for (i = 0; i < ns; i++) {
1295 if (opt->ionoopt >= IONOOPT_EST) {
1296 im[i] = (ionmapf(posu, azel + iu[i] * 2) + ionmapf(posr, azel + ir[i] * 2)) / 2.0;
1297 }
1298 if (opt->tropopt >= TROPOPT_EST) {
1299 tropu[i] = prectrop(rtk->sol.time, posu, 0, azel + iu[i] * 2, opt, x, dtdxu + i * 3);
1300 tropr[i] = prectrop(rtk->sol.time, posr, 1, azel + ir[i] * 2, opt, x, dtdxr + i * 3);
1301 }
1302 }
1303 for (m = 0; m < 6; m++) { /* m=0:GPS/SBS,1:GLO,2:GAL,3:BDS,4:QZS,5:IRN */
1304
1305 for (f = opt->mode > PMODE_DGPS ? 0 : nf; f < nf * 2; f++) {
1306 /* search reference satellite with highest elevation */
1307 for (i = -1, j = 0; j < ns; j++) {
1308 sysi = rtk->ssat[sat[j] - 1].sys;
1309 if (!test_sys(sysi, m)) {
1310 continue;
1311 }
1312 if (!validobs(iu[j], ir[j], f, nf, y)) {
1313 continue;
1314 }
1315 if (i < 0 || azel[1 + iu[j] * 2] >= azel[1 + iu[i] * 2]) {
1316 i = j;
1317 }
1318 }
1319 if (i < 0) {
1320 continue;
1321 }
1322
1323 /* make DD (double difference) */
1324 for (j = 0; j < ns; j++) {
1325 if (i == j) {
1326 continue;
1327 }
1328 sysi = rtk->ssat[sat[i] - 1].sys;
1329 sysj = rtk->ssat[sat[j] - 1].sys;
1330 freqi = freq[f % nf + iu[i] * nf];
1331 freqj = freq[f % nf + iu[j] * nf];
1332 if (!test_sys(sysj, m)) {
1333 continue;
1334 }
1335 if (!validobs(iu[j], ir[j], f, nf, y)) {
1336 continue;
1337 }
1338
1339 if (H) {
1340 Hi = H + nv * rtk->nx;
1341 for (k = 0; k < rtk->nx; k++) {
1342 Hi[k] = 0.0;
1343 }
1344 }
1345 /* DD residual */
1346 v[nv] =
1347 (y[f + iu[i] * nf * 2] - y[f + ir[i] * nf * 2]) - (y[f + iu[j] * nf * 2] - y[f + ir[j] * nf * 2]);
1348
1349 /* partial derivatives by rover position */
1350 if (H) {
1351 for (k = 0; k < 3; k++) {
1352 Hi[k] = -e[k + iu[i] * 3] + e[k + iu[j] * 3];
1353 }
1354 }
1355 /* DD ionospheric delay term */
1356 if (opt->ionoopt == IONOOPT_EST) {
1357 didxi = (f < nf ? -1.0 : 1.0) * im[i] * SQR(FREQ1 / freqi);
1358 didxj = (f < nf ? -1.0 : 1.0) * im[j] * SQR(FREQ1 / freqj);
1359 v[nv] -= didxi * x[II(sat[i], opt)] - didxj * x[II(sat[j], opt)];
1360 if (H) {
1361 Hi[II(sat[i], opt)] = didxi;
1362 Hi[II(sat[j], opt)] = -didxj;
1363 }
1364 }
1365 /* DD tropospheric delay term */
1366 if (opt->tropopt == TROPOPT_EST || opt->tropopt == TROPOPT_ESTG) {
1367 v[nv] -= (tropu[i] - tropu[j]) - (tropr[i] - tropr[j]);
1368 for (k = 0; k < (opt->tropopt < TROPOPT_ESTG ? 1 : 3); k++) {
1369 if (!H) {
1370 continue;
1371 }
1372 Hi[IT(0, opt) + k] = (dtdxu[k + i * 3] - dtdxu[k + j * 3]);
1373 Hi[IT(1, opt) + k] = -(dtdxr[k + i * 3] - dtdxr[k + j * 3]);
1374 }
1375 }
1376 /* DD phase-bias term */
1377 if (f < nf) {
1378 if (opt->ionoopt != IONOOPT_IFLC) {
1379 v[nv] -= CLIGHT / freqi * x[IB(sat[i], f, opt)] - CLIGHT / freqj * x[IB(sat[j], f, opt)];
1380 if (H) {
1381 Hi[IB(sat[i], f, opt)] = CLIGHT / freqi;
1382 Hi[IB(sat[j], f, opt)] = -CLIGHT / freqj;
1383 }
1384 } else {
1385 v[nv] -= x[IB(sat[i], f, opt)] - x[IB(sat[j], f, opt)];
1386 if (H) {
1387 Hi[IB(sat[i], f, opt)] = 1.0;
1388 Hi[IB(sat[j], f, opt)] = -1.0;
1389 }
1390 }
1391 }
1392 if (f < nf) {
1393 rtk->ssat[sat[j] - 1].resc[f] = v[nv];
1394 } else {
1395 rtk->ssat[sat[j] - 1].resp[f - nf] = v[nv];
1396 }
1397
1398 /* test innovation */
1399 if (opt->maxinno > 0.0) {
1400 /* open threshold 10x if either phase-bias just initialized */
1401 double threshadj = 1.0;
1402 if (f < nf && opt->mode > PMODE_DGPS && P) {
1403 int ii = IB(sat[i], f, opt), jj = IB(sat[j], f, opt);
1404 if (P[ii + rtk->nx * ii] == SQR(opt->std[0]) || P[jj + rtk->nx * jj] == SQR(opt->std[0])) {
1405 threshadj = 10.0;
1406 }
1407 }
1408 if (fabs(v[nv]) > opt->maxinno * threshadj) {
1409 /* only count phase outliers: rejc[] drives phase-bias resets
1410 in udbias(); code outliers must not trigger phase resets */
1411 if (f < nf) {
1412 rtk->ssat[sat[j] - 1].rejc[f]++;
1413 }
1414 errmsg(rtk, "outlier rejected (sat=%3d-%3d %s%d v=%.3f)\n", sat[i], sat[j], f < nf ? "L" : "P",
1415 f % nf + 1, v[nv]);
1416 continue;
1417 }
1418 }
1419 /* SD (single-differenced) measurement error variances */
1420 {
1421 int frq = f % nf, code = (f < nf) ? 0 : 1;
1422 Ri[nv] = varerr(sat[i], sysi, azel[1 + iu[i] * 2], obs[iu[i]].SNR[frq] * SNR_UNIT,
1423 obs[ir[i]].SNR[frq] * SNR_UNIT, bl, dt, f, opt);
1424 Rj[nv] = varerr(sat[j], sysj, azel[1 + iu[j] * 2], obs[iu[j]].SNR[frq] * SNR_UNIT,
1425 obs[ir[j]].SNR[frq] * SNR_UNIT, bl, dt, f, opt);
1426 /* inflate variance for half-cycle uncertain phase observations */
1427 if (!code && (obs[iu[i]].LLI[frq] & LLI_HALFC)) {
1428 Ri[nv] += 0.01;
1429 }
1430 if (!code && (obs[iu[j]].LLI[frq] & LLI_HALFC)) {
1431 Rj[nv] += 0.01;
1432 }
1433
1434 /* set valid data flags: use code DD to flag valid (E5 reverted:
1435 demo5 phase-only vsat caused nfix=0 deadlock in urban canyons) */
1436 if (code) {
1437 rtk->ssat[sat[i] - 1].vsat[frq] = rtk->ssat[sat[j] - 1].vsat[frq] = 1;
1438 }
1439 }
1440 trace(NULL, 4, "sat=%3d-%3d %s%d v=%13.3f R=%8.6f %8.6f\n", sat[i], sat[j], f < nf ? "L" : "P",
1441 f % nf + 1, v[nv], Ri[nv], Rj[nv]);
1442
1443 vflg[nv++] = (sat[i] << 16) | (sat[j] << 8) | ((f < nf ? 0 : 1) << 4) | (f % nf);
1444 nb[b]++;
1445 }
1446 b++;
1447 }
1448 }
1449 /* end of system loop */
1450
1451 /* baseline length constraint for moving baseline */
1452 if (opt->mode == PMODE_MOVEB && constbl(rtk, x, P, v, H, Ri, Rj, nv)) {
1453 vflg[nv++] = 3 << 4;
1454 nb[b++]++;
1455 }
1456 if (H) {
1457 trace(NULL, 5, "H=\n");
1458 tracemat(NULL, 5, H, rtk->nx, nv, 7, 4);
1459 }
1460
1461 /* DD measurement error covariance */
1462 ddcov(nb, b, Ri, Rj, nv, R);
1463
1464 free(Ri);
1465 free(Rj);
1466 free(im);
1467 free(tropu);
1468 free(tropr);
1469 free(dtdxu);
1470 free(dtdxr);
1471
1472 return nv;
1473}
1474/* time-interpolation of residuals (for post-processing) ---------------------*/
1475static double intpres(gtime_t time, const obsd_t* obs, int n, const nav_t* nav, rtk_t* rtk, double* y) {
1476 static obsd_t obsb[MAXOBS];
1477 static double yb[MAXOBS * NFREQ * 2], rs[MAXOBS * 6], dts[MAXOBS * 2], var[MAXOBS];
1478 static double e[MAXOBS * 3], azel[MAXOBS * 2], freq[MAXOBS * NFREQ];
1479 static int nb = 0, svh[MAXOBS * 2];
1480 prcopt_t* opt = &rtk->opt;
1481 double tt = timediff(time, obs[0].time), ttb, *p, *q;
1482 int i, j, k, nf = NF(opt);
1483
1484 trace(NULL, 3, "intpres : n=%d tt=%.1f\n", n, tt);
1485
1486 if (nb == 0 || fabs(tt) < DTTOL) {
1487 nb = n;
1488 for (i = 0; i < n; i++) {
1489 obsb[i] = obs[i];
1490 }
1491 return tt;
1492 }
1493 ttb = timediff(time, obsb[0].time);
1494 if (fabs(ttb) > opt->maxtdiff * 2.0 || ttb == tt) {
1495 return tt;
1496 }
1497
1498 satposs(time, obsb, nb, nav, opt->sateph, rs, dts, var, svh);
1499
1500 if (!zdres(1, obsb, nb, rs, dts, var, svh, nav, rtk->rb, opt, 1, yb, e, azel, freq)) {
1501 return tt;
1502 }
1503 for (i = 0; i < n; i++) {
1504 for (j = 0; j < nb; j++) {
1505 if (obsb[j].sat == obs[i].sat) {
1506 break;
1507 }
1508 }
1509 if (j >= nb) {
1510 continue;
1511 }
1512 for (k = 0, p = y + i * nf * 2, q = yb + j * nf * 2; k < nf * 2; k++, p++, q++) {
1513 if (*p == 0.0 || *q == 0.0) {
1514 *p = 0.0;
1515 } else {
1516 *p = (ttb * (*p) - tt * (*q)) / (ttb - tt);
1517 }
1518 }
1519 }
1520 return fabs(ttb) > fabs(tt) ? ttb : tt;
1521}
1522/* index for SD to DD transformation matrix D --------------------------------*/
1523static int ddidx(rtk_t* rtk, int* ix, int gps, int glo, int sbs) {
1524 int i, j, k, m, f, n, nb = 0, na = rtk->na, nf = NF(&rtk->opt), nofix;
1525
1526 trace(NULL, 3, "ddidx : gps=%d glo=%d sbs=%d\n", gps, glo, sbs);
1527
1528 /* clear fix flag for all sats (1=float, 2=fix) */
1529 for (i = 0; i < MAXSAT; i++) {
1530 for (j = 0; j < NFREQ; j++) {
1531 rtk->ssat[i].fix[j] = 0;
1532 }
1533 }
1534 for (m = 0; m < 6; m++) { /* m=0:GPS/SBS,1:GLO,2:GAL,3:BDS,4:QZS,5:IRN */
1535
1536 /* skip if ambiguity resolution turned off for this constellation */
1537 nofix = (m == 0 && gps == 0) || (m == 1 && glo == 0) || (m == 3 && rtk->opt.bdsmodear == 0);
1538
1539 for (f = 0, k = na; f < nf; f++, k += MAXSAT) {
1540 /* find reference satellite: must be valid, locked, no half-cycle,
1541 above elevation mask, and not in a no-fix constellation */
1542 for (i = k; i < k + MAXSAT; i++) {
1543 if (rtk->x[i] == 0.0 || !test_sys(rtk->ssat[i - k].sys, m) || !rtk->ssat[i - k].vsat[f]) {
1544 continue;
1545 }
1546 if (rtk->ssat[i - k].lock[f] >= 0 && !(rtk->ssat[i - k].slip[f] & 2) &&
1547 rtk->ssat[i - k].azel[1] >= rtk->opt.elmaskar && !nofix) {
1548 rtk->ssat[i - k].fix[f] = 2; /* fix */
1549 break;
1550 } else {
1551 rtk->ssat[i - k].fix[f] = 1;
1552 }
1553 }
1554 /* no valid reference satellite found for this freq/constellation */
1555 if (i >= k + MAXSAT || rtk->ssat[i - k].fix[f] != 2) {
1556 continue;
1557 }
1558
1559 /* build DD pairs using found reference satellite */
1560 for (n = 0, j = k; j < k + MAXSAT; j++) {
1561 if (i == j || rtk->x[j] == 0.0 || !test_sys(rtk->ssat[j - k].sys, m) || !rtk->ssat[j - k].vsat[f]) {
1562 continue;
1563 }
1564 /* skip SBAS sats if sbs flag not set */
1565 if (sbs == 0 && satsys(j - k + 1, NULL) == SYS_SBS) {
1566 continue;
1567 }
1568 if (rtk->ssat[j - k].lock[f] >= 0 && !(rtk->ssat[j - k].slip[f] & 2) && rtk->ssat[j - k].vsat[f] &&
1569 rtk->ssat[j - k].azel[1] >= rtk->opt.elmaskar && !nofix) {
1570 ix[nb * 2] = i; /* state index of ref bias */
1571 ix[nb * 2 + 1] = j; /* state index of target bias */
1572 nb++;
1573 rtk->ssat[j - k].fix[f] = 2; /* fix */
1574 n++;
1575 } else {
1576 rtk->ssat[j - k].fix[f] = 1;
1577 }
1578 }
1579 /* don't use ref sat if no valid pairs found */
1580 if (n == 0) {
1581 rtk->ssat[i - k].fix[f] = 1;
1582 }
1583 }
1584 }
1585 return nb;
1586}
1587/* restore SD (single-differenced) ambiguity ---------------------------------*/
1588static void restamb(rtk_t* rtk, const double* bias, double* xa) {
1589 int i, n, m, f, index[MAXSAT], nv = 0, nf = NF(&rtk->opt);
1590
1591 trace(NULL, 3, "restamb :\n");
1592
1593 for (i = 0; i < rtk->nx; i++) {
1594 xa[i] = rtk->x[i];
1595 }
1596 for (i = 0; i < rtk->na; i++) {
1597 xa[i] = rtk->xa[i];
1598 }
1599
1600 for (m = 0; m < 5; m++) {
1601 for (f = 0; f < nf; f++) {
1602 for (n = i = 0; i < MAXSAT; i++) {
1603 if (!test_sys(rtk->ssat[i].sys, m) || rtk->ssat[i].fix[f] != 2) {
1604 continue;
1605 }
1606 index[n++] = IB(i + 1, f, &rtk->opt);
1607 }
1608 if (n < 2) {
1609 continue;
1610 }
1611
1612 xa[index[0]] = rtk->x[index[0]];
1613
1614 for (i = 1; i < n; i++) {
1615 xa[index[i]] = xa[index[0]] - bias[nv++];
1616 }
1617 }
1618 }
1619}
1620/* hold integer ambiguity ----------------------------------------------------*/
1621static void holdamb(rtk_t* rtk, const double* xa) {
1622 double *v, *H, *R;
1623 int i, n, m, f, info, index[MAXSAT], nb = rtk->nx - rtk->na, nv = 0, nf = NF(&rtk->opt);
1624
1625 trace(NULL, 3, "holdamb :\n");
1626
1627 v = mat(nb, 1);
1628 H = zeros(rtk->nx, nb);
1629
1630 for (m = 0; m < 6; m++) {
1631 for (f = 0; f < nf; f++) {
1632 for (n = i = 0; i < MAXSAT; i++) {
1633 if (!test_sys(rtk->ssat[i].sys, m) || rtk->ssat[i].fix[f] != 2 ||
1634 rtk->ssat[i].azel[1] < rtk->opt.elmaskhold) {
1635 continue;
1636 }
1637 index[n++] = IB(i + 1, f, &rtk->opt);
1638 rtk->ssat[i].fix[f] = 3; /* hold */
1639 }
1640 /* constraint to fixed ambiguity */
1641 for (i = 1; i < n; i++) {
1642 v[nv] = (xa[index[0]] - xa[index[i]]) - (rtk->x[index[0]] - rtk->x[index[i]]);
1643
1644 H[index[0] + nv * rtk->nx] = 1.0;
1645 H[index[i] + nv * rtk->nx] = -1.0;
1646 nv++;
1647 }
1648 }
1649 }
1650 /* skip hold if not enough sats (only enforced in fix-and-hold mode) */
1651 if (rtk->opt.modear == ARMODE_FIXHOLD && rtk->opt.minholdsats > 0 && nv < rtk->opt.minholdsats) {
1652 trace(NULL, 3, "holdamb: not enough sats to hold ambiguity (nv=%d)\n", nv);
1653 free(v);
1654 free(H);
1655 return;
1656 }
1657 rtk->holdamb = 1; /* record that hold has occurred at least once */
1658
1659 if (nv > 0) {
1660 R = zeros(nv, nv);
1661 for (i = 0; i < nv; i++) {
1662 R[i + i * nv] = VAR_HOLDAMB;
1663 }
1664
1665 /* update states with constraints */
1666 if ((info = filter(rtk->x, rtk->P, H, v, R, rtk->nx, nv))) {
1667 errmsg(rtk, "filter error (info=%d)\n", info);
1668 }
1669 free(R);
1670 }
1671 free(v);
1672 free(H);
1673}
1674/* resolve integer ambiguity by LAMBDA for one constellation set --------------*/
1675static int resamb_LAMBDA(rtk_t* rtk, double* bias, double* xa, int gps, int glo, int sbs) {
1676 prcopt_t* opt = &rtk->opt;
1677 int i, j, nb, info, nx = rtk->nx, na = rtk->na;
1678 double *DP, *y, *b, *db, *Qb, *Qab, *QQ, s[2];
1679 int* ix;
1680
1681 trace(NULL, 3, "resamb_LAMBDA : nx=%d\n", nx);
1682
1683 rtk->sol.ratio = 0.0;
1684 rtk->nb_ar = 0;
1685
1686 /* create index of SD→DD transformation matrix */
1687 ix = imat(nx, 2);
1688 if ((nb = ddidx(rtk, ix, gps, glo, sbs)) < (opt->minfixsats > 0 ? opt->minfixsats - 1 : 0)) {
1689 errmsg(rtk, "not enough valid double-differences\n");
1690 free(ix);
1691 return -1; /* flag abort — insufficient pairs */
1692 }
1693 rtk->nb_ar = nb;
1694
1695 y = mat(nb, 1);
1696 DP = mat(nb, nx - na);
1697 b = mat(nb, 2);
1698 db = mat(nb, 1);
1699 Qb = mat(nb, nb);
1700 Qab = mat(na, nb);
1701 QQ = mat(na, nb);
1702
1703 /* y=D*xc, Qb=D*Qc*D', Qab=Qac*D' */
1704 for (i = 0; i < nb; i++) {
1705 y[i] = rtk->x[ix[i * 2]] - rtk->x[ix[i * 2 + 1]];
1706 }
1707 for (j = 0; j < nx - na; j++) {
1708 for (i = 0; i < nb; i++) {
1709 DP[i + j * nb] = rtk->P[ix[i * 2] + (na + j) * nx] - rtk->P[ix[i * 2 + 1] + (na + j) * nx];
1710 }
1711 }
1712 for (j = 0; j < nb; j++) {
1713 for (i = 0; i < nb; i++) {
1714 Qb[i + j * nb] = DP[i + (ix[j * 2] - na) * nb] - DP[i + (ix[j * 2 + 1] - na) * nb];
1715 }
1716 }
1717 for (j = 0; j < nb; j++) {
1718 for (i = 0; i < na; i++) {
1719 Qab[i + j * na] = rtk->P[i + ix[j * 2] * nx] - rtk->P[i + ix[j * 2 + 1] * nx];
1720 }
1721 }
1722 /* LAMBDA/MLAMBDA ILS (integer least-square) estimation */
1723 if (!(info = lambda(nb, 2, y, Qb, b, s))) {
1724 trace(NULL, 4, "N(1)=");
1725 tracemat(NULL, 4, b, 1, nb, 10, 3);
1726 trace(NULL, 4, "N(2)=");
1727 tracemat(NULL, 4, b + nb, 1, nb, 10, 3);
1728
1729 rtk->sol.ratio = s[0] > 0 ? (float)(s[1] / s[0]) : 0.0f;
1730 if (rtk->sol.ratio > 999.9) {
1731 rtk->sol.ratio = 999.9f;
1732 }
1733
1734 /* adaptive AR ratio threshold: polynomial in 1/(nb+1) if min≠max */
1735 if (opt->thresar[5] != opt->thresar[6]) {
1736 int nb1 = nb < 50 ? nb : 50; /* poly fitted for up to 50 sat pairs */
1737 double coeff[3];
1738 /* evaluate each polynomial coefficient at thresar[0] */
1739 for (int ki = 0; ki < 3; ki++) {
1740 coeff[ki] = ar_poly_coeffs[ki][0];
1741 for (int kj = 1; kj < 5; kj++) {
1742 coeff[ki] = coeff[ki] * opt->thresar[0] + ar_poly_coeffs[ki][kj];
1743 }
1744 }
1745 /* evaluate threshold polynomial in 1/(nb1+1) */
1746 rtk->sol.thres = (float)coeff[0];
1747 for (int ki = 1; ki < 3; ki++) {
1748 rtk->sol.thres = rtk->sol.thres * 1.0f / (float)(nb1 + 1) + (float)coeff[ki];
1749 }
1750 /* clamp to [thresar[5], thresar[6]] */
1751 if (rtk->sol.thres < (float)opt->thresar[5]) {
1752 rtk->sol.thres = (float)opt->thresar[5];
1753 }
1754 if (rtk->sol.thres > (float)opt->thresar[6]) {
1755 rtk->sol.thres = (float)opt->thresar[6];
1756 }
1757 } else {
1758 rtk->sol.thres = (float)opt->thresar[0];
1759 }
1760
1761 /* validation by popular ratio-test */
1762 if (s[0] <= 0.0 || s[1] / s[0] >= rtk->sol.thres) {
1763 /* transform float to fixed solution (xa=x-Qab*Qb\‍(b0-b)) */
1764 for (i = 0; i < na; i++) {
1765 rtk->xa[i] = rtk->x[i];
1766 for (j = 0; j < na; j++) {
1767 rtk->Pa[i + j * na] = rtk->P[i + j * nx];
1768 }
1769 }
1770 for (i = 0; i < nb; i++) {
1771 bias[i] = b[i];
1772 y[i] -= b[i];
1773 }
1774 if (!matinv(Qb, nb)) {
1775 matmul("NN", nb, 1, nb, 1.0, Qb, y, 0.0, db);
1776 matmul("NN", na, 1, nb, -1.0, Qab, db, 1.0, rtk->xa);
1777
1778 /* covariance of fixed solution (Qa=Qa-Qab*Qb^-1*Qab') */
1779 matmul("NN", na, nb, nb, 1.0, Qab, Qb, 0.0, QQ);
1780 matmul("NT", na, na, nb, -1.0, QQ, Qab, 1.0, rtk->Pa);
1781
1782 trace(NULL, 3, "resamb : validation ok (nb=%d ratio=%.2f thresh=%.2f s=%.2f/%.2f)\n", nb,
1783 s[0] == 0.0 ? 0.0 : s[1] / s[0], rtk->sol.thres, s[0], s[1]);
1784
1785 /* restore SD ambiguity */
1786 restamb(rtk, bias, xa);
1787 } else {
1788 nb = 0;
1789 }
1790 } else { /* validation failed */
1791 errmsg(rtk, "ambiguity validation failed (nb=%d ratio=%.2f thresh=%.2f s=%.2f/%.2f)\n", nb, s[1] / s[0],
1792 rtk->sol.thres, s[0], s[1]);
1793 nb = 0;
1794 }
1795 } else {
1796 errmsg(rtk, "lambda error (info=%d)\n", info);
1797 nb = 0;
1798 }
1799 free(ix);
1800 free(y);
1801 free(DP);
1802 free(b);
1803 free(db);
1804 free(Qb);
1805 free(Qab);
1806 free(QQ);
1807
1808 return nb; /* number of ambiguities (0=failed, -1=aborted) */
1809}
1810/* resolve integer ambiguity with partial AR and multi-attempt strategy -------*/
1811static int manage_amb_LAMBDA(rtk_t* rtk, double* bias, double* xa, const int* sat, int nf, int ns) {
1812 int gps1 = -1, glo1 = -1, sbas1 = -1, gps2, glo2, sbas2, nb, rerun, dly;
1813 float ratio1, posvar = 0.0f;
1814
1815 /* compute mean diagonal position variance to guard against false fix */
1816 for (int i = 0; i < 3; i++) {
1817 posvar += rtk->P[i + i * rtk->nx];
1818 }
1819 posvar /= 3.0f;
1820
1821 trace(NULL, 3, "manage_amb_LAMBDA: posvar=%.6f prevR=%.3f/%.3f nb_ar=%d\n", posvar, rtk->sol.prev_ratio1,
1822 rtk->sol.prev_ratio2, rtk->nb_ar);
1823
1824 /* skip AR if basic criteria not met */
1825 if (rtk->opt.mode <= PMODE_DGPS || rtk->opt.modear == ARMODE_OFF || rtk->opt.thresar[0] < 1.0 ||
1826 (rtk->opt.thresar[1] > 0.0 && posvar > rtk->opt.thresar[1])) {
1827 trace(NULL, 3, "manage_amb_LAMBDA: skip AR\n");
1828 rtk->sol.ratio = 0.0f;
1829 rtk->sol.prev_ratio1 = rtk->sol.prev_ratio2 = 0.0f;
1830 rtk->nb_ar = 0;
1831 return 0;
1832 }
1833
1834 /* partial AR: exclude next weakest satellite if still not fixed */
1835 int lockc[NFREQ], excsat = 0;
1836 if (rtk->opt.mindropsats > 0 && rtk->sol.prev_ratio2 < rtk->sol.thres && rtk->nb_ar >= rtk->opt.mindropsats) {
1837 /* find position of last excluded sat, then step to next */
1838 int i = 0;
1839 if (rtk->excsat != 0) {
1840 for (; i < ns; i++) {
1841 if (rtk->excsat == sat[i]) {
1842 i++;
1843 break;
1844 }
1845 }
1846 if (i >= ns) {
1847 i = 0; /* wrap if not found */
1848 }
1849 }
1850 /* find next sat currently used in AR */
1851 for (; i < ns; i++) {
1852 for (int f = 0; f < nf; f++) {
1853 if (rtk->ssat[sat[i] - 1].vsat[f] && rtk->ssat[sat[i] - 1].lock[f] >= 0 &&
1854 rtk->ssat[sat[i] - 1].azel[1] >= rtk->opt.elmin) {
1855 excsat = sat[i];
1856 break;
1857 }
1858 }
1859 if (excsat) {
1860 break;
1861 }
1862 }
1863 if (excsat) {
1864 for (int f = 0; f < nf; f++) {
1865 lockc[f] = rtk->ssat[excsat - 1].lock[f];
1866 /* temporarily mark sat as newly-locked to exclude from AR */
1867 rtk->ssat[excsat - 1].lock[f] = -rtk->nb_ar;
1868 }
1869 trace(NULL, 3, "manage_amb_LAMBDA: exclude sat %d\n", excsat);
1870 }
1871 rtk->excsat = excsat;
1872 }
1873
1874 /* initial AR pass: always enable GPS; GLO depends on fix-and-hold state */
1875 gps1 = 1;
1876 glo1 = (rtk->opt.navsys & SYS_GLO) ? ((rtk->opt.glomodear == GLO_ARMODE_FIXHOLD && !rtk->holdamb) ? 0 : 1) : 0;
1877 sbas1 = (rtk->opt.navsys & SYS_SBS) ? 1 : 0; /* independent of GLO state */
1878 nb = resamb_LAMBDA(rtk, bias, xa, gps1, glo1, sbas1);
1879 ratio1 = rtk->sol.ratio;
1880
1881 /* AR filter: if ratio dropped badly compared to previous epoch,
1882 remove any newly-locked satellites and rerun (Phase 2A: arfilter=1) */
1883 if (rtk->opt.arfilter) {
1884 rerun = 0;
1885 if (nb >= 0 && rtk->sol.prev_ratio2 >= rtk->sol.thres &&
1886 ((rtk->sol.ratio < rtk->sol.thres) ||
1887 (rtk->sol.ratio < rtk->opt.thresar[0] * 1.1 && rtk->sol.ratio < rtk->sol.prev_ratio1 / 2.0))) {
1888 trace(NULL, 3, "manage_amb_LAMBDA: low ratio, check new sats\n");
1889 dly = 2;
1890 for (int i = 0; i < ns; i++) {
1891 for (int f = 0; f < nf; f++) {
1892 if (rtk->ssat[sat[i] - 1].fix[f] != 2) {
1893 continue;
1894 }
1895 if (rtk->ssat[sat[i] - 1].lock[f] == 0) {
1896 trace(NULL, 3, "manage_amb_LAMBDA: remove sat %d:%d\n", sat[i], f);
1897 rtk->ssat[sat[i] - 1].lock[f] = -rtk->opt.minlock - dly;
1898 dly += 2;
1899 rerun = 1;
1900 }
1901 }
1902 }
1903 }
1904 if (rerun) {
1905 trace(NULL, 3, "manage_amb_LAMBDA: rerun with new sats removed\n");
1906 nb = resamb_LAMBDA(rtk, bias, xa, gps1, glo1, sbas1);
1907 }
1908 }
1909 rtk->sol.prev_ratio1 = ratio1;
1910
1911 /* GLO fix-and-hold second pass: re-run with final GLO settings if needed */
1912 if ((rtk->opt.navsys & SYS_GLO) && rtk->opt.glomodear == GLO_ARMODE_FIXHOLD && rtk->sol.ratio < rtk->sol.thres) {
1913 glo2 = sbas2 = 0;
1914 /* turn off GPS if gpsmodear==0 and already have good fix (debug use) */
1915 gps2 = (rtk->opt.gpsmodear == 0 && rtk->sol.ratio >= rtk->sol.thres) ? 0 : 1;
1916 if (glo1 != glo2 || gps1 != gps2) {
1917 nb = resamb_LAMBDA(rtk, bias, xa, gps2, glo2, sbas2);
1918 }
1919 }
1920
1921 /* restore excluded satellite if still no fix or ratio didn't improve */
1922 if (excsat && (rtk->sol.ratio < rtk->sol.thres) && (rtk->sol.ratio < 1.5f * rtk->sol.prev_ratio2)) {
1923 for (int f = 0; f < nf; f++) {
1924 rtk->ssat[excsat - 1].lock[f] = lockc[f];
1925 }
1926 trace(NULL, 3, "manage_amb_LAMBDA: restore sat %d\n", excsat);
1927 }
1928
1929 rtk->sol.prev_ratio1 = ratio1 > 0.0f ? ratio1 : rtk->sol.ratio;
1930 rtk->sol.prev_ratio2 = rtk->sol.ratio;
1931
1932 return nb;
1933}
1934/* validation of solution ----------------------------------------------------*/
1935static int valpos(rtk_t* rtk, const double* v, const double* R, const int* vflg, int nv, double thres) {
1936 double fact = thres * thres;
1937 int i, stat = 1, sat1, sat2, type, freq;
1938 char* stype;
1939
1940 trace(NULL, 3, "valpos : nv=%d thres=%.1f\n", nv, thres);
1941
1942 /* post-fit residual test */
1943 for (i = 0; i < nv; i++) {
1944 if (v[i] * v[i] <= fact * R[i + i * nv]) {
1945 continue;
1946 }
1947 sat1 = (vflg[i] >> 16) & 0xFF;
1948 sat2 = (vflg[i] >> 8) & 0xFF;
1949 type = (vflg[i] >> 4) & 0xF;
1950 freq = vflg[i] & 0xF;
1951 stype = type == 0 ? "L" : (type == 1 ? "L" : "C");
1952 errmsg(rtk, "large residual (sat=%2d-%2d %s%d v=%6.3f sig=%.3f)\n", sat1, sat2, stype, freq + 1, v[i],
1953 SQRT(R[i + i * nv]));
1954 }
1955 return stat;
1956}
1957/* relative positioning ------------------------------------------------------*/
1958static int relpos(rtk_t* rtk, const obsd_t* obs, int nu, int nr, const nav_t* nav) {
1959 prcopt_t* opt = &rtk->opt;
1960 gtime_t time = obs[0].time;
1961 double *rs, *dts, *var, *y, *e, *azel, *freq, *v, *H, *R, *xp, *Pp, *xa, *bias, dt;
1962 int i, j, f, n = nu + nr, ns, ny, nv, sat[MAXSAT], iu[MAXSAT], ir[MAXSAT], niter;
1963 int info, vflg[MAXOBS * NFREQ * 2 + 1], svh[MAXOBS * 2];
1964 int stat = rtk->opt.mode <= PMODE_DGPS ? SOLQ_DGPS : SOLQ_FLOAT;
1965 int nf = opt->ionoopt == IONOOPT_IFLC ? 1 : opt->nf;
1966
1967 trace(NULL, 3, "relpos : nx=%d nu=%d nr=%d\n", rtk->nx, nu, nr);
1968
1969 dt = timediff(time, obs[nu].time);
1970
1971 rs = mat(6, n);
1972 dts = mat(2, n);
1973 var = mat(1, n);
1974 y = mat(nf * 2, n);
1975 e = mat(3, n);
1976 azel = zeros(2, n);
1977 freq = zeros(nf, n);
1978
1979 for (i = 0; i < MAXSAT; i++) {
1980 rtk->ssat[i].sys = satsys(i + 1, NULL);
1981 for (j = 0; j < NFREQ; j++) {
1982 rtk->ssat[i].vsat[j] = 0;
1983 }
1984 for (j = 1; j < NFREQ; j++) {
1985 rtk->ssat[i].snr[j] = 0;
1986 }
1987 }
1988 /* satellite positions/clocks */
1989 satposs(time, obs, n, nav, opt->sateph, rs, dts, var, svh);
1990
1991 /* UD (undifferenced) residuals for base station */
1992 if (!zdres(1, obs + nu, nr, rs + nu * 6, dts + nu * 2, var + nu, svh + nu, nav, rtk->rb, opt, 1, y + nu * nf * 2,
1993 e + nu * 3, azel + nu * 2, freq + nu * nf)) {
1994 errmsg(rtk, "initial base station position error\n");
1995
1996 free(rs);
1997 free(dts);
1998 free(var);
1999 free(y);
2000 free(e);
2001 free(azel);
2002 free(freq);
2003 return 0;
2004 }
2005 /* time-interpolation of residuals (for post-processing) */
2006 if (opt->intpref) {
2007 dt = intpres(time, obs + nu, nr, nav, rtk, y + nu * nf * 2);
2008 }
2009 /* select common satellites between rover and base-station */
2010 if ((ns = selsat(obs, azel, nu, nr, opt, sat, iu, ir)) <= 0) {
2011 errmsg(rtk, "no common satellite\n");
2012
2013 free(rs);
2014 free(dts);
2015 free(var);
2016 free(y);
2017 free(e);
2018 free(azel);
2019 free(freq);
2020 return 0;
2021 }
2022 /* temporal update of states */
2023 udstate(rtk, obs, sat, iu, ir, ns, nav);
2024
2025 trace(NULL, 4, "x(0)=");
2026 tracemat(NULL, 4, rtk->x, 1, NR(opt), 13, 4);
2027
2028 xp = mat(rtk->nx, 1);
2029 Pp = zeros(rtk->nx, rtk->nx);
2030 xa = mat(rtk->nx, 1);
2031 matcpy(xp, rtk->x, rtk->nx, 1);
2032
2033 ny = ns * nf * 2 + 2;
2034 v = mat(ny, 1);
2035 H = zeros(rtk->nx, ny);
2036 R = mat(ny, ny);
2037 bias = mat(rtk->nx, 1);
2038
2039 /* add 2 iterations for baseline-constraint moving-base */
2040 niter = opt->niter + (opt->mode == PMODE_MOVEB && opt->baseline[0] > 0.0 ? 2 : 0);
2041
2042 for (i = 0; i < niter; i++) {
2043 /* UD (undifferenced) residuals for rover */
2044 if (!zdres(0, obs, nu, rs, dts, var, svh, nav, xp, opt, 0, y, e, azel, freq)) {
2045 errmsg(rtk, "rover initial position error\n");
2046 stat = SOLQ_NONE;
2047 break;
2048 }
2049 /* DD (double-differenced) residuals and partial derivatives */
2050 if ((nv = ddres(rtk, nav, obs, dt, xp, Pp, sat, y, e, azel, freq, iu, ir, ns, v, H, R, vflg)) < 4) {
2051 errmsg(rtk, "insufficient double-differenced residuals (%d)\n", nv);
2052 stat = SOLQ_NONE;
2053 break;
2054 }
2055 /* Kalman filter measurement update */
2056 matcpy(Pp, rtk->P, rtk->nx, rtk->nx);
2057 if ((info = filter(xp, Pp, H, v, R, rtk->nx, nv))) {
2058 errmsg(rtk, "filter error (info=%d)\n", info);
2059 stat = SOLQ_NONE;
2060 break;
2061 }
2062 trace(NULL, 4, "x(%d)=", i + 1);
2063 tracemat(NULL, 4, xp, 1, NR(opt), 13, 4);
2064 }
2065 if (stat != SOLQ_NONE && zdres(0, obs, nu, rs, dts, var, svh, nav, xp, opt, 0, y, e, azel, freq)) {
2066 /* post-fit residuals for float solution */
2067 nv = ddres(rtk, nav, obs, dt, xp, Pp, sat, y, e, azel, freq, iu, ir, ns, v, NULL, R, vflg);
2068
2069 /* validation of float solution */
2070 if (valpos(rtk, v, R, vflg, nv, 4.0)) {
2071 /* update state and covariance matrix */
2072 matcpy(rtk->x, xp, rtk->nx, 1);
2073 matcpy(rtk->P, Pp, rtk->nx, rtk->nx);
2074
2075 /* update ambiguity control struct */
2076 rtk->sol.ns = 0;
2077 for (i = 0; i < ns; i++) {
2078 for (f = 0; f < nf; f++) {
2079 if (!rtk->ssat[sat[i] - 1].vsat[f]) {
2080 continue;
2081 }
2082 /* increment lock unconditionally — conditional (demo5 Phase 4D) froze
2083 lock counts on nfix=0 epochs, preventing satellite-set diversification
2084 needed to escape false-fix/holdamb cycles in urban canyons */
2085 rtk->ssat[sat[i] - 1].lock[f]++;
2086 rtk->ssat[sat[i] - 1].outc[f] = 0;
2087 if (f == 0) {
2088 rtk->sol.ns++; /* valid satellite count by L1 */
2089 }
2090 }
2091 }
2092 /* lack of valid satellites */
2093 if (rtk->sol.ns < 4) {
2094 stat = SOLQ_DGPS; /* degrade, not invalidate */
2095 }
2096 } else {
2097 stat = SOLQ_NONE;
2098 }
2099 }
2100 /* resolve integer ambiguity by LAMBDA (partial AR + multi-attempt) */
2101 if (stat == SOLQ_FLOAT && manage_amb_LAMBDA(rtk, bias, xa, sat, nf, ns) > 1) {
2102 if (zdres(0, obs, nu, rs, dts, var, svh, nav, xa, opt, 0, y, e, azel, freq)) {
2103 /* post-fit residuals for fixed solution */
2104 nv = ddres(rtk, nav, obs, dt, xa, NULL, sat, y, e, azel, freq, iu, ir, ns, v, NULL, R, vflg);
2105
2106 /* validation of fixed solution */
2107 if (valpos(rtk, v, R, vflg, nv, 4.0)) {
2108 /* hold integer ambiguity if meet minfix count */
2109 if (++rtk->nfix >= rtk->opt.minfix && rtk->opt.modear == ARMODE_FIXHOLD) {
2110 holdamb(rtk, xa);
2111 }
2112 stat = SOLQ_FIX;
2113 }
2114 }
2115 }
2116 /* save solution status */
2117 if (stat == SOLQ_FIX) {
2118 for (i = 0; i < 3; i++) {
2119 rtk->sol.rr[i] = rtk->xa[i];
2120 rtk->sol.qr[i] = (float)rtk->Pa[i + i * rtk->na];
2121 }
2122 rtk->sol.qr[3] = (float)rtk->Pa[1];
2123 rtk->sol.qr[4] = (float)rtk->Pa[1 + 2 * rtk->na];
2124 rtk->sol.qr[5] = (float)rtk->Pa[2];
2125
2126 if (rtk->opt.dynamics) { /* velocity and covariance */
2127 for (i = 3; i < 6; i++) {
2128 rtk->sol.rr[i] = rtk->xa[i];
2129 rtk->sol.qv[i - 3] = (float)rtk->Pa[i + i * rtk->na];
2130 }
2131 rtk->sol.qv[3] = (float)rtk->Pa[4 + 3 * rtk->na];
2132 rtk->sol.qv[4] = (float)rtk->Pa[5 + 4 * rtk->na];
2133 rtk->sol.qv[5] = (float)rtk->Pa[5 + 3 * rtk->na];
2134 }
2135 } else {
2136 for (i = 0; i < 3; i++) {
2137 rtk->sol.rr[i] = rtk->x[i];
2138 rtk->sol.qr[i] = (float)rtk->P[i + i * rtk->nx];
2139 }
2140 rtk->sol.qr[3] = (float)rtk->P[1];
2141 rtk->sol.qr[4] = (float)rtk->P[1 + 2 * rtk->nx];
2142 rtk->sol.qr[5] = (float)rtk->P[2];
2143
2144 if (rtk->opt.dynamics) { /* velocity and covariance */
2145 for (i = 3; i < 6; i++) {
2146 rtk->sol.rr[i] = rtk->x[i];
2147 rtk->sol.qv[i - 3] = (float)rtk->P[i + i * rtk->nx];
2148 }
2149 rtk->sol.qv[3] = (float)rtk->P[4 + 3 * rtk->nx];
2150 rtk->sol.qv[4] = (float)rtk->P[5 + 4 * rtk->nx];
2151 rtk->sol.qv[5] = (float)rtk->P[5 + 3 * rtk->nx];
2152 }
2153 rtk->nfix = 0;
2154 }
2155 for (i = 0; i < n; i++) {
2156 for (j = 0; j < nf; j++) {
2157 if (obs[i].L[j] == 0.0) {
2158 continue;
2159 }
2160 rtk->ssat[obs[i].sat - 1].pt[obs[i].rcv - 1][j] = obs[i].time;
2161 rtk->ssat[obs[i].sat - 1].ph[obs[i].rcv - 1][j] = obs[i].L[j];
2162 }
2163 }
2164 for (i = 0; i < ns; i++) {
2165 for (j = 0; j < nf; j++) {
2166 /* output snr of rover receiver */
2167 rtk->ssat[sat[i] - 1].snr[j] = obs[iu[i]].SNR[j];
2168 }
2169 }
2170 for (i = 0; i < MAXSAT; i++) {
2171 for (j = 0; j < nf; j++) {
2172 /* retain fix[j]=2 across epochs to preserve AR candidates (demo5) */
2173 /* if (rtk->ssat[i].fix[j]==2&&stat!=SOLQ_FIX) rtk->ssat[i].fix[j]=1; */
2174 if (rtk->ssat[i].slip[j] & 1) {
2175 rtk->ssat[i].slipc[j]++;
2176 }
2177 }
2178 }
2179 free(rs);
2180 free(dts);
2181 free(var);
2182 free(y);
2183 free(e);
2184 free(azel);
2185 free(freq);
2186 free(xp);
2187 free(Pp);
2188 free(xa);
2189 free(v);
2190 free(H);
2191 free(R);
2192 free(bias);
2193
2194 if (stat != SOLQ_NONE) {
2195 rtk->sol.stat = stat;
2196 }
2197
2198 return stat != SOLQ_NONE;
2199}
2200/* initialize RTK control ------------------------------------------------------
2201 * initialize RTK control struct
2202 * args : rtk_t *rtk IO TKk control/result struct
2203 * prcopt_t *opt I positioning options (see rtklib.h)
2204 * return : none
2205 *-----------------------------------------------------------------------------*/
2206extern void rtkinit(rtk_t* rtk, const prcopt_t* opt) {
2207 sol_t sol0 = {{0}};
2208 ambc_t ambc0 = {{{0}}};
2209 ssat_t ssat0 = {0};
2210 int i;
2211
2212 trace(NULL, 3, "rtkinit :\n");
2213
2214 rtk->sol = sol0;
2215 for (i = 0; i < 6; i++) {
2216 rtk->rb[i] = 0.0;
2217 }
2218 if (opt->mode == PMODE_PPP_RTK) {
2219 rtk->nx = ppp_rtk_nx(opt);
2220 rtk->na = ppp_rtk_na(opt);
2221 } else if (opt->mode <= PMODE_FIXED || opt->mode == PMODE_VRS_RTK) {
2222 rtk->nx = NX(opt);
2223 rtk->na = NR(opt);
2224 } else {
2225 rtk->nx = pppnx(opt);
2226 rtk->na = pppnx(opt);
2227 }
2228 rtk->tt = 0.0;
2229 rtk->x = zeros(rtk->nx, 1);
2230 rtk->P = zeros(rtk->nx, rtk->nx);
2231 rtk->xa = zeros(rtk->na, 1);
2232 rtk->Pa = zeros(rtk->na, rtk->na);
2233 rtk->Q = zeros(rtk->nx, rtk->nx);
2234 rtk->nfix = rtk->neb = 0;
2235 for (i = 0; i < MAXSAT; i++) {
2236 rtk->ambc[i] = ambc0;
2237 rtk->ssat[i] = ssat0;
2238 }
2239 for (i = 0; i < MAXERRMSG; i++) {
2240 rtk->errbuf[i] = 0;
2241 }
2242 for (i = 0; i < 6; i++) {
2243 rtk->prev_qr[i] = 0.0;
2244 }
2245 rtk->excsat = 0;
2246 rtk->nb_ar = 0;
2247 rtk->holdamb = 0;
2248 rtk->opt = *opt;
2249 rtk->sol.thres = (float)opt->thresar[0];
2250
2251 /* BeiDou generation of PRN 1-18 and BDS-3 IGSO (#401) */
2252 if (!set_bdsopt(opt, 0)) {
2253 errmsg(rtk, "invalid bds2_sats=%s or bds3_igso_sats=%s\n", opt->bds2prn, opt->bds3igso);
2254 }
2255}
2256/* free rtk control ------------------------------------------------------------
2257 * free memory for rtk control struct
2258 * args : rtk_t *rtk IO rtk control/result struct
2259 * return : none
2260 *-----------------------------------------------------------------------------*/
2261extern void rtkfree(rtk_t* rtk) {
2262 trace(NULL, 3, "rtkfree :\n");
2263
2264 rtk->nx = rtk->na = 0;
2265 free(rtk->x);
2266 rtk->x = NULL;
2267 free(rtk->P);
2268 rtk->P = NULL;
2269 free(rtk->Q);
2270 rtk->Q = NULL;
2271 free(rtk->xa);
2272 rtk->xa = NULL;
2273 free(rtk->Pa);
2274 rtk->Pa = NULL;
2275}
2276
2277/* update satellite code bias ------------------------------------------------*/
2278static void udsatcb(gtime_t gt, nav_t* nav, osb_t* biaosb, const prcopt_t* popt) {
2279 int i, j, sys, udcnt = 0, ssrcode = CODE_NONE;
2280 double vp = MAXAGESSRL6;
2281 char* p;
2282
2283 /* check the vendor and apply specific settings */
2284 for (i = 0; i < MAXSAT; i++) {
2285 if (nav->ssr_ch[0][i].vendor == SSR_VENDOR_RTCM) {
2286 if ((p = strstr(popt->rtcmopt, "-RTCM_CB_VALID=")) && sscanf(p, "-RTCM_CB_VALID=%lf", &vp) == 1) {
2287 }
2288 break;
2289 } else if (nav->ssr_ch[0][i].vendor == SSR_VENDOR_L6) {
2290 break;
2291 }
2292 }
2293
2294 /* auto or ssr */
2295 if (popt->pppsatcb == 0 || popt->pppsatcb == 1) {
2296 for (i = 0; i < MAXSAT; i++) {
2297 sys = satsys(i + 1, NULL);
2298 if (timediff(gt, nav->ssr_ch[0][i].t0[4]) > vp) {
2299 continue;
2300 }
2301 for (j = 0; j < MAXCODE; j++) {
2302 ssrcode = mcssr_sel_biascode(sys, j + 1);
2303 if (ssrcode == CODE_NONE) {
2304 continue;
2305 }
2306 nav->osb.vscb[i][j] = 1;
2307 nav->osb.scb[i][j] = nav->ssr_ch[0][i].cbias[ssrcode - 1];
2308 udcnt++;
2309 }
2310 if (0 < udcnt) {
2311 nav->osb.gt[0] = nav->ssr_ch[0][i].t0[4];
2312 }
2313 }
2314 }
2315 if (0 < udcnt || popt->pppsatcb == 1) {
2316 trace(NULL, 4, "%s ssr update satellite code bias cnt=%d vp=%.f\n", time_str(gt, 0), udcnt, vp);
2317 return;
2318 }
2319
2320 /* auto or bia */
2321 if ((popt->pppsatcb == 0 || popt->pppsatcb == 2) && timediff(gt, biaosb->gt[0]) <= popt->maxbiasdt) {
2322 for (i = 0; i < MAXSAT; i++) {
2323 for (j = 0; j < MAXCODE; j++) {
2324 if (biaosb->vscb[i][j] != 0) {
2325 udcnt++;
2326 }
2327 }
2328 }
2329 if (0 < udcnt) {
2330 memcpy(nav->osb.vscb, biaosb->vscb, sizeof(nav->osb.vscb));
2331 memcpy(nav->osb.scb, biaosb->scb, sizeof(nav->osb.scb));
2332 for (i = 0; i < MAXSAT; i++) {
2333 for (j = 0; j < MAXCODE; j++) {
2334 nav->osb.scb[i][j] *= -1;
2335 }
2336 }
2337 nav->osb.gt[0] = biaosb->gt[0];
2338 }
2339 }
2340
2341 if (0 < udcnt || popt->pppsatcb == 2) {
2342 trace(NULL, 4, "%s bia update satellite code bias cnt=%d dt=%.f\n", time_str(gt, 0), udcnt, popt->maxbiasdt);
2343 return;
2344 }
2345
2346 /* auto or dcb */
2347 if (popt->pppsatcb == 0 || popt->pppsatcb == 3) {
2348 for (i = 0; i < MAXSAT; i++) {
2349 /* P1->P2 is not used, so it is excluded from processing */
2350
2351 /* C1->P1 */
2352 if (nav->cbias[i][1] != 0.0) {
2353 nav->osb.vscb[i][CODE_L1C - 1] = 1;
2354 nav->osb.scb[i][CODE_L1C - 1] = nav->cbias[i][1];
2355 udcnt++;
2356 }
2357 /* C2->P2 */
2358 if (nav->cbias[i][2] != 0.0) {
2359 nav->osb.vscb[i][CODE_L2C - 1] = 1;
2360 nav->osb.scb[i][CODE_L2C - 1] = nav->cbias[i][2];
2361 udcnt++;
2362 }
2363 }
2364 }
2365 if (0 < udcnt) {
2366 trace(NULL, 4, "%s dcb update satellite code bias cnt=%d\n", time_str(gt, 0), udcnt);
2367 return;
2368 }
2369 trace(NULL, 3, "%s no update satellite code bias \n", time_str(gt, 0));
2370 return;
2371}
2372
2373/* update satellite phase bias ------------------------------------------------*/
2374static void udsatpb(gtime_t gt, nav_t* nav, osb_t* biaosb, osb_t* fcbosb, const prcopt_t* popt) {
2375 int i, j, sys, udcnt = 0, ssrcode = CODE_NONE;
2376 double vp = MAXAGESSRL6;
2377 char* p;
2378
2379 /* check the vendor and apply specific settings */
2380 for (i = 0; i < MAXSAT; i++) {
2381 if (nav->ssr_ch[0][i].vendor == SSR_VENDOR_RTCM) {
2382 if ((p = strstr(popt->rtcmopt, "-RTCM_PB_VALID=")) && sscanf(p, "-RTCM_PB_VALID=%lf", &vp) == 1) {
2383 }
2384 break;
2385 }
2386 }
2387
2388 /* auto or ssr */
2389 if (popt->pppsatpb == 0 || popt->pppsatpb == 1) {
2390 for (i = 0; i < MAXSAT; i++) {
2391 sys = satsys(i + 1, NULL);
2392 if (timediff(gt, nav->ssr_ch[0][i].t0[5]) > vp) {
2393 continue;
2394 }
2395 for (j = 0; j < MAXCODE; j++) {
2396 ssrcode = mcssr_sel_biascode(sys, j + 1);
2397 if (ssrcode == CODE_NONE) {
2398 continue;
2399 }
2400 nav->osb.vspb[i][j] = 1;
2401 nav->osb.spb[i][j] = nav->ssr_ch[0][i].pbias[ssrcode - 1];
2402 if (nav->ssr_ch[0][i].vendor == SSR_VENDOR_RTCM) {
2403 nav->osb.spb[i][j] *= -1;
2404 }
2405 udcnt++;
2406 }
2407 if (0 < udcnt) {
2408 nav->osb.gt[0] = nav->ssr_ch[0][i].t0[5];
2409 }
2410 }
2411 }
2412 if (0 < udcnt || popt->pppsatpb == 1) {
2413 trace(NULL, 4, "%s ssr update satellite phase bias cnt=%d vp=%.f\n", time_str(gt, 0), udcnt, vp);
2414 return;
2415 }
2416
2417 /* auto or fcb */
2418 if ((popt->pppsatpb == 0 || popt->pppsatpb == 3) && timediff(gt, fcbosb->gt[0]) <= popt->maxbiasdt) {
2419 for (i = 0; i < MAXSAT; i++) {
2420 for (j = 0; j < MAXCODE; j++) {
2421 if (fcbosb->vspb[i][j] != 0) {
2422 udcnt++;
2423 }
2424 }
2425 }
2426 if (0 < udcnt) {
2427 memcpy(nav->osb.vspb, fcbosb->vspb, sizeof(nav->osb.vspb));
2428 memcpy(nav->osb.spb, fcbosb->spb, sizeof(nav->osb.spb));
2429 for (i = 0; i < MAXSAT; i++) {
2430 for (j = 0; j < MAXCODE; j++) {
2431 nav->osb.spb[i][j] *= -1;
2432 }
2433 }
2434 nav->osb.gt[0] = fcbosb->gt[0];
2435 }
2436 }
2437 if (0 < udcnt) {
2438 trace(NULL, 4, "%s fcb update satellite phase bias cnt=%d\n", time_str(gt, 0), udcnt);
2439 return;
2440 }
2441
2442 /* auto or bia (#142: satellite phase bias from the Bias-SINEX OSB product,
2443 * e.g. COD MGEX IAR OSB.BIA). Mirrors the bia branch of udsatcb(): only
2444 * reached when no SSR/FCB phase bias was found, so SSR-based engines
2445 * (MADOCA/CLAS/VRS) and FCB users are unaffected. Stored negated (*=-1) to
2446 * match the consumer-adds convention used for the code bias. */
2447 if ((popt->pppsatpb == 0 || popt->pppsatpb == 2) && timediff(gt, biaosb->gt[0]) <= popt->maxbiasdt) {
2448 for (i = 0; i < MAXSAT; i++) {
2449 for (j = 0; j < MAXCODE; j++) {
2450 if (biaosb->vspb[i][j] != 0) {
2451 udcnt++;
2452 }
2453 }
2454 }
2455 if (0 < udcnt) {
2456 memcpy(nav->osb.vspb, biaosb->vspb, sizeof(nav->osb.vspb));
2457 memcpy(nav->osb.spb, biaosb->spb, sizeof(nav->osb.spb));
2458 for (i = 0; i < MAXSAT; i++) {
2459 for (j = 0; j < MAXCODE; j++) {
2460 nav->osb.spb[i][j] *= -1;
2461 }
2462 }
2463 nav->osb.gt[0] = biaosb->gt[0];
2464 }
2465 }
2466 if (0 < udcnt || popt->pppsatpb == 2) {
2467 if (udcnt == 0) {
2468 /* forced bia but no valid OSB phase bias (e.g. expired past
2469 * maxbiasdt): drop any stale phase biases so corr_meas() does not
2470 * keep applying them. */
2471 memset(nav->osb.vspb, 0x00, sizeof(nav->osb.vspb));
2472 memset(nav->osb.spb, 0x00, sizeof(nav->osb.spb));
2473 }
2474 trace(NULL, 4, "%s bia update satellite phase bias cnt=%d dt=%.f\n", time_str(gt, 0), udcnt, popt->maxbiasdt);
2475 return;
2476 }
2477
2478 trace(NULL, 3, "%s no update satellite phase bias \n", time_str(gt, 0));
2479 return;
2480}
2481
2482/* update station code bias ------------------------------------------------*/
2483static void udstacb(gtime_t gt, nav_t* nav, osb_t* biaosb, const prcopt_t* popt) {
2484 int i, j, udcnt = 0;
2485
2486 /* auto or bia */
2487 if (timediff(gt, biaosb->gt[0]) <= popt->maxbiasdt) {
2488 for (i = 0; i < MAXBSNXSYS; i++) {
2489 for (j = 0; j < MAXCODE; j++) {
2490 if (biaosb->vrsyscb[i][j] != 0) {
2491 udcnt++;
2492 }
2493 }
2494 }
2495 if (0 < udcnt) {
2496 memcpy(nav->osb.vrsyscb, biaosb->vrsyscb, sizeof(nav->osb.vrsyscb));
2497 memcpy(nav->osb.rsyscb, biaosb->rsyscb, sizeof(nav->osb.rsyscb));
2498 for (i = 0; i < MAXBSNXSYS; i++) {
2499 for (j = 0; j < MAXCODE; j++) {
2500 nav->osb.rsyscb[i][j] *= -1;
2501 }
2502 }
2503 nav->osb.gt[1] = biaosb->gt[1];
2504 }
2505 for (i = 0; i < MAXSAT; i++) {
2506 for (j = 0; j < MAXCODE; j++) {
2507 if (biaosb->vrsatcb[i][j] != 0) {
2508 udcnt++;
2509 }
2510 }
2511 }
2512 if (0 < udcnt) {
2513 memcpy(nav->osb.vrsatcb, biaosb->vrsatcb, sizeof(nav->osb.vrsatcb));
2514 memcpy(nav->osb.rsatcb, biaosb->rsatcb, sizeof(nav->osb.rsatcb));
2515 for (i = 0; i < MAXSAT; i++) {
2516 for (j = 0; j < MAXCODE; j++) {
2517 nav->osb.rsatcb[i][j] *= -1;
2518 }
2519 }
2520 nav->osb.gt[1] = biaosb->gt[1];
2521 }
2522 }
2523
2524 if (0 < udcnt) {
2525 trace(NULL, 4, "%s bia update station code bias cnt=%d dt=%.f\n", time_str(gt, 0), udcnt, popt->maxbiasdt);
2526 return;
2527 }
2528
2529 trace(NULL, 3, "%s no update station code bias \n", time_str(gt, 0));
2530 return;
2531}
2532
2533/* update biass --------------------------------------------------------------*/
2534static void udbiass(gtime_t gt, prcopt_t* popt, nav_t* nav) {
2535 char path[1024];
2536 osb_t *biaosb = NULL, *fcbosb = NULL;
2537 int ret;
2538
2539 /* heap-allocate osb_t (~700KB each) to avoid thread stack overflow */
2540 biaosb = (osb_t*)calloc(1, sizeof(osb_t));
2541 fcbosb = (osb_t*)calloc(1, sizeof(osb_t));
2542 if (!biaosb || !fcbosb) {
2543 free(biaosb);
2544 free(fcbosb);
2545 return;
2546 }
2547
2548 /* update bia file data */
2549 if (strlen(nav->biapath) > 0) {
2550 reppath(nav->biapath, path, gt, "", "");
2551 if (strcmp(nav->pr_biapath, path)) {
2552 if (readbsnx(path) != 0) {
2553 trace(NULL, 3, "update bias sinex file data %s -> %s\n", nav->pr_biapath, path);
2554 strncpy(nav->pr_biapath, path, sizeof(nav->pr_biapath));
2555 }
2556 }
2557 /* update satellite osb */
2558 if ((ret = udosb_sat(biaosb, gt, 0)) == 0) {
2559 if (popt->maxbiasdt > 0) {
2560 if ((ret = udosb_sat(biaosb, gt, 1)) != 0) {
2561 if (timediff(gt, biaosb->gt[0]) > popt->maxbiasdt) {
2562 ret = 0;
2563 }
2564 }
2565 }
2566 }
2567 if (ret == 0) {
2568 memset(biaosb->scb, 0x00, sizeof(biaosb->scb));
2569 memset(biaosb->spb, 0x00, sizeof(biaosb->spb));
2570 memset(biaosb->vscb, 0x00, sizeof(biaosb->vscb));
2571 memset(biaosb->vspb, 0x00, sizeof(biaosb->vspb));
2572 }
2573 trace(NULL, 4, "update satellite osb: %s num=%d\n", time_str(gt, 3), ret);
2574
2575 /* update station osb */
2576 if ((ret = udosb_station(biaosb, gt, 0, popt->staname)) == 0) {
2577 if (popt->maxbiasdt > 0) {
2578 if ((ret = udosb_station(biaosb, gt, 1, popt->staname)) != 0) {
2579 if (timediff(gt, biaosb->gt[1]) > popt->maxbiasdt) {
2580 ret = 0;
2581 }
2582 }
2583 }
2584 }
2585 if (ret == 0) {
2586 memset(biaosb->rsyscb, 0x00, sizeof(biaosb->rsyscb));
2587 memset(biaosb->rsatcb, 0x00, sizeof(biaosb->rsatcb));
2588 memset(biaosb->vrsyscb, 0x00, sizeof(biaosb->vrsyscb));
2589 memset(biaosb->vrsatcb, 0x00, sizeof(biaosb->vrsatcb));
2590 }
2591 trace(NULL, 4, "update station osb: %s num=%d\n", time_str(gt, 3), ret);
2592 }
2593 /* update fcb file data */
2594 if (strlen(nav->fcbpath) > 0) {
2595 reppath(nav->fcbpath, path, gt, "", "");
2596 if (strcmp(nav->pr_fcbpath, path)) {
2597 if (readfcb(path) != 0) {
2598 trace(NULL, 3, "update fcb product file data %s -> %s\n", nav->pr_fcbpath, path);
2599 strncpy(nav->pr_fcbpath, path, sizeof(nav->pr_fcbpath));
2600 }
2601 }
2602 /* update satellite fcb */
2603 if ((ret = udfcb_sat(fcbosb, gt, 0)) == 0) {
2604 if (popt->maxbiasdt > 0) {
2605 if ((ret = udfcb_sat(fcbosb, gt, 1)) != 0) {
2606 if (timediff(gt, fcbosb->gt[0]) > popt->maxbiasdt) {
2607 ret = 0;
2608 }
2609 }
2610 }
2611 }
2612 if (ret == 0) {
2613 memset(fcbosb->spb, 0x00, sizeof(fcbosb->scb));
2614 memset(fcbosb->vspb, 0x00, sizeof(fcbosb->vscb));
2615 }
2616 trace(NULL, 4, "update satellite fcb: %s num=%d\n", time_str(gt, 3), ret);
2617 }
2618
2619 /* update satellite code bias and station codebias correction */
2620 udsatcb(gt, nav, biaosb, popt);
2621 udsatpb(gt, nav, biaosb, fcbosb, popt);
2622 udstacb(gt, nav, biaosb, popt);
2623
2624 free(biaosb);
2625 free(fcbosb);
2626}
2627
2628/* precise positioning ---------------------------------------------------------
2629 * input observation data and navigation message, compute rover position by
2630 * precise positioning
2631 * args : rtk_t *rtk IO RTK control/result struct
2632 * rtk->sol IO solution
2633 * .time O solution time
2634 * .rr[] IO rover position/velocity
2635 * (I:fixed mode,O:single mode)
2636 * .dtr[0] O receiver clock bias (s)
2637 * .dtr[1-5] O receiver GLO/GAL/BDS/IRN/QZS-GPS time offset (s)
2638 * .Qr[] O rover position covarinace
2639 * .stat O solution status (SOLQ_???)
2640 * .ns O number of valid satellites
2641 * .age O age of differential (s)
2642 * .ratio O ratio factor for ambiguity validation
2643 * rtk->rb[] IO base station position/velocity
2644 * (I:relative mode,O:moving-base mode)
2645 * rtk->nx I number of all states
2646 * rtk->na I number of integer states
2647 * rtk->ns O number of valid satellites in use
2648 * rtk->tt O time difference between current and previous (s)
2649 * rtk->x[] IO float states pre-filter and post-filter
2650 * rtk->P[] IO float covariance pre-filter and post-filter
2651 * rtk->xa[] O fixed states after AR
2652 * rtk->Pa[] O fixed covariance after AR
2653 * rtk->ssat[s] IO satellite {s+1} status
2654 * .sys O system (SYS_???)
2655 * .az [r] O azimuth angle (rad) (r=0:rover,1:base)
2656 * .el [r] O elevation angle (rad) (r=0:rover,1:base)
2657 * .vs [r] O data valid single (r=0:rover,1:base)
2658 * .resp [f] O freq(f+1) pseudorange residual (m)
2659 * .resc [f] O freq(f+1) carrier-phase residual (m)
2660 * .vsat [f] O freq(f+1) data vaild (0:invalid,1:valid)
2661 * .fix [f] O freq(f+1) ambiguity flag
2662 * (0:nodata,1:float,2:fix,3:hold)
2663 * .slip [f] O freq(f+1) cycle slip flag
2664 * (bit8-7:rcv1 LLI, bit6-5:rcv2 LLI,
2665 * bit2:parity unknown, bit1:slip)
2666 * .lock [f] IO freq(f+1) carrier lock count
2667 * .outc [f] IO freq(f+1) carrier outage count
2668 * .slipc[f] IO freq(f+1) cycle slip count
2669 * .rejc [f] IO freq(f+1) data reject count
2670 * .gf IO geometry-free phase (L1-L2) (m)
2671 * .gf2 IO geometry-free phase (L1-L5) (m)
2672 * rtk->nfix IO number of continuous fixes of ambiguity
2673 * rtk->neb IO bytes of error message buffer
2674 * rtk->errbuf IO error message buffer
2675 * rtk->tstr O time string for debug
2676 * rtk->opt I processing options
2677 * obsd_t *obs I observation data for an epoch
2678 * obs[i].rcv=1:rover,2:reference
2679 * sorted by receiver and satellte
2680 * int n I number of observation data
2681 * nav_t *nav I navigation messages
2682 * return : status (0:no solution,1:valid solution)
2683 * notes : before calling function, base station position rtk->sol.rb[] should
2684 * be properly set for relative mode except for moving-baseline
2685 *-----------------------------------------------------------------------------*/
2686extern int rtkpos(mrtk_ctx_t* ctx, rtk_t* rtk, const obsd_t* obs, int n, nav_t* nav) {
2687 prcopt_t* opt = &rtk->opt;
2688 sol_t solb = {{0}};
2689 gtime_t time;
2690 int i, nu, nr[2] = {0};
2691 char msg[128] = "";
2692
2693 trace(ctx, 3, "rtkpos : time=%s n=%d\n", time_str(obs[0].time, 3), n);
2694
2695 /* set base station position */
2696 if (opt->refpos <= POSOPT_RINEX && opt->mode != PMODE_SINGLE && opt->mode != PMODE_MOVEB) {
2697 for (i = 0; i < 6; i++) {
2698 rtk->rb[i] = i < 3 ? opt->rb[i] : 0.0;
2699 }
2700 }
2701 /* count rover/base station observations */
2702 for (nu = 0; nu < n && obs[nu].rcv == 1; nu++) {
2703 /* count via loop header */
2704 }
2705 for (nr[0] = 0; nu + nr[0] < n && obs[nu + nr[0]].rcv == 2; nr[0]++) {
2706 /* count via loop header */
2707 }
2708 if (opt->l6mrg) {
2709 for (nr[1] = 0; nu + nr[0] + nr[1] < n && obs[nu + nr[0] + nr[1]].rcv == 3; nr[1]++) {
2710 /* count via loop header */
2711 }
2712 }
2713
2714 time = rtk->sol.time; /* previous epoch */
2715
2716 /* update code/phase bias */
2717 udbiass(obs[0].time, opt, nav);
2718
2719 /* rover position by single point positioning */
2720 {
2721 prcopt_t sppopt = rtk->opt;
2722 /* PPP-RTK/VRS: force broadcast ephemeris for SPP initial position */
2723 if (sppopt.mode == PMODE_PPP_RTK || sppopt.mode == PMODE_VRS_RTK) {
2724 sppopt.sateph = EPHOPT_BRDC;
2725 /* Upstream CLAS applies SPP pseudorange large-residual exclusion with
2726 * defaults of 20 m and 7 satellites. Keep standalone SPP bit-compatible
2727 * by applying the auto-default only to this private CLAS seed copy.
2728 * Set pos2-rejethres < 0 to force it off for CLAS diagnostics. */
2729 if (sppopt.rejethres == 0.0) {
2730 sppopt.rejethres = 20.0;
2731 }
2732 if (sppopt.rejeminsat <= 0) {
2733 sppopt.rejeminsat = 7;
2734 }
2735 /* Enhanced seed (default SEEDENH_BASE = cn0+tdcp; set OFF to disable):
2736 * apply the proven v0.6.10 SPP error model to this PRIVATE copy only.
2737 * err[5]/err[6] are the C/N0 snr_max/snr_error coefficients here; the
2738 * CLAS measurement model uses its dedicated atmosphere-error fields.
2739 * When the profile is OFF the seed is bit-identical to prior behaviour. */
2740 if (rtk->opt.enhanced_spp_seed >= SEEDENH_BASE) {
2741 sppopt.err[5] = 50.0; /* C/N0 reference snr_max (dB-Hz) */
2742 sppopt.err[6] = 0.5; /* C/N0 weighting coefficient (m) */
2743 sppopt.tdcp = 1; /* TDCP velocity + jump-rejection QC */
2744 sppopt.tdcpjump = 5.0; /* reject seed epoch on code-vs-TDCP jump > 5 m */
2745 sppopt.thresdop = 1.0; /* TDCP cycle-slip Doppler threshold (cyc/s) */
2746 /* IGG-III robust is an explicit opt-in (SEEDENH_ROBUST), not the
2747 * recommended default: it improves average seed accuracy but in deep
2748 * urban canyons perturbs the seed enough to trip the filter into
2749 * fix-and-hold mis-fixes (tokyo_run3 misfix 5->71). Auto-gating it
2750 * from seed-local diagnostics proved infeasible (the harm is a filter
2751 * coupling effect, invisible to the robust residual statistics); see
2752 * docs/design/spp-accuracy.md. So it is left to the operator, who
2753 * knows whether the environment is open-sky (use it) or urban (don't). */
2754 if (rtk->opt.enhanced_spp_seed == SEEDENH_ROBUST) {
2755 sppopt.robust = 1;
2756 sppopt.robustk[0] = 1.5;
2757 sppopt.robustk[1] = 4.0;
2758 }
2759 }
2760 }
2761 if (!pntpos(ctx, obs, nu, nav, &sppopt, &rtk->sol, NULL, rtk->ssat, msg)) {
2762 errmsg(rtk, "point pos error (%s)\n", msg);
2763
2764 if (!rtk->opt.dynamics) {
2765 trace(ctx, 4, "obs=\n");
2766 traceobs(ctx, 4, obs, n);
2767 outsolstat(rtk);
2768 return 0;
2769 }
2770 }
2771 }
2772 if (time.time != 0) {
2773 rtk->tt = timediff(rtk->sol.time, time);
2774 }
2775
2776 /* single point positioning */
2777 if (opt->mode == PMODE_SINGLE) {
2778 trace(ctx, 4, "obs=\n");
2779 traceobs(ctx, 4, obs, n);
2780 outsolstat(rtk);
2781 return 1;
2782 }
2783 /* suppress output of single solution */
2784 if (!opt->outsingle) {
2785 rtk->sol.stat = SOLQ_NONE;
2786 }
2787 /* SSR2OSR conversion mode (upstream compat: dispatched via rtkpos) */
2788 if (opt->mode >= PMODE_SSR2OSR && opt->mode <= PMODE_SSR2OSR_FIXED) {
2789 clas_ctx_t* clas = (clas_ctx_t*)nav->clas_ctx;
2790 static clas_osrd_t osr_buf[MAXOBS];
2791 if (clas) {
2792 int nobs = clas_ssr2osr(rtk, obs, nu, nav, osr_buf, 0, clas);
2793 if (nobs > 0) {
2794 /* set solution for outsol() output */
2795 for (i = 0; i < 3; i++) {
2796 rtk->sol.rr[i] = rtk->x[i];
2797 }
2798 rtk->sol.stat = SOLQ_PPP;
2799 rtk->sol.time = obs[0].time;
2800 rtk->sol.ns = (uint8_t)nobs;
2801 }
2802 }
2803 outsolstat(rtk);
2804 return 1;
2805 }
2806 /* PPP-RTK positioning (CLAS) */
2807 if (opt->mode == PMODE_PPP_RTK) {
2808 trace(ctx, 4, "obs=\n");
2809 traceobs(ctx, 4, obs, n);
2810 ppp_rtk_pos(rtk, obs, nu, nav);
2811 outsolstat(rtk);
2812 return 1;
2813 }
2814 /* VRS-RTK positioning (CLAS) */
2815 if (opt->mode == PMODE_VRS_RTK) {
2816 static int floatcnt = 0;
2817 trace(ctx, 4, "obs=\n");
2818 traceobs(ctx, 4, obs, n);
2819 relposvrs(rtk, obs, nu, nr, nav);
2820 if (opt->floatcnt > 0 && rtk->sol.stat == SOLQ_FLOAT) {
2821 if (++floatcnt > opt->floatcnt) {
2822 trace(ctx, 1, "float status continued than %d epoch: filter reset.\n", opt->floatcnt);
2823 nav->filreset = 1;
2824 floatcnt = 0;
2825 }
2826 } else {
2827 floatcnt = 0;
2828 }
2829 outsolstat(rtk);
2830 return 1;
2831 }
2832 /* precise point positioning */
2833 if (opt->mode >= PMODE_PPP_KINEMA) {
2834 miono_get_corr(rtk->sol.rr, nav);
2835 trace(ctx, 4, "obs=\n");
2836 traceobs(ctx, 4, obs, n);
2837 pppos(ctx, rtk, obs, nu, nav);
2838 outsolstat(rtk);
2839 return 1;
2840 }
2841 trace(ctx, 4, "obs=\n");
2842 traceobs(ctx, 4, obs, n);
2843 /* check number of data of base station and age of differential */
2844 if (nr[0] == 0) {
2845 errmsg(rtk, "no base station observation data for rtk\n");
2846 outsolstat(rtk);
2847 return 1;
2848 }
2849 if (opt->mode == PMODE_MOVEB) { /* moving baseline */
2850
2851 /* estimate position/velocity of base station */
2852 if (!pntpos(ctx, obs + nu, nr[0], nav, &rtk->opt, &solb, NULL, NULL, msg)) {
2853 errmsg(rtk, "base station position error (%s)\n", msg);
2854 return 0;
2855 }
2856 rtk->sol.age = (float)timediff(rtk->sol.time, solb.time);
2857
2858 if (fabs(rtk->sol.age) > TTOL_MOVEB) {
2859 errmsg(rtk, "time sync error for moving-base (age=%.1f)\n", rtk->sol.age);
2860 return 0;
2861 }
2862 for (i = 0; i < 6; i++) {
2863 rtk->rb[i] = solb.rr[i];
2864 }
2865
2866 /* time-synchronized position of base station */
2867 for (i = 0; i < 3; i++) {
2868 rtk->rb[i] += rtk->rb[i + 3] * rtk->sol.age;
2869 }
2870 } else {
2871 rtk->sol.age = (float)timediff(obs[0].time, obs[nu].time);
2872
2873 if (fabs(rtk->sol.age) > opt->maxtdiff) {
2874 errmsg(rtk, "age of differential error (age=%.1f)\n", rtk->sol.age);
2875 outsolstat(rtk);
2876 return 1;
2877 }
2878 }
2879 /* relative positioning */
2880 relpos(rtk, obs, nu, nr[0], nav);
2881 outsolstat(rtk);
2882
2883 return 1;
2884}
antmodel
void antmodel(const pcv_t *pcv, const double *del, const double *azel, int opt, double *dant)
Compute receiver antenna offset by antenna phase center parameters.
Definition mrtk_antenna.c:481
mrtk_antenna.h
MRTKLIB Antenna Module — Antenna phase center model and parameter I/O.
tropmodel
double tropmodel(gtime_t time, const double *pos, const double *azel, double humi)
Compute tropospheric delay by standard atmosphere and Saastamoinen.
Definition mrtk_atmos.c:162
tropmapf
double tropmapf(gtime_t time, const double pos[], const double azel[], double *mapfw)
Definition mrtk_atmos.c:256
ionmapf
double ionmapf(const double *pos, const double *azel)
Compute ionospheric delay mapping function by single layer model.
Definition mrtk_atmos.c:114
mrtk_atmos.h
MRTKLIB Atmosphere Module — Ionospheric and tropospheric models.
satno2id
void satno2id(int sat, char *id)
Definition mrtk_nav.c:817
udosb_station
int udosb_station(osb_t *osb, gtime_t gt, int mode, char *name)
Update station OSB structure with current biases.
Definition mrtk_bias_sinex.c:578
readbsnx
int readbsnx(const char *file)
Read BIAS-SINEX file.
Definition mrtk_bias_sinex.c:134
code2obs
char * code2obs(uint8_t code)
Definition mrtk_obs.c:190
udosb_sat
int udosb_sat(osb_t *osb, gtime_t gt, int mode)
Update satellite OSB structure with current biases.
Definition mrtk_bias_sinex.c:528
mrtk_bias_sinex.h
MRTKLIB Bias SINEX Module — BIAS-SINEX data types and I/O functions.
setbitu
void setbitu(uint8_t *buff, int pos, int len, uint32_t data)
Set unsigned bits in byte data.
Definition mrtk_bits.c:122
getbitu
uint32_t getbitu(const uint8_t *buff, int pos, int len)
Extract unsigned bits from byte data.
Definition mrtk_bits.c:99
mrtk_bits.h
MRTKLIB Bits Module — Binary bit operations and CRC functions.
SQR
#define SQR(x)
Definition mrtk_clas.c:46
mrtk_clas.h
CLAS (QZSS L6D) CSSR decoder and grid correction type definitions.
clas_ssr2osr
int clas_ssr2osr(rtk_t *rtk, obsd_t *obs, int n, nav_t *nav, clas_osrd_t *osr, int mode, clas_ctx_t *clas)
Convert compact SSR corrections to observation-space representations.
Definition mrtk_clas_osr.c:1493
sat2freq
double sat2freq(int sat, uint8_t code, const nav_t *nav)
Definition mrtk_obs.c:244
mrtk_const.h
MRTKLIB Constants — Physical, navigation, and system constants.
EFACT_IRN
#define EFACT_IRN
Definition mrtk_const.h:178
EFACT_CMP
#define EFACT_CMP
Definition mrtk_const.h:175
EFACT_GAL
#define EFACT_GAL
Definition mrtk_const.h:169
EFACT_QZS
#define EFACT_QZS
Definition mrtk_const.h:172
geodist
double geodist(const double *rs, const double *rr, double *e)
Compute geometric distance and receiver-to-satellite unit vector.
Definition mrtk_coords.c:493
D2R
#define D2R
Definition mrtk_coords.c:36
satazel
double satazel(const double *pos, const double *e, double *azel)
Compute satellite azimuth/elevation angle.
Definition mrtk_coords.c:517
CLIGHT
#define CLIGHT
Definition mrtk_coords.c:41
covecef
void covecef(const double *pos, const double *Q, double *P)
Transform local ENU covariance to ECEF coordinate.
Definition mrtk_coords.c:202
MAXSAT
#define MAXSAT
Definition mrtk_coords.c:45
ecef2pos
void ecef2pos(const double *r, double *pos)
Transform ECEF position to geodetic position.
Definition mrtk_coords.c:102
ecef2enu
void ecef2enu(const double *pos, const double *r, double *e)
Transform ECEF vector to local tangential coordinate.
Definition mrtk_coords.c:162
mrtk_coords.h
MRTKLIB Coordinates Module — Coordinate transformations and geometry.
satsys
int satsys(int sat, int *prn)
Definition mrtk_nav.c:117
SYS_GPS
#define SYS_GPS
Definition mrtk_eph.c:65
SYS_SBS
#define SYS_SBS
Definition mrtk_eph.c:66
SYS_QZS
#define SYS_QZS
Definition mrtk_eph.c:69
satposs
void satposs(gtime_t teph, const obsd_t *obs, int n, const nav_t *nav, int ephopt, double *rs, double *dts, double *var, int *svh)
Compute satellite positions, velocities and clocks.
Definition mrtk_eph.c:940
SYS_IRN
#define SYS_IRN
Definition mrtk_eph.c:71
SYS_CMP
#define SYS_CMP
Definition mrtk_eph.c:70
SYS_GAL
#define SYS_GAL
Definition mrtk_eph.c:68
EPHOPT_BRDC
#define EPHOPT_BRDC
Definition mrtk_eph.c:83
SYS_GLO
#define SYS_GLO
Definition mrtk_eph.c:67
nav_t
struct nav_s nav_t
Definition mrtk_eph.h:40
udfcb_sat
int udfcb_sat(osb_t *osb, gtime_t gt, int mode)
Update satellite FCB structure with current biases.
Definition mrtk_fcb.c:169
readfcb
int readfcb(const char *file)
Read FCB (Fractional Cycle Bias) file.
Definition mrtk_fcb.c:136
mrtk_fcb.h
MRTKLIB FCB Module — Fractional Cycle Bias data types and functions.
SNR_UNIT
#define SNR_UNIT
Definition mrtk_foundation.h:54
NFREQGLO
#define NFREQGLO
Definition mrtk_foundation.h:47
MAXOBS
#define MAXOBS
Definition mrtk_foundation.h:151
MAXCODE
#define MAXCODE
Definition mrtk_foundation.h:256
CODE_L1C
#define CODE_L1C
Definition mrtk_foundation.h:186
CODE_L2C
#define CODE_L2C
Definition mrtk_foundation.h:199
DTTOL
#define DTTOL
Definition mrtk_foundation.h:158
CODE_NONE
#define CODE_NONE
Definition mrtk_foundation.h:185
MAXERRMSG
#define MAXERRMSG
Definition mrtk_foundation.h:172
NFREQPCV
#define NFREQPCV
Definition mrtk_foundation.h:48
MAXSOLMSG
#define MAXSOLMSG
Definition mrtk_foundation.h:170
MAXBSNXSYS
#define MAXBSNXSYS
Definition mrtk_foundation.h:307
lambda
int lambda(int n, int m, const double *a, const double *Q, double *F, double *s)
Integer least-square estimation (LAMBDA/MLAMBDA).
Definition mrtk_lambda.c:225
mrtk_lambda.h
MRTKLIB Integer Ambiguity Resolution Module — LAMBDA method.
mcssr_sel_biascode
int mcssr_sel_biascode(const int sys, const int code)
Select MADOCA-PPP code/phase bias code from observation code.
Definition mrtk_madoca.c:988
mrtk_madoca.h
MRTKLIB MADOCA-PPP Module — QZSS L6E CSSR message decode functions.
mat
double * mat(int n, int m)
Allocate memory for a double matrix.
Definition mrtk_mat.c:76
matinv
int matinv(double *A, int n)
Inverse of matrix (A = A^-1).
Definition mrtk_mat.c:421
zeros
double * zeros(int n, int m)
Allocate a new zero matrix.
Definition mrtk_mat.c:108
matcpy
void matcpy(double *A, const double *B, int n, int m)
Copy matrix.
Definition mrtk_mat.c:205
imat
int * imat(int n, int m)
Allocate memory for an integer matrix.
Definition mrtk_mat.c:92
filter
int filter(double *x, double *P, const double *H, const double *v, const double *R, int n, int m)
Kalman filter state update.
Definition mrtk_mat.c:524
norm
double norm(const double *a, int n)
Euclidean norm of vector.
Definition mrtk_mat.c:165
matmul
void matmul(const char *tr, int n, int k, int m, double alpha, const double *A, const double *B, double beta, double *C)
Multiply matrix by matrix (C = alpha*A*B + beta*C).
Definition mrtk_mat.c:310
eye
double * eye(int n)
Generate a new identity matrix.
Definition mrtk_mat.c:129
mrtk_mat.h
MRTKLIB Matrix Module — Matrix/vector operations and estimation.
satexclude
int satexclude(int sat, double var, int svh, const struct prcopt_t *opt)
Test excluded satellite.
Definition mrtk_nav.c:1169
set_bdsopt
int set_bdsopt(const prcopt_t *opt, int warn)
Apply the BeiDou PRN-reassignment options process-wide.
Definition mrtk_nav.c:282
R2D
#define R2D
Definition mrtk_obs.c:51
FREQ1
#define FREQ1
Definition mrtk_obs.c:67
testelmask
int testelmask(const double *azel, const int16_t *elmask)
Test elevation mask.
Definition mrtk_obs.c:160
testsnr
int testsnr(int base, int idx, double el, double snr, const snrmask_t *mask)
Test SNR mask.
Definition mrtk_obs.c:132
PMODE_PPP_RTK
#define PMODE_PPP_RTK
Definition mrtk_opt.h:54
GLO_ARMODE_FIXHOLD
#define GLO_ARMODE_FIXHOLD
Definition mrtk_opt.h:138
POSOPT_RINEX
#define POSOPT_RINEX
Definition mrtk_opt.h:152
PMODE_MOVEB
#define PMODE_MOVEB
Definition mrtk_opt.h:49
SOLQ_FLOAT
#define SOLQ_FLOAT
Definition mrtk_opt.h:72
IONOOPT_EST
#define IONOOPT_EST
Definition mrtk_opt.h:96
SEEDENH_BASE
#define SEEDENH_BASE
Definition mrtk_opt.h:344
ARMODE_FIXHOLD
#define ARMODE_FIXHOLD
Definition mrtk_opt.h:131
SOLQ_DGPS
#define SOLQ_DGPS
Definition mrtk_opt.h:74
PMODE_PPP_KINEMA
#define PMODE_PPP_KINEMA
Definition mrtk_opt.h:51
PMODE_DGPS
#define PMODE_DGPS
Definition mrtk_opt.h:46
TROPOPT_EST
#define TROPOPT_EST
Definition mrtk_opt.h:105
PMODE_SSR2OSR_FIXED
#define PMODE_SSR2OSR_FIXED
Definition mrtk_opt.h:56
PMODE_FIXED
#define PMODE_FIXED
Definition mrtk_opt.h:50
TROPOPT_ESTG
#define TROPOPT_ESTG
Definition mrtk_opt.h:106
SOLQ_NONE
#define SOLQ_NONE
Definition mrtk_opt.h:70
ARMODE_OFF
#define ARMODE_OFF
Definition mrtk_opt.h:128
SOLQ_PPP
#define SOLQ_PPP
Definition mrtk_opt.h:76
PMODE_STATIC
#define PMODE_STATIC
Definition mrtk_opt.h:48
PMODE_VRS_RTK
#define PMODE_VRS_RTK
Definition mrtk_opt.h:57
PMODE_SSR2OSR
#define PMODE_SSR2OSR
Definition mrtk_opt.h:55
SEEDENH_ROBUST
#define SEEDENH_ROBUST
Definition mrtk_opt.h:345
IONOOPT_IFLC
#define IONOOPT_IFLC
Definition mrtk_opt.h:95
ARMODE_INST
#define ARMODE_INST
Definition mrtk_opt.h:130
SOLQ_FIX
#define SOLQ_FIX
Definition mrtk_opt.h:71
PMODE_SINGLE
#define PMODE_SINGLE
Definition mrtk_opt.h:45
mrtk_peph.h
MRTKLIB Precise Ephemeris Module — Earth rotation parameter types and ERP I/O functions.
solb
static sol_t * solb
Definition mrtk_postpos.c:92
fp_stat
static FILE * fp_stat
Definition mrtk_postpos.c:117
NR
#define NR(opt)
Definition mrtk_ppp.c:222
VAR_VEL
#define VAR_VEL
Definition mrtk_ppp.c:195
pppoutstat
int pppoutstat(rtk_t *rtk, char *buff)
Write PPP solution status to buffer.
Definition mrtk_ppp.c:240
GAP_RESION
#define GAP_RESION
Definition mrtk_ppp.c:208
pppos
void pppos(mrtk_ctx_t *ctx, rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
Precise point positioning.
Definition mrtk_ppp.c:1783
pppnx
int pppnx(const prcopt_t *opt)
Number of estimated states for PPP.
Definition mrtk_ppp.c:1681
VAR_GRA
#define VAR_GRA
Definition mrtk_ppp.c:199
IB
#define IB(s, f, opt)
Definition mrtk_ppp.c:230
VAR_POS
#define VAR_POS
Definition mrtk_ppp.c:194
VAR_ACC
#define VAR_ACC
Definition mrtk_ppp.c:196
mrtk_ppp.h
MRTKLIB Precise Point Positioning (PPP) Module.
ppp_rtk_pos
void ppp_rtk_pos(rtk_t *rtk, const obsd_t *obs, int n, nav_t *nav)
PPP-RTK positioning for one epoch (CLAS).
Definition mrtk_ppp_rtk.c:2540
test_sys
static int test_sys(int sys, int m, int qzsmodear, int code)
Test if satellite system matches group index.
Definition mrtk_ppp_rtk.c:1089
ppp_rtk_na
int ppp_rtk_na(const prcopt_t *opt)
Compute number of non-ambiguity states for PPP-RTK.
Definition mrtk_ppp_rtk.c:2498
INIT_ZWD
#define INIT_ZWD
Definition mrtk_ppp_rtk.c:69
ppp_rtk_nx
int ppp_rtk_nx(const prcopt_t *opt)
Compute number of estimated states for PPP-RTK.
Definition mrtk_ppp_rtk.c:2496
mrtk_ppp_rtk.h
PPP-RTK positioning engine declarations.
LLI_HALFC
#define LLI_HALFC
Definition mrtk_rinex.c:57
errmsg
static void errmsg(rtk_t *rtk, const char *format,...)
Definition mrtk_rtkpos.c:382
selsat
static int selsat(const obsd_t *obs, double *azel, int nu, int nr, const prcopt_t *opt, int *sat, int *iu, int *ir)
Definition mrtk_rtkpos.c:495
rtkpos
int rtkpos(mrtk_ctx_t *ctx, rtk_t *rtk, const obsd_t *obs, int n, nav_t *nav)
Compute rover position by precise positioning.
Definition mrtk_rtkpos.c:2686
valpos
static int valpos(rtk_t *rtk, const double *v, const double *R, const int *vflg, int nv, double thres)
Definition mrtk_rtkpos.c:1935
INT_SWAP_STAT
#define INT_SWAP_STAT
Definition mrtk_rtkpos.c:64
NX
#define NX(opt)
Definition mrtk_rtkpos.c:103
sdobs
static double sdobs(const obsd_t *obs, int i, int j, int k)
Definition mrtk_rtkpos.c:397
udsatcb
static void udsatcb(gtime_t gt, nav_t *nav, osb_t *biaosb, const prcopt_t *popt)
Definition mrtk_rtkpos.c:2278
rtkoutstat
int rtkoutstat(rtk_t *rtk, char *buff)
Write solution status to buffer.
Definition mrtk_rtkpos.c:230
test_sys
static int test_sys(int sys, int m)
Definition mrtk_rtkpos.c:1246
NR
#define NR(opt)
Definition mrtk_rtkpos.c:102
holdamb
static void holdamb(rtk_t *rtk, const double *xa)
Definition mrtk_rtkpos.c:1621
EFACT_GPS
#define EFACT_GPS
Definition mrtk_rtkpos.c:60
miono_get_corr
int miono_get_corr(const double *rr, nav_t *nav)
Definition mrtk_ppp_iono.c:174
ddidx
static int ddidx(rtk_t *rtk, int *ix, int gps, int glo, int sbs)
Definition mrtk_rtkpos.c:1523
udbias
static void udbias(rtk_t *rtk, double tt, const obsd_t *obs, const int *sat, const int *iu, const int *ir, int ns, const nav_t *nav)
Definition mrtk_rtkpos.c:877
zdres_sat
static void zdres_sat(int base, double r, const obsd_t *obs, const nav_t *nav, const double *azel, const double *dant, const prcopt_t *opt, double *y, double *freq)
Definition mrtk_rtkpos.c:1027
EFACT_GLO
#define EFACT_GLO
Definition mrtk_rtkpos.c:61
udpos
static void udpos(rtk_t *rtk, double tt)
Definition mrtk_rtkpos.c:515
initx
static void initx(rtk_t *rtk, double xi, double var, int i)
Definition mrtk_rtkpos.c:487
baseline
static double baseline(const double *ru, const double *rb, double *dr)
Definition mrtk_rtkpos.c:479
detslp_code
static void detslp_code(rtk_t *rtk, const obsd_t *obs, int i, int rcv)
Definition mrtk_rtkpos.c:713
IL
#define IL(f, opt)
Definition mrtk_rtkpos.c:108
SSR_VENDOR_L6
#define SSR_VENDOR_L6
Definition mrtk_rtkpos.c:65
EFACT_SBS
#define EFACT_SBS
Definition mrtk_rtkpos.c:62
restamb
static void restamb(rtk_t *rtk, const double *bias, double *xa)
Definition mrtk_rtkpos.c:1588
udtrop
static void udtrop(rtk_t *rtk, double tt, double bl)
Definition mrtk_rtkpos.c:664
VAR_HOLDAMB
#define VAR_HOLDAMB
Definition mrtk_rtkpos.c:90
prectrop
static double prectrop(gtime_t time, const double *pos, int r, const double *azel, const prcopt_t *opt, const double *x, double *dtdx)
Definition mrtk_rtkpos.c:1223
rtkfree
void rtkfree(rtk_t *rtk)
Free memory for RTK control struct.
Definition mrtk_rtkpos.c:2261
udstate
static void udstate(rtk_t *rtk, const obsd_t *obs, const int *sat, const int *iu, const int *ir, int ns, const nav_t *nav)
Definition mrtk_rtkpos.c:999
intpres
static double intpres(gtime_t time, const obsd_t *obs, int n, const nav_t *nav, rtk_t *rtk, double *y)
Definition mrtk_rtkpos.c:1475
manage_amb_LAMBDA
static int manage_amb_LAMBDA(rtk_t *rtk, double *bias, double *xa, const int *sat, int nf, int ns)
Definition mrtk_rtkpos.c:1811
PRN_HWBIAS
#define PRN_HWBIAS
Definition mrtk_rtkpos.c:86
swapsolstat
static void swapsolstat(void)
Definition mrtk_rtkpos.c:315
rtkopenstat
int rtkopenstat(const char *file, int level)
Open solution status file and set output level.
Definition mrtk_rtkpos.c:193
MAXAGESSRL6
#define MAXAGESSRL6
Definition mrtk_rtkpos.c:67
varerr
static double varerr(int sat, int sys, double el, double snr_rover, double snr_base, double bl, double dt, int f, const prcopt_t *opt)
Definition mrtk_rtkpos.c:418
MIN
#define MIN(x, y)
Definition mrtk_rtkpos.c:76
VAR_HWBIAS
#define VAR_HWBIAS
Definition mrtk_rtkpos.c:82
udstacb
static void udstacb(gtime_t gt, nav_t *nav, osb_t *biaosb, const prcopt_t *popt)
Definition mrtk_rtkpos.c:2483
IT
#define IT(r, opt)
Definition mrtk_rtkpos.c:107
rtkclosestat
void rtkclosestat(void)
Close solution status file.
Definition mrtk_rtkpos.c:219
NF
#define NF(opt)
Definition mrtk_rtkpos.c:96
SSR_VENDOR_RTCM
#define SSR_VENDOR_RTCM
Definition mrtk_rtkpos.c:66
udbiass
static void udbiass(gtime_t gt, prcopt_t *popt, nav_t *nav)
Definition mrtk_rtkpos.c:2534
udsatpb
static void udsatpb(gtime_t gt, nav_t *nav, osb_t *biaosb, osb_t *fcbosb, const prcopt_t *popt)
Definition mrtk_rtkpos.c:2374
resamb_LAMBDA
static int resamb_LAMBDA(rtk_t *rtk, double *bias, double *xa, int gps, int glo, int sbs)
Definition mrtk_rtkpos.c:1675
statlevel
static int statlevel
Definition mrtk_rtkpos.c:119
detslp_dop
static void detslp_dop(rtk_t *rtk, const obsd_t *obs, const int *ix, int ns, int rcv, const nav_t *nav)
Definition mrtk_rtkpos.c:819
detslp_gf
static void detslp_gf(rtk_t *rtk, const obsd_t *obs, int i, int j, const nav_t *nav)
Definition mrtk_rtkpos.c:783
ddcov
static void ddcov(const int *nb, int n, const double *Ri, const double *Rj, int nv, double *R)
Definition mrtk_rtkpos.c:1157
zdres
static int zdres(int base, const obsd_t *obs, int n, const double *rs, const double *dts, const double *var, const int *svh, const nav_t *nav, const double *rr, const prcopt_t *opt, int index, double *y, double *e, double *azel, double *freq)
Definition mrtk_rtkpos.c:1077
SQR
#define SQR(x)
Definition mrtk_rtkpos.c:74
udrcvbias
static void udrcvbias(rtk_t *rtk, double tt)
Definition mrtk_rtkpos.c:692
seliflc
static int seliflc(int optnf, int sys)
Definition mrtk_rtkpos.c:417
ddres
static int ddres(rtk_t *rtk, const nav_t *nav, const obsd_t *obs, double dt, const double *x, const double *P, const int *sat, double *y, double *e, double *azel, double *freq, const int *iu, const int *ir, int ns, double *v, double *H, double *R, int *vflg)
Definition mrtk_rtkpos.c:1266
validobs
static int validobs(int i, int j, int f, int nf, double *y)
Definition mrtk_rtkpos.c:1151
IB
#define IB(s, f, opt)
Definition mrtk_rtkpos.c:109
outsolstat
static void outsolstat(rtk_t *rtk)
Definition mrtk_rtkpos.c:338
gfobs
static double gfobs(const obsd_t *obs, int i, int j, int k, const nav_t *nav)
Definition mrtk_rtkpos.c:403
file_stat
static char file_stat[1024]
Definition mrtk_rtkpos.c:121
SQRT
#define SQRT(x)
Definition mrtk_rtkpos.c:75
udion
static void udion(rtk_t *rtk, double tt, double bl, const int *sat, int ns)
Definition mrtk_rtkpos.c:638
rtkinit
void rtkinit(rtk_t *rtk, const prcopt_t *opt)
Initialize RTK control struct.
Definition mrtk_rtkpos.c:2206
relpos
static int relpos(rtk_t *rtk, const obsd_t *obs, int nu, int nr, const nav_t *nav)
Definition mrtk_rtkpos.c:1958
ar_poly_coeffs
static const double ar_poly_coeffs[3][5]
Definition mrtk_rtkpos.c:113
detslp_ll
static void detslp_ll(rtk_t *rtk, const obsd_t *obs, int i, int rcv)
Definition mrtk_rtkpos.c:735
constbl
static int constbl(rtk_t *rtk, const double *x, const double *P, double *v, double *H, double *Ri, double *Rj, int index)
Definition mrtk_rtkpos.c:1176
time_stat
static gtime_t time_stat
Definition mrtk_rtkpos.c:122
II
#define II(s, opt)
Definition mrtk_rtkpos.c:106
TTOL_MOVEB
#define TTOL_MOVEB
Definition mrtk_rtkpos.c:92
mrtk_rtkpos.h
MRTKLIB RTK Positioning Module — relative and precise positioning.
pntpos
int pntpos(mrtk_ctx_t *ctx, const obsd_t *obs, int n, const nav_t *nav, const prcopt_t *opt, sol_t *sol, double *azel, ssat_t *ssat, char *msg)
Compute receiver position, velocity, clock bias by single-point positioning with pseudorange and dopp...
Definition mrtk_spp.c:865
mrtk_spp.h
MRTKLIB Single-Point Positioning Module — SPP functions.
reppath
int reppath(const char *path, char *rpath, gtime_t time, const char *rov, const char *base)
Replace keywords in file path with date, time, rover and base station id.
Definition mrtk_sys.c:248
mrtk_sys.h
MRTKLIB System Module — OS/filesystem utility functions.
tidedisp
void tidedisp(gtime_t tutc, const double *rr, int opt, const erp_t *erp, const double *odisp, double *dr)
Compute displacements by earth tides.
Definition mrtk_tides.c:261
mrtk_tides.h
MRTKLIB Tides Module — Tidal displacement corrections.
time2str
void time2str(gtime_t t, char *s, int n)
Convert gtime_t struct to string.
Definition mrtk_time.c:507
gpst2utc
gtime_t gpst2utc(gtime_t t)
Convert GPS time to UTC considering leap seconds.
Definition mrtk_time.c:422
timediff
double timediff(gtime_t t1, gtime_t t2)
Compute time difference (t1 - t2).
Definition mrtk_time.c:279
time2gpst
double time2gpst(gtime_t t, int *week)
Convert gtime_t to GPS week/tow.
Definition mrtk_time.c:181
time_str
char * time_str(gtime_t t, int n)
Get time string (not reentrant).
Definition mrtk_time.c:530
utc2gpst
gtime_t utc2gpst(gtime_t t)
Convert UTC to GPS time considering leap seconds.
Definition mrtk_time.c:440
timeget
gtime_t timeget(void)
Get current time in UTC.
Definition mrtk_time.c:289
traceobs
void traceobs(mrtk_ctx_t *ctx, int level, const obsd_t *obs, int n)
Output trace observation data.
Definition mrtk_trace.c:122
tracemat
void tracemat(mrtk_ctx_t *ctx, int level, const double *A, int n, int m, int p, int q)
Output trace matrix.
Definition mrtk_trace.c:106
trace
void trace(mrtk_ctx_t *ctx, int level, const char *format,...)
Output trace message.
Definition mrtk_trace.c:78
mrtk_trace.h
MRTKLIB Trace Module — debug trace and logging functions.
relposvrs
int relposvrs(rtk_t *rtk, const obsd_t *obs, int nu, int *nr, nav_t *nav)
VRS-RTK relative positioning.
Definition mrtk_vrs.c:2163
mrtk_vrs.h
MRTKLIB VRS-RTK Module — double-differenced RTK using CLAS VRS.
ambc_t
Definition mrtk_rtkpos.h:40
clas_ctx_t
Definition mrtk_clas.h:580
clas_osrd_t
Definition mrtk_clas.h:477
gtime_t
GNSS time representation.
Definition mrtk_time.h:49
gtime_t::time
time_t time
Definition mrtk_time.h:50
mrtk_ctx_t
Runtime context for MRTKLIB processing.
Definition mrtk_context.h:73
nav_s::filreset
int filreset
Definition mrtk_nav.h:436
nav_s::fcbpath
char fcbpath[MAXSTRPATH]
Definition mrtk_nav.h:443
nav_s::erp
erp_t erp
Definition mrtk_nav.h:413
nav_s::clas_ctx
void * clas_ctx
Definition mrtk_nav.h:439
nav_s::cbias
double cbias[MAXSAT][3]
Definition mrtk_nav.h:427
nav_s::pr_biapath
char pr_biapath[MAXSTRPATH]
Definition mrtk_nav.h:444
nav_s::biapath
char biapath[MAXSTRPATH]
Definition mrtk_nav.h:442
nav_s::ssr_ch
ssr_t ssr_ch[SSR_CH_NUM][MAXSAT]
Definition mrtk_nav.h:434
nav_s::osb
osb_t osb
Definition mrtk_nav.h:441
nav_s::pr_fcbpath
char pr_fcbpath[MAXSTRPATH]
Definition mrtk_nav.h:445
obsd_t
Observation data record.
Definition mrtk_obs.h:44
obsd_t::D
float D[5+5]
Definition mrtk_obs.h:52
obsd_t::time
gtime_t time
Definition mrtk_obs.h:45
obsd_t::sat
uint8_t sat
Definition mrtk_obs.h:46
obsd_t::LLI
uint8_t LLI[5+5]
Definition mrtk_obs.h:48
obsd_t::L
double L[5+5]
Definition mrtk_obs.h:50
obsd_t::rcv
uint8_t rcv
Definition mrtk_obs.h:46
obsd_t::P
double P[5+5]
Definition mrtk_obs.h:51
obsd_t::SNR
uint16_t SNR[5+5]
Definition mrtk_obs.h:47
obsd_t::code
uint8_t code[5+5]
Definition mrtk_obs.h:49
osb_t
Definition mrtk_nav.h:198
osb_t::gt
gtime_t gt[2]
Definition mrtk_nav.h:199
osb_t::rsyscb
double rsyscb[MAXBSNXSYS][MAXCODE]
Definition mrtk_nav.h:203
osb_t::scb
double scb[MAXSAT][MAXCODE]
Definition mrtk_nav.h:201
osb_t::vrsatcb
int vrsatcb[MAXSAT][MAXCODE]
Definition mrtk_nav.h:208
osb_t::vspb
int vspb[MAXSAT][MAXCODE]
Definition mrtk_nav.h:206
osb_t::vscb
int vscb[MAXSAT][MAXCODE]
Definition mrtk_nav.h:205
osb_t::vrsyscb
int vrsyscb[MAXBSNXSYS][MAXCODE]
Definition mrtk_nav.h:207
osb_t::spb
double spb[MAXSAT][MAXCODE]
Definition mrtk_nav.h:202
osb_t::rsatcb
double rsatcb[MAXSAT][MAXCODE]
Definition mrtk_nav.h:204
prcopt_t
Definition mrtk_opt.h:171
prcopt_t::thresdop
double thresdop
Definition mrtk_opt.h:215
prcopt_t::snrmask
snrmask_t snrmask
Definition mrtk_opt.h:177
prcopt_t::minfixsats
int minfixsats
Definition mrtk_opt.h:289
prcopt_t::outsingle
int outsingle
Definition mrtk_opt.h:235
prcopt_t::niter
int niter
Definition mrtk_opt.h:191
prcopt_t::minfix
int minfix
Definition mrtk_opt.h:185
prcopt_t::elmin
double elmin
Definition mrtk_opt.h:176
prcopt_t::nf
int nf
Definition mrtk_opt.h:174
prcopt_t::bdsmodear
int bdsmodear
Definition mrtk_opt.h:181
prcopt_t::robustk
double robustk[2]
Definition mrtk_opt.h:304
prcopt_t::mindropsats
int mindropsats
Definition mrtk_opt.h:290
prcopt_t::ionoopt
int ionoopt
Definition mrtk_opt.h:187
prcopt_t::pcvr
pcv_t pcvr[2]
Definition mrtk_opt.h:224
prcopt_t::elmaskhold
double elmaskhold
Definition mrtk_opt.h:213
prcopt_t::thresslip
double thresslip
Definition mrtk_opt.h:214
prcopt_t::tdcpjump
double tdcpjump
Definition mrtk_opt.h:308
prcopt_t::prn
double prn[7]
Definition mrtk_opt.h:209
prcopt_t::minlock
int minlock
Definition mrtk_opt.h:184
prcopt_t::dynamics
int dynamics
Definition mrtk_opt.h:189
prcopt_t::err
double err[8]
Definition mrtk_opt.h:201
prcopt_t::floatcnt
int floatcnt
Definition mrtk_opt.h:257
prcopt_t::intpref
int intpref
Definition mrtk_opt.h:193
prcopt_t::eratio
double eratio[5]
Definition mrtk_opt.h:200
prcopt_t::tropopt
int tropopt
Definition mrtk_opt.h:188
prcopt_t::tdcp
int tdcp
Definition mrtk_opt.h:307
prcopt_t::pppsatcb
int pppsatcb
Definition mrtk_opt.h:229
prcopt_t::baseline
double baseline[2]
Definition mrtk_opt.h:219
prcopt_t::odisp
double odisp[2][6 *11]
Definition mrtk_opt.h:239
prcopt_t::arfilter
int arfilter
Definition mrtk_opt.h:292
prcopt_t::sclkstab
double sclkstab
Definition mrtk_opt.h:210
prcopt_t::glomodear
int glomodear
Definition mrtk_opt.h:180
prcopt_t::minholdsats
int minholdsats
Definition mrtk_opt.h:291
prcopt_t::rejeminsat
int rejeminsat
Definition mrtk_opt.h:321
prcopt_t::posopt
int posopt[13]
Definition mrtk_opt.h:237
prcopt_t::maxbiasdt
double maxbiasdt
Definition mrtk_opt.h:232
prcopt_t::thresar
double thresar[8]
Definition mrtk_opt.h:211
prcopt_t::elmaskopt
int16_t elmaskopt[360]
Definition mrtk_opt.h:242
prcopt_t::maxout
int maxout
Definition mrtk_opt.h:183
prcopt_t::rejethres
double rejethres
Definition mrtk_opt.h:320
prcopt_t::gpsmodear
int gpsmodear
Definition mrtk_opt.h:293
prcopt_t::elmaskar
double elmaskar
Definition mrtk_opt.h:212
prcopt_t::modear
int modear
Definition mrtk_opt.h:179
prcopt_t::mode
int mode
Definition mrtk_opt.h:172
prcopt_t::tidecorr
int tidecorr
Definition mrtk_opt.h:190
prcopt_t::enhanced_spp_seed
int enhanced_spp_seed
Definition mrtk_opt.h:317
prcopt_t::antdel
double antdel[2][3]
Definition mrtk_opt.h:223
prcopt_t::robust
int robust
Definition mrtk_opt.h:303
prcopt_t::ru
double ru[3]
Definition mrtk_opt.h:220
prcopt_t::std
double std[3]
Definition mrtk_opt.h:208
prcopt_t::navsys
int navsys
Definition mrtk_opt.h:175
prcopt_t::maxtdiff
double maxtdiff
Definition mrtk_opt.h:216
prcopt_t::sateph
int sateph
Definition mrtk_opt.h:178
prcopt_t::bds3igso
char bds3igso[256]
Definition mrtk_opt.h:334
prcopt_t::bds2prn
char bds2prn[256]
Definition mrtk_opt.h:333
prcopt_t::pppsatpb
int pppsatpb
Definition mrtk_opt.h:230
prcopt_t::rb
double rb[3]
Definition mrtk_opt.h:221
prcopt_t::maxinno
double maxinno
Definition mrtk_opt.h:217
prcopt_t::staname
char staname[32]
Definition mrtk_opt.h:246
prcopt_t::l6mrg
int l6mrg
Definition mrtk_opt.h:275
prcopt_t::refpos
int refpos
Definition mrtk_opt.h:197
prcopt_t::rtcmopt
char rtcmopt[256]
Definition mrtk_opt.h:244
rtk_t
Definition mrtk_rtkpos.h:49
rtk_t::Pa
double * Pa
Definition mrtk_rtkpos.h:56
rtk_t::tt
double tt
Definition mrtk_rtkpos.h:53
rtk_t::nx
int nx
Definition mrtk_rtkpos.h:52
rtk_t::prev_qr
double prev_qr[6]
Definition mrtk_rtkpos.h:65
rtk_t::nb_ar
int nb_ar
Definition mrtk_rtkpos.h:67
rtk_t::x
double * x
Definition mrtk_rtkpos.h:54
rtk_t::opt
prcopt_t opt
Definition mrtk_rtkpos.h:63
rtk_t::ambc
ambc_t ambc[MAXSAT]
Definition mrtk_rtkpos.h:59
rtk_t::excsat
int excsat
Definition mrtk_rtkpos.h:66
rtk_t::na
int na
Definition mrtk_rtkpos.h:52
rtk_t::P
double * P
Definition mrtk_rtkpos.h:54
rtk_t::neb
int neb
Definition mrtk_rtkpos.h:61
rtk_t::errbuf
char errbuf[MAXERRMSG]
Definition mrtk_rtkpos.h:62
rtk_t::xa
double * xa
Definition mrtk_rtkpos.h:56
rtk_t::nfix
int nfix
Definition mrtk_rtkpos.h:57
rtk_t::rb
double rb[6]
Definition mrtk_rtkpos.h:51
rtk_t::holdamb
char holdamb
Definition mrtk_rtkpos.h:68
rtk_t::ssat
ssat_t ssat[MAXSAT]
Definition mrtk_rtkpos.h:60
rtk_t::sol
sol_t sol
Definition mrtk_rtkpos.h:50
rtk_t::Q
double * Q
Definition mrtk_rtkpos.h:55
sol_t
Definition mrtk_sol.h:44
sol_t::ratio
float ratio
Definition mrtk_sol.h:57
sol_t::qr
float qr[6]
Definition mrtk_sol.h:48
sol_t::rr
double rr[6]
Definition mrtk_sol.h:46
sol_t::thres
float thres
Definition mrtk_sol.h:58
sol_t::prev_ratio2
float prev_ratio2
Definition mrtk_sol.h:60
sol_t::stat
uint8_t stat
Definition mrtk_sol.h:54
sol_t::dtr
double dtr[NSYS+1]
Definition mrtk_sol.h:52
sol_t::prev_ratio1
float prev_ratio1
Definition mrtk_sol.h:59
sol_t::age
float age
Definition mrtk_sol.h:56
sol_t::ns
uint8_t ns
Definition mrtk_sol.h:55
sol_t::qv
float qv[6]
Definition mrtk_sol.h:51
sol_t::time
gtime_t time
Definition mrtk_sol.h:45
ssat_t
Definition mrtk_sol.h:98
ssat_t::pt
gtime_t pt[2][5]
Definition mrtk_sol.h:116
ssat_t::vsat
uint8_t vsat[5]
Definition mrtk_sol.h:104
ssat_t::slipc
uint32_t slipc[5]
Definition mrtk_sol.h:111
ssat_t::vs
uint8_t vs
Definition mrtk_sol.h:100
ssat_t::ph
double ph[2][5]
Definition mrtk_sol.h:117
ssat_t::outc
uint32_t outc[5]
Definition mrtk_sol.h:110
ssat_t::rejc
uint32_t rejc[5]
Definition mrtk_sol.h:112
ssat_t::codeprev
uint8_t codeprev[5][2]
Definition mrtk_sol.h:126
ssat_t::azel
double azel[2]
Definition mrtk_sol.h:101
ssat_t::resp
double resp[5]
Definition mrtk_sol.h:102
ssat_t::lock
int lock[5]
Definition mrtk_sol.h:109
ssat_t::slip
uint8_t slip[5]
Definition mrtk_sol.h:107
ssat_t::gf
double gf[5]
Definition mrtk_sol.h:113
ssat_t::snr
uint16_t snr[5]
Definition mrtk_sol.h:105
ssat_t::half
uint8_t half[5]
Definition mrtk_sol.h:108
ssat_t::sys
uint8_t sys
Definition mrtk_sol.h:99
ssat_t::fix
uint8_t fix[5]
Definition mrtk_sol.h:106
ssat_t::resc
double resc[5]
Definition mrtk_sol.h:103
ssr_t::t0
gtime_t t0[6]
Definition mrtk_nav.h:176
ssr_t::pbias
double pbias[MAXCODE]
Definition mrtk_nav.h:188
ssr_t::cbias
float cbias[MAXCODE]
Definition mrtk_nav.h:187
ssr_t::vendor
int vendor
Definition mrtk_nav.h:195
calloc
#define calloc(x, y)
Definition toml.c:54
free
#define free(x)
Definition toml.c:53