MRTKLIB 0.7.14
Modernized GNSS Positioning Library
Loading...
Searching...
No Matches
mrtk_vrs.c
Go to the documentation of this file.
1/*------------------------------------------------------------------------------
2 * mrtk_vrs.c : VRS-RTK positioning engine
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) 2007-2023 T.TAKASU
12 *
13 * SPDX-License-Identifier: BSD-2-Clause
14 *----------------------------------------------------------------------------*/
15/**
16 * @file mrtk_vrs.c
17 * @brief VRS-RTK positioning engine.
18 *
19 * Notes:
20 * - Ported from upstream claslib/src/rtkvrs.c with zero algorithmic changes.
21 * - VRS uses double-differenced RTK with LAMBDA PAR and fix-and-hold.
22 * - VRS observations are pre-generated by ssr2osr utility.
23 */
24#include "mrtklib/mrtk_vrs.h"
25
26#include <math.h>
27#include <stdarg.h>
28#include <stdlib.h>
29#include <string.h>
30
31#include "mrtklib/mrtk_antenna.h"
32#include "mrtklib/mrtk_atmos.h"
33#include "mrtklib/mrtk_clas.h"
34#include "mrtklib/mrtk_const.h"
35#include "mrtklib/mrtk_context.h"
36#include "mrtklib/mrtk_coords.h"
37#include "mrtklib/mrtk_eph.h"
38#include "mrtklib/mrtk_foundation.h"
39#include "mrtklib/mrtk_lambda.h"
40#include "mrtklib/mrtk_mat.h"
41#include "mrtklib/mrtk_nav.h"
42#include "mrtklib/mrtk_obs.h"
43#include "mrtklib/mrtk_opt.h"
44#include "mrtklib/mrtk_rtkpos.h"
45#include "mrtklib/mrtk_sol.h"
46#include "mrtklib/mrtk_tides.h"
47#include "mrtklib/mrtk_time.h"
48#include "mrtklib/mrtk_trace.h"
49
50/* constants/macros ----------------------------------------------------------*/
51#define MAXREF 2 /* max iteration of searching reference sat */
52
53#define SQR(x) ((x) * (x))
54#define MAXITER 20 /* max iteration of PAR */
55#define SQRT(x) ((x) <= 0.0 ? 0.0 : sqrt(x))
56#define MIN(x, y) ((x) <= (y) ? (x) : (y))
57#define ROUND(x) (int)floor((x) + 0.5)
58
59#define VAR_POS SQR(30.0) /* initial variance of receiver pos ((m)^2) */
60#define VAR_VEL SQR(1.0) /* initial variance of receiver acc ((m/s)^2) */
61#define VAR_ACC SQR(1.0) /* initial variance of receiver acc ((m/ss)^2) */
62#define VAR_HWBIAS SQR(1.0) /* initial variance of h/w bias ((m/MHz)^2) */
63#define VAR_GRA SQR(0.001) /* initial variance of gradient (m^2) */
64#define INIT_ZWD 0.15 /* initial zwd (m) */
65
66#define PRN_HWBIAS 1E-6 /* process noise of h/w bias (m/MHz/sqrt(s)) */
67#define GAP_RESION 120 /* gap to reset ionosphere parameters (epochs) */
68#define MAXACC 30.0 /* max accel for doppler slip detection (m/s^2) */
69
70#define VAR_HOLDAMB 0.001 /* constraint to hold ambiguity (cycle^2) */
71
72#define TTOL_MOVEB (1.0 + 2 * DTTOL)
73/* time sync tolerance for moving-baseline (s) */
74
75/* frequency pair selection constants */
76#define POSL1 1 /* L1 single freq positioning */
77#define POSL1L2 2 /* L1+L2 dual freq positioning */
78#define POSL1L5 3 /* L1+L5 dual freq positioning */
79#define POSL1L2L5 4 /* L1+L2+L5 triple freq positioning */
80#define POSL1L5_L2 5 /* L1+L5(or L2) dual freq positioning */
81
82/* carrier wavelength constants */
83#define LAM_CARR_L1 (CLIGHT / FREQ1)
84#define LAM_CARR_L2 (CLIGHT / FREQ2)
85#define LAM_CARR_L5 (CLIGHT / FREQ5)
86
87static const double lam_carr[3] = {LAM_CARR_L1, LAM_CARR_L2, LAM_CARR_L5};
88
89/* number of parameters (pos,ionos,tropos,hw-bias,phase-bias,real,estimated) */
90#define NF(opt) ((opt)->ionoopt == IONOOPT_IFLC ? 1 : (opt)->nf)
91#define NP(opt) ((opt)->dynamics == 0 ? 3 : 9)
92#define NI(opt) ((((opt)->ionoopt != IONOOPT_EST) && ((opt)->ionoopt != IONOOPT_EST_ADPT)) ? 0 : MAXSAT)
93#define NT(opt) ((opt)->tropopt < TROPOPT_EST ? 0 : ((opt)->tropopt < TROPOPT_ESTG ? 2 : 6))
94#define NL(opt) ((opt)->glomodear != 2 ? 0 : NFREQGLO)
95#define NB(opt) ((opt)->mode <= PMODE_DGPS ? 0 : MAXSAT * NF(opt))
96#define NR(opt) (NP(opt) + NI(opt) + NT(opt) + NL(opt))
97#define NX(opt) (NR(opt) + NB(opt))
98
99/* state variable index */
100#define II(s, opt) (NP(opt) + (s) - 1) /* ionos (s:satellite no) */
101#define IT(r, opt) (NP(opt) + NI(opt) + NT(opt) / 2 * (r)) /* tropos (r:0=rov,1:ref) */
102#define IL(f, opt) (NP(opt) + NI(opt) + NT(opt) + (f)) /* receiver h/w bias */
103#define IB(s, f, opt) (NR(opt) + MAXSAT * (f) + (s) - 1) /* phase bias (s:satno,f:freq) */
104#define ITT(opt) (NP(opt) + NI(opt)) /* tropos (r:0=rov,1:ref) */
105
106enum ddres_order { FST = 1, SND, TRD };
107
108/* chi-squared table (alpha=0.001) */
109static const double chisqr[100] = {
110 10.8, 13.8, 16.3, 18.5, 20.5, 22.5, 24.3, 26.1, 27.9, 29.6, 31.3, 32.9, 34.5, 36.1, 37.7, 39.3, 40.8,
111 42.3, 43.8, 45.3, 46.8, 48.3, 49.7, 51.2, 52.6, 54.1, 55.5, 56.9, 58.3, 59.7, 61.1, 62.5, 63.9, 65.2,
112 66.6, 68.0, 69.3, 70.7, 72.1, 73.4, 74.7, 76.0, 77.3, 78.6, 80.0, 81.4, 82.7, 84.0, 85.4, 86.7, 88.0,
113 89.3, 90.6, 91.9, 93.3, 94.7, 96.0, 97.4, 98.7, 100, 101, 102, 103, 104, 105, 107, 108, 109,
114 110, 111, 112, 113, 114, 115, 116, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127,
115 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142};
116
117/* F-distribution table for ratio test */
118static const double qfrtk[6][60] = {
119 /* significance level, alpha=0.001(0.1%) */
120 {405284.07, 999.00, 141.11, 53.44, 29.75, 20.03, 15.02, 12.05, 10.11, 8.75, 7.76, 7.00, 6.41, 5.93, 5.54,
121 5.20, 4.92, 4.68, 4.47, 4.29, 4.13, 3.98, 3.85, 3.74, 3.63, 3.53, 3.44, 3.36, 3.29, 3.22,
122 3.15, 3.09, 3.04, 2.98, 2.93, 2.89, 2.84, 2.80, 2.76, 2.73, 2.69, 2.66, 2.63, 2.60, 2.57,
123 2.54, 2.51, 2.49, 2.46, 2.44, 2.42, 2.40, 2.38, 2.36, 2.34, 2.32, 2.30, 2.28, 2.27, 2.25},
124 /* significance level, alpha=0.005(0.5%) */
125 {16210.72, 199.00, 47.47, 23.15, 14.94, 11.07, 8.89, 7.50, 6.54, 5.85, 5.32, 4.91, 4.57, 4.30, 4.07,
126 3.87, 3.71, 3.56, 3.43, 3.32, 3.22, 3.12, 3.04, 2.97, 2.90, 2.84, 2.78, 2.72, 2.67, 2.63,
127 2.58, 2.54, 2.51, 2.47, 2.44, 2.41, 2.38, 2.35, 2.32, 2.30, 2.27, 2.25, 2.23, 2.21, 2.19,
128 2.17, 2.15, 2.13, 2.11, 2.10, 2.08, 2.07, 2.05, 2.04, 2.02, 2.01, 2.00, 1.99, 1.97, 1.96},
129 /* significance level, alpha=0.010(1.0%) */
130 {4052.18, 99.00, 29.46, 15.98, 10.97, 8.47, 6.99, 6.03, 5.35, 4.85, 4.46, 4.16, 3.91, 3.70, 3.52,
131 3.37, 3.24, 3.13, 3.03, 2.94, 2.86, 2.78, 2.72, 2.66, 2.60, 2.55, 2.51, 2.46, 2.42, 2.39,
132 2.35, 2.32, 2.29, 2.26, 2.23, 2.21, 2.18, 2.16, 2.14, 2.11, 2.09, 2.08, 2.06, 2.04, 2.02,
133 2.01, 1.99, 1.98, 1.96, 1.95, 1.94, 1.92, 1.91, 1.90, 1.89, 1.88, 1.87, 1.86, 1.85, 1.84},
134 /* significance level, alpha=0.050(5.0%) */
135 {161.45, 19.00, 9.28, 6.39, 5.05, 4.28, 3.79, 3.44, 3.18, 2.98, 2.82, 2.69, 2.58, 2.48, 2.40,
136 2.33, 2.27, 2.22, 2.17, 2.12, 2.08, 2.05, 2.01, 1.98, 1.96, 1.93, 1.90, 1.88, 1.86, 1.84,
137 1.82, 1.80, 1.79, 1.77, 1.76, 1.74, 1.73, 1.72, 1.70, 1.69, 1.68, 1.67, 1.66, 1.65, 1.64,
138 1.63, 1.62, 1.62, 1.61, 1.60, 1.59, 1.58, 1.58, 1.57, 1.56, 1.56, 1.55, 1.55, 1.54, 1.53},
139 /* significance level, alpha=0.100(10.0%) */
140 {39.86, 9.00, 5.39, 4.11, 3.45, 3.05, 2.78, 2.59, 2.44, 2.32, 2.23, 2.15, 2.08, 2.02, 1.97,
141 1.93, 1.89, 1.85, 1.82, 1.79, 1.77, 1.74, 1.72, 1.70, 1.68, 1.67, 1.65, 1.63, 1.62, 1.61,
142 1.59, 1.58, 1.57, 1.56, 1.55, 1.54, 1.53, 1.52, 1.51, 1.51, 1.50, 1.49, 1.48, 1.48, 1.47,
143 1.46, 1.46, 1.45, 1.45, 1.44, 1.44, 1.43, 1.43, 1.42, 1.42, 1.41, 1.41, 1.40, 1.40, 1.40},
144 /* significance level, alpha=0.200(20.0%) */
145 {9.47, 4.00, 2.94, 2.48, 2.23, 2.06, 1.94, 1.86, 1.79, 1.73, 1.69, 1.65, 1.61, 1.58, 1.56,
146 1.54, 1.52, 1.50, 1.48, 1.47, 1.45, 1.44, 1.43, 1.42, 1.41, 1.40, 1.39, 1.38, 1.37, 1.36,
147 1.36, 1.35, 1.34, 1.34, 1.33, 1.33, 1.32, 1.32, 1.31, 1.31, 1.30, 1.30, 1.30, 1.29, 1.29,
148 1.28, 1.28, 1.28, 1.27, 1.27, 1.27, 1.26, 1.26, 1.26, 1.26, 1.25, 1.25, 1.25, 1.25, 1.24}};
149
150/*============================================================================
151 * Satellite wavelength helper (MRTKLIB uses sat2freq, not nav->lam)
152 *===========================================================================*/
153static uint8_t obs_code_for_freq_vrs(int sat, int f) {
154 int sys = satsys(sat, NULL);
155 switch (f) {
156 case 0:
157 return (sys == SYS_CMP) ? CODE_L2I : CODE_L1C;
158 case 1:
159 if (sys == SYS_GPS) {
160 return CODE_L2W;
161 }
162 if (sys == SYS_QZS) {
163 return CODE_L2S;
164 }
165 if (sys == SYS_GAL) {
166 return CODE_L5Q;
167 }
168 if (sys == SYS_CMP) {
169 return CODE_L7I;
170 }
171 return CODE_L2W;
172 case 2:
173 if (sys == SYS_CMP) {
174 return CODE_L6I;
175 }
176 return CODE_L5Q;
177 }
178 return CODE_L1C;
179}
180
181static double sat_lambda(int sat, int f, const nav_t* nav) {
182 double freq = sat2freq(sat, obs_code_for_freq_vrs(sat, f), nav);
183 return freq > 0.0 ? CLIGHT / freq : 0.0;
184}
185
186/* module-level chi-squared validation value (sol_t has no chisq field) */
187static double vrs_chisq = 0.0;
188
189/* clearsatcorr stub (no-op in MRTKLIB) --------------------------------------*/
190static void clearsatcorr(void) { /* no-op */ }
191
192/* save error message --------------------------------------------------------*/
193static void errmsg(rtk_t* rtk, const char* format, ...) {
194 char buff[256], tstr[32];
195 int n;
196 va_list ap;
197 time2str(rtk->sol.time, tstr, 2);
198 n = sprintf(buff, "%s: ", tstr + 11);
199 va_start(ap, format);
200 n += vsprintf(buff + n, format, ap);
201 va_end(ap);
202 n = n < MAXERRMSG - rtk->neb ? n : MAXERRMSG - rtk->neb;
203 memcpy(rtk->errbuf + rtk->neb, buff, n);
204 rtk->neb += n;
205 trace(NULL, 2, "%s", buff);
206}
207/* single-differenced measurement error variance -----------------------------*/
208static double varerr(int sat, int sys, double el, double bl, double dt, int f, const prcopt_t* opt) {
209 double a, b, c = opt->err[3] * bl / 1E4, d = CLIGHT * opt->sclkstab * dt, fact = 1.0;
210 double sinel = sin(el);
211 int i = sys == SYS_GLO ? 1 : (sys == SYS_GAL ? 2 : 0), nf = NF(opt);
212
213 /* extended error model */
214 if (f >= nf && opt->exterr.ena[0]) { /* code */
215 a = opt->exterr.cerr[i][(f - nf) * 2];
216 b = opt->exterr.cerr[i][1 + (f - nf) * 2];
217 if (sys == SYS_SBS) {
218 a *= EFACT_SBS;
219 b *= EFACT_SBS;
220 }
221 } else if (f < nf && opt->exterr.ena[1]) { /* phase */
222 a = opt->exterr.perr[i][f * 2];
223 b = opt->exterr.perr[i][1 + f * 2];
224 if (sys == SYS_SBS) {
225 a *= EFACT_SBS;
226 b *= EFACT_SBS;
227 }
228 } else { /* normal error model */
229 if (f >= nf) {
230 fact = opt->eratio[f - nf];
231 }
232 if (fact <= 0.0) {
233 fact = opt->eratio[0];
234 }
235 fact *= sys == SYS_GLO ? EFACT_GLO : (sys == SYS_SBS ? EFACT_SBS : EFACT_GPS);
236 a = fact * opt->err[1];
237 b = fact * opt->err[2];
238 c = d = 0.0;
239 }
240
241 return 1.0 * (opt->ionoopt == IONOOPT_IFLC ? 3.0 : 1.0) * (a * a + b * b / sinel / sinel + c * c) + d * d;
242}
243/* baseline length -----------------------------------------------------------*/
244static double baseline(const double* ru, const double* rb, double* dr) {
245 int i;
246 for (i = 0; i < 3; i++) {
247 dr[i] = ru[i] - rb[i];
248 }
249 return norm(dr, 3);
250}
251/* initialize state and covariance -------------------------------------------*/
252static void initx(rtk_t* rtk, double xi, double var, int i) {
253 int j;
254 rtk->x[i] = xi;
255 for (j = 0; j < rtk->nx; j++) {
256 rtk->P[i + j * rtk->nx] = rtk->P[j + i * rtk->nx] = i == j ? var : 0.0;
257 }
258}
259/* select common satellites between rover and reference station --------------*/
260static int selsat(const obsd_t* obs, double* azel, int nu, int* nr, const prcopt_t* opt, int* sat, int* iu, int* ir,
261 int ch) {
262 int i, j, k = 0, st, end;
263
264 trace(NULL, 3, "selsat : nu=%d nr[0]=%d nr[1]=%d ch=%d\n", nu, nr[0], nr[1], ch);
265 for (i = 0; i < MAXSAT; i++) {
266 sat[i] = -1;
267 ir[i] = -1;
268 iu[i] = -1;
269 }
270
271 st = (ch == 0) ? nu : nu + nr[0];
272 end = (ch == 0) ? nu + nr[0] : nu + nr[0] + nr[1];
273 for (i = 0, j = st; i < nu && j < end; i++, j++) {
274 if (obs[i].sat < obs[j].sat) {
275 j--;
276 } else if (obs[i].sat > obs[j].sat) {
277 i--;
278 } else if (azel[1 + j * 2] >= opt->elmin) { /* elevation at base station */
279 if (timediff(obs[i].time, obs[j].time) < -1.0) {
280 trace(NULL, 5, "ch:%d obs[i].sat:%d dt:%f\n", ch, obs[i].sat, timediff(obs[i].time, obs[j].time));
281 continue;
282 }
283 sat[k] = obs[i].sat;
284 iu[k] = i;
285 ir[k++] = j;
286 trace(NULL, 4, "(%2d) sat=%3d iu=%2d ir=%2d\n", k - 1, obs[i].sat, i, j);
287 }
288 }
289 return k;
290}
291static int selfreqpair(const int sat, const prcopt_t* opt, const obsd_t* obs) {
292 int optf = opt->posopt[10];
293 if (NFREQ == 1 || optf == POSL1) {
294 return 0;
295 }
296 if (NFREQ >= 3) {
297 if ((satsys(sat, NULL)) & (SYS_GAL)) {
298 return 2; /* GAL */
299 }
300 if (optf == POSL1L2L5) {
301 return 1 + 2;
302 }
303 if (optf == POSL1L5) {
304 return 2;
305 }
306 if (optf == POSL1L5_L2 && obs->L[2] != 0.0 && obs->P[2] != 0.0) {
307 return 2;
308 }
309 }
310 return 1;
311}
312/* temporal update of position/velocity/acceleration -------------------------*/
313static void udpos(rtk_t* rtk, const obsd_t* obs, double tt) {
314 double *F, *FP, *xp, pos[3], Q[9] = {0}, Qv[9], var = 0.0;
315 double prnpos_h, prnpos_v;
316 int i, j, na, nb, nx, ni, flag_init = 0;
317 double *Fx, *Pxx, *Pxi, *Pxb, *Pxi_, *Pxb_, *xp_, ki;
318
319 trace(NULL, 3, "udpos : tt=%.3f\n", tt);
320
321 nx = rtk->nx;
322
323 /* fixed mode */
324 if (rtk->opt.mode == PMODE_FIXED) {
325 for (i = 0; i < 3; i++) {
326 initx(rtk, rtk->opt.ru[i], 1E-8, i);
327 }
328 return;
329 }
330 /* initialize position for first epoch */
331 if (norm(rtk->x, 3) <= 0.0) {
332 for (i = 0; i < 3; i++) {
333 initx(rtk, rtk->sol.rr[i], VAR_POS, i);
334 }
335 if (rtk->opt.dynamics) {
336 for (i = 3; i < 6; i++) {
337 initx(rtk, rtk->sol.rr[i], VAR_VEL, i);
338 }
339 for (i = 6; i < 9; i++) {
340 initx(rtk, 1E-6, VAR_ACC, i);
341 }
342 }
343 flag_init = 1;
344 }
345
346 /* check variance of estimated postion */
347 for (i = 0; i < 3; i++) {
348 var += rtk->P[i + i * nx];
349 }
350 var /= 3.0;
351
352 if (var > VAR_POS) {
353 /* reset position with large variance */
354 for (i = 0; i < 3; i++) {
355 initx(rtk, rtk->sol.rr[i], VAR_POS, i);
356 }
357 for (i = 3; i < 6; i++) {
358 initx(rtk, rtk->sol.rr[i], VAR_VEL, i);
359 }
360 for (i = 6; i < 9; i++) {
361 initx(rtk, 1E-6, VAR_ACC, i);
362 }
363 trace(NULL, 2, "reset rtk position due to large variance: var=%.3f\n", var);
364 return;
365 }
366
367 /* state transition of position/velocity/acceleration */
368 if (rtk->opt.ionoopt >= IONOOPT_EST) {
369 F = eye(nx);
370
371 /* dynamics on */
372 if (rtk->opt.dynamics) {
373 for (i = 0; i < 6; i++) {
374 F[i + (i + 3) * nx] = tt;
375 }
376 for (i = 0; i < 3; i++) {
377 F[i + (i + 6) * nx] = SQR(tt) / 2.0;
378 }
379 }
380
381 na = NP(&rtk->opt);
382 ni = MAXSAT;
383 nb = nx - na - ni;
384 ki = rtk->opt.beta > 0.0 ? exp(-tt / rtk->opt.beta) : 1.0;
385 Fx = mat(na, na);
386 xp = mat(na, 1);
387 xp_ = mat(na, 1);
388 FP = mat(na, na);
389 Pxx = mat(na, na);
390 Pxi = mat(na, ni);
391 Pxi_ = mat(na, ni);
392 Pxb = mat(na, nb);
393 Pxb_ = mat(na, nb);
394
395 for (i = 0; i < na; i++) {
396 for (j = 0; j < na; j++) {
397 Pxx[i + j * na] = rtk->P[i + j * nx];
398 }
399 }
400 for (i = 0; i < na; i++) {
401 for (j = 0; j < ni; j++) {
402 Pxi_[i + j * na] = 0.5 * (rtk->P[i + (na + j) * nx] + rtk->P[na + j + (i)*nx]);
403 }
404 }
405 for (i = 0; i < na; i++) {
406 for (j = 0; j < nb; j++) {
407 Pxb_[i + j * na] = 0.5 * (rtk->P[i + (na + ni + j) * nx] + rtk->P[na + ni + j + (i)*nx]);
408 }
409 }
410 for (i = 0; i < na; i++) {
411 for (j = 0; j < na; j++) {
412 Fx[i + j * na] = F[i + j * nx];
413 }
414 }
415 for (i = 0; i < na; i++) {
416 xp_[i] = rtk->x[i];
417 }
418
419 /* x=F*x */
420 matmul("NN", na, 1, na, 1.0, Fx, xp_, 0.0, xp);
421 for (i = 0; i < na; i++) {
422 rtk->x[i] = xp[i];
423 }
424 for (i = 0; i < ni; i++) {
425 rtk->x[i + na] *= ki;
426 }
427
428 /* Pxx=Fx*Pxx*Fx */
429 matmul("NN", na, na, na, 1.0, Fx, Pxx, 0.0, FP);
430 matmul("NT", na, na, na, 1.0, FP, Fx, 0.0, Pxx);
431 for (i = 0; i < na; i++) {
432 for (j = 0; j < na; j++) {
433 rtk->P[i + j * nx] = Pxx[i + j * na];
434 }
435 }
436 /* Pxi=k*Fx*Pxi*/
437 matmul("NN", na, ni, na, ki, Fx, Pxi_, 0.0, Pxi);
438 for (i = 0; i < na; i++) {
439 for (j = 0; j < ni; j++) {
440 rtk->P[i + (na + j) * nx] = rtk->P[(na + j) + i * nx] = Pxi[i + j * na];
441 }
442 }
443 /* Pxb=Fx*Pxb*/
444 matmul("NN", na, nb, na, 1.0, Fx, Pxb_, 0.0, Pxb);
445 for (i = 0; i < na; i++) {
446 for (j = 0; j < nb; j++) {
447 rtk->P[i + (na + ni + j) * nx] = rtk->P[(na + ni + j) + i * nx] = Pxb[i + j * na];
448 }
449 }
450
451 for (i = 0; i < ni; i++) {
452 for (j = 0; j < ni; j++) {
453 rtk->P[(i + na) + (j + na) * nx] *= (ki * ki);
454 }
455 }
456 for (i = 0; i < ni; i++) {
457 for (j = 0; j < nb; j++) {
458 rtk->P[(i + na) + (j + na + ni) * nx] *= ki;
459 }
460 }
461 for (i = 0; i < nb; i++) {
462 for (j = 0; j < ni; j++) {
463 rtk->P[(i + na + ni) + (j + na) * nx] *= ki;
464 }
465 }
466
467 free(Pxx);
468 free(Pxi);
469 free(Pxb);
470 free(Pxi_);
471 free(Pxb_);
472 free(Fx);
473 free(xp_);
474 free(F);
475 free(FP);
476 free(xp);
477 } else if (rtk->opt.dynamics) {
478 /* calculate state transition of position/velocity/acceleration */
479 /* x=Fx, P=FPF' */
480 F = eye(9);
481 xp = mat(9, 1);
482 FP = mat(nx, nx);
483 /* generate F matrix to update position and velocity */
484 for (i = 0; i < 6; i++) {
485 F[i + (i + 3) * 9] = tt;
486 }
487 for (i = 0; i < 3; i++) {
488 F[i + (i + 6) * 9] = SQR(tt) / 2.0;
489 }
490
491 /* x=F*x, only calculate pos/vel/acc states to save time, the rest are unchanged */
492 matmul("NN", 9, 1, 9, 1.0, F, rtk->x, 0.0, xp);
493 matcpy(rtk->x, xp, 9, 1);
494 /* P=F*P */
495 matcpy(FP, rtk->P, nx, nx);
496 for (j = 0; j < nx; j++) {
497 for (i = 0; i < 6; i++) {
498 FP[i + j * nx] += rtk->P[i + 3 + j * nx] * tt;
499 }
500 }
501 /* P=FP*F'*/
502 matcpy(rtk->P, FP, nx, nx);
503 for (j = 0; j < nx; j++) {
504 for (i = 0; i < 6; i++) {
505 rtk->P[j + i * nx] += rtk->P[j + (i + 3) * nx] * tt;
506 }
507 }
508 free(F);
509 free(FP);
510 free(xp);
511 }
512
513 if (rtk->opt.dynamics) {
514 /* process noise added to only acceleration */
515 Q[0] = Q[4] = SQR(rtk->opt.prn[3]) * fabs(tt);
516 Q[8] = SQR(rtk->opt.prn[4]) * fabs(tt);
517 ecef2pos(rtk->x, pos);
518 covecef(pos, Q, Qv);
519 for (i = 0; i < 3; i++) {
520 for (j = 0; j < 3; j++) {
521 rtk->P[i + 6 + (j + 6) * nx] += Qv[i + j * 3];
522 }
523 }
524 if (rtk->opt.prnadpt && flag_init == 0) {
525 for (i = 0; i < 9; i++) {
526 for (j = 0; j < 9; j++) {
527 rtk->P[i + j * nx] += rtk->Q[i + j * nx] * fabs(tt);
528 }
529 }
530 }
531 } else {
532 /* process noise added to only position */
533 /* Keep CLAS position noise independent of PPP's IFB slot prn[6]. */
534 prnpos_h = rtk->opt.stats_prnposith;
535 prnpos_v = rtk->opt.stats_prnpositv;
536 Q[0] = Q[4] = SQR(prnpos_h) * fabs(tt);
537 Q[8] = SQR(prnpos_v) * fabs(tt);
538 ecef2pos(rtk->x, pos);
539 covecef(pos, Q, Qv);
540 if (rtk->opt.prnadpt == 0 || flag_init) {
541 for (i = 0; i < 3; i++) {
542 for (j = 0; j < 3; j++) {
543 rtk->P[i + j * nx] += Qv[i + j * 3];
544 }
545 }
546 } else {
547 for (i = 0; i < 3; i++) {
548 for (j = 0; j < 3; j++) {
549 rtk->P[i + j * nx] += rtk->Q[i + j * nx] * fabs(tt);
550 }
551 }
552 }
553 }
554}
555/* temporal update of ionospheric parameters ---------------------------------*/
556static void udion(rtk_t* rtk, double tt, double bl, const obsd_t* obs, int ns) {
557 double el, fact, qi;
558 int i, j, f, k, m;
559
560 trace(NULL, 3, "udion : tt=%.1f bl=%.0f ns=%d\n", tt, bl, ns);
561
562 for (i = 1; i <= MAXSAT; i++) {
563 j = II(i, &rtk->opt);
564 if (rtk->x[j] == 0.0) {
565 continue;
566 }
567 for (k = 0; k < ns; k++) {
568 if (obs[k].sat == i) {
569 break;
570 }
571 }
572 m = selfreqpair(i, &rtk->opt, obs + k);
573 for (f = 0; f < rtk->opt.nf; f++) {
574 if (f > 0 && !(f & m)) {
575 continue;
576 }
577 if ((rtk->ssat[i - 1].outc[f] > (unsigned int)rtk->opt.maxout)) {
578 trace(NULL, 3, "reset iono estimation, sat=%2d outc=%6d \n", i, rtk->ssat[i - 1].outc[f]);
579 rtk->x[j] = 0.0;
580 }
581 }
582 }
583 for (i = 0; i < ns; i++) {
584 j = II(obs[i].sat, &rtk->opt);
585 if (rtk->x[j] == 0.0) {
586 initx(rtk, 1E-6, SQR(rtk->opt.std[1]), j);
587 rtk->Q[j + j * rtk->nx] = 0.0;
588 } else {
589 /* elevation dependent factor of process noise */
590 el = rtk->ssat[obs[i].sat - 1].azel[1];
591 fact = cos(el);
592 if (rtk->opt.ionoopt == IONOOPT_EST_ADPT) {
593 if (rtk->Q[j + j * rtk->nx] == 0.0) {
594 rtk->Q[j + j * rtk->nx] = SQR(rtk->opt.prn[1] * fact);
595 } else {
596 if (rtk->Q[j + j * rtk->nx] > SQR(rtk->opt.prnionomax)) {
597 rtk->Q[j + j * rtk->nx] = SQR(rtk->opt.prnionomax);
598 } else if (rtk->Q[j + j * rtk->nx] < SQR(rtk->opt.prn[1])) {
599 rtk->Q[j + j * rtk->nx] = SQR(rtk->opt.prn[1]);
600 }
601 }
602 qi = rtk->Q[j + j * rtk->nx];
603 } else {
604 qi = SQR(rtk->opt.prn[1] * fact);
605 rtk->Q[j + j * rtk->nx] = qi;
606 }
607
608 rtk->P[j + j * rtk->nx] += qi * tt;
609 }
610 }
611}
612/* temporal update of tropospheric parameters --------------------------------*/
613static void udtrop(rtk_t* rtk, double tt, double bl) {
614 int i, j, k;
615
616 trace(NULL, 3, "udtrop : tt=%.1f\n", tt);
617
618 for (i = 0; i < 2; i++) {
619 j = IT(i, &rtk->opt);
620
621 if (rtk->x[j] == 0.0) {
622 initx(rtk, INIT_ZWD, SQR(rtk->opt.std[2]), j); /* initial zwd */
623
624 if (rtk->opt.tropopt >= TROPOPT_ESTG) {
625 for (k = 0; k < 2; k++) {
626 initx(rtk, 1E-6, VAR_GRA, ++j);
627 }
628 }
629 } else {
630 rtk->P[j + j * rtk->nx] += SQR(rtk->opt.prn[2]) * tt;
631
632 if (rtk->opt.tropopt >= TROPOPT_ESTG) {
633 for (k = 0; k < 2; k++) {
634 rtk->P[++j * (1 + rtk->nx)] += SQR(rtk->opt.prn[2] * 0.3) * fabs(rtk->tt);
635 }
636 }
637 }
638 }
639}
640/* temporal update of receiver h/w biases ------------------------------------*/
641static void udrcvbias(rtk_t* rtk, double tt) {
642 int i, j;
643
644 trace(NULL, 3, "udrcvbias: tt=%.1f\n", tt);
645
646 for (i = 0; i < NFREQGLO; i++) {
647 j = IL(i, &rtk->opt);
648
649 if (rtk->x[j] == 0.0) {
650 initx(rtk, 1E-6, VAR_HWBIAS, j);
651 }
652 /* hold to fixed solution */
653 else if (rtk->nfix >= rtk->opt.minfix && rtk->sol.ratio > rtk->opt.thresar[0]) {
654 initx(rtk, rtk->xa[j], rtk->Pa[j + j * rtk->na], j);
655 } else {
656 rtk->P[j + j * rtk->nx] += SQR(PRN_HWBIAS) * tt;
657 }
658 }
659}
660/* L1/L2 geometry-free phase measurement -------------------------------------*/
661static double gfmeasL1L2(const obsd_t* obs, const nav_t* nav) {
662 double lam0 = sat_lambda(obs->sat, 0, nav);
663 double lam1 = sat_lambda(obs->sat, 1, nav);
664 if (lam0 == 0.0 || lam1 == 0.0 || obs->L[0] == 0.0 || obs->L[1] == 0.0) {
665 return 0.0;
666 }
667 return lam0 * obs->L[0] - lam1 * obs->L[1];
668}
669/* L1/L5 geometry-free phase measurement -------------------------------------*/
670static double gfmeasL1L5(const obsd_t* obs, const nav_t* nav) {
671 double lam0 = sat_lambda(obs->sat, 0, nav);
672 double lam2 = sat_lambda(obs->sat, 2, nav);
673 if (lam0 == 0.0 || lam2 == 0.0 || obs->L[0] == 0.0 || obs->L[2] == 0.0) {
674 return 0.0;
675 }
676 return lam0 * obs->L[0] - lam2 * obs->L[2];
677}
678
679/* detect cycle slip by LLI --------------------------------------------------*/
680static void detslp_ll(rtk_t* rtk, const obsd_t* obs, int n) {
681 int i, j;
682
683 for (i = 0; i < n && i < MAXOBS; i++) {
684 for (j = 0; j < rtk->opt.nf; j++) {
685 if (obs[i].L[j] == 0.0 || !(obs[i].LLI[j] & 3)) {
686 continue;
687 }
688
689 trace(NULL, 2, "detslp_ll: slip detected sat=%2d f=%d\n", obs[i].sat, j + 1);
690
691 rtk->ssat[obs[i].sat - 1].slip[j] = 1;
692 }
693 }
694}
695
696/* detect cycle slip by geometry free phase jump -----------------------------*/
697static void detslp_gf(rtk_t* rtk, const obsd_t* obs, int n, const nav_t* nav) {
698 double g0, g1;
699 int i;
700
701 trace(NULL, 3, "detslip_gf_L1_L2\n");
702
703 /* slip detection for L1-L2 gf */
704 for (i = 0; i < n && i < MAXOBS; i++) {
705 if ((g1 = gfmeasL1L2(obs + i, nav)) == 0.0) {
706 continue;
707 }
708 g0 = rtk->ssat[obs[i].sat - 1].gf[0];
709 rtk->ssat[obs[i].sat - 1].gf[0] = g1;
710 if (g0 != 0.0 && fabs(g1 - g0) > rtk->opt.thresslip) {
711 trace(NULL, 2, "detslip_gf_L1_L2: slip detected sat=%2d gf=%8.3f->%8.3f\n", obs[i].sat, g0, g1);
712 rtk->ssat[obs[i].sat - 1].slip[0] |= 1;
713 rtk->ssat[obs[i].sat - 1].slip[1] |= 1;
714 }
715 }
716
717 trace(NULL, 3, "detslip_gf_L1_L5\n");
718
719 /* slip detection for L1-L5 gf */
720 if (rtk->opt.nf <= 2) {
721 return;
722 }
723 for (i = 0; i < n && i < MAXOBS; i++) {
724 if ((g1 = gfmeasL1L5(obs + i, nav)) == 0.0) {
725 continue;
726 }
727 g0 = rtk->ssat[obs[i].sat - 1].gf[1];
728 rtk->ssat[obs[i].sat - 1].gf[1] = g1;
729 if (g0 != 0.0 && fabs(g1 - g0) > rtk->opt.thresslip) {
730 trace(NULL, 2, "detslip_gf_L1_L5: slip detected sat=%2d gf=%8.3f->%8.3f\n", obs[i].sat, g0, g1);
731 rtk->ssat[obs[i].sat - 1].slip[0] |= 1;
732 rtk->ssat[obs[i].sat - 1].slip[2] |= 1;
733 }
734 }
735}
736static void udbias_ppp(rtk_t* rtk, double tt, const obsd_t* obs, int n, const nav_t* nav) {
737 double offset, lami, cp, pr;
738 double* bias;
739 int i, j, f, slip, nf = NF(&rtk->opt), reset, isat;
740
741 trace(NULL, 3, "udbias : tt=%.1f n=%d\n", tt, n);
742
743 for (i = 0; i < MAXSAT; i++) {
744 for (j = 0; j < rtk->opt.nf; j++) {
745 rtk->ssat[i].slip[j] = 0;
746 rtk->ssat[i].pbreset[j] = 0;
747 }
748 }
749
750 /* detect cycle slip by LLI */
751 detslp_ll(rtk, obs, n);
752
753 /* detect cycle slip by geometry-free phase jump */
754 detslp_gf(rtk, obs, n, nav);
755
756 for (f = 0; f < nf; f++) {
757 /* reset phase-bias if expire obs outage counter */
758 for (i = 0; i < MAXSAT; i++) {
759 reset = ++rtk->ssat[i].outc[f] > (unsigned int)rtk->opt.maxout;
760
761 if (rtk->opt.modear == ARMODE_INST && rtk->x[IB(i + 1, f, &rtk->opt)] != 0.0) {
762 initx(rtk, 0.0, 0.0, IB(i + 1, f, &rtk->opt));
763 } else if (reset && rtk->x[IB(i + 1, f, &rtk->opt)] != 0.0) {
764 /* reset udbias */
765 initx(rtk, 0.0, 0.0, IB(i + 1, f, &rtk->opt));
766 trace(NULL, 3, "udbias : obs outage counter overflow (sat=%3d L%d n=%5d)\n", i + 1, f + 1,
767 rtk->ssat[i].outc[f]);
768 }
769 if (rtk->opt.modear != ARMODE_INST && reset) {
770 rtk->ssat[i].lock[f] = -rtk->opt.minlock;
771 }
772 }
773
774 /* reset phase-bias if detecting cycle slip */
775 for (i = 0; i < n && i < MAXOBS; i++) {
776 j = IB(obs[i].sat, f, &rtk->opt);
777 rtk->P[j + j * rtk->nx] += rtk->opt.prn[0] * rtk->opt.prn[0] * tt;
778 slip = rtk->ssat[obs[i].sat - 1].slip[f];
779 if (rtk->opt.ionoopt == IONOOPT_IFLC) {
780 slip |= rtk->ssat[obs[i].sat - 1].slip[1];
781 }
782 if (rtk->opt.modear == ARMODE_INST || !(slip & 1)) {
783 continue;
784 }
785 if (!(slip & 1)) {
786 continue;
787 }
788 rtk->x[j] = 0.0;
789 rtk->ssat[obs[i].sat - 1].lock[f] = -rtk->opt.minlock;
790 }
791 bias = zeros(n, 1);
792
793 /* estimate approximate phase-bias by phase - code */
794 for (i = j = 0, offset = 0.0; i < n; i++) {
795 cp = obs[i].L[f];
796 pr = obs[i].P[f];
797 lami = sat_lambda(obs[i].sat, f, nav);
798 if (cp == 0.0 || pr == 0.0 || lami == 0.0) {
799 continue;
800 }
801
802 bias[i] = cp * lami - pr;
803
804 /* offset = sum of (bias - phase-bias) for all valid sats in meters */
805 if (rtk->x[IB(obs[i].sat, f, &rtk->opt)] != 0.0) {
806 lami = sat_lambda(obs[i].sat, f, nav);
807 offset += bias[i] - rtk->x[IB(obs[i].sat, f, &rtk->opt)] * lami;
808 j++;
809 }
810 }
811
812 rtk->com_bias = j > 0 ? offset / j : 0; /* save offset for initialization below */
813 /* set initial states of phase-bias */
814 for (i = 0; i < n; i++) {
815 isat = obs[i].sat;
816 if (bias[i] == 0.0 || rtk->x[IB(isat, f, &rtk->opt)] != 0.0) {
817 continue;
818 }
819 trace(NULL, 3, "set initial states of phase-bias, sat=%2d f=%2d\n", obs[i].sat, f + 1);
820 lami = sat_lambda(isat, f, nav);
821 initx(rtk, (bias[i] - rtk->com_bias) / lami, SQR(rtk->opt.std[0]), IB(obs[i].sat, f, &rtk->opt));
822 rtk->ssat[isat - 1].lock[f] = -rtk->opt.minlock;
823 rtk->ssat[isat - 1].pbreset[f] = 1;
824 }
825
826 free(bias);
827 }
828}
829/* temporal update of states --------------------------------------------------*/
830static void udstate(rtk_t* rtk, const obsd_t* obs, int nu, const nav_t* nav) {
831 double tt = fabs(rtk->tt), bl = 0.0, dr[3];
832
833 trace(NULL, 3, "udstate : nu=%d\n", nu);
834
835 /* temporal update of position/velocity/acceleration */
836 udpos(rtk, obs, tt);
837
838 /* temporal update of tropospheric parameters */
839 if (rtk->opt.tropopt >= TROPOPT_EST) {
840 udtrop(rtk, tt, bl);
841 }
842 /* temporal update of phase-bias */
843 if (rtk->opt.mode > PMODE_DGPS) {
844 udbias_ppp(rtk, tt, obs, nu, nav);
845 }
846 /* temporal update of ionospheric parameters */
847 if (rtk->opt.ionoopt >= IONOOPT_EST) {
848 bl = baseline(rtk->x, rtk->rb, dr);
849 udion(rtk, tt, bl, obs, nu);
850 }
851 /* temporal update of eceiver h/w bias */
852 if (rtk->opt.glomodear == 2 && (rtk->opt.navsys & SYS_GLO)) {
853 udrcvbias(rtk, tt);
854 }
855}
856/* undifferenced phase/code residual for satellite ---------------------------*/
857static void zdres_sat(int base, double r, const obsd_t* obs, const nav_t* nav, const double* azel, const double* dant,
858 const prcopt_t* opt, double* y) {
859 double lam_arr[NFREQ];
860 int i, nf = NF(opt), qj;
861 for (i = 0; i < NFREQ; i++) {
862 lam_arr[i] = sat_lambda(obs->sat, i, nav);
863 }
864
865 qj = selfreqpair(obs->sat, opt, obs);
866
867 /* fallback: if preferred secondary freq has no data, try index 1
868 * (e.g. Galileo E5b selected but only E5a available from cssr2rtcm3) */
869 if (qj > 1 && obs->L[qj] == 0.0 && obs->P[qj] == 0.0 && (obs->L[1] != 0.0 || obs->P[1] != 0.0)) {
870 qj = 1;
871 }
872
873 for (i = 0; i < nf; i++) {
874 if (lam_arr[i] == 0.0) {
875 continue;
876 }
877
878 /* check snr mask */
879 if (testsnr(base, i, azel[1], obs->SNR[i] * 0.25, &opt->snrmask)) {
880 continue;
881 }
882 if (i > 0 && !(i & qj)) {
883 continue;
884 }
885 /* residuals = observable - pseudorange */
886 if (obs->L[i] != 0.0) {
887 y[i] = obs->L[i] * lam_arr[i] - r - dant[i];
888 }
889 if (obs->P[i] != 0.0) {
890 y[i + nf] = obs->P[i] - r - dant[i];
891 }
892 }
893}
894/* undifferenced phase/code residuals ----------------------------------------*/
895static int zdres(rtk_t* rtk, int base, const obsd_t* obs, int n, const double* rs, const double* dts, const int* svh,
896 const nav_t* nav, const double* rr, const prcopt_t* opt, int index, double* y, double* e,
897 double* azel) {
898 double r, rr_[3], pos[3], dant[NFREQPCV] = {0}, disp[3];
899 double zhd, zazel[] = {0.0, 90.0 * D2R}, ydif = 0.0;
900 int i, j, nf = NF(opt);
901
902 trace(NULL, 3, "zdres : n=%d\n", n);
903
904 for (i = 0; i < n * nf * 2; i++) {
905 y[i] = 0.0;
906 }
907
908 if (norm(rr, 3) <= 0.0) {
909 return 0; /* no receiver position */
910 }
911
912 for (i = 0; i < 3; i++) {
913 rr_[i] = rr[i];
914 }
915
916 /* earth tide correction */
917 if (opt->tidecorr) {
918 tidedisp(gpst2utc(obs[0].time), rr_, opt->tidecorr, &nav->erp, opt->odisp[base], disp);
919 for (i = 0; i < 3; i++) {
920 rr_[i] += disp[i];
921 }
922 }
923 ecef2pos(rr_, pos);
924
925 for (i = 0; i < n; i++) {
926 /* compute geometric-range and azimuth/elevation angle */
927 if ((r = geodist(rs + i * 6, rr_, e + i * 3)) <= 0.0) {
928 continue;
929 }
930 if (satazel(pos, e + i * 3, azel + i * 2) < opt->elmin) {
931 continue;
932 }
933
934 /* excluded satellite? */
935 if (satexclude(obs[i].sat, 0.0, svh[i], opt)) {
936 continue;
937 }
938
939 /* satellite clock-bias */
940 r += -CLIGHT * dts[i * 2];
941
942 /* troposphere delay model (hydrostatic) */
943 zhd = tropmodel(obs[0].time, pos, zazel, 0.0);
944 r += tropmapf(obs[i].time, pos, azel + i * 2, NULL) * zhd;
945
946 /* receiver antenna phase center correction */
947 antmodel(opt->pcvr + index, opt->antdel[index], azel + i * 2, opt->posopt[1], dant);
948
949 /* undifferenced phase/code residual for satellite */
950 zdres_sat(base, r, obs + i, nav, azel + i * 2, dant, opt, y + i * nf * 2);
951
952 /* 1/4 cycle phase shift correction in case of L2C */
953 if (opt->phasshft == ISBOPT_TABLE && isL2C(obs[i].code[1])) {
954 for (j = 0; j < nav->sfts.n; j++) {
955 if (strcmp(nav->sfts.rectyp[j], opt->rectype[0]) == 0) {
956 ydif = nav->sfts.bias[j];
957 }
958 }
959 y[1 + i * nf * 2] += ydif * sat_lambda(obs[i].sat, 1, nav);
960 }
961 }
962 trace(NULL, 4, "rr_=%.3f %.3f %.3f\n", rr_[0], rr_[1], rr_[2]);
963 trace(NULL, 4, "pos=%.9f %.9f %.3f\n", pos[0] * R2D, pos[1] * R2D, pos[2]);
964 for (i = 0; i < n; i++) {
965 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],
966 rs[2 + i * 6], dts[i * 2], azel[i * 2] * R2D, azel[1 + i * 2] * R2D);
967 }
968 trace(NULL, 4, "y=\n");
969 tracemat(NULL, 4, y, nf * 2, n, 13, 3);
970
971 return 1;
972}
973/* test valid observation data -----------------------------------------------*/
974static int validobs(int* iu, int* ir, int jj, int f, int nf, double* y) {
975 int i, j;
976 if (jj < 0) {
977 return 0;
978 }
979 i = iu[jj];
980 j = ir[jj];
981 if (i < 0 || j < 0) {
982 return 0;
983 }
984 /* if no phase observable, psudorange is also unusable */
985 return y[f + i * nf * 2] != 0.0 && y[f + j * nf * 2] != 0.0 &&
986 (f < nf || (y[f - nf + i * nf * 2] != 0.0 && y[f - nf + j * nf * 2] != 0.0));
987}
988/* double-differenced measurement error covariance ---------------------------*/
989static void ddcov(const int* nb, int n, const double* Ri, const double* Rj, int nv, double* R) {
990 int i, j, k = 0, b;
991
992 trace(NULL, 3, "ddcov : n=%d\n", n);
993
994 for (i = 0; i < nv * nv; i++) {
995 R[i] = 0.0;
996 }
997 for (b = 0; b < n; k += nb[b++]) {
998 for (i = 0; i < nb[b]; i++) {
999 for (j = 0; j < nb[b]; j++) {
1000 R[k + i + (k + j) * nv] = Ri[k + i] + (i == j ? Rj[k + i] : 0.0);
1001 }
1002 }
1003 }
1004 trace(NULL, 5, "R=\n");
1005 tracemat(NULL, 5, R, nv, nv, 8, 6);
1006}
1007/* precise tropspheric model -------------------------------------------------*/
1008static double prectrop(gtime_t time, const double* pos, int r, const double* azel, const prcopt_t* opt, const double* x,
1009 double* dtdx) {
1010 double m_w = 0.0, cotz, grad_n, grad_e;
1011 int i = IT(r, opt);
1012
1013 /* wet mapping function */
1014 tropmapf(time, pos, azel, &m_w);
1015
1016 if (opt->tropopt >= TROPOPT_ESTG && azel[1] > 0.0) {
1017 /* m_w=m_0+m_0*cot(el)*(Gn*cos(az)+Ge*sin(az)): ref [6] */
1018 cotz = 1.0 / tan(azel[1]);
1019 grad_n = m_w * cotz * cos(azel[0]);
1020 grad_e = m_w * cotz * sin(azel[0]);
1021 m_w += grad_n * x[i + 1] + grad_e * x[i + 2];
1022 dtdx[1] = grad_n * x[i];
1023 dtdx[2] = grad_e * x[i];
1024 } else {
1025 dtdx[1] = dtdx[2] = 0.0;
1026 }
1027 dtdx[0] = m_w;
1028 return m_w * x[i];
1029}
1030/* glonass inter-channel bias correction -------------------------------------*/
1031static double gloicbcorr(int sat1, int sat2, const prcopt_t* opt, double lam1, double lam2, int f) {
1032 double dfreq;
1033
1034 if (f >= NFREQGLO || f >= opt->nf || !opt->exterr.ena[2]) {
1035 return 0.0;
1036 }
1037
1038 dfreq = (CLIGHT / lam1 - CLIGHT / lam2) / (f == 0 ? DFRQ1_GLO : DFRQ2_GLO);
1039
1040 return opt->exterr.gloicb[f] * 0.01 * dfreq; /* (m) */
1041}
1042/* test navi system (m=0:gps/qzs/sbs,1:glo,2:gal,3:bds,4:qzs) ----------------*/
1043static int test_sys(int sys, int m, int qzsmodear, int code) {
1044 switch (sys) {
1045 case SYS_GPS:
1046 return m == (code == CODE_L2X ? 5 : 0);
1047 case SYS_QZS:
1048 return m == (qzsmodear == 2 ? 0 : 4);
1049 case SYS_SBS:
1050 return m == 0;
1051 case SYS_GLO:
1052 return m == 1;
1053 case SYS_GAL:
1054 return m == 2;
1055 case SYS_CMP:
1056 return m == 3;
1057 }
1058 return 0;
1059}
1060
1061static int check_same_fac(const obsd_t* obs, int ir[][MAXSAT]) {
1062 return (obs[ir[0][0]].facility == obs[ir[1][0]].facility);
1063}
1064
1065static int get_iriu_index(int sati, int ns, const int* sat) {
1066 int i;
1067 for (i = 0; i < ns; i++) {
1068 if (sati == sat[i]) {
1069 return i;
1070 }
1071 }
1072 return -1;
1073}
1074
1075static int valid_zd_sig(rtk_t* rtk, const obsd_t* obs, const nav_t* nav, double* y, int nf, int f, int j, int m, int ch,
1076 int sati, int iu[][MAXSAT], int ir[][MAXSAT], int sat[][MAXSAT], int* ff, int* ns, int* satj,
1077 int* jj, int* sysi, int* sysj, double* lami, double* lamj) {
1078 prcopt_t* opt = &rtk->opt;
1079 int code;
1080
1081 *satj = obs[j].sat;
1082 if ((*jj = get_iriu_index(*satj, ns[ch], sat[ch])) == -1) {
1083 return 0;
1084 }
1085 if (sati == *satj) {
1086 return 0;
1087 }
1088 *sysi = rtk->ssat[sati - 1].sys;
1089 *sysj = rtk->ssat[*satj - 1].sys;
1090 code = rtk->ssat[*satj - 1].code[f % nf];
1091 if (!test_sys(*sysj, m, opt->qzsmodear, code)) {
1092 return 0;
1093 }
1094 if (!validobs(iu[ch], ir[ch], *jj, f, nf, y)) {
1095 return 0;
1096 }
1097
1098 *ff = f % nf;
1099 *lami = sat_lambda(sati, *ff, nav);
1100 *lamj = sat_lambda(*satj, *ff, nav);
1101 if (*lami <= 0.0 || *lamj <= 0.0) {
1102 return 0;
1103 }
1104
1105 return 1;
1106}
1107
1108/* set_init_pb() is created by copying the second half of the udbias_ppp() */
1109static void set_init_pb(rtk_t* rtk, const nav_t* nav, const obsd_t* obs, int n, int f, int satj) {
1110 int i, j, isat;
1111 double offset, lami, cp, pr;
1112 double* bias;
1113 bias = zeros(n, 1);
1114
1115 /* estimate approximate phase-bias by phase - code */
1116 for (i = j = 0, offset = 0.0; i < n; i++) {
1117 cp = obs[i].L[f];
1118 pr = obs[i].P[f];
1119 lami = sat_lambda(obs[i].sat, f, nav);
1120 if (cp == 0.0 || pr == 0.0 || lami == 0.0) {
1121 continue;
1122 }
1123
1124 bias[i] = cp * lami - pr;
1125
1126 /* offset = sum of (bias - phase-bias) for all valid sats in meters */
1127 if (rtk->x[IB(obs[i].sat, f, &rtk->opt)] != 0.0) {
1128 lami = sat_lambda(obs[i].sat, f, nav);
1129 offset += bias[i] - rtk->x[IB(obs[i].sat, f, &rtk->opt)] * lami;
1130 j++;
1131 }
1132 }
1133 rtk->com_bias = j > 0 ? offset / j : 0; /* save offset for initialization below */
1134 /* set initial states of phase-bias */
1135 for (i = 0; i < n; i++) {
1136 isat = obs[i].sat;
1137 if (isat != satj) {
1138 continue;
1139 }
1140 if (bias[i] == 0.0) {
1141 continue;
1142 }
1143 trace(NULL, 2, "set initial states of phase-bias(set_init_pb), sat=%2d f=%2d\n", obs[i].sat, f + 1);
1144 lami = sat_lambda(isat, f, nav);
1145 initx(rtk, (bias[i] - rtk->com_bias) / lami, SQR(rtk->opt.std[0]), IB(obs[i].sat, f, &rtk->opt));
1146 rtk->ssat[isat - 1].lock[f] = -rtk->opt.minlock;
1147 rtk->ssat[isat - 1].pbreset[f] = 1;
1148 }
1149 free(bias);
1150}
1151/* double-differenced phase/code residuals -----------------------------------*/
1152static int ddres(rtk_t* rtk, const nav_t* nav, double* dt, const double* x, const double* P, const obsd_t* obs, int n,
1153 int sat[][MAXSAT], double* y, double* e, double* azel, int iu[][MAXSAT], int ir[][MAXSAT], int* ns,
1154 double* v, double* H, double* R, int* vflg, int niter, int order) {
1155 prcopt_t* opt = &rtk->opt;
1156 double bl, dr[3], posu[3], posr[3], didxi = 0.0, didxj = 0.0, *im[SSR_CH_NUM];
1157 double *tropr[SSR_CH_NUM], *tropu[SSR_CH_NUM], *dtdxr, *dtdxu, *Ri, *Rj, lami, lamj, fi, fj, df, *Hi = NULL;
1158 int h, i, j, k, m, f, ff, nv = 0, nb[NFREQ * 6 * 2 + 1] = {0}, b = 0, sysi, sysj, nf = NF(opt);
1159 static int refsat[NFREQ * 2 * MAXREF * 6 + 1] = {0};
1160 static int refsat2[NFREQ * 2 * 6 + 1] = {0};
1161 int flg = 0, code, ch, sch, pch = 0, samefac = 0, refchgflg, ns_, sati, satj, refsatch, rch;
1162 int v0, v1, j0, j1, jj, vidx;
1163 double azeli = 0.0;
1164 static int use_ch[NFREQ][MAXSAT], pre_use_ch[NFREQ][MAXSAT];
1165 double tow = time2gpst(obs[0].time, NULL);
1166
1167 ns_ = opt->l6mrg ? ns[0] + ns[1] : ns[0];
1168 trace(NULL, 2, "ddres : tow=%.1f dt[0]=%.1f dt[1]=%.1f nx=%d ns=%d\n", tow, dt[0], dt[1], rtk->nx, ns_);
1169
1170 bl = baseline(x, rtk->rb, dr);
1171 ecef2pos(x, posu);
1172 ecef2pos(rtk->rb, posr);
1173 memcpy(pre_use_ch, use_ch, sizeof(use_ch));
1174 for (f = 0; f < NFREQ; f++) {
1175 for (i = 0; i < MAXSAT; i++) {
1176 use_ch[f][i] = -1;
1177 }
1178 }
1179
1180 Ri = mat(ns_ * nf * 2 + 2, 1);
1181 Rj = mat(ns_ * nf * 2 + 2, 1);
1182 dtdxu = mat(ns_, 3);
1183 dtdxr = mat(ns_, 3);
1184
1185 for (ch = 0; ch < (opt->l6mrg ? SSR_CH_NUM : 1); ch++) {
1186 tropu[ch] = mat(ns[ch], 1);
1187 tropr[ch] = mat(ns[ch], 1);
1188 im[ch] = mat(ns[ch], 1);
1189 }
1190
1191 for (i = 0; i < MAXSAT; i++) {
1192 for (j = 0; j < NFREQ; j++) {
1193 rtk->ssat[i].resp[j] = rtk->ssat[i].resc[j] = 0.0;
1194 }
1195 }
1196 /* compute factors of ionospheric and tropospheric delay */
1197 for (ch = 0; ch < (opt->l6mrg ? SSR_CH_NUM : 1); ch++) {
1198 for (i = 0; i < ns[ch]; i++) {
1199 if (opt->ionoopt >= IONOOPT_EST) {
1200 im[ch][i] = (ionmapf(posu, azel + iu[ch][i] * 2) + ionmapf(posr, azel + ir[ch][i] * 2)) / 2.0;
1201 }
1202 if (opt->tropopt >= TROPOPT_EST) {
1203 tropu[ch][i] = prectrop(rtk->sol.time, posu, 0, azel + iu[ch][i] * 2, opt, x, dtdxu + i * 3);
1204 tropr[ch][i] = prectrop(rtk->sol.time, posr, 1, azel + ir[ch][i] * 2, opt, x, dtdxr + i * 3);
1205 }
1206 }
1207 }
1208
1209 /* check if the same facility is used */
1210 if (opt->l6mrg) {
1211 samefac = check_same_fac(obs, ir);
1212 }
1213
1214 for (m = 0; m < 6; m++) { /* m=0:gps/qzs/sbs,1:glo,2:gal,3:bds,4:qzs */
1215
1216 for (f = opt->mode > PMODE_DGPS ? 0 : nf; f < nf * 2; f++) {
1217 /* search reference satellite with highest elevation */
1218 for (h = 0, sch = -1, refchgflg = 0, refsatch = -1; h < ((opt->l6mrg && !samefac) ? SSR_CH_NUM : 1); h++) {
1219 azeli = 0.0;
1220 j0 = -1;
1221 j1 = -1;
1222 for (i = -1, j = 0, flg = 0; j < n; j++) {
1223 sati = obs[j].sat;
1224 sysi = rtk->ssat[sati - 1].sys;
1225 code = rtk->ssat[obs[j].sat - 1].code[f % nf];
1226 if (!test_sys(sysi, m, opt->qzsmodear, code)) {
1227 continue;
1228 }
1229 if (niter > 0) {
1230 for (k = 0; k < niter; k++) {
1231 if (refsat[NFREQ * 2 * MAXREF * m + NFREQ * 2 * k + f] == sati) {
1232 flg = 1;
1233 }
1234 }
1235 if (flg == 1) {
1236 continue;
1237 }
1238 }
1239 j0 = get_iriu_index(sati, ns[0], sat[0]);
1240 j1 = get_iriu_index(sati, ns[1], sat[1]); /* only valid for 2ch */
1241 if (j0 == -1 && j1 == -1) {
1242 continue;
1243 }
1244
1245 v0 = 0;
1246 v1 = 0;
1247 if (samefac || h == 1) {
1248 /* refsat is selected from one side ch */
1249 v0 = validobs(iu[0], ir[0], j0, f, nf, y);
1250 v1 = validobs(iu[1], ir[1], j1, f, nf, y);
1251 if (!v0 && !v1) {
1252 continue;
1253 }
1254 trace(NULL, 5, "h:%d sati:%d v0:%d v1:%d ir[0][j0]:%d\n", h, sati, v0, v1,
1255 ir[0][j0 == -1 ? j1 : j0]);
1256 ch = (!v0) ? 1 : 0;
1257 } else { /* h=0 || samefac==0 */
1258 /* select common refsat between 2ch */
1259 if (!(v0 = validobs(iu[0], ir[0], j0, f, nf, y))) {
1260 continue;
1261 }
1262 if (opt->l6mrg && !(v1 = validobs(iu[1], ir[1], j1, f, nf, y))) {
1263 continue;
1264 }
1265 ch = (!v0) ? 1 : 0;
1266 }
1267 if (opt->qzsmodear == 2 && opt->posopt[8] && sysi == SYS_QZS) {
1268 continue;
1269 }
1270
1271 vidx = v0 ? j0 : j1;
1272 if (i < 0 || azel[1 + iu[ch][vidx] * 2] >= azeli) {
1273 i = j;
1274 azeli = azel[1 + iu[ch][vidx] * 2];
1275 if (samefac) {
1276 refsatch = ch;
1277 }
1278 if (h == 1) {
1279 sch = ch;
1280 }
1281 trace(NULL, 5, "refsatroop m:%d i:%d sat:%d azel:%f\n", m, i, sati, azeli);
1282 }
1283 }
1284 if (i >= 0) {
1285 break;
1286 }
1287 }
1288 if (i < 0) {
1289 continue;
1290 }
1291 sati = obs[i].sat;
1292 if (f < nf) {
1293 use_ch[f][sati] = 2; /* use_ch=2 means ref sat */
1294 }
1295 if (f < nf && refsat2[NFREQ * 2 * m + NFREQ * 2 + f] != sati &&
1296 refsat2[NFREQ * 2 * m + NFREQ * 2 + f] != 0) {
1297 trace(NULL, 2, "refsat changed2 %s L%d sat=%2d -> %2d\n", time_str(rtk->sol.time, 0), f % nf + 1,
1298 refsat2[NFREQ * 2 * m + NFREQ * 2 + f], sati);
1299 refchgflg = 1;
1300 }
1301 refsat[NFREQ * 2 * MAXREF * m + NFREQ * 2 * niter + f] = sati;
1302 refsat2[NFREQ * 2 * m + NFREQ * 2 + f] = sati;
1303 if (niter > 0) {
1304 trace(NULL, 2, "refsat changed %s L%d sat=%2d -> %2d\n", time_str(rtk->sol.time, 0), f % nf + 1,
1305 refsat[NFREQ * MAXREF * 2 * m + NFREQ * 2 * (niter - 1) + f],
1306 refsat[NFREQ * MAXREF * 2 * m + NFREQ * 2 * niter + f]);
1307 }
1308 if (f < nf) {
1309 trace(NULL, 4, "refsat m:%d L%d %s: %d niter:%d refsatch:%d\n", m, f + 1, time_str(rtk->sol.time, 0),
1310 sati, niter, refsatch);
1311 }
1312
1313 /* make double difference */
1314 if (opt->l6mrg) {
1315 pch = (ns[0] >= ns[1]) ? 0 : 1;
1316 }
1317 for (j = 0; j < n; j++) {
1318 if (samefac) {
1319 rch = refsatch;
1320 } else {
1321 rch = (sch != -1) ? sch : pch;
1322 }
1323 i = get_iriu_index(sati, ns[rch], sat[rch]);
1324 if (!opt->l6mrg) {
1325 ch = 0;
1326 if (!valid_zd_sig(rtk, obs, nav, y, nf, f, j, m, ch, sati, iu, ir, sat, &ff, ns, &satj, &jj, &sysi,
1327 &sysj, &lami, &lamj)) {
1328 continue;
1329 }
1330 } else if (sch == -1) {
1331 /* priority channel selection */
1332 ch = pch;
1333 if (!valid_zd_sig(rtk, obs, nav, y, nf, f, j, m, ch, sati, iu, ir, sat, &ff, ns, &satj, &jj, &sysi,
1334 &sysj, &lami, &lamj)) {
1335 ch = (ch == 1) ? 0 : 1;
1336 if (!samefac) {
1337 i = get_iriu_index(sati, ns[ch], sat[ch]);
1338 rch = ch;
1339 }
1340 if (!valid_zd_sig(rtk, obs, nav, y, nf, f, j, m, ch, sati, iu, ir, sat, &ff, ns, &satj, &jj,
1341 &sysi, &sysj, &lami, &lamj)) {
1342 continue;
1343 }
1344 }
1345 } else {
1346 trace(NULL, 3, "even though it is a clas 2ch case, refsat is selected from one side ch.(ch=%d).\n",
1347 sch);
1348 ch = sch;
1349 if (!valid_zd_sig(rtk, obs, nav, y, nf, f, j, m, ch, sati, iu, ir, sat, &ff, ns, &satj, &jj, &sysi,
1350 &sysj, &lami, &lamj)) {
1351 continue;
1352 }
1353 }
1354 if (f < nf) {
1355 use_ch[f][satj] = ch;
1356 }
1357 if (order == FST && f < nf && opt->l6mrg && !samefac) {
1358 if ((pre_use_ch[f][sati] != ch && refchgflg) || (pre_use_ch[f][satj] != ch)) {
1359 set_init_pb(rtk, nav, obs, n, f, satj);
1360 }
1361 }
1362
1363 if (H) {
1364 Hi = H + nv * rtk->nx;
1365 }
1366
1367 /* double-differenced residual */
1368 v[nv] = (y[f + iu[rch][i] * nf * 2] - y[f + ir[rch][i] * nf * 2]) -
1369 (y[f + iu[ch][jj] * nf * 2] - y[f + ir[ch][jj] * nf * 2]);
1370
1371 if (order == FST) {
1372 trace(NULL, 5, "dd_res sati:%d satj:%d rch:%d ch:%d v:%f\n", sati, satj, rch, ch, v[nv]);
1373 }
1374
1375 /* partial derivatives by rover position */
1376 if (H) {
1377 for (k = 0; k < 3; k++) {
1378 Hi[k] = -e[k + iu[rch][i] * 3] + e[k + iu[ch][jj] * 3];
1379 }
1380 }
1381 /* double-differenced ionospheric delay term */
1382 if (opt->ionoopt == IONOOPT_EST || opt->ionoopt == IONOOPT_EST_ADPT) {
1383 fi = lami / lam_carr[0];
1384 fj = lamj / lam_carr[0];
1385 didxi = (f < nf ? -1.0 : 1.0) * fi * fi * im[rch][i];
1386 didxj = (f < nf ? -1.0 : 1.0) * fj * fj * im[ch][jj];
1387 v[nv] -= didxi * x[II(sati, opt)] - didxj * x[II(satj, opt)];
1388 if (H) {
1389 Hi[II(sati, opt)] = didxi;
1390 Hi[II(satj, opt)] = -didxj;
1391 }
1392 }
1393 /* double-differenced tropospheric delay term */
1394 if (opt->tropopt == TROPOPT_EST || opt->tropopt == TROPOPT_ESTG) {
1395 v[nv] -= (tropu[rch][i] - tropu[ch][jj]) - (tropr[rch][i] - tropr[ch][jj]);
1396 for (k = 0; k < (opt->tropopt < TROPOPT_ESTG ? 1 : 3); k++) {
1397 if (!H) {
1398 continue;
1399 }
1400 Hi[IT(0, opt) + k] = (dtdxu[k + i * 3] - dtdxu[k + jj * 3]);
1401 Hi[IT(1, opt) + k] = -(dtdxr[k + i * 3] - dtdxr[k + jj * 3]);
1402 }
1403 }
1404 /* double-differenced phase-bias term */
1405 if (f < nf) {
1406 if (opt->ionoopt != IONOOPT_IFLC) {
1407 v[nv] -= lami * x[IB(sati, f, opt)] - lamj * x[IB(satj, f, opt)];
1408 if (H) {
1409 Hi[IB(sati, f, opt)] = lami;
1410 Hi[IB(satj, f, opt)] = -lamj;
1411 }
1412 } else {
1413 v[nv] -= x[IB(sati, f, opt)] - x[IB(satj, f, opt)];
1414 if (H) {
1415 Hi[IB(sati, f, opt)] = 1.0;
1416 Hi[IB(satj, f, opt)] = -1.0;
1417 }
1418 }
1419 }
1420 /* glonass receiver h/w bias term */
1421 if (rtk->opt.glomodear == 2 && sysi == SYS_GLO && sysj == SYS_GLO && ff < NFREQGLO) {
1422 df = (CLIGHT / lami - CLIGHT / lamj) / 1E6; /* freq-difference (MHz) */
1423 v[nv] -= df * x[IL(ff, opt)];
1424 if (H) {
1425 Hi[IL(ff, opt)] = df;
1426 }
1427 }
1428 /* glonass interchannel bias correction */
1429 else if (sysi == SYS_GLO && sysj == SYS_GLO) {
1430 v[nv] -= gloicbcorr(sati, satj, &rtk->opt, lami, lamj, f);
1431 }
1432
1433 if (f < nf) {
1434 rtk->ssat[satj - 1].resc[f] = v[nv];
1435 } else {
1436 rtk->ssat[satj - 1].resp[f - nf] = v[nv];
1437 }
1438
1439 /* single-differenced measurement error variances */
1440 Ri[nv] = varerr(sati, sysi, azel[1 + iu[rch][i] * 2], bl, dt[rch], f, opt);
1441 Rj[nv] = varerr(satj, sysj, azel[1 + iu[ch][jj] * 2], bl, dt[ch], f, opt);
1442 if (f == 1) {
1443 Ri[nv] *= SQR(2.55 / 1.55);
1444 Rj[nv] *= SQR(2.55 / 1.55);
1445 }
1446 /* set valid data flags */
1447 if (f < nf) {
1448 rtk->ssat[sati - 1].vsat[f] = rtk->ssat[satj - 1].vsat[f] = 1;
1449 }
1450
1451 trace(NULL, 4, "sat=%3d-%3d %s%d v=%13.3f R=%8.6f %8.6f\n", sati, satj, f < nf ? "L" : "P", f % nf + 1,
1452 v[nv], Ri[nv], Rj[nv]);
1453
1454 vflg[nv++] = (sati << 16) | (satj << 8) | ((f < nf ? 0 : 1) << 4) | (f % nf);
1455 nb[b]++;
1456 }
1457 b++;
1458 }
1459 }
1460 /* end of system loop */
1461
1462 /* double-differenced measurement error covariance */
1463 ddcov(nb, b, Ri, Rj, nv, R);
1464
1465 free(Ri);
1466 free(Rj);
1467 free(dtdxu);
1468 free(dtdxr);
1469 for (ch = 0; ch < (opt->l6mrg ? SSR_CH_NUM : 1); ch++) {
1470 free(tropu[ch]);
1471 free(tropr[ch]);
1472 free(im[ch]);
1473 }
1474
1475 return nv;
1476}
1477/* time-interpolation of residuals (for post-mission) ------------------------*/
1478static double intpres(gtime_t time, const obsd_t* obs, int n, const nav_t* nav, rtk_t* rtk, int ch, double* y) {
1479 static obsd_t obsb[SSR_CH_NUM][MAXOBS];
1480 static double yb[SSR_CH_NUM][MAXOBS * NFREQ * 2], rs[SSR_CH_NUM][MAXOBS * 6], dts[SSR_CH_NUM][MAXOBS * 2],
1481 var[SSR_CH_NUM][MAXOBS];
1482 static double e[SSR_CH_NUM][MAXOBS * 3], azel[SSR_CH_NUM][MAXOBS * 2];
1483 static int nb[SSR_CH_NUM] = {0}, svh[SSR_CH_NUM][MAXOBS * 2];
1484 prcopt_t* opt = &rtk->opt;
1485 double tt = timediff(time, obs[0].time), ttb, *p, *q;
1486 int i, j, k, nf = NF(opt);
1487
1488 trace(NULL, 3, "intpres : n=%d tt=%.1f\n", n, tt);
1489
1490 if (nb[ch] == 0 || fabs(tt) < DTTOL) {
1491 nb[ch] = n;
1492 for (i = 0; i < n; i++) {
1493 obsb[ch][i] = obs[i];
1494 }
1495 return tt;
1496 }
1497 ttb = timediff(time, obsb[ch][0].time);
1498 if (fabs(ttb) > opt->maxtdiff * 2.0 || ttb == tt) {
1499 return tt;
1500 }
1501
1502 satposs(time, obsb[ch], nb[ch], nav, opt->sateph, rs[ch], dts[ch], var[ch], svh[ch]);
1503
1504 if (!zdres(rtk, 1, obsb[ch], nb[ch], rs[ch], dts[ch], svh[ch], nav, rtk->rb, opt, 1, yb[ch], e[ch], azel[ch])) {
1505 return tt;
1506 }
1507 for (i = 0; i < n; i++) {
1508 for (j = 0; j < nb[ch]; j++) {
1509 if (obsb[ch][j].sat == obs[i].sat) {
1510 break;
1511 }
1512 }
1513 if (j >= nb[ch]) {
1514 continue;
1515 }
1516 for (k = 0, p = y + i * nf * 2, q = yb[ch] + j * nf * 2; k < nf * 2; k++, p++, q++) {
1517 if (*p == 0.0 || *q == 0.0) {
1518 *p = 0.0;
1519 } else {
1520 *p = (ttb * (*p) - tt * (*q)) / (ttb - tt);
1521 }
1522 }
1523 }
1524 return fabs(ttb) > fabs(tt) ? ttb : tt;
1525}
1526/* single to double-difference transformation matrix (D') --------------------*/
1527static int ddmat(rtk_t* rtk, double* D) {
1528 int i, j, k, m, f, nb = 0, nx = rtk->nx, na = rtk->na, nf = NF(&rtk->opt), code;
1529
1530 trace(NULL, 3, "ddmat :\n");
1531
1532 for (i = 0; i < MAXSAT; i++) {
1533 for (j = 0; j < NFREQ; j++) {
1534 rtk->ssat[i].fix[j] = 0;
1535 }
1536 }
1537 for (i = 0; i < na; i++) {
1538 D[i + i * nx] = 1.0;
1539 }
1540
1541 for (m = 0; m < 6; m++) { /* m=0:gps/qzs/sbs,1:glo,2:gal,3:bds,4:qzs */
1542
1543 if (m == 1 && rtk->opt.glomodear == 0) {
1544 continue;
1545 }
1546 if (m == 3 && rtk->opt.bdsmodear == 0) {
1547 continue;
1548 }
1549 if (m == 4 && rtk->opt.qzsmodear != 1) {
1550 continue;
1551 }
1552
1553 for (f = 0, k = na; f < nf; f++, k += MAXSAT) {
1554 for (i = k; i < k + MAXSAT; i++) {
1555 code = rtk->ssat[i - k].code[f];
1556 if (rtk->x[i] == 0.0 || !test_sys(rtk->ssat[i - k].sys, m, rtk->opt.qzsmodear, code) ||
1557 !rtk->ssat[i - k].vsat[f]) {
1558 continue;
1559 }
1560 if (rtk->ssat[i - k].lock[f] > 0 && !(rtk->ssat[i - k].slip[f] & 2) &&
1561 rtk->ssat[i - k].azel[1] >= rtk->opt.elmaskar) {
1562 rtk->ssat[i - k].fix[f] = 2; /* fix */
1563 break;
1564 } else {
1565 rtk->ssat[i - k].fix[f] = 1;
1566 }
1567 }
1568 for (j = k; j < k + MAXSAT; j++) {
1569 code = rtk->ssat[j - k].code[f];
1570 if (i == j || rtk->x[j] == 0.0 || !test_sys(rtk->ssat[j - k].sys, m, rtk->opt.qzsmodear, code) ||
1571 !rtk->ssat[j - k].vsat[f]) {
1572 continue;
1573 }
1574 if (rtk->ssat[j - k].lock[f] > 0 && !(rtk->ssat[j - k].slip[f] & 2) &&
1575 rtk->ssat[j - k].azel[1] >= rtk->opt.elmaskar) {
1576 D[i + (na + nb) * nx] = 1.0;
1577 D[j + (na + nb) * nx] = -1.0;
1578 nb++;
1579 rtk->ssat[j - k].fix[f] = 2; /* fix */
1580 } else {
1581 rtk->ssat[j - k].fix[f] = 1;
1582 }
1583 }
1584 }
1585 }
1586 trace(NULL, 5, "D=\n");
1587 tracemat(NULL, 5, D, nx, na + nb, 2, 0);
1588 return nb;
1589}
1590/* restore single-differenced ambiguity --------------------------------------*/
1591static void restamb(rtk_t* rtk, const double* bias, int nb, double* xa) {
1592 int i, n, m, f, index[MAXSAT], nv = 0, nf = NF(&rtk->opt), code;
1593
1594 trace(NULL, 3, "restamb :\n");
1595
1596 for (i = 0; i < rtk->nx; i++) {
1597 xa[i] = rtk->x[i];
1598 }
1599 for (i = 0; i < rtk->na; i++) {
1600 xa[i] = rtk->xa[i];
1601 }
1602
1603 for (m = 0; m < 6; m++) {
1604 for (f = 0; f < nf; f++) {
1605 for (n = i = 0; i < MAXSAT; i++) {
1606 code = rtk->ssat[i].code[f];
1607 if (!test_sys(rtk->ssat[i].sys, m, rtk->opt.qzsmodear, code) || rtk->ssat[i].fix[f] != 2) {
1608 continue;
1609 }
1610 index[n++] = IB(i + 1, f, &rtk->opt);
1611 }
1612 if (n < 2) {
1613 continue;
1614 }
1615
1616 xa[index[0]] = rtk->x[index[0]];
1617
1618 for (i = 1; i < n; i++) {
1619 xa[index[i]] = xa[index[0]] - bias[nv++];
1620 }
1621 }
1622 }
1623}
1624/* hold integer ambiguity ----------------------------------------------------*/
1625static void holdamb(rtk_t* rtk, const double* xa) {
1626 double *v, *H, *R;
1627 int i, n, m, f, info, index[MAXSAT], nb = rtk->nx - rtk->na, nv = 0, nf = NF(&rtk->opt);
1628 int code;
1629
1630 trace(NULL, 3, "holdamb :\n");
1631
1632 v = mat(nb, 1);
1633 H = zeros(nb, rtk->nx);
1634
1635 for (m = 0; m < 6; m++) {
1636 for (f = 0; f < nf; f++) {
1637 for (n = i = 0; i < MAXSAT; i++) {
1638 code = rtk->ssat[i].code[f];
1639 if (!test_sys(rtk->ssat[i].sys, m, rtk->opt.qzsmodear, code) || rtk->ssat[i].fix[f] != 2 ||
1640 rtk->ssat[i].azel[1] < rtk->opt.elmaskhold || !rtk->ssat[i].vsat[f]) {
1641 continue;
1642 }
1643 index[n++] = IB(i + 1, f, &rtk->opt);
1644 rtk->ssat[i].fix[f] = 3; /* hold */
1645 }
1646 /* constraint to fixed ambiguity */
1647 for (i = 1; i < n; i++) {
1648 v[nv] = (xa[index[0]] - xa[index[i]]) - (rtk->x[index[0]] - rtk->x[index[i]]);
1649
1650 H[index[0] + nv * rtk->nx] = 1.0;
1651 H[index[i] + nv * rtk->nx] = -1.0;
1652 nv++;
1653 }
1654 }
1655 }
1656 if (nv > 0) {
1657 R = zeros(nv, nv);
1658 for (i = 0; i < nv; i++) {
1659 R[i + i * nv] = rtk->opt.varholdamb;
1660 }
1661
1662 /* update states with constraints */
1663 if ((info = filter(rtk->x, rtk->P, H, v, R, rtk->nx, nv))) {
1664 errmsg(rtk, "filter error (info=%d)\n", info);
1665 }
1666 free(R);
1667 }
1668 free(v);
1669 free(H);
1670}
1671/* resolve integer ambiguity by LAMBDA ---------------------------------------*/
1672static int resamb_LAMBDA(rtk_t* rtk, double* bias, double* xa) {
1673 prcopt_t* opt = &rtk->opt;
1674 int i, j, ny, nb, info, nx = rtk->nx, na = rtk->na;
1675 double *D, *DP, *y, *Qy, *b, *db, *Qb, *Qab, *QQ, s[2];
1676 double *Ds, *Pac, *Pcc, *D_P;
1677 int nc;
1678
1679 trace(NULL, 3, "resamb_LAMBDA : nx=%d\n", nx);
1680
1681 rtk->sol.ratio = 0.0;
1682
1683 if (rtk->opt.mode <= PMODE_DGPS || rtk->opt.modear == ARMODE_OFF || rtk->opt.thresar[0] < 1.0) {
1684 return 0;
1685 }
1686 /* single to double-difference transformation matrix (D') */
1687 D = zeros(nx, nx);
1688 if ((nb = ddmat(rtk, D)) <= 0) {
1689 errmsg(rtk, "no valid double-difference\n");
1690 free(D);
1691 return 0;
1692 }
1693 ny = na + nb;
1694 y = mat(ny, 1);
1695 Qy = mat(ny, ny);
1696 DP = mat(ny, nx);
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 nc = nx - na;
1704 Ds = mat(nc, nb);
1705 Pac = mat(na, nc);
1706 Pcc = mat(nc, nc);
1707 D_P = mat(nb, nc);
1708 for (i = 0; i < nc; i++) {
1709 for (j = 0; j < nb; j++) {
1710 Ds[i + j * nc] = D[na + i + (na + j) * nx];
1711 }
1712 }
1713 for (i = 0; i < na; i++) {
1714 for (j = 0; j < nc; j++) {
1715 Pac[i + j * na] = rtk->P[i + (na + j) * nx];
1716 }
1717 }
1718 for (i = 0; i < nc; i++) {
1719 for (j = 0; j < nc; j++) {
1720 Pcc[i + j * nc] = rtk->P[na + i + (na + j) * nx];
1721 }
1722 }
1723 matmul("TN", ny, 1, nx, 1.0, D, rtk->x, 0.0, y);
1724 matmul("NN", na, nb, nc, 1.0, Pac, Ds, 0.0, Qab);
1725 matmul("TN", nb, nc, nc, 1.0, Ds, Pcc, 0.0, D_P);
1726 matmul("NN", nb, nb, nc, 1.0, D_P, Ds, 0.0, Qb);
1727 free(Ds);
1728 free(Pac);
1729 free(Pcc);
1730 free(D_P);
1731
1732 trace(NULL, 4, "N(0)=");
1733 tracemat(NULL, 4, y + na, 1, nb, 10, 3);
1734
1735 /* lambda/mlambda integer least-square estimation */
1736 if (!(info = lambda(nb, 2, y + na, Qb, b, s))) {
1737 trace(NULL, 4, "N(1)=");
1738 tracemat(NULL, 4, b, 1, nb, 10, 3);
1739 trace(NULL, 4, "N(2)=");
1740 tracemat(NULL, 4, b + nb, 1, nb, 10, 3);
1741
1742 rtk->sol.ratio = s[0] > 0 ? (float)(s[1] / s[0]) : 0.0f;
1743 if (rtk->sol.ratio > 999.9) {
1744 rtk->sol.ratio = 999.9f;
1745 }
1746
1747 /* validation by popular ratio-test */
1748 if (s[0] <= 0.0 || s[1] / s[0] >= qfrtk[opt->alphaar][nb - 1]) {
1749 /* transform float to fixed solution (xa=xa-Qab*Qb\‍(b0-b)) */
1750 for (i = 0; i < na; i++) {
1751 rtk->xa[i] = rtk->x[i];
1752 for (j = 0; j < na; j++) {
1753 rtk->Pa[i + j * na] = rtk->P[i + j * nx];
1754 }
1755 }
1756 for (i = 0; i < nb; i++) {
1757 bias[i] = b[i];
1758 y[na + i] -= b[i];
1759 }
1760 if (!matinv(Qb, nb)) {
1761 matmul("NN", nb, 1, nb, 1.0, Qb, y + na, 0.0, db);
1762 matmul("NN", na, 1, nb, -1.0, Qab, db, 1.0, rtk->xa);
1763
1764 /* covariance of fixed solution (Qa=Qa-Qab*Qb^-1*Qab') */
1765 matmul("NN", na, nb, nb, 1.0, Qab, Qb, 0.0, QQ);
1766 matmul("NT", na, na, nb, -1.0, QQ, Qab, 1.0, rtk->Pa);
1767
1768 trace(NULL, 3, "resamb : validation ok (nb=%d ratio=%.2f s=%.2f/%.2f)\n", nb,
1769 s[0] == 0.0 ? 0.0 : s[1] / s[0], s[0], s[1]);
1770
1771 /* restore single-differenced ambiguity */
1772 restamb(rtk, bias, nb, xa);
1773 } else {
1774 nb = 0;
1775 }
1776 } else { /* validation failed */
1777 nb = 0;
1778 }
1779 } else {
1780 errmsg(rtk, "lambda error (info=%d)\n", info);
1781 }
1782 free(D);
1783 free(y);
1784 free(Qy);
1785 free(DP);
1786 free(b);
1787 free(db);
1788 free(Qb);
1789 free(Qab);
1790 free(QQ);
1791
1792 return nb; /* number of ambiguities */
1793}
1794static void elsort(rtk_t* rtk, const obsd_t* obs, int* isat, double* el, const int n) {
1795 double el_;
1796 int i, j, k, m = 0, sati;
1797
1798 for (i = 0; i < n; i++) {
1799 sati = obs[i].sat;
1800 el_ = rtk->ssat[sati - 1].azel[1];
1801 if (m <= 0) {
1802 isat[m] = sati;
1803 el[m++] = el_;
1804 } else {
1805 for (j = 0; j < m; j++) {
1806 if (el_ < el[j]) {
1807 break;
1808 }
1809 }
1810 if (j >= n) {
1811 continue;
1812 }
1813 for (k = m; k > j; k--) {
1814 isat[k] = isat[k - 1];
1815 el[k] = el[k - 1];
1816 }
1817 isat[j] = sati;
1818 el[j] = el_;
1819 if (m < n) {
1820 m++;
1821 }
1822 }
1823 }
1824}
1825
1826/* compute single difference dop */
1827static double sddop(const rtk_t* rtk, const double* H, const int* vflg, const double elmask, const int nv) {
1828 int i, cnt, freq, type, sat;
1829 double *H_, Q[9], dop = 1000.0;
1830
1831 H_ = mat(3, nv);
1832
1833 for (i = cnt = 0; i < nv; i++) {
1834 sat = (vflg[i] >> 8) & 0xFF;
1835 type = (vflg[i] >> 4) & 0xF;
1836 freq = (vflg[i]) & 0xF;
1837 if (freq > 0 || type > 0 || rtk->ssat[sat - 1].vsat[0] == 0 || rtk->ssat[sat - 1].azel[1] <= elmask) {
1838 continue;
1839 }
1840 H_[cnt * 3] = H[i * rtk->nx];
1841 H_[cnt * 3 + 1] = H[i * rtk->nx + 1];
1842 H_[cnt * 3 + 2] = H[i * rtk->nx + 2];
1843 cnt++;
1844 }
1845 matmul("NT", 3, 3, cnt, 1.0, H_, H_, 0.0, Q);
1846 if (cnt > 0 && !matinv(Q, 3)) {
1847 dop = sqrt(Q[0] + Q[4] + Q[8]); /*PDOP*/
1848 }
1849
1850 free(H_);
1851
1852 return dop;
1853}
1854
1855/* compute zero difference (conventional) dop */
1856static double zddop(const rtk_t* rtk, const double* azel, const obsd_t* obs, const double elmask, const int n) {
1857 int i, j;
1858 double *azel_, dop[4];
1859
1860 azel_ = mat(2, n);
1861
1862 for (i = j = 0; i < n; i++) {
1863 if (rtk->ssat[obs[i].sat - 1].vsat[0] == 0) {
1864 continue;
1865 }
1866 azel_[j * 2] = azel[i * 2];
1867 azel_[j * 2 + 1] = azel[i * 2 + 1];
1868 j++;
1869 }
1870
1871 dops(j, azel_, elmask, dop);
1872
1873 free(azel_);
1874
1875 return dop[1]; /* PDOP */
1876}
1877
1878/* compute dop */
1879static double calpdop(const rtk_t* rtk, const double* H, const double* azel, const int* vflg, const obsd_t* obs,
1880 const double elmask, const int n, const int nv) {
1881 switch (rtk->opt.refdop) {
1882 case 0:
1883 return zddop(rtk, azel, obs, elmask, n);
1884 break;
1885 case 1:
1886 return sddop(rtk, H, vflg, elmask, nv);
1887 break;
1888 }
1889
1890 return 0.0;
1891}
1892static void resetfilterflag(nav_t* nav) {
1893 nav->filreset = 0;
1894 clearsatcorr();
1895}
1896
1897static void check_vrs_facility(nav_t* nav, const obsd_t* obs, int nu, int* nr, int l6mrg) {
1898 static int savefacility[SSR_CH_NUM] = {-1, -1};
1899 int ch, idx;
1900
1901 for (ch = 0; ch < (l6mrg ? SSR_CH_NUM : 1); ch++) {
1902 idx = (ch == 0) ? nu : nu + nr[0];
1903 if (savefacility[ch] != obs[idx].facility) {
1904 if (savefacility[ch] != -1) {
1905 trace(NULL, 1, "VRS facility changed[ch:%d], %d ---> %d\n", ch, savefacility[ch], obs[idx].facility);
1906 nav->filreset = 1;
1907 } else {
1908 trace(NULL, 1, "VRS facility changed[ch:%d], ---> %d\n", ch, obs[idx].facility);
1909 }
1910 }
1911 savefacility[ch] = obs[idx].facility;
1912 }
1913}
1914
1915/*============================================================================
1916 * Residual test (ported from upstream claslib/src/rtkcmn.c)
1917 *===========================================================================*/
1918static int residual_test(rtk_t* rtk, const int* vflg, const double* v, const double* Q, int* iv, const int m, int* ns,
1919 const int flg) {
1920 int i, k, sat1, sat2, type, freq, freq_, nf = rtk->opt.nf;
1921 char* stype;
1922 double *inno = rtk->opt.maxinno_ext, gamma;
1923 double *resc, *resp, sig2, vv0_, vvf_, vv = 0.0;
1924 int qi, qj, cnt, cnt2;
1925
1926 resc = zeros(MAXSAT * nf, 1);
1927 resp = zeros(MAXSAT * nf, 1);
1928
1929 /* residuals validation */
1930 for (i = cnt2 = 0; i < m; i++) {
1931 sat2 = (vflg[i] >> 8) & 0xFF;
1932 type = (vflg[i] >> 4) & 0xF;
1933 freq = (vflg[i]) & 0xF;
1934 if (type == 0) {
1935 stype = "L";
1936 } else if (type == 1) {
1937 stype = "P";
1938 } else {
1939 stype = "C";
1940 }
1941 (void)stype;
1942 if (type == 0) {
1943 resc[(sat2 - 1) * nf + freq] = v[i];
1944 resp[(sat2 - 1) * nf + freq] = Q[i + i * m];
1945 cnt2++;
1946 }
1947 }
1948
1949 for (i = 0; i < MAXSAT; i++) {
1950 for (k = 1; k < nf; k++) {
1951 qi = 0;
1952 qj = k;
1953 gamma = SQR(lam_carr[qj]) / SQR(lam_carr[qi]);
1954 if (resp[i * nf + qi] == 0.0 || resp[i * nf + qj] == 0.0) {
1955 continue;
1956 }
1957 vvf_ = FREQ1 / FREQ2 * (resc[i * nf + qi] - resc[i * nf + qj]) / (1.0 - gamma);
1958 vv0_ = (gamma * resc[i * nf + qi] - resc[i * nf + qj]) / (gamma - 1.0);
1959 sig2 = resp[i * nf + qi] > resp[i * nf + qj] ? resp[i * nf + qi] : resp[i * nf + qj];
1960
1961 /* residuals validation for dispersive term */
1962 if ((inno[1] > 0.0) && (vvf_ * vvf_ > inno[1] * inno[1] * sig2)) {
1963 rtk->ssat[i].vsat[qi] = rtk->ssat[i].vsat[qj] = 0;
1964 rtk->ssat[i].rejc[qi]++;
1965 rtk->ssat[i].rejc[qj]++;
1966 trace(NULL, 2,
1967 "Residuals validation:freq outlier detected (stage=%1d sat=%2d f=%d %d v=%6.3f sig=%.3f)\n", flg,
1968 i + 1, qi, qj, vvf_, sqrt(sig2));
1969 continue;
1970 }
1971 /* residuals validation for non-dispersive(L0) term */
1972 if ((inno[2] > 0.0) && (vv0_ * vv0_ > inno[2] * inno[2] * sig2)) {
1973 rtk->ssat[i].vsat[qi] = rtk->ssat[i].vsat[qj] = 0;
1974 rtk->ssat[i].rejc[qi]++;
1975 rtk->ssat[i].rejc[qj]++;
1976 trace(NULL, 2, "Residuals validation:L0 outlier detected (stage=%1d sat=%2d v=%6.3f sig=%.3f)\n", flg,
1977 i + 1, vv0_, sqrt(sig2));
1978 continue;
1979 }
1980 }
1981 }
1982
1983 /* pre/post residuals test */
1984 for (i = k = cnt = 0; i < m; i++) {
1985 sat1 = (vflg[i] >> 16) & 0xFF;
1986 sat2 = (vflg[i] >> 8) & 0xFF;
1987 type = (vflg[i] >> 4) & 0xF;
1988 freq = (vflg[i]) & 0xF;
1989 if (type == 0) {
1990 stype = "L";
1991 } else if (type == 1) {
1992 stype = "P";
1993 } else {
1994 stype = "C";
1995 }
1996 (void)stype;
1997 freq_ = type < 1 ? freq : freq * type + nf;
1998
1999 /* validation for L1/L2 residuals */
2000 if (!Q[i + i * m] || inno[0] == 0 || (v[i] * v[i] < Q[i + i * m] * inno[0] * inno[0])) {
2001 if (type == 0) {
2002 if (!rtk->ssat[sat2 - 1].vsat[freq]) {
2003 continue;
2004 }
2005 vv += v[i] * v[i] / Q[i + i * m];
2006 cnt++;
2007 }
2008 ns[freq_]++;
2009 iv[k++] = i;
2010 } else {
2011 if (type == 0) {
2012 rtk->ssat[sat2 - 1].vsat[freq] = 0;
2013 }
2014 trace(NULL, 2, "meas update outlier detected (sat=%2d-%2d %s%d v=%6.3f sig=%6.3f outc=%2d)\n", sat1, sat2,
2015 stype, freq + 1, v[i], sqrt(Q[i + i * m]), rtk->ssat[sat2 - 1].outc[freq]);
2016 }
2017 }
2018
2019 /* chi-square validation for L1/L2 carrier residuals */
2020 if (cnt > NP(&rtk->opt)) {
2021 if (vv > 1.0 * chisqr[cnt - NP(&rtk->opt) - 1]) {
2022 trace(NULL, 2, "L1/L2 chi-square validation failed (stage=%d cnt=%d np=%d vv=%.2f cs=%5.3f)\n", flg, cnt,
2023 NP(&rtk->opt), vv, chisqr[cnt - NP(&rtk->opt) - 1]);
2024 }
2025 vrs_chisq = vv / chisqr[cnt - NP(&rtk->opt) - 1];
2026 } else {
2027 vrs_chisq = ((double)cnt / cnt2) < 0.5 ? 100.0 : 0.0;
2028 }
2029
2030 free(resc);
2031 free(resp);
2032
2033 return k;
2034}
2035
2036/*============================================================================
2037 * Kalman filter with state compaction (ported from upstream rtkcmn.c)
2038 *===========================================================================*/
2039static int filter2_(rtk_t* rtk, const double* x, const double* P, const double* H, const double* v, const double* R,
2040 const int n, const int m, double* xp, double* Pp, double* Qp, const int* vflg, const int flg) {
2041 double *K, *v_, *Q_, *F_, *H_, *I, *F, *Q, *vd;
2042 int i, j, k, info = 0, *iv, nn[3 * 2];
2043
2044 Q = mat(m, m);
2045 F = mat(n, m);
2046 iv = imat(m, 1);
2047 nn[0] = nn[1] = nn[2] = nn[3] = 0;
2048 matcpy(Q, R, m, m);
2049 matcpy(xp, x, n, 1);
2050 matmul("NN", n, m, n, 1.0, P, H, 0.0, F); /* Q=H'*P*H+R */
2051 matmul("TN", m, m, n, 1.0, H, F, 1.0, Q);
2052
2053 /* prefit/postfit residuals test */
2054 if (!(k = residual_test(rtk, vflg, v, Q, iv, m, nn, flg))) {
2055 trace(NULL, 2, "filter status (all measurements are rejected): stage=%1d\n", flg);
2056 } else {
2057 trace(NULL, 2, "filter status : stage=%1d refsat=%2d states=%2d obs=%2d,vobs=%2d,vsat(L1)=%2d,chi-val=%5.3f\n",
2058 flg, (vflg[iv[0]] >> 16) & 0xFF, n, m, k, nn[0], vrs_chisq);
2059 }
2060
2061 if (flg > 0 || k == 0) {
2062 matcpy(Pp, P, n, n);
2063 free(Q);
2064 free(F);
2065 free(iv);
2066 return info;
2067 }
2068
2069 v_ = mat(k, 1);
2070 Q_ = mat(k, k);
2071 F_ = mat(n, k);
2072 H_ = mat(n, k);
2073 K = mat(n, k);
2074 vd = mat(n, 1);
2075 for (i = 0; i < k; i++) {
2076 v_[i] = v[iv[i]];
2077 for (j = 0; j < k; j++) {
2078 Q_[i + j * k] = Q[iv[i] + iv[j] * m];
2079 }
2080 for (j = 0; j < n; j++) {
2081 F_[j + i * n] = F[j + iv[i] * n];
2082 }
2083 for (j = 0; j < n; j++) {
2084 H_[j + i * n] = H[j + iv[i] * n];
2085 }
2086 }
2087
2088 /* -------LU decomposition -------------------------*/
2089 I = eye(n);
2090 if (!(info = matinv(Q_, k))) {
2091 matmul("NN", n, k, k, 1.0, F_, Q_, 0.0, K); /* K=P*H*Q^-1 */
2092 matmul("NN", n, 1, k, 1.0, K, v_, 1.0, xp); /* xp=x+K*v */
2093 matmul("NT", n, n, k, -1.0, K, H_, 1.0, I);
2094 matmul("NN", n, n, n, 1.0, I, P, 0.0, Pp); /* Pp=(I-K*H')*P */
2095 if (flg == 0 && (rtk->opt.ionoopt == IONOOPT_EST_ADPT || rtk->opt.prnadpt)) {
2096 matmul("NN", n, 1, k, 1.0, K, v_, 0.0, vd); /* vd=K*v */
2097 matmul("NT", n, n, 1, 1.0, vd, vd, 0.0, Qp); /* Q=vd(vd)' */
2098 }
2099 }
2100 free(I);
2101
2102 free(iv);
2103 free(v_);
2104 free(Q_);
2105 free(F_);
2106 free(H_);
2107 free(K);
2108 free(Q);
2109 free(F);
2110 free(vd);
2111
2112 return info;
2113}
2114
2115static int filter2(rtk_t* rtk, double* x, double* P, double* Q, const double* H, const double* v, const double* R,
2116 const int n, const int m, const int* vflg, const int flg) {
2117 double *x_, *xp_, *P_, *Pp_, *Qp_, *H_;
2118 int i, j, k, info, *ix;
2119
2120 ix = imat(n, 1);
2121 for (i = k = 0; i < n; i++) {
2122 if (x[i] != 0.0 && P[i + i * n] > 0.0) {
2123 ix[k++] = i;
2124 }
2125 }
2126 x_ = mat(k, 1);
2127 xp_ = mat(k, 1);
2128 P_ = mat(k, k);
2129 Pp_ = mat(k, k);
2130 H_ = mat(k, m);
2131 Qp_ = mat(k, k);
2132 for (i = 0; i < k; i++) {
2133 x_[i] = x[ix[i]];
2134 for (j = 0; j < k; j++) {
2135 P_[i + j * k] = P[ix[i] + ix[j] * n];
2136 }
2137 for (j = 0; j < m; j++) {
2138 H_[i + j * k] = H[ix[i] + j * n];
2139 }
2140 }
2141 info = filter2_(rtk, x_, P_, H_, v, R, k, m, xp_, Pp_, Qp_, vflg, flg);
2142 for (i = 0; i < k; i++) {
2143 x[ix[i]] = xp_[i];
2144 for (j = 0; j < k; j++) {
2145 P[ix[i] + ix[j] * n] = Pp_[i + j * k];
2146 }
2147 if (flg == 0 && (rtk->opt.ionoopt == IONOOPT_EST_ADPT || rtk->opt.prnadpt)) {
2148 for (j = 0; j < k; j++) {
2149 Q[ix[i] + ix[j] * n] = Qp_[i + j * k];
2150 }
2151 }
2152 }
2153 free(ix);
2154 free(x_);
2155 free(xp_);
2156 free(P_);
2157 free(Pp_);
2158 free(H_);
2159 free(Qp_);
2160 return info;
2161}
2162/* relative positioning ------------------------------------------------------*/
2163extern int relposvrs(rtk_t* rtk, const obsd_t* obs, int nu, int* nr, nav_t* nav) {
2164 static gtime_t regularly = {-1, 0.0};
2165 prcopt_t* opt = &rtk->opt;
2166 gtime_t time = obs[0].time;
2167 double *rs, *dts, *var, *y, *e, *azel, *v, *H, *R, *xp, *Pp, *Qp, *xa, *bias, dt[SSR_CH_NUM];
2168 int i, j, k = 0, f, l, n, ns[SSR_CH_NUM] = {0}, ny, nv = 0, sat[SSR_CH_NUM][MAXSAT], iu[SSR_CH_NUM][MAXSAT],
2169 ir[SSR_CH_NUM][MAXSAT], sati;
2170 int info, vflg[MAXOBS * NFREQ * 2 + 1], svh[MAXOBS * 2];
2171 int stat = rtk->opt.mode <= PMODE_DGPS ? SOLQ_DGPS : SOLQ_FLOAT;
2172 int nf = opt->ionoopt == IONOOPT_IFLC ? 1 : opt->nf;
2173 int ch, ofst;
2174 double pdop;
2175 int nb, isat[MAXSAT];
2176 double el[MAXSAT], dist;
2177 static int resetcnt = 0, cntdiffp = 0, np;
2178 double tow = time2gpst(obs[0].time, NULL);
2179
2180 trace(NULL, 2, "relposvrs : time=%s nx=%d nu=%d nr[0]=%d nr[1]=%d\n", time_str(obs[0].time, 0), rtk->nx, nu, nr[0],
2181 nr[1]);
2182
2183 n = opt->l6mrg ? nu + nr[0] + nr[1] : nu + nr[0];
2184
2185 check_vrs_facility(nav, obs, nu, nr, opt->l6mrg);
2186
2187 regularly = (regularly.time == -1 ? obs[0].time : regularly);
2188 if (opt->regularly != 0 && timediff(obs[0].time, regularly) >= (double)opt->regularly) {
2189 trace(NULL, 1, "relposvrs(): regularly reset filter, tow=%.1f\n", tow);
2190 for (i = 0; i < rtk->nx; i++) {
2191 rtk->x[i] = 0.0;
2192 }
2193 rtk->sol.stat = SOLQ_SINGLE;
2194 regularly = obs[0].time;
2195 resetfilterflag(nav);
2196 resetcnt = 0;
2197 return 0;
2198 }
2199
2200 if (nav->filreset == 1) {
2201 trace(NULL, 1, "relposvrs(): reset filter, tow=%.1f\n", time2gpst(obs[0].time, NULL));
2202 for (i = NP(opt); i < rtk->nx; i++) {
2203 rtk->x[i] = 0.0;
2204 }
2205 rtk->sol.stat = SOLQ_FLOAT;
2206 resetfilterflag(nav);
2207 stat = SOLQ_FLOAT;
2208 resetcnt = 0;
2209 }
2210 ++resetcnt;
2211
2212 /* temporal update of states */
2213 udstate(rtk, obs, nu, nav);
2214
2215 dt[0] = timediff(time, obs[nu].time);
2216 if (opt->l6mrg) {
2217 dt[1] = timediff(time, obs[nu + nr[0]].time);
2218 }
2219
2220 rs = mat(6, n);
2221 dts = mat(2, n);
2222 var = mat(1, n);
2223 y = mat(nf * 2, n);
2224 e = mat(3, n);
2225 azel = zeros(2, n);
2226
2227 for (i = 0; i < MAXSAT; i++) {
2228 rtk->ssat[i].sys = satsys(i + 1, NULL);
2229 for (j = 0; j < NFREQ; j++) {
2230 rtk->ssat[i].vsat[j] = rtk->ssat[i].snr[j] = 0;
2231 }
2232 }
2233 /* satellite positions/clocks */
2234 satposs(time, obs, n, nav, opt->sateph, rs, dts, var, svh);
2235
2236 for (ch = 0; ch < (opt->l6mrg ? SSR_CH_NUM : 1); ch++) {
2237 ofst = ch == 1 ? nu + nr[0] : nu;
2238 /* undifferenced residuals for base station */
2239 if (!zdres(rtk, 1, obs + ofst, nr[ch], rs + ofst * 6, dts + ofst * 2, svh + ofst, nav, rtk->rb, opt, 1,
2240 y + ofst * nf * 2, e + ofst * 3, azel + ofst * 2)) {
2241 errmsg(rtk, "initial base station position error\n");
2242
2243 free(rs);
2244 free(dts);
2245 free(var);
2246 free(y);
2247 free(e);
2248 free(azel);
2249 return 0;
2250 }
2251 /* time-interpolation of residuals (for post-processing) */
2252 if (opt->intpref) {
2253 dt[ch] = intpres(time, obs + ofst, nr[ch], nav, rtk, ch, y + ofst * nf * 2);
2254 }
2255 /* select common satellites between rover and base-station */
2256 ns[ch] = selsat(obs, azel, nu, nr, opt, sat[ch], iu[ch], ir[ch], ch);
2257 }
2258 if (ns[0] + ns[1] == 0) {
2259 errmsg(rtk, "no common satellite\n");
2260 free(rs);
2261 free(dts);
2262 free(var);
2263 free(y);
2264 free(e);
2265 free(azel);
2266 return 0;
2267 }
2268
2269 trace(NULL, 4, "x(0)=");
2270 tracemat(NULL, 4, rtk->x, 1, NR(opt), 13, 4);
2271
2272 xp = mat(rtk->nx, 1);
2273 Pp = zeros(rtk->nx, rtk->nx);
2274 xa = mat(rtk->nx, 1);
2275 matcpy(xp, rtk->x, rtk->nx, 1);
2276
2277 ny = (ns[0] + ns[1]) * nf * 2 + 2;
2278 v = mat(ny, 1);
2279 H = zeros(rtk->nx, ny);
2280 R = mat(ny, ny);
2281 bias = mat(rtk->nx, 1);
2282 Qp = zeros(rtk->nx, rtk->nx);
2283
2284 ch = 0;
2285 for (i = 0; i < rtk->opt.niter; i++) {
2286 for (k = 0; k < MAXREF; k++) {
2287 nv = zdres(rtk, 0, obs, nu, rs, dts, svh, nav, xp, opt, 0, y, e, azel);
2288 /* trop estimation params is set to be zero when trop ssr is valid */
2289 if ((resetcnt == 1) && rtk->opt.tropopt >= TROPOPT_EST && nav->invtrop == 0) {
2290 j = ITT(&rtk->opt);
2291 initx(rtk, INIT_ZWD, 0.0, j);
2292 rtk->opt.std[2] = rtk->opt.prn[2] = 0.0;
2293 if (rtk->opt.tropopt == TROPOPT_ESTG) {
2294 for (l = 0; l < 2; l++) {
2295 initx(rtk, 1E-6, 0.0, ++j);
2296 }
2297 }
2298 }
2299 if (!nv) {
2300 trace(NULL, 2, "rover initial position error\n");
2301 stat = SOLQ_NONE;
2302 break;
2303 }
2304 /* double-differenced residuals and partial derivatives */
2305 if ((nv = ddres(rtk, nav, dt, xp, Pp, obs, nu, sat, y, e, azel, iu, ir, ns, v, H, R, vflg, k, FST)) <= 0) {
2306 trace(NULL, 2, "no double-differenced residual\n");
2307 stat = SOLQ_NONE;
2308 break;
2309 }
2310 /* Kalman filter measurement update */
2311 matcpy(Pp, rtk->P, rtk->nx, rtk->nx);
2312 if ((info = filter2(rtk, xp, Pp, Qp, H, v, R, rtk->nx, nv, vflg, 0)) != 0) {
2313 trace(NULL, 2, "ppp-rtk filter error %s info=%d\n", time_str(rtk->sol.time, 0), info);
2314 continue;
2315 }
2316 /* postfit residuals of float solution */
2317 if (zdres(rtk, 0, obs, nu, rs, dts, svh, nav, xp, opt, 0, y, e, azel)) {
2318 nv = ddres(rtk, nav, dt, xp, Pp, obs, nu, sat, y, e, azel, iu, ir, ns, v, H, R, vflg, k, SND);
2319 /* validation of float solution */
2320 filter2(rtk, xp, Pp, NULL, H, v, R, rtk->nx, nv, vflg, 1);
2321 if (vrs_chisq < 100.0) {
2322 break;
2323 }
2324 }
2325 matcpy(xp, rtk->x, rtk->nx, 1);
2326 if (k == (MAXREF - 1)) {
2327 stat = SOLQ_NONE;
2328 trace(NULL, 2, "measurement update failed \n");
2329 }
2330 }
2331 }
2332
2333 /* update process noise of ionoshepre delay */
2334 if (rtk->opt.ionoopt == IONOOPT_EST_ADPT) {
2335 for (i = 0; i < nu; i++) {
2336 j = II(obs[i].sat, &rtk->opt);
2337 rtk->Q[j + j * rtk->nx] = rtk->opt.forgetion * rtk->Q[j + j * rtk->nx] +
2338 (1.0 - rtk->opt.forgetion) * SQR(rtk->opt.afgainion) * Qp[j + j * rtk->nx];
2339 }
2340 }
2341 /* update process noise of pos, vel, and acc */
2342 if (rtk->opt.prnadpt) {
2343 np = NP(opt);
2344 for (i = 0; i < np; i++) {
2345 for (j = 0; j < np; j++) {
2346 rtk->Q[i + j * rtk->nx] = rtk->opt.forgetpva * rtk->Q[i + j * rtk->nx] +
2347 (1.0 - rtk->opt.forgetpva) * SQR(rtk->opt.afgainpva) * Qp[i + j * rtk->nx];
2348 }
2349 }
2350 }
2351
2352 if (stat != SOLQ_NONE) {
2353 /* upate state and covariance matrix */
2354 matcpy(rtk->x, xp, rtk->nx, 1);
2355 matcpy(rtk->P, Pp, rtk->nx, rtk->nx);
2356 /* update ambiguity control struct */
2357 rtk->sol.ns = 0;
2358 for (i = 0; i < nu; i++) {
2359 sati = obs[i].sat;
2360 for (f = 0; f < nf; f++) {
2361 if (!rtk->ssat[sati - 1].vsat[f]) {
2362 continue;
2363 }
2364 rtk->ssat[sati - 1].lock[f]++;
2365 rtk->ssat[sati - 1].outc[f] = 0;
2366 if (f == 0) {
2367 rtk->sol.ns++; /* valid satellite count by L1 */
2368 }
2369 }
2370 }
2371 /* lack of valid satellites */
2372 if (rtk->sol.ns < 4) {
2373 stat = SOLQ_NONE;
2374 }
2375 }
2376
2377 /* keep symmetry of covariance matrix */
2378 for (i = 0; i < rtk->nx; i++) {
2379 for (j = 0; j < i; j++) {
2380 rtk->P[i + j * rtk->nx] = rtk->P[j + i * rtk->nx] =
2381 (rtk->P[i + j * rtk->nx] + rtk->P[j + i * rtk->nx]) * 0.5;
2382 }
2383 }
2384
2385 /* resolve integer ambiguity by LAMBDA */
2386 pdop = calpdop(rtk, H, azel, vflg, obs, rtk->opt.elmaskar, nu, nv);
2387 if (stat != SOLQ_NONE && (rtk->opt.maxpdopar == 0.0 || pdop <= rtk->opt.maxpdopar)) {
2388 if ((nb = resamb_LAMBDA(rtk, bias, xa)) <= 1 && rtk->opt.posopt[6]) {
2389 int exsat = 0;
2390 double ratio = 0.0;
2391 elsort(rtk, obs, isat, el, nu);
2392 for (l = 0; l < rtk->opt.armaxdelsat; l++) {
2393 for (i = 0; i < nu; i++) {
2394 sati = isat[i];
2395 if (el[i] <= rtk->opt.elmaskar) {
2396 continue;
2397 }
2398 for (f = j = 0; f < nf; f++) {
2399 if (!rtk->ssat[sati - 1].vsat[f] || rtk->ssat[sati - 1].lock[f] <= 0) {
2400 j++;
2401 }
2402 }
2403 if (j > 1) {
2404 continue;
2405 }
2406 for (f = 0; f < nf; f++) {
2407 rtk->ssat[sati - 1].vsat[f] = 0;
2408 }
2409 pdop = calpdop(rtk, H, azel, vflg, obs, rtk->opt.elmaskar, nu, nv);
2410 if ((rtk->opt.maxpdopar == 0.0 || pdop <= rtk->opt.maxpdopar) &&
2411 (nb = resamb_LAMBDA(rtk, bias, xa)) > 1) {
2412 trace(NULL, 2, "PAR OK l=%d i=%2d nb=%2d excluded sat=%2d el=%4.1f ratio=%5.2f\n", l, i, nb,
2413 sati, el[i] * R2D, rtk->sol.ratio);
2414 break;
2415 } else {
2416 for (f = 0; f < nf; f++) {
2417 rtk->ssat[sati - 1].vsat[f] = 1;
2418 }
2419 if (ratio < rtk->sol.ratio) {
2420 exsat = sati;
2421 ratio = rtk->sol.ratio;
2422 }
2423 }
2424 }
2425 if (nb > 1) {
2426 break;
2427 } else {
2428 for (f = 0; f < nf; f++) {
2429 rtk->ssat[exsat - 1].vsat[f] = 0;
2430 }
2431 }
2432 }
2433 }
2434
2435 if (nb > 1 && zdres(rtk, 0, obs, nu, rs, dts, svh, nav, rtk->opt.mode == PMODE_FIXED ? rtk->opt.ru : xa, opt, 0,
2436 y, e, azel)) {
2437 /* post-fix residuals for fixed solution */
2438 nv = ddres(rtk, nav, dt, xa, NULL, obs, nu, sat, y, e, azel, iu, ir, ns, v, NULL, R, vflg, k, TRD);
2439 /* validation of fixed solution */
2440 filter2(rtk, xp, Pp, NULL, H, v, R, rtk->nx, nv, vflg, 2);
2441
2442 if (vrs_chisq < opt->maxinno_ext[3]) {
2443 /* hold integer ambiguity */
2444 pdop = calpdop(rtk, H, azel, vflg, obs, rtk->opt.elmaskhold, nu, nv);
2445 if (++rtk->nfix >= rtk->opt.minfix && rtk->opt.modear == ARMODE_FIXHOLD &&
2446 (rtk->opt.maxpdophold == 0.0 || pdop <= rtk->opt.maxpdophold)) {
2447 holdamb(rtk, xa);
2448 }
2449 }
2450 if (vrs_chisq < opt->maxinno_ext[4]) {
2451 stat = SOLQ_FIX;
2452 }
2453 }
2454 }
2455
2456 /* reset all states due to large differece in solution between point pos and pprtk solutions */
2457 if (rtk->opt.maxdiffp > 0.0 && norm(rtk->x, 3) > 0.0 && norm(rtk->sol.rr, 3) > 0.0) {
2458 dist = SQR(rtk->x[0] - rtk->sol.rr[0]) + SQR(rtk->x[1] - rtk->sol.rr[1]) + SQR(rtk->x[2] - rtk->sol.rr[2]);
2459 if (dist > SQR(rtk->opt.maxdiffp)) {
2460 ++cntdiffp;
2461 if (cntdiffp > rtk->opt.poserrcnt) {
2462 trace(NULL, 2, "reset all states due to large pos error dist = %.2f\n", sqrt(dist));
2463 stat = SOLQ_NONE;
2464 for (i = 0; i < 3; i++) {
2465 rtk->x[i] = rtk->sol.rr[i];
2466 rtk->P[i + i * rtk->nx] = (double)rtk->sol.qr[i];
2467 }
2468 cntdiffp = 0;
2469 }
2470 } else {
2471 cntdiffp = 0;
2472 }
2473 }
2474
2475 /* save solution status */
2476 if (stat == SOLQ_FIX) {
2477 for (i = 0; i < 3; i++) {
2478 rtk->sol.rr[i] = rtk->xa[i];
2479 rtk->sol.qr[i] = (float)rtk->Pa[i + i * rtk->na];
2480 }
2481 rtk->sol.qr[3] = (float)rtk->Pa[1];
2482 rtk->sol.qr[4] = (float)rtk->Pa[1 + 2 * rtk->na];
2483 rtk->sol.qr[5] = (float)rtk->Pa[2];
2484 } else {
2485 for (i = 0; i < 3; i++) {
2486 rtk->sol.rr[i] = rtk->x[i];
2487 rtk->sol.qr[i] = (float)rtk->P[i + i * rtk->nx];
2488 }
2489 rtk->sol.qr[3] = (float)rtk->P[1];
2490 rtk->sol.qr[4] = (float)rtk->P[1 + 2 * rtk->nx];
2491 rtk->sol.qr[5] = (float)rtk->P[2];
2492 rtk->nfix = 0;
2493 }
2494 for (i = 0; i < n; i++) {
2495 for (j = 0; j < nf; j++) {
2496 if (obs[i].L[j] == 0.0) {
2497 continue;
2498 }
2499 rtk->ssat[obs[i].sat - 1].pt[obs[i].rcv - 1][j] = obs[i].time;
2500 rtk->ssat[obs[i].sat - 1].ph[obs[i].rcv - 1][j] = obs[i].L[j];
2501 }
2502 }
2503 for (i = 0; i < ns[ch]; i++) {
2504 for (j = 0; j < nf; j++) {
2505 /* output snr of rover receiver */
2506 rtk->ssat[sat[ch][i] - 1].snr[j] = obs[iu[ch][i]].SNR[j];
2507 }
2508 }
2509 for (i = 0; i < MAXSAT; i++) {
2510 for (j = 0; j < nf; j++) {
2511 if (rtk->ssat[i].fix[j] == 2 && stat != SOLQ_FIX) {
2512 rtk->ssat[i].fix[j] = 1;
2513 }
2514 if (rtk->ssat[i].slip[j] & 1) {
2515 rtk->ssat[i].slipc[j]++;
2516 }
2517 }
2518 }
2519
2520 free(rs);
2521 free(dts);
2522 free(var);
2523 free(y);
2524 free(e);
2525 free(azel);
2526 free(xp);
2527 free(Pp);
2528 free(xa);
2529 free(v);
2530 free(H);
2531 free(R);
2532 free(bias);
2533 free(Qp);
2534
2535 if (stat != SOLQ_NONE) {
2536 rtk->sol.stat = stat;
2537 }
2538
2539 if (stat == SOLQ_NONE) {
2540 for (i = 0; i < rtk->nx; i++) {
2541 rtk->x[i] = 0.0;
2542 }
2543 }
2544
2545 return stat != SOLQ_NONE;
2546}
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.
SQR
#define SQR(x)
Definition mrtk_clas.c:46
mrtk_clas.h
CLAS (QZSS L6D) CSSR decoder and grid correction type definitions.
isL2C
int isL2C(int code)
Check if observation code is L2C.
Definition mrtk_clas_isb.c:546
POSL1L2L5
#define POSL1L2L5
Definition mrtk_clas_osr.c:116
POSL1
#define POSL1
Definition mrtk_clas_osr.c:113
POSL1L5_L2
#define POSL1L5_L2
Definition mrtk_clas_osr.c:117
POSL1L5
#define POSL1L5
Definition mrtk_clas_osr.c:115
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.
mrtk_context.h
MRTKLIB Context — runtime context for processing, trace, and error handling.
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
dops
void dops(int ns, const double *azel, double elmin, double *dop)
Compute DOP (dilution of precision).
Definition mrtk_coords.c:545
ecef2pos
void ecef2pos(const double *r, double *pos)
Transform ECEF position to geodetic position.
Definition mrtk_coords.c:102
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_CMP
#define SYS_CMP
Definition mrtk_eph.c:70
SYS_GAL
#define SYS_GAL
Definition mrtk_eph.c:68
SYS_GLO
#define SYS_GLO
Definition mrtk_eph.c:67
mrtk_eph.h
MRTKLIB Broadcast Ephemeris Module — Satellite orbit/clock types and Kepler/GLONASS/SBAS orbit comput...
nav_t
struct nav_s nav_t
Definition mrtk_eph.h:40
mrtk_foundation.h
MRTKLIB Foundation — Common types, standard includes, and system constants (MAXSAT,...
CODE_L5Q
#define CODE_L5Q
Definition mrtk_foundation.h:210
NFREQGLO
#define NFREQGLO
Definition mrtk_foundation.h:47
CODE_L2X
#define CODE_L2X
Definition mrtk_foundation.h:203
MAXOBS
#define MAXOBS
Definition mrtk_foundation.h:151
CODE_L2I
#define CODE_L2I
Definition mrtk_foundation.h:225
CODE_L1C
#define CODE_L1C
Definition mrtk_foundation.h:186
CODE_L6I
#define CODE_L6I
Definition mrtk_foundation.h:227
CODE_L7I
#define CODE_L7I
Definition mrtk_foundation.h:212
DTTOL
#define DTTOL
Definition mrtk_foundation.h:158
MAXERRMSG
#define MAXERRMSG
Definition mrtk_foundation.h:172
CODE_L2S
#define CODE_L2S
Definition mrtk_foundation.h:201
NFREQPCV
#define NFREQPCV
Definition mrtk_foundation.h:48
SSR_CH_NUM
#define SSR_CH_NUM
Definition mrtk_foundation.h:161
CODE_L2W
#define CODE_L2W
Definition mrtk_foundation.h:205
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.
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
mrtk_nav.h
MRTKLIB Navigation Module — Navigation data types and satellite/nav utility functions.
ISBOPT_TABLE
#define ISBOPT_TABLE
Definition mrtk_nav.h:349
DFRQ2_GLO
#define DFRQ2_GLO
Definition mrtk_obs.c:77
FREQ2
#define FREQ2
Definition mrtk_obs.c:68
R2D
#define R2D
Definition mrtk_obs.c:51
FREQ1
#define FREQ1
Definition mrtk_obs.c:67
testsnr
int testsnr(int base, int idx, double el, double snr, const snrmask_t *mask)
Test SNR mask.
Definition mrtk_obs.c:132
DFRQ1_GLO
#define DFRQ1_GLO
Definition mrtk_obs.c:75
mrtk_obs.h
MRTKLIB Observation Module — Observation data types and utility functions.
mrtk_opt.h
MRTKLIB Options Module — Processing and solution option types.
SOLQ_SINGLE
#define SOLQ_SINGLE
Definition mrtk_opt.h:75
SOLQ_FLOAT
#define SOLQ_FLOAT
Definition mrtk_opt.h:72
IONOOPT_EST
#define IONOOPT_EST
Definition mrtk_opt.h:96
ARMODE_FIXHOLD
#define ARMODE_FIXHOLD
Definition mrtk_opt.h:131
SOLQ_DGPS
#define SOLQ_DGPS
Definition mrtk_opt.h:74
PMODE_DGPS
#define PMODE_DGPS
Definition mrtk_opt.h:46
TROPOPT_EST
#define TROPOPT_EST
Definition mrtk_opt.h:105
IONOOPT_EST_ADPT
#define IONOOPT_EST_ADPT
Definition mrtk_opt.h:100
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
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
NR
#define NR(opt)
Definition mrtk_ppp.c:222
VAR_VEL
#define VAR_VEL
Definition mrtk_ppp.c:195
NP
#define NP(opt)
Definition mrtk_ppp.c:212
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
resamb_LAMBDA
static int resamb_LAMBDA(rtk_t *rtk, double *bias, double *xa)
Resolve integer ambiguity by LAMBDA.
Definition mrtk_ppp_rtk.c:2280
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
refsat2
int refsat2[5 *2 *6+1]
Definition mrtk_ppp_rtk.c:163
cntdiffp
int cntdiffp
Definition mrtk_ppp_rtk.c:191
LAM_CARR_L5
#define LAM_CARR_L5
Definition mrtk_ppp_rtk.c:84
ddres_order
ddres_order
Definition mrtk_ppp_rtk.c:1439
TRD
@ TRD
Definition mrtk_ppp_rtk.c:1439
FST
@ FST
Definition mrtk_ppp_rtk.c:1439
SND
@ SND
Definition mrtk_ppp_rtk.c:1439
MAXREF
#define MAXREF
Definition mrtk_ppp_rtk.c:62
LAM_CARR_L2
#define LAM_CARR_L2
Definition mrtk_ppp_rtk.c:83
refsat
int refsat[5 *2 *2 *6+1]
Definition mrtk_ppp_rtk.c:162
chisqr
static const double chisqr[100]
Definition mrtk_ppp_rtk.c:114
lam_carr
static const double lam_carr[3]
Definition mrtk_ppp_rtk.c:86
Qp
double * Qp
Definition mrtk_ppp_rtk.c:184
INIT_ZWD
#define INIT_ZWD
Definition mrtk_ppp_rtk.c:69
use_ch
int use_ch[5][MAXSAT+1]
Definition mrtk_ppp_rtk.c:174
validobs
static int validobs(int i, int f, int nf, double *y)
Check if zero-difference observation is valid.
Definition mrtk_ppp_rtk.c:1045
LAM_CARR_L1
#define LAM_CARR_L1
Definition mrtk_ppp_rtk.c:82
pre_use_ch
int pre_use_ch[5][MAXSAT+1]
Definition mrtk_ppp_rtk.c:175
EFACT_GPS
#define EFACT_GPS
Definition mrtk_rtkpos.c:60
EFACT_GLO
#define EFACT_GLO
Definition mrtk_rtkpos.c:61
EFACT_SBS
#define EFACT_SBS
Definition mrtk_rtkpos.c:62
PRN_HWBIAS
#define PRN_HWBIAS
Definition mrtk_rtkpos.c:86
VAR_HWBIAS
#define VAR_HWBIAS
Definition mrtk_rtkpos.c:82
mrtk_rtkpos.h
MRTKLIB RTK Positioning Module — relative and precise positioning.
mrtk_sol.h
MRTKLIB Solution Module — Solution data types and I/O 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
mrtk_time.h
MRTKLIB Time Module — GNSS time types and conversion functions.
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.
clearsatcorr
static void clearsatcorr(void)
Definition mrtk_vrs.c:190
errmsg
static void errmsg(rtk_t *rtk, const char *format,...)
Definition mrtk_vrs.c:193
detslp_ll
static void detslp_ll(rtk_t *rtk, const obsd_t *obs, int n)
Definition mrtk_vrs.c:680
validobs
static int validobs(int *iu, int *ir, int jj, int f, int nf, double *y)
Definition mrtk_vrs.c:974
resamb_LAMBDA
static int resamb_LAMBDA(rtk_t *rtk, double *bias, double *xa)
Definition mrtk_vrs.c:1672
get_iriu_index
static int get_iriu_index(int sati, int ns, const int *sat)
Definition mrtk_vrs.c:1065
zdres
static int zdres(rtk_t *rtk, int base, const obsd_t *obs, int n, const double *rs, const double *dts, const int *svh, const nav_t *nav, const double *rr, const prcopt_t *opt, int index, double *y, double *e, double *azel)
Definition mrtk_vrs.c:895
zddop
static double zddop(const rtk_t *rtk, const double *azel, const obsd_t *obs, const double elmask, const int n)
Definition mrtk_vrs.c:1856
test_sys
static int test_sys(int sys, int m, int qzsmodear, int code)
Definition mrtk_vrs.c:1043
holdamb
static void holdamb(rtk_t *rtk, const double *xa)
Definition mrtk_vrs.c:1625
valid_zd_sig
static int valid_zd_sig(rtk_t *rtk, const obsd_t *obs, const nav_t *nav, double *y, int nf, int f, int j, int m, int ch, int sati, int iu[][MAXSAT], int ir[][MAXSAT], int sat[][MAXSAT], int *ff, int *ns, int *satj, int *jj, int *sysi, int *sysj, double *lami, double *lamj)
Definition mrtk_vrs.c:1075
ddmat
static int ddmat(rtk_t *rtk, double *D)
Definition mrtk_vrs.c:1527
ddres
static int ddres(rtk_t *rtk, const nav_t *nav, double *dt, const double *x, const double *P, const obsd_t *obs, int n, int sat[][MAXSAT], double *y, double *e, double *azel, int iu[][MAXSAT], int ir[][MAXSAT], int *ns, double *v, double *H, double *R, int *vflg, int niter, int order)
Definition mrtk_vrs.c:1152
initx
static void initx(rtk_t *rtk, double xi, double var, int i)
Definition mrtk_vrs.c:252
sddop
static double sddop(const rtk_t *rtk, const double *H, const int *vflg, const double elmask, const int nv)
Definition mrtk_vrs.c:1827
baseline
static double baseline(const double *ru, const double *rb, double *dr)
Definition mrtk_vrs.c:244
varerr
static double varerr(int sat, int sys, double el, double bl, double dt, int f, const prcopt_t *opt)
Definition mrtk_vrs.c:208
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)
Definition mrtk_vrs.c:857
IL
#define IL(f, opt)
Definition mrtk_vrs.c:102
filter2_
static int filter2_(rtk_t *rtk, const double *x, const double *P, const double *H, const double *v, const double *R, const int n, const int m, double *xp, double *Pp, double *Qp, const int *vflg, const int flg)
Definition mrtk_vrs.c:2039
selsat
static int selsat(const obsd_t *obs, double *azel, int nu, int *nr, const prcopt_t *opt, int *sat, int *iu, int *ir, int ch)
Definition mrtk_vrs.c:260
udtrop
static void udtrop(rtk_t *rtk, double tt, double bl)
Definition mrtk_vrs.c:613
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_vrs.c:1008
residual_test
static int residual_test(rtk_t *rtk, const int *vflg, const double *v, const double *Q, int *iv, const int m, int *ns, const int flg)
Definition mrtk_vrs.c:1918
check_same_fac
static int check_same_fac(const obsd_t *obs, int ir[][MAXSAT])
Definition mrtk_vrs.c:1061
NP
#define NP(opt)
Definition mrtk_vrs.c:91
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
ITT
#define ITT(opt)
Definition mrtk_vrs.c:104
selfreqpair
static int selfreqpair(const int sat, const prcopt_t *opt, const obsd_t *obs)
Definition mrtk_vrs.c:291
gfmeasL1L2
static double gfmeasL1L2(const obsd_t *obs, const nav_t *nav)
Definition mrtk_vrs.c:661
resetfilterflag
static void resetfilterflag(nav_t *nav)
Definition mrtk_vrs.c:1892
IT
#define IT(r, opt)
Definition mrtk_vrs.c:101
NF
#define NF(opt)
Definition mrtk_vrs.c:90
qfrtk
static const double qfrtk[6][60]
Definition mrtk_vrs.c:118
gloicbcorr
static double gloicbcorr(int sat1, int sat2, const prcopt_t *opt, double lam1, double lam2, int f)
Definition mrtk_vrs.c:1031
obs_code_for_freq_vrs
static uint8_t obs_code_for_freq_vrs(int sat, int f)
Definition mrtk_vrs.c:153
ddcov
static void ddcov(const int *nb, int n, const double *Ri, const double *Rj, int nv, double *R)
Definition mrtk_vrs.c:989
check_vrs_facility
static void check_vrs_facility(nav_t *nav, const obsd_t *obs, int nu, int *nr, int l6mrg)
Definition mrtk_vrs.c:1897
udpos
static void udpos(rtk_t *rtk, const obsd_t *obs, double tt)
Definition mrtk_vrs.c:313
SQR
#define SQR(x)
Definition mrtk_vrs.c:53
restamb
static void restamb(rtk_t *rtk, const double *bias, int nb, double *xa)
Definition mrtk_vrs.c:1591
intpres
static double intpres(gtime_t time, const obsd_t *obs, int n, const nav_t *nav, rtk_t *rtk, int ch, double *y)
Definition mrtk_vrs.c:1478
elsort
static void elsort(rtk_t *rtk, const obsd_t *obs, int *isat, double *el, const int n)
Definition mrtk_vrs.c:1794
vrs_chisq
static double vrs_chisq
Definition mrtk_vrs.c:187
udstate
static void udstate(rtk_t *rtk, const obsd_t *obs, int nu, const nav_t *nav)
Definition mrtk_vrs.c:830
detslp_gf
static void detslp_gf(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
Definition mrtk_vrs.c:697
udrcvbias
static void udrcvbias(rtk_t *rtk, double tt)
Definition mrtk_vrs.c:641
udion
static void udion(rtk_t *rtk, double tt, double bl, const obsd_t *obs, int ns)
Definition mrtk_vrs.c:556
calpdop
static double calpdop(const rtk_t *rtk, const double *H, const double *azel, const int *vflg, const obsd_t *obs, const double elmask, const int n, const int nv)
Definition mrtk_vrs.c:1879
IB
#define IB(s, f, opt)
Definition mrtk_vrs.c:103
udbias_ppp
static void udbias_ppp(rtk_t *rtk, double tt, const obsd_t *obs, int n, const nav_t *nav)
Definition mrtk_vrs.c:736
filter2
static int filter2(rtk_t *rtk, double *x, double *P, double *Q, const double *H, const double *v, const double *R, const int n, const int m, const int *vflg, const int flg)
Definition mrtk_vrs.c:2115
gfmeasL1L5
static double gfmeasL1L5(const obsd_t *obs, const nav_t *nav)
Definition mrtk_vrs.c:670
sat_lambda
static double sat_lambda(int sat, int f, const nav_t *nav)
Definition mrtk_vrs.c:181
set_init_pb
static void set_init_pb(rtk_t *rtk, const nav_t *nav, const obsd_t *obs, int n, int f, int satj)
Definition mrtk_vrs.c:1109
II
#define II(s, opt)
Definition mrtk_vrs.c:100
mrtk_vrs.h
MRTKLIB VRS-RTK Module — double-differenced RTK using CLAS VRS.
exterr_t::cerr
double cerr[4][5 *2]
Definition mrtk_opt.h:161
exterr_t::ena
int ena[4]
Definition mrtk_opt.h:160
exterr_t::perr
double perr[4][5 *2]
Definition mrtk_opt.h:162
exterr_t::gloicb
double gloicb[5]
Definition mrtk_opt.h:164
gtime_t
GNSS time representation.
Definition mrtk_time.h:49
gtime_t::time
time_t time
Definition mrtk_time.h:50
nav_s::filreset
int filreset
Definition mrtk_nav.h:436
nav_s::sfts
sft_t sfts
Definition mrtk_nav.h:450
nav_s::invtrop
int invtrop
Definition mrtk_nav.h:437
nav_s::erp
erp_t erp
Definition mrtk_nav.h:413
obsd_t
Observation data record.
Definition mrtk_obs.h:44
obsd_t::time
gtime_t time
Definition mrtk_obs.h:45
obsd_t::sat
uint8_t sat
Definition mrtk_obs.h:46
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::facility
int facility
Definition mrtk_obs.h:53
prcopt_t
Definition mrtk_opt.h:171
prcopt_t::prnadpt
int prnadpt
Definition mrtk_opt.h:261
prcopt_t::snrmask
snrmask_t snrmask
Definition mrtk_opt.h:177
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::qzsmodear
int qzsmodear
Definition mrtk_opt.h:254
prcopt_t::nf
int nf
Definition mrtk_opt.h:174
prcopt_t::bdsmodear
int bdsmodear
Definition mrtk_opt.h:181
prcopt_t::afgainion
double afgainion
Definition mrtk_opt.h:267
prcopt_t::poserrcnt
int poserrcnt
Definition mrtk_opt.h:258
prcopt_t::maxdiffp
double maxdiffp
Definition mrtk_opt.h:264
prcopt_t::phasshft
int phasshft
Definition mrtk_opt.h:259
prcopt_t::ionoopt
int ionoopt
Definition mrtk_opt.h:187
prcopt_t::armaxdelsat
int armaxdelsat
Definition mrtk_opt.h:256
prcopt_t::pcvr
pcv_t pcvr[2]
Definition mrtk_opt.h:224
prcopt_t::stats_prnpositv
double stats_prnpositv
Definition mrtk_opt.h:272
prcopt_t::elmaskhold
double elmaskhold
Definition mrtk_opt.h:213
prcopt_t::thresslip
double thresslip
Definition mrtk_opt.h:214
prcopt_t::maxpdopar
double maxpdopar
Definition mrtk_opt.h:279
prcopt_t::forgetpva
double forgetpva
Definition mrtk_opt.h:268
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::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::rectype
char rectype[2][MAXANT]
Definition mrtk_opt.h:274
prcopt_t::odisp
double odisp[2][6 *11]
Definition mrtk_opt.h:239
prcopt_t::sclkstab
double sclkstab
Definition mrtk_opt.h:210
prcopt_t::glomodear
int glomodear
Definition mrtk_opt.h:180
prcopt_t::maxpdophold
double maxpdophold
Definition mrtk_opt.h:280
prcopt_t::afgainpva
double afgainpva
Definition mrtk_opt.h:269
prcopt_t::regularly
int regularly
Definition mrtk_opt.h:276
prcopt_t::stats_prnposith
double stats_prnposith
Definition mrtk_opt.h:271
prcopt_t::posopt
int posopt[13]
Definition mrtk_opt.h:237
prcopt_t::thresar
double thresar[8]
Definition mrtk_opt.h:211
prcopt_t::beta
double beta
Definition mrtk_opt.h:282
prcopt_t::maxout
int maxout
Definition mrtk_opt.h:183
prcopt_t::elmaskar
double elmaskar
Definition mrtk_opt.h:212
prcopt_t::modear
int modear
Definition mrtk_opt.h:179
prcopt_t::refdop
int refdop
Definition mrtk_opt.h:281
prcopt_t::mode
int mode
Definition mrtk_opt.h:172
prcopt_t::tidecorr
int tidecorr
Definition mrtk_opt.h:190
prcopt_t::exterr
exterr_t exterr
Definition mrtk_opt.h:240
prcopt_t::varholdamb
double varholdamb
Definition mrtk_opt.h:263
prcopt_t::forgetion
double forgetion
Definition mrtk_opt.h:266
prcopt_t::antdel
double antdel[2][3]
Definition mrtk_opt.h:223
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::prnionomax
double prnionomax
Definition mrtk_opt.h:270
prcopt_t::sateph
int sateph
Definition mrtk_opt.h:178
prcopt_t::maxinno_ext
double maxinno_ext[5]
Definition mrtk_opt.h:262
prcopt_t::alphaar
int alphaar
Definition mrtk_opt.h:253
prcopt_t::l6mrg
int l6mrg
Definition mrtk_opt.h:275
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::x
double * x
Definition mrtk_rtkpos.h:54
rtk_t::opt
prcopt_t opt
Definition mrtk_rtkpos.h:63
rtk_t::com_bias
double com_bias
Definition mrtk_rtkpos.h:58
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::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
sft_t::rectyp
char rectyp[100][MAXANT]
Definition mrtk_nav.h:361
sft_t::n
int n
Definition mrtk_nav.h:360
sft_t::bias
double bias[100]
Definition mrtk_nav.h:362
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::stat
uint8_t stat
Definition mrtk_sol.h:54
sol_t::ns
uint8_t ns
Definition mrtk_sol.h:55
sol_t::time
gtime_t time
Definition mrtk_sol.h:45
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::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::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::sys
uint8_t sys
Definition mrtk_sol.h:99
ssat_t::pbreset
int pbreset[5]
Definition mrtk_sol.h:125
ssat_t::fix
uint8_t fix[5]
Definition mrtk_sol.h:106
ssat_t::code
uint8_t code[5]
Definition mrtk_sol.h:124
ssat_t::resc
double resc[5]
Definition mrtk_sol.h:103
free
#define free(x)
Definition toml.c:53