MRTKLIB 0.7.14
Modernized GNSS Positioning Library
Loading...
Searching...
No Matches
mrtk_ppp_rtk.c
Go to the documentation of this file.
1/*------------------------------------------------------------------------------
2 * mrtk_ppp_rtk.c : PPP-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) 2014 T.SUZUKI
12 * Copyright (C) 2007-2023 T.TAKASU
13 *
14 * SPDX-License-Identifier: BSD-2-Clause
15 *----------------------------------------------------------------------------*/
16/**
17 * @file mrtk_ppp_rtk.c
18 * @brief PPP-RTK positioning engine.
19 *
20 * References:
21 * [1] claslib ppprtk.c — upstream PPP-RTK implementation
22 * [2] claslib rtkcmn.c — upstream filter2/residual_test
23 * [6] MacMillan et al., Atmospheric gradients and the VLBI terrestrial and
24 * celestial reference frames, Geophys. Res. Let., 1997
25 *
26 * Notes:
27 * - Ported from upstream claslib ppprtk.c as an independent positioning
28 * engine. Uses double-differencing (DD) with CLAS grid corrections,
29 * unlike the undifferenced PPP engine in mrtk_ppp.c.
30 * - Dual-channel (l6mrg) processing merges the two L6 transmit patterns:
31 * residuals are computed per channel and combined in ddres(), which picks
32 * a channel per satellite and, when the channels carry different
33 * facilities, resets the phase bias of a satellite that switches channel.
34 * With l6mrg off the engine runs on channel 0 only.
35 */
36#include "mrtklib/mrtk_ppp_rtk.h"
37
38#include <math.h>
39#include <stdio.h>
40#include <stdlib.h>
41#include <string.h>
42
43#include "mrtklib/mrtk_atmos.h"
44#include "mrtklib/mrtk_clas.h"
45#include "mrtklib/mrtk_const.h"
46#include "mrtklib/mrtk_coords.h"
47#include "mrtklib/mrtk_eph.h"
48#include "mrtklib/mrtk_foundation.h"
49#include "mrtklib/mrtk_lambda.h"
50#include "mrtklib/mrtk_mat.h"
51#include "mrtklib/mrtk_nav.h"
52#include "mrtklib/mrtk_time.h"
53#include "mrtklib/mrtk_trace.h"
54
55/*============================================================================
56 * Local Macros
57 *===========================================================================*/
58#define SQR(x) ((x) * (x))
59#define SQRT(x) ((x) <= 0.0 || (x) != (x) ? 0.0 : sqrt(x))
60#define MIN(x, y) ((x) <= (y) ? (x) : (y))
61
62#define MAXREF 2 /* max iterations searching reference sat */
63
64/* initial variances */
65#define VAR_POS SQR(30.0) /* initial variance of receiver pos (m^2) */
66#define VAR_VEL SQR(1.0) /* initial variance of receiver vel (m/s)^2 */
67#define VAR_ACC SQR(1.0) /* initial variance of receiver acc (m/ss)^2 */
68#define VAR_GRA SQR(0.001) /* gradient variance (m^2) */
69#define INIT_ZWD 0.001 /* initial ZWD (m) */
70
71#define BASELEN 1000.0 /* baseline length for iono update (m) */
72
73#define MAXPBCORSSR 20.0 /* max phase bias correction of SSR (m) */
74
75/* compile-time defaults for missing prcopt fields */
76#define VARHOLDAMB 0.001 /* hold-ambiguity constraint variance */
77#define ALPHAAR_DEFAULT 3 /* F-test significance index (5%) */
78#define MINAMB_DEFAULT 2 /* minimum ambiguities for AR */
79#define MAXOBSLOSS_DEFAULT 300.0 /* max obs gap to reset states (s) */
80
81/* carrier wavelengths */
82#define LAM_CARR_L1 (CLIGHT / FREQ1)
83#define LAM_CARR_L2 (CLIGHT / FREQ2)
84#define LAM_CARR_L5 (CLIGHT / FREQ5)
85
86static const double lam_carr[3] = {LAM_CARR_L1, LAM_CARR_L2, LAM_CARR_L5};
87
88/*============================================================================
89 * State Index Macros (from upstream ppprtk.c)
90 *
91 * State vector layout:
92 * [0..NP-1] : position (3 or 9 with dynamics)
93 * [NP..NP+NI-1] : ionosphere per satellite (MAXSAT)
94 * [NP+NI..NR-1] : troposphere (1 or 3)
95 * [NR..NR+NB-1] : phase bias per satellite per frequency
96 *===========================================================================*/
97
98#define NF_RTK(opt) ((opt)->nf)
99#define NP_RTK(opt) ((opt)->dynamics == 0 ? 3 : 9)
100#define NI_RTK(opt) (((opt)->ionoopt == IONOOPT_EST || (opt)->ionoopt == IONOOPT_EST_ADPT) ? MAXSAT : 0)
101#define NT_RTK(opt) ((opt)->tropopt < TROPOPT_EST ? 0 : ((opt)->tropopt < TROPOPT_ESTG ? 1 : 3))
102#define NR_RTK(opt) (NP_RTK(opt) + NI_RTK(opt) + NT_RTK(opt))
103#define NB_RTK(opt) (MAXSAT * NF_RTK(opt))
104#define NX_RTK(opt) (NR_RTK(opt) + NB_RTK(opt))
105
106/* state index accessors */
107#define II_RTK(s, opt) (NP_RTK(opt) + (s) - 1) /* iono state index */
108#define IT_RTK(opt) (NP_RTK(opt) + NI_RTK(opt)) /* trop state index */
109#define IB_RTK(s, f, opt) (NR_RTK(opt) + MAXSAT * (f) + (s) - 1) /* bias state index */
110
111/*============================================================================
112 * Chi-squared Table (alpha=0.001)
113 *===========================================================================*/
114static const double chisqr[100] = {
115 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,
116 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,
117 66.6, 68.0, 69.3, 70.7, 72.1, 73.4, 74.7, 76.0, 77.3, 78.6, 80.0, 81.3, 82.6, 84.0, 85.4, 86.7, 88.0,
118 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,
119 110, 112, 113, 114, 115, 116, 118, 119, 120, 122, 123, 125, 126, 127, 128, 129, 131,
120 132, 133, 134, 135, 137, 138, 139, 140, 142, 143, 144, 145, 147, 148, 149};
121
122/*============================================================================
123 * F Inverse Distribution Table
124 *===========================================================================*/
125static const double qf[6][60] = {
126 /* significance level, alpha=0.001(0.1%) */
127 {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,
128 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,
129 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,
130 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},
131 /* significance level, alpha=0.005(0.5%) */
132 {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,
133 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,
134 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,
135 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},
136 /* significance level, alpha=0.010(1.0%) */
137 {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,
138 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,
139 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,
140 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},
141 /* significance level, alpha=0.050(5.0%) */
142 {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,
143 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,
144 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,
145 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},
146 /* significance level, alpha=0.100(10.0%) */
147 {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,
148 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,
149 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,
150 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},
151 /* significance level, alpha=0.200(20.0%) */
152 {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,
153 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,
154 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,
155 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}};
156
157/*============================================================================
158 * File-scope persistent state (replaces static variables from upstream)
159 *===========================================================================*/
160static struct {
161 /* DD reference satellite tracking */
162 int refsat[NFREQ * 2 * MAXREF * 6 + 1];
163 int refsat2[NFREQ * 2 * 6 + 1];
164
165 /* CPC/pt0 persistent state from zdres, per L6 channel
166 * (CPC is freq-major: f * MAXSAT + sat-1) */
167 double cpc[SSR_CH_NUM][MAXSAT * NFREQ];
168 gtime_t pt0[SSR_CH_NUM][MAXSAT];
169
170 /* l6mrg: per-frequency channel selected for each satellite in ddres(),
171 * and the selection of the previous epoch. Indexed by satellite NUMBER
172 * (1..MAXSAT) as upstream does, hence the MAXSAT+1 width.
173 * -1 = unselected, 0/1 = L6 channel, 2 = reference satellite. */
174 int use_ch[NFREQ][MAXSAT + 1];
175 int pre_use_ch[NFREQ][MAXSAT + 1];
176
177 /* OSR context */
178 clas_osr_ctx_t osr_ctx;
179
180 /* chi-squared validation value (sol_t has no chisq field) */
181 double chisq;
182
183 /* innovation covariance Qp = K*v*(K*v)' from filter2_ (adaptive noise) */
184 double* Qp;
185 int Qp_nx;
186
187 /* consecutive FLOAT epoch counter for filter reset (floatcnt mechanism) */
188 int float_count;
189
190 /* consecutive position error counter for maxdiffp reset */
191 int cntdiffp;
192
193 /* flag: floatcnt just triggered a mass reset (suppress pbreset) */
194 int floatcnt_reset;
195
196 int initialized;
197} prtk_ctx;
198
199static void load_ppp_osr_cpc(clas_osr_ctx_t* ctx, const double* cpc, const gtime_t* pt0) {
200 int f, i;
201
202 for (f = 0; f < NFREQ; f++) {
203 for (i = 0; i < MAXSAT; i++) {
204 ctx->cpctmp[f * MAXSAT + i] = cpc[f * MAXSAT + i];
205 }
206 }
207 for (i = 0; i < MAXSAT; i++) {
208 ctx->pt0tmp[i] = pt0[i];
209 }
210}
211
212static void save_ppp_osr_cpc(const clas_osr_ctx_t* ctx, double* cpc, gtime_t* pt0) {
213 int f, i;
214
215 for (f = 0; f < NFREQ; f++) {
216 for (i = 0; i < MAXSAT; i++) {
217 cpc[f * MAXSAT + i] = ctx->cpctmp[f * MAXSAT + i];
218 }
219 }
220 for (i = 0; i < MAXSAT; i++) {
221 pt0[i] = ctx->pt0tmp[i];
222 }
223}
224
225/*============================================================================
226 * Helper: satellite-specific wavelength
227 *===========================================================================*/
228
229/**
230 * @brief Return default observation code for frequency index and system.
231 * @param[in] sat Satellite number
232 * @param[in] f Frequency index (0=L1, 1=L2, 2=L5)
233 * @return Observation code (CODE_L*)
234 */
235static uint8_t obs_code_for_freq(int sat, int f) {
236 int sys = satsys(sat, NULL);
237 switch (f) {
238 case 0:
239 return (sys == SYS_CMP) ? CODE_L2I : CODE_L1C;
240 case 1:
241 if (sys == SYS_GPS) {
242 return CODE_L2W;
243 }
244 if (sys == SYS_QZS) {
245 return CODE_L5Q; /* obsdef: QZS slot1 = L5 */
246 }
247 if (sys == SYS_GAL) {
248 return CODE_L5Q;
249 }
250 if (sys == SYS_CMP) {
251 return CODE_L7I;
252 }
253 return CODE_L2W;
254 case 2:
255 if (sys == SYS_CMP) {
256 return CODE_L6I;
257 }
258 return CODE_L5Q;
259 }
260 return CODE_L1C;
261}
262
263/**
264 * @brief Compute satellite-specific wavelength for frequency index.
265 *
266 * Uses direct obs code to frequency mapping, bypassing obsdef tables which
267 * may be rearranged by apply_pppsig() for MADOCA PPP signal selection.
268 *
269 * @param[in] sat Satellite number
270 * @param[in] f Frequency index (0,1,2)
271 * @param[in] nav Navigation data (unused, kept for API compatibility)
272 * @return Wavelength in meters, or 0.0 if invalid
273 */
274static double sat_lambda(int sat, int f, const nav_t* nav) {
275 uint8_t code = obs_code_for_freq(sat, f);
276 char* obs = code2obs(code);
277 int prn = 0, sys = satsys(sat, &prn);
278 int freq_num = freq_idx2freq_num(sys, f);
279 int fcn = 0;
280 double freq = 0.0;
281
282 if (freq_num > 0) {
283 if (sys == SYS_GLO && nav && prn >= 1 && prn <= 32 && nav->glo_fcn[prn - 1] > 0) {
284 fcn = nav->glo_fcn[prn - 1] - 8;
285 }
286 freq = freq_num2freq(sys, freq_num, fcn);
287 if (freq > 0.0) {
288 return CLIGHT / freq;
289 }
290 }
291
292 if (!obs || strlen(obs) < 2) {
293 return 0.0;
294 }
295
296 switch (obs[0]) {
297 case '1':
298 freq = FREQ1;
299 break;
300 case '2':
301 freq = FREQ2;
302 break;
303 case '5':
304 freq = FREQ5;
305 break;
306 case '6':
307 freq = FREQ6;
308 break;
309 case '7':
310 freq = FREQ7;
311 break;
312 case '8':
313 freq = FREQ8;
314 break;
315 default:
316 return 0.0;
317 }
318 return CLIGHT / freq;
319}
320
321/*============================================================================
322 * Measurement error variance
323 *===========================================================================*/
324
325/**
326 * @brief Compute measurement error variance.
327 * @param[in] sat Satellite number
328 * @param[in] sys Satellite system
329 * @param[in] el Elevation angle (rad)
330 * @param[in] type Measurement type (0:phase, 1:code, 2:iono, 3:trop)
331 * @param[in] opt Processing options
332 * @return Variance (m^2)
333 */
334static double varerr(int sat, int sys, double el, int type, const prcopt_t* opt) {
335 double a, b, c = 0.0, d = 0.0, fact = 1.0;
336 double sinel = sin(el);
337
338 (void)sat;
339
340 if (type == 2) { /* iono */
341 return 0.0;
342 } else if (type == 3) { /* trop */
343 return 0.0;
344 } else {
345 /* normal error model */
346 if (type == 1) {
347 fact *= opt->eratio[0];
348 }
349 fact *= sys == SYS_GLO
350 ? EFACT_GLO
351 : (sys == SYS_GAL
352 ? EFACT_GAL
353 : (sys == SYS_QZS
354 ? EFACT_QZS
355 : (sys == SYS_CMP
356 ? EFACT_CMP
357 : (sys == SYS_IRN ? EFACT_IRN : (sys == SYS_SBS ? EFACT_SBS : EFACT_GPS)))));
358 if (opt->ionoopt == IONOOPT_IFLC) {
359 fact *= 3.0;
360 }
361 a = fact * opt->err[1];
362 b = fact * opt->err[2];
363 if (opt->ionoopt == IONOOPT_EST || opt->ionoopt == IONOOPT_EST_ADPT) {
364 c = opt->stats_erriono;
365 }
366 if (opt->tropopt >= TROPOPT_EST) {
367 d = opt->stats_errtrop;
368 }
369 }
370 return a * a + b * b / sinel / sinel + c * c + d * d;
371}
372
373/*============================================================================
374 * State initialization
375 *===========================================================================*/
376
377/**
378 * @brief Initialize a state and its covariance.
379 * @param[in,out] rtk RTK control struct
380 * @param[in] xi Initial state value
381 * @param[in] var Initial variance
382 * @param[in] i State index
383 */
384static void initx(rtk_t* rtk, double xi, double var, int i) {
385 int j;
386 rtk->x[i] = xi;
387 for (j = 0; j < rtk->nx; j++) {
388 rtk->P[i + j * rtk->nx] = rtk->P[j + i * rtk->nx] = i == j ? var : 0.0;
389 }
390}
391
392/*============================================================================
393 * Position time update
394 *===========================================================================*/
395
396/**
397 * @brief Temporal update of position states.
398 * @param[in,out] rtk RTK control struct
399 * @param[in] obs Observation data
400 * @param[in] n Number of observations
401 * @param[in] tt Time step (s)
402 */
403static void udpos_ppp(rtk_t* rtk, const obsd_t* obs, int n, double tt) {
404 double *F, *FP, *xp, pos[3], Q[9] = {0}, Qv[9], var = 0.0;
405 double *Fx, *Pxx, *Pxi, *Pxb, *Pxi_, *Pxb_, *xp_;
406 double prnpos_h, prnpos_v;
407 double ki;
408 int i, j, nx, na, ni, nb;
409
410 (void)obs;
411 (void)n;
412
413 trace(NULL, 3, "udpos_ppp:\n");
414
415 nx = rtk->nx;
416
417 /* fixed mode */
418 if (rtk->opt.mode == PMODE_PPP_FIXED) {
419 for (i = 0; i < 3; i++) {
420 initx(rtk, rtk->opt.ru[i], 1E-8, i);
421 }
422 return;
423 }
424
425 /* initialize position for first epoch */
426 if (norm(rtk->x, 3) <= 0.0) {
427 for (i = 0; i < 3; i++) {
428 initx(rtk, rtk->sol.rr[i], VAR_POS, i);
429 }
430 if (rtk->opt.dynamics) {
431 for (i = 3; i < 6; i++) {
432 initx(rtk, rtk->sol.rr[i], VAR_VEL, i);
433 }
434 for (i = 6; i < 9; i++) {
435 initx(rtk, 1E-6, VAR_ACC, i);
436 }
437 }
438 return;
439 }
440
441 /* check variance of estimated position */
442 for (i = 0; i < 3; i++) {
443 var += rtk->P[i + i * nx];
444 }
445 var /= 3.0;
446
447 if (var > VAR_POS) {
448 /* reset position with large variance */
449 for (i = 0; i < 3; i++) {
450 initx(rtk, rtk->sol.rr[i], VAR_POS, i);
451 }
452 for (i = 3; i < 6; i++) {
453 initx(rtk, rtk->sol.rr[i], VAR_VEL, i);
454 }
455 for (i = 6; i < 9; i++) {
456 initx(rtk, 1E-6, VAR_ACC, i);
457 }
458 trace(NULL, 2, "reset rtk position due to large variance: var=%.3f\n", var);
459 return;
460 }
461
462 /* state transition of position/velocity/acceleration. For PPP-RTK iono
463 * states, keep the same joint pos/iono/bias covariance propagation as
464 * claslib: pos evolves by F, iono decays by ki, and cross-covariances are
465 * transformed consistently. */
466 if (NI_RTK(&rtk->opt) > 0) {
467 F = eye(nx);
468
469 if (rtk->opt.dynamics) {
470 for (i = 0; i < 6; i++) {
471 F[i + (i + 3) * nx] = tt;
472 }
473 for (i = 0; i < 3; i++) {
474 F[i + (i + 6) * nx] = SQR(tt) / 2.0;
475 }
476 }
477
478 na = NP_RTK(&rtk->opt);
479 ni = NI_RTK(&rtk->opt);
480 nb = nx - na - ni;
481 ki = rtk->opt.stats_tconstiono > 0.0 ? exp(-tt / rtk->opt.stats_tconstiono) : 1.0;
482
483 Fx = mat(na, na);
484 xp = mat(na, 1);
485 xp_ = mat(na, 1);
486 FP = mat(na, na);
487 Pxx = mat(na, na);
488 Pxi = mat(na, ni);
489 Pxi_ = mat(na, ni);
490 Pxb = mat(na, nb);
491 Pxb_ = mat(na, nb);
492
493 for (i = 0; i < na; i++) {
494 for (j = 0; j < na; j++) {
495 Pxx[i + j * na] = rtk->P[i + j * nx];
496 }
497 }
498 for (i = 0; i < na; i++) {
499 for (j = 0; j < ni; j++) {
500 Pxi_[i + j * na] = 0.5 * (rtk->P[i + (na + j) * nx] + rtk->P[na + j + i * nx]);
501 }
502 }
503 for (i = 0; i < na; i++) {
504 for (j = 0; j < nb; j++) {
505 Pxb_[i + j * na] = 0.5 * (rtk->P[i + (na + ni + j) * nx] + rtk->P[na + ni + j + i * nx]);
506 }
507 }
508 for (i = 0; i < na; i++) {
509 for (j = 0; j < na; j++) {
510 Fx[i + j * na] = F[i + j * nx];
511 }
512 xp_[i] = rtk->x[i];
513 }
514
515 matmul("NN", na, 1, na, 1.0, Fx, xp_, 0.0, xp);
516 for (i = 0; i < na; i++) {
517 rtk->x[i] = xp[i];
518 }
519 for (i = 0; i < ni; i++) {
520 rtk->x[na + i] *= ki;
521 }
522
523 matmul("NN", na, na, na, 1.0, Fx, Pxx, 0.0, FP);
524 matmul("NT", na, na, na, 1.0, FP, Fx, 0.0, Pxx);
525 for (i = 0; i < na; i++) {
526 for (j = 0; j < na; j++) {
527 rtk->P[i + j * nx] = Pxx[i + j * na];
528 }
529 }
530
531 matmul("NN", na, ni, na, ki, Fx, Pxi_, 0.0, Pxi);
532 for (i = 0; i < na; i++) {
533 for (j = 0; j < ni; j++) {
534 rtk->P[i + (na + j) * nx] = rtk->P[(na + j) + i * nx] = Pxi[i + j * na];
535 }
536 }
537
538 matmul("NN", na, nb, na, 1.0, Fx, Pxb_, 0.0, Pxb);
539 for (i = 0; i < na; i++) {
540 for (j = 0; j < nb; j++) {
541 rtk->P[i + (na + ni + j) * nx] = rtk->P[(na + ni + j) + i * nx] = Pxb[i + j * na];
542 }
543 }
544
545 for (i = 0; i < ni; i++) {
546 for (j = 0; j < ni; j++) {
547 rtk->P[(na + i) + (na + j) * nx] *= ki * ki;
548 }
549 }
550 for (i = 0; i < ni; i++) {
551 for (j = 0; j < nb; j++) {
552 rtk->P[(na + i) + (na + ni + j) * nx] *= ki;
553 }
554 }
555 for (i = 0; i < nb; i++) {
556 for (j = 0; j < ni; j++) {
557 rtk->P[(na + ni + i) + (na + j) * nx] *= ki;
558 }
559 }
560
561 free(Pxx);
562 free(Pxi);
563 free(Pxb);
564 free(Pxi_);
565 free(Pxb_);
566 free(Fx);
567 free(xp_);
568 free(F);
569 free(FP);
570 free(xp);
571 } else if (rtk->opt.dynamics) {
572 F = eye(9);
573 xp = mat(9, 1);
574 FP = mat(nx, nx);
575
576 for (i = 0; i < 6; i++) {
577 F[i + (i + 3) * 9] = tt;
578 }
579 for (i = 0; i < 3; i++) {
580 F[i + (i + 6) * 9] = SQR(tt) / 2.0;
581 }
582
583 matmul("NN", 9, 1, 9, 1.0, F, rtk->x, 0.0, xp);
584 matcpy(rtk->x, xp, 9, 1);
585
586 matcpy(FP, rtk->P, nx, nx);
587 for (j = 0; j < nx; j++) {
588 for (i = 0; i < 6; i++) {
589 FP[i + j * nx] += rtk->P[i + 3 + j * nx] * tt;
590 }
591 }
592
593 matcpy(rtk->P, FP, nx, nx);
594 for (j = 0; j < nx; j++) {
595 for (i = 0; i < 6; i++) {
596 rtk->P[j + i * nx] += rtk->P[j + (i + 3) * nx] * tt;
597 }
598 }
599
600 free(F);
601 free(FP);
602 free(xp);
603 }
604
605 if (rtk->opt.dynamics) {
606 /* process noise added to acceleration */
607 Q[0] = Q[4] = SQR(rtk->opt.prn[3]) * fabs(tt);
608 Q[8] = SQR(rtk->opt.prn[4]) * fabs(tt);
609 ecef2pos(rtk->x, pos);
610 covecef(pos, Q, Qv);
611 for (i = 0; i < 3; i++) {
612 for (j = 0; j < 3; j++) {
613 rtk->P[i + 6 + (j + 6) * nx] += Qv[i + j * 3];
614 }
615 }
616 } else {
617 /* process noise added to position */
618 /* Do not fall back to prn[6]: PPP uses that slot for IFB process
619 * noise. CLAS position process noise has dedicated fields. */
620 prnpos_h = rtk->opt.stats_prnposith;
621 prnpos_v = rtk->opt.stats_prnpositv;
622 Q[0] = Q[4] = SQR(prnpos_h) * fabs(tt);
623 Q[8] = SQR(prnpos_v) * fabs(tt);
624 ecef2pos(rtk->x, pos);
625 covecef(pos, Q, Qv);
626 for (i = 0; i < 3; i++) {
627 for (j = 0; j < 3; j++) {
628 rtk->P[i + j * nx] += Qv[i + j * 3];
629 }
630 }
631 }
632}
633
634/*============================================================================
635 * Troposphere time update
636 *===========================================================================*/
637
638/**
639 * @brief Temporal update of tropospheric parameters.
640 * @param[in,out] rtk RTK control struct
641 * @param[in] tt Time step (s)
642 */
643static void udtrop(rtk_t* rtk, double tt) {
644 int j, k;
645
646 trace(NULL, 3, "udtrop : tt=%.1f\n", tt);
647
648 j = IT_RTK(&rtk->opt);
649
650 if (rtk->x[j] == 0.0) {
651 initx(rtk, INIT_ZWD, SQR(rtk->opt.std[2]), j); /* initial ZWD */
652 if (rtk->opt.tropopt >= TROPOPT_ESTG) {
653 for (k = 0; k < 2; k++) {
654 initx(rtk, 1E-6, VAR_GRA, ++j);
655 }
656 }
657 } else {
658 rtk->P[j + j * rtk->nx] += SQR(rtk->opt.prn[2]) * tt;
659 if (rtk->opt.tropopt >= TROPOPT_ESTG) {
660 for (k = 0; k < 2; k++) {
661 rtk->P[++j * (1 + rtk->nx)] += SQR(rtk->opt.prn[2] * 0.3) * fabs(rtk->tt);
662 }
663 }
664 }
665}
666
667/*============================================================================
668 * Ionosphere time update
669 *===========================================================================*/
670
671/**
672 * @brief Temporal update of ionospheric parameters.
673 * @param[in,out] rtk RTK control struct
674 * @param[in] tt Time step (s)
675 * @param[in] bl Baseline length (m)
676 * @param[in] obs Observation data
677 * @param[in] ns Number of observations
678 */
679static void udion(rtk_t* rtk, double tt, double bl, const obsd_t* obs, int ns) {
680 double el, fact, qi;
681 int i, j, k, f, m;
682
683 (void)bl;
684
685 trace(NULL, 3, "udion : tt=%.1f ns=%d\n", tt, ns);
686
687 for (i = 1; i <= MAXSAT; i++) {
688 j = II_RTK(i, &rtk->opt);
689 if (rtk->x[j] == 0.0) {
690 continue;
691 }
692 for (k = 0; k < ns; k++) {
693 if (obs[k].sat == i) {
694 break;
695 }
696 }
697 if (ns == k) {
698 continue;
699 }
700 m = clas_osr_selfreqpair(i, &rtk->opt, obs + k);
701 for (f = 0; f < rtk->opt.nf; f++) {
702 if (f > 0 && !(f & m)) {
703 continue;
704 }
705 if ((rtk->ssat[i - 1].outc[f] > (unsigned int)rtk->opt.maxout)) {
706 trace(NULL, 3, "reset iono estimation, sat=%2d outc=%6d\n", i, rtk->ssat[i - 1].outc[f]);
707 rtk->x[j] = 0.0;
708 }
709 }
710 }
711 for (i = 0; i < ns; i++) {
712 j = II_RTK(obs[i].sat, &rtk->opt);
713 if (rtk->x[j] == 0.0) {
714 initx(rtk, 1E-6, SQR(rtk->opt.std[1]), j);
715 rtk->Q[j + j * rtk->nx] = 0.0;
716 } else {
717 /* elevation dependent factor of process noise */
718 el = rtk->ssat[obs[i].sat - 1].azel[1];
719 fact = cos(el);
720 if (rtk->opt.ionoopt == IONOOPT_EST_ADPT) {
721 /* adaptive: clamp Q within [prn[1]^2, prnionomax^2] */
722 if (rtk->Q[j + j * rtk->nx] == 0.0) {
723 rtk->Q[j + j * rtk->nx] = SQR(rtk->opt.prn[1] * fact);
724 } else {
725 if (rtk->Q[j + j * rtk->nx] > SQR(rtk->opt.prnionomax)) {
726 rtk->Q[j + j * rtk->nx] = SQR(rtk->opt.prnionomax);
727 } else if (rtk->Q[j + j * rtk->nx] < SQR(rtk->opt.prn[1])) {
728 rtk->Q[j + j * rtk->nx] = SQR(rtk->opt.prn[1]);
729 }
730 }
731 qi = rtk->Q[j + j * rtk->nx];
732 } else {
733 qi = SQR(rtk->opt.prn[1] * fact);
734 rtk->Q[j + j * rtk->nx] = qi;
735 }
736 rtk->P[j + j * rtk->nx] += qi * tt;
737 }
738 }
739}
740
741/*============================================================================
742 * Cycle slip detection: LLI
743 *===========================================================================*/
744
745/**
746 * @brief Detect cycle slip by LLI flag.
747 * @param[in,out] rtk RTK control struct
748 * @param[in] obs Observation data
749 * @param[in] n Number of observations
750 */
751static void detslp_ll(rtk_t* rtk, const obsd_t* obs, int n) {
752 int i, j;
753
754 for (i = 0; i < n && i < MAXOBS; i++) {
755 for (j = 0; j < rtk->opt.nf; j++) {
756 if (obs[i].L[j] == 0.0 || !(obs[i].LLI[j] & 3)) {
757 continue;
758 }
759 trace(NULL, 2, "detslp_ll: slip detected sat=%2d f=%d\n", obs[i].sat, j + 1);
760 rtk->ssat[obs[i].sat - 1].slip[j] = 1;
761 }
762 }
763}
764
765/*============================================================================
766 * Geometry-free phase measurements
767 *===========================================================================*/
768
769/**
770 * @brief L1/L2 geometry-free phase measurement.
771 */
772static double gfmeas_L1L2(const obsd_t* obs, const nav_t* nav) {
773 double lam0 = sat_lambda(obs->sat, 0, nav);
774 double lam1 = sat_lambda(obs->sat, 1, nav);
775 if (lam0 == 0.0 || lam1 == 0.0 || obs->L[0] == 0.0 || obs->L[1] == 0.0) {
776 return 0.0;
777 }
778 return lam0 * obs->L[0] - lam1 * obs->L[1];
779}
780
781/**
782 * @brief L1/L5 geometry-free phase measurement.
783 */
784static double gfmeas_L1L5(const obsd_t* obs, const nav_t* nav) {
785 double lam0 = sat_lambda(obs->sat, 0, nav);
786 double lam2 = sat_lambda(obs->sat, 2, nav);
787 if (lam0 == 0.0 || lam2 == 0.0 || obs->L[0] == 0.0 || obs->L[2] == 0.0) {
788 return 0.0;
789 }
790 return lam0 * obs->L[0] - lam2 * obs->L[2];
791}
792
793/*============================================================================
794 * Cycle slip detection: geometry-free
795 *===========================================================================*/
796
797/**
798 * @brief Detect cycle slip by geometry-free phase jump.
799 * @param[in,out] rtk RTK control struct
800 * @param[in] obs Observation data
801 * @param[in] n Number of observations
802 * @param[in] nav Navigation data
803 */
804static void detslp_gf(rtk_t* rtk, const obsd_t* obs, int n, const nav_t* nav) {
805 double g0, g1;
806 int i;
807
808 trace(NULL, 3, "detslp_gf_L1_L2\n");
809
810 /* slip detection for L1-L2 GF */
811 for (i = 0; i < n && i < MAXOBS; i++) {
812 if ((g1 = gfmeas_L1L2(obs + i, nav)) == 0.0) {
813 continue;
814 }
815 g0 = rtk->ssat[obs[i].sat - 1].gf[0]; /* gf[0] = L1-L2 */
816 rtk->ssat[obs[i].sat - 1].gf[0] = g1;
817 if (g0 != 0.0 && fabs(g1 - g0) > rtk->opt.thresslip) {
818 trace(NULL, 2, "detslp_gf_L1_L2: slip detected sat=%2d gf=%8.3f->%8.3f\n", obs[i].sat, g0, g1);
819 rtk->ssat[obs[i].sat - 1].slip[0] |= 1;
820 rtk->ssat[obs[i].sat - 1].slip[1] |= 1;
821 }
822 }
823
824 trace(NULL, 3, "detslp_gf_L1_L5\n");
825
826 /* slip detection for L1-L5 GF */
827 if (rtk->opt.nf <= 2) {
828 return;
829 }
830 for (i = 0; i < n && i < MAXOBS; i++) {
831 if ((g1 = gfmeas_L1L5(obs + i, nav)) == 0.0) {
832 continue;
833 }
834 g0 = rtk->ssat[obs[i].sat - 1].gf[1]; /* gf[1] = L1-L5 */
835 rtk->ssat[obs[i].sat - 1].gf[1] = g1;
836 if (g0 != 0.0 && fabs(g1 - g0) > rtk->opt.thresslip) {
837 trace(NULL, 2, "detslp_gf_L1_L5: slip detected sat=%2d gf=%8.3f->%8.3f\n", obs[i].sat, g0, g1);
838 rtk->ssat[obs[i].sat - 1].slip[0] |= 1;
839 rtk->ssat[obs[i].sat - 1].slip[2] |= 1;
840 }
841 }
842}
843
844/*============================================================================
845 * Phase bias time update
846 *===========================================================================*/
847
848/**
849 * @brief Temporal update of phase biases.
850 * @param[in,out] rtk RTK control struct
851 * @param[in] tt Time step (s)
852 * @param[in] obs Observation data
853 * @param[in] n Number of observations
854 * @param[in] nav Navigation data
855 */
856static void udbias_ppp(rtk_t* rtk, double tt, const obsd_t* obs, int n, const nav_t* nav) {
857 double offset, lami, cp, pr, com_bias;
858 double* bias;
859 int i, j, f, slip, nf = NF_RTK(&rtk->opt), reset, isat;
860
861 trace(NULL, 3, "udbias : tt=%.1f n=%d\n", tt, n);
862
863 for (i = 0; i < MAXSAT; i++) {
864 for (j = 0; j < rtk->opt.nf; j++) {
865 rtk->ssat[i].slip[j] = 0;
866 rtk->ssat[i].pbreset[j] = 0;
867 }
868 }
869
870 /* detect cycle slip by LLI */
871 detslp_ll(rtk, obs, n);
872
873 /* detect cycle slip by geometry-free phase jump */
874 detslp_gf(rtk, obs, n, nav);
875
876 for (f = 0; f < nf; f++) {
877 /* reset phase-bias if expire obs outage counter */
878 for (i = 0; i < MAXSAT; i++) {
879 reset = ++rtk->ssat[i].outc[f] > (unsigned int)rtk->opt.maxout;
880
881 if (rtk->opt.modear == ARMODE_INST && rtk->x[IB_RTK(i + 1, f, &rtk->opt)] != 0.0) {
882 initx(rtk, 0.0, 0.0, IB_RTK(i + 1, f, &rtk->opt));
883 } else if (reset && rtk->x[IB_RTK(i + 1, f, &rtk->opt)] != 0.0) {
884 initx(rtk, 0.0, 0.0, IB_RTK(i + 1, f, &rtk->opt));
885 trace(NULL, 3, "udbias : obs outage counter overflow (sat=%3d L%d n=%5d)\n", i + 1, f + 1,
886 rtk->ssat[i].outc[f]);
887 }
888 if (rtk->opt.modear != ARMODE_INST && reset) {
889 rtk->ssat[i].lock[f] = -rtk->opt.minlock;
890 }
891 }
892
893 /* reset phase-bias if detecting cycle slip */
894 for (i = 0; i < n && i < MAXOBS; i++) {
895 j = IB_RTK(obs[i].sat, f, &rtk->opt);
896 rtk->P[j + j * rtk->nx] += rtk->opt.prn[0] * rtk->opt.prn[0] * tt;
897 slip = rtk->ssat[obs[i].sat - 1].slip[f];
898 if (rtk->opt.ionoopt == IONOOPT_IFLC) {
899 slip |= rtk->ssat[obs[i].sat - 1].slip[1];
900 }
901 if (rtk->opt.modear == ARMODE_INST || !(slip & 1)) {
902 continue;
903 }
904 if (!(slip & 1)) {
905 continue;
906 }
907 rtk->x[j] = 0.0;
908 rtk->ssat[obs[i].sat - 1].lock[f] = -rtk->opt.minlock;
909 }
910
911 bias = zeros(n, 1);
912
913 /* estimate approximate phase-bias by phase - code */
914 for (i = j = 0, offset = 0.0; i < n; i++) {
915 cp = obs[i].L[f];
916 pr = obs[i].P[f];
917 lami = sat_lambda(obs[i].sat, f, nav);
918 if (cp == 0.0 || pr == 0.0 || lami == 0.0) {
919 continue;
920 }
921
922 bias[i] = cp * lami - pr;
923
924 /* offset = sum of (bias - phase-bias) for all valid sats in meters */
925 if (rtk->x[IB_RTK(obs[i].sat, f, &rtk->opt)] != 0.0) {
926 lami = sat_lambda(obs[i].sat, f, nav);
927 offset += bias[i] - rtk->x[IB_RTK(obs[i].sat, f, &rtk->opt)] * lami;
928 j++;
929 }
930 }
931
932 com_bias = j > 0 ? offset / j : 0.0;
933
934 /* set initial states of phase-bias */
935 for (i = 0; i < n; i++) {
936 isat = obs[i].sat;
937 if (bias[i] == 0.0 || rtk->x[IB_RTK(isat, f, &rtk->opt)] != 0.0) {
938 continue;
939 }
940 trace(NULL, 3, "set initial states of phase-bias, sat=%2d f=%2d\n", obs[i].sat, f + 1);
941 lami = sat_lambda(isat, f, nav);
942 if (lami == 0.0) {
943 continue;
944 }
945 initx(rtk, (bias[i] - com_bias) / lami, SQR(rtk->opt.std[0]), IB_RTK(obs[i].sat, f, &rtk->opt));
946 rtk->ssat[isat - 1].lock[f] = -rtk->opt.minlock;
947 /* suppress pbreset after floatcnt mass reset: not an actual
948 * phase bias change, so compensatedisp should continue normally */
949 if (!prtk_ctx.floatcnt_reset) {
950 rtk->ssat[isat - 1].pbreset[f] = 1;
951 }
952 }
953
954 free(bias);
955 }
956 /* clear floatcnt_reset flag after bias update */
957 prtk_ctx.floatcnt_reset = 0;
958}
959
960/*============================================================================
961 * State update dispatcher
962 *===========================================================================*/
963
964/**
965 * @brief Temporal update of all PPP-RTK states.
966 * @param[in,out] rtk RTK control struct
967 * @param[in] obs Observation data
968 * @param[in] n Number of observations
969 * @param[in] nav Navigation data
970 */
971static void udstate_ppp(rtk_t* rtk, const obsd_t* obs, int n, const nav_t* nav) {
972 double tt = fabs(rtk->tt);
973
974 trace(NULL, 3, "udstate_ppp: n=%d tt=%f\n", n, tt);
975
976 /* temporal update of position */
977 udpos_ppp(rtk, obs, n, tt);
978
979 /* temporal update of tropospheric parameters */
980 if (rtk->opt.tropopt >= TROPOPT_EST) {
981 udtrop(rtk, tt);
982 }
983
984 /* temporal update of phase-bias */
985 udbias_ppp(rtk, tt, obs, n, nav);
986
987 /* temporal update of ionospheric parameters */
988 if (rtk->opt.ionoopt >= IONOOPT_EST) {
989 udion(rtk, tt, BASELEN, obs, n);
990 }
991}
992
993/*============================================================================
994 * Precise troposphere model for KF estimation
995 *===========================================================================*/
996
997/**
998 * @brief Compute precise tropospheric delay for KF.
999 * @param[in] time Observation time
1000 * @param[in] pos Receiver position (lat/lon/hgt)
1001 * @param[in] r Unused (receiver index, always 0)
1002 * @param[in] azel Azimuth/elevation (rad)
1003 * @param[in] opt Processing options
1004 * @param[in] x State vector
1005 * @param[out] dtdx Partial derivatives [3]
1006 * @return Tropospheric delay (m)
1007 */
1008static double prectrop2(gtime_t time, const double* pos, int r, const double* azel, const prcopt_t* opt,
1009 const double* x, double* dtdx) {
1010 double m_w = 0.0, cotz, grad_n, grad_e;
1011 int i = IT_RTK(opt);
1012
1013 (void)r;
1014
1015 /* wet mapping function */
1016 tropmapf(time, pos, azel, &m_w);
1017
1018 if (opt->tropopt >= TROPOPT_ESTG && azel[1] > 0.0) {
1019 /* m_w=m_0+m_0*cot(el)*(Gn*cos(az)+Ge*sin(az)): ref [6] */
1020 cotz = 1.0 / tan(azel[1]);
1021 grad_n = m_w * cotz * cos(azel[0]);
1022 grad_e = m_w * cotz * sin(azel[0]);
1023 m_w += grad_n * x[i + 1] + grad_e * x[i + 2];
1024 dtdx[1] = grad_n * x[i];
1025 dtdx[2] = grad_e * x[i];
1026 } else {
1027 dtdx[1] = dtdx[2] = 0.0;
1028 }
1029 dtdx[0] = m_w;
1030 return m_w * x[i];
1031}
1032
1033/*============================================================================
1034 * Validity check for ZD observation
1035 *===========================================================================*/
1036
1037/**
1038 * @brief Check if zero-difference observation is valid.
1039 * @param[in] i Observation index
1040 * @param[in] f Frequency index (0..nf*2-1)
1041 * @param[in] nf Number of frequencies
1042 * @param[in] y ZD residual array
1043 * @return 1 if valid, 0 otherwise
1044 */
1045static int validobs(int i, int f, int nf, double* y) {
1046 return y[f + i * nf * 2] != 0.0 && (f < nf || (y[f - nf + i * nf * 2] != 0.0));
1047}
1048
1049/*============================================================================
1050 * DD covariance matrix
1051 *===========================================================================*/
1052
1053/**
1054 * @brief Build DD covariance from SD variances.
1055 * @param[in] nb Block sizes array
1056 * @param[in] n Number of blocks
1057 * @param[in] Ri Reference satellite variances
1058 * @param[in] Rj Non-reference satellite variances
1059 * @param[in] nv Total number of DD observations
1060 * @param[out] R DD covariance matrix (nv x nv)
1061 */
1062static void ddcov(const int* nb, int n, const double* Ri, const double* Rj, int nv, double* R) {
1063 int i, j, k = 0, b;
1064
1065 for (i = 0; i < nv * nv; i++) {
1066 R[i] = 0.0;
1067 }
1068 for (b = 0; b < n; k += nb[b++]) {
1069 for (i = 0; i < nb[b]; i++) {
1070 for (j = 0; j < nb[b]; j++) {
1071 R[k + i + (k + j) * nv] = Ri[k + i] + ((i == j) ? Rj[k + i] : 0.0);
1072 }
1073 }
1074 }
1075}
1076
1077/*============================================================================
1078 * Simplified system test
1079 *===========================================================================*/
1080
1081/**
1082 * @brief Test if satellite system matches group index.
1083 * @param[in] sys Satellite system (SYS_*)
1084 * @param[in] m Group index (0:GPS/SBS, 1:GLO, 2:GAL, 3:BDS, 4:QZS, 5:GPS-L2C)
1085 * @param[in] qzsmodear QZSS AR mode (0:off, 1:on, 2:gps-qzs)
1086 * @param[in] code Observation code indicator (CODE_*)
1087 * @return 1 if match, 0 otherwise
1088 */
1089static int test_sys(int sys, int m, int qzsmodear, int code) {
1090 switch (sys) {
1091 case SYS_GPS:
1092 return m == (code == CODE_L2X ? 5 : 0);
1093 case SYS_QZS:
1094 return m == (qzsmodear == 2 ? 0 : 4);
1095 case SYS_SBS:
1096 return m == 0;
1097 case SYS_GLO:
1098 return m == 1;
1099 case SYS_GAL:
1100 return m == 2;
1101 case SYS_CMP:
1102 return m == 3;
1103 }
1104 return 0;
1105}
1106
1107/*============================================================================
1108 * Residual test (chi-squared outlier rejection)
1109 *===========================================================================*/
1110
1111/**
1112 * @brief Validate residuals and reject outliers.
1113 * @param[in,out] rtk RTK control struct
1114 * @param[in] vflg Residual flags (sat/type/freq encoded)
1115 * @param[in] v Residual vector
1116 * @param[in] Q Innovation covariance matrix (m x m)
1117 * @param[out] iv Indices of valid residuals
1118 * @param[in] m Total number of residuals
1119 * @param[out] ns Count per frequency/type
1120 * @param[in] flg Stage flag (0:prefit, 1:postfit float, 2:postfit fix)
1121 * @return Number of valid residuals
1122 */
1123static int residual_test(rtk_t* rtk, const int* vflg, const double* v, const double* Q, int* iv, const int m, int* ns,
1124 const int flg) {
1125 int i, k, sat2, type, freq, nf = rtk->opt.nf;
1126 double inno0, inno1, inno2, gamma;
1127 double *resc, *resp, sig2, vv0_, vvf_, vv = 0.0;
1128 int cnt, cnt2;
1129 char* stype;
1130
1131 /* per-component rejection thresholds (rejionno1-3) or fallback to maxinno */
1132 inno0 = rtk->opt.maxinno_ext[0] > 0.0 ? rtk->opt.maxinno_ext[0] : rtk->opt.maxinno;
1133 inno1 = rtk->opt.maxinno_ext[1] > 0.0 ? rtk->opt.maxinno_ext[1] : rtk->opt.maxinno;
1134 inno2 = rtk->opt.maxinno_ext[2] > 0.0 ? rtk->opt.maxinno_ext[2] : rtk->opt.maxinno;
1135
1136 resc = zeros(MAXSAT * nf, 1);
1137 resp = zeros(MAXSAT * nf, 1);
1138
1139 /* collect carrier phase residuals */
1140 for (i = cnt2 = 0; i < m; i++) {
1141 sat2 = (vflg[i] >> 8) & 0xFF;
1142 type = (vflg[i] >> 4) & 0xF;
1143 freq = (vflg[i]) & 0xF;
1144 if (type == 0) {
1145 resc[(sat2 - 1) * nf + freq] = v[i];
1146 resp[(sat2 - 1) * nf + freq] = Q[i + i * m];
1147 cnt2++;
1148 }
1149 }
1150
1151 /* multi-frequency dispersive/non-dispersive validation */
1152 for (i = 0; i < MAXSAT; i++) {
1153 for (k = 1; k < nf; k++) {
1154 int qi = 0, qj = k;
1155 double q0, q1;
1156
1157 if (resp[i * nf + qi] == 0.0 || resp[i * nf + qj] == 0.0) {
1158 continue;
1159 }
1160 gamma = SQR(lam_carr[qj]) / SQR(lam_carr[qi]);
1161 q0 = resp[i * nf + qi];
1162 q1 = resp[i * nf + qj];
1163 vvf_ = FREQ1 / FREQ2 * (resc[i * nf + qi] - resc[i * nf + qj]) / (1.0 - gamma);
1164 vv0_ = (gamma * resc[i * nf + qi] - resc[i * nf + qj]) / (gamma - 1.0);
1165 sig2 = q0 > q1 ? q0 : q1;
1166
1167 /* dispersive term validation */
1168 if ((inno1 > 0.0) && (vvf_ * vvf_ > inno1 * inno1 * sig2)) {
1169 rtk->ssat[i].vsat[qi] = rtk->ssat[i].vsat[qj] = 0;
1170 rtk->ssat[i].rejc[qi]++;
1171 rtk->ssat[i].rejc[qj]++;
1172 trace(NULL, 2,
1173 "Residuals validation:freq outlier detected "
1174 "(stage=%1d sat=%2d f=%d %d v=%6.3f sig=%.3f)\n",
1175 flg, i + 1, qi, qj, vvf_, sqrt(sig2));
1176 continue;
1177 }
1178 /* non-dispersive (L0) term validation */
1179 if ((inno2 > 0.0) && (vv0_ * vv0_ > inno2 * inno2 * sig2)) {
1180 rtk->ssat[i].vsat[qi] = rtk->ssat[i].vsat[qj] = 0;
1181 rtk->ssat[i].rejc[qi]++;
1182 rtk->ssat[i].rejc[qj]++;
1183 trace(NULL, 2,
1184 "Residuals validation:L0 outlier detected "
1185 "(stage=%1d sat=%2d v=%6.3f sig=%.3f)\n",
1186 flg, i + 1, vv0_, sqrt(sig2));
1187 continue;
1188 }
1189 }
1190 }
1191
1192 /* pre/post residuals test */
1193 for (i = k = cnt = 0; i < m; i++) {
1194 sat2 = (vflg[i] >> 8) & 0xFF;
1195 type = (vflg[i] >> 4) & 0xF;
1196 freq = (vflg[i]) & 0xF;
1197 if (type == 0) {
1198 stype = "L";
1199 } else if (type == 1) {
1200 stype = "P";
1201 } else {
1202 stype = "C";
1203 }
1204 (void)stype;
1205
1206 /* validate individual residuals */
1207 double qdiag = Q[i + i * m];
1208 if (!qdiag || inno0 == 0 || (v[i] * v[i] < qdiag * inno0 * inno0)) {
1209 if (type == 0) {
1210 if (!rtk->ssat[sat2 - 1].vsat[freq]) {
1211 continue;
1212 }
1213 vv += v[i] * v[i] / qdiag;
1214 cnt++;
1215 }
1216 ns[type < 1 ? freq : freq * type + nf]++;
1217 iv[k++] = i;
1218 } else {
1219 if (type == 0) {
1220 rtk->ssat[sat2 - 1].vsat[freq] = 0;
1221 }
1222 trace(NULL, 2,
1223 "meas update outlier detected "
1224 "(sat=%2d %s%d v=%6.3f sig=%6.3f outc=%2d)\n",
1225 sat2, stype, freq + 1, v[i], sqrt(qdiag), rtk->ssat[sat2 - 1].outc[freq]);
1226 }
1227 }
1228
1229 /* chi-square validation for carrier residuals */
1230 if (cnt > NP_RTK(&rtk->opt)) {
1231 if (vv > 1.0 * chisqr[cnt - NP_RTK(&rtk->opt) - 1]) {
1232 trace(NULL, 2,
1233 "L1/L2 chi-square validation failed "
1234 "(stage=%d cnt=%d np=%d vv=%.2f cs=%5.3f)\n",
1235 flg, cnt, NP_RTK(&rtk->opt), vv, chisqr[cnt - NP_RTK(&rtk->opt) - 1]);
1236 }
1237 prtk_ctx.chisq = vv / chisqr[cnt - NP_RTK(&rtk->opt) - 1];
1238 } else {
1239 prtk_ctx.chisq = ((double)cnt / (cnt2 > 0 ? cnt2 : 1)) < 0.5 ? 100.0 : 0.0;
1240 }
1241
1242 free(resc);
1243 free(resp);
1244
1245 return k;
1246}
1247
1248/*============================================================================
1249 * Inner Kalman filter update with outlier rejection
1250 *===========================================================================*/
1251
1252/**
1253 * @brief Inner KF measurement update with residual test.
1254 * @param[in,out] rtk RTK control struct
1255 * @param[in] x State vector (n)
1256 * @param[in] P Covariance matrix (n x n)
1257 * @param[in] H Observation matrix (n x m)
1258 * @param[in] v Residual vector (m)
1259 * @param[in] R Measurement noise (m x m)
1260 * @param[in] n Number of states
1261 * @param[in] m Number of measurements
1262 * @param[out] xp Updated state vector (n)
1263 * @param[out] Pp Updated covariance (n x n)
1264 * @param[in] vflg Residual flags
1265 * @param[in] flg Stage flag (0:prefit, 1+:postfit)
1266 * @return 0 on success, non-zero on error
1267 */
1268static int filter2_(rtk_t* rtk, const double* x, const double* P, const double* H, const double* v, const double* R,
1269 const int n, const int m, double* xp, double* Pp, double* Qp, const int* vflg, const int flg) {
1270 double *K, *v_, *Q_, *F_, *H_, *I, *F, *Q;
1271 int i, j, k, info = 0, *iv, nn[6];
1272
1273 Q = mat(m, m);
1274 F = mat(n, m);
1275 iv = imat(m, 1);
1276 memset(nn, 0, sizeof(nn));
1277 matcpy(Q, R, m, m);
1278 matcpy(xp, x, n, 1);
1279 matmul("NN", n, m, n, 1.0, P, H, 0.0, F); /* F=P*H */
1280 matmul("TN", m, m, n, 1.0, H, F, 1.0, Q); /* Q=H'*P*H+R */
1281
1282 /* prefit/postfit residuals test */
1283 if (!(k = residual_test(rtk, vflg, v, Q, iv, m, nn, flg))) {
1284 trace(NULL, 2, "filter status (all measurements rejected): stage=%1d\n", flg);
1285 } else {
1286 trace(NULL, 2,
1287 "filter status: stage=%1d refsat=%2d states=%2d obs=%2d "
1288 "vobs=%2d vsat(L1)=%2d chi-val=%5.3f\n",
1289 flg, (vflg[iv[0]] >> 16) & 0xFF, n, m, k, nn[0], prtk_ctx.chisq);
1290 }
1291
1292 if (flg > 0 || k == 0) {
1293 matcpy(Pp, P, n, n);
1294 free(Q);
1295 free(F);
1296 free(iv);
1297 return info;
1298 }
1299
1300 v_ = mat(k, 1);
1301 Q_ = mat(k, k);
1302 F_ = mat(n, k);
1303 H_ = mat(n, k);
1304 K = mat(n, k);
1305 for (i = 0; i < k; i++) {
1306 v_[i] = v[iv[i]];
1307 for (j = 0; j < k; j++) {
1308 Q_[i + j * k] = Q[iv[i] + iv[j] * m];
1309 }
1310 for (j = 0; j < n; j++) {
1311 F_[j + i * n] = F[j + iv[i] * n];
1312 }
1313 for (j = 0; j < n; j++) {
1314 H_[j + i * n] = H[j + iv[i] * n];
1315 }
1316 }
1317
1318 I = eye(n);
1319 if (!(info = matinv(Q_, k))) {
1320 matmul("NN", n, k, k, 1.0, F_, Q_, 0.0, K); /* K=P*H*Q^-1 */
1321 matmul("NN", n, 1, k, 1.0, K, v_, 1.0, xp); /* xp=x+K*v */
1322 matmul("NT", n, n, k, -1.0, K, H_, 1.0, I); /* I=I-K*H' */
1323 matmul("NN", n, n, n, 1.0, I, P, 0.0, Pp); /* Pp=(I-K*H')*P */
1324
1325 /* innovation covariance Qp = K*v * (K*v)^T for adaptive noise, in the
1326 compacted (k) state space; filter2() scatters it back to full nx. */
1327 if (rtk->opt.ionoopt == IONOOPT_EST_ADPT || rtk->opt.prnadpt) {
1328 double* vd = mat(n, 1);
1329 matmul("NN", n, 1, k, 1.0, K, v_, 0.0, vd); /* vd=K*v */
1330 matmul("NT", n, n, 1, 1.0, vd, vd, 0.0, Qp); /* Qp=vd*vd' */
1331 free(vd);
1332 }
1333 }
1334 free(I);
1335
1336 free(iv);
1337 free(v_);
1338 free(Q_);
1339 free(F_);
1340 free(H_);
1341 free(K);
1342 free(Q);
1343 free(F);
1344
1345 return info;
1346}
1347
1348/*============================================================================
1349 * Outer KF wrapper (state compaction)
1350 *===========================================================================*/
1351
1352/**
1353 * @brief KF measurement update with state compaction.
1354 *
1355 * Compacts the state vector to exclude zero/invalid states before calling
1356 * the inner filter, then maps results back.
1357 *
1358 * @param[in,out] rtk RTK control struct
1359 * @param[in,out] x State vector (n)
1360 * @param[in,out] P Covariance matrix (n x n)
1361 * @param[in] H Observation matrix (n x m)
1362 * @param[in] v Residual vector (m)
1363 * @param[in] R Measurement noise (m x m)
1364 * @param[in] n Number of states
1365 * @param[in] m Number of measurements
1366 * @param[in] vflg Residual flags
1367 * @param[in] flg Stage flag
1368 * @return 0 on success, non-zero on error
1369 */
1370static int filter2(rtk_t* rtk, double* x, double* P, const double* H, const double* v, const double* R, const int n,
1371 const int m, const int* vflg, const int flg) {
1372 double *x_, *xp_, *P_, *Pp_, *Qp_, *H_;
1373 int i, j, k, info, *ix;
1374 int adapt = (flg == 0 && (rtk->opt.ionoopt == IONOOPT_EST_ADPT || rtk->opt.prnadpt));
1375
1376 ix = imat(n, 1);
1377 for (i = k = 0; i < n; i++) {
1378 if (x[i] != 0.0 && P[i + i * n] > 0.0) {
1379 ix[k++] = i;
1380 }
1381 }
1382 x_ = mat(k, 1);
1383 xp_ = mat(k, 1);
1384 P_ = mat(k, k);
1385 Pp_ = mat(k, k);
1386 Qp_ = mat(k, k);
1387 H_ = mat(k, m);
1388
1389 for (i = 0; i < k; i++) {
1390 x_[i] = x[ix[i]];
1391 for (j = 0; j < k; j++) {
1392 P_[i + j * k] = P[ix[i] + ix[j] * n];
1393 }
1394 for (j = 0; j < m; j++) {
1395 H_[i + j * k] = H[ix[i] + j * n];
1396 }
1397 }
1398
1399 info = filter2_(rtk, x_, P_, H_, v, R, k, m, xp_, Pp_, Qp_, vflg, flg);
1400
1401 /* scatter the compacted innovation covariance back into the full-nx Qp
1402 store, so the adaptive process-noise update can index it by II_RTK(). */
1403 if (adapt && prtk_ctx.Qp_nx != n) {
1404 free(prtk_ctx.Qp);
1405 prtk_ctx.Qp = zeros(n, n);
1406 prtk_ctx.Qp_nx = n;
1407 }
1408
1409 /* Only scatter on a successful update: filter2_() leaves Qp_ untouched
1410 (and mat() does not zero it) when matinv() fails, so scattering it then
1411 would seed the persistent Qp store with garbage. */
1412 for (i = 0; i < k; i++) {
1413 x[ix[i]] = xp_[i];
1414 for (j = 0; j < k; j++) {
1415 P[ix[i] + ix[j] * n] = Pp_[i + j * k];
1416 }
1417 if (adapt && info == 0) {
1418 for (j = 0; j < k; j++) {
1419 prtk_ctx.Qp[ix[i] + ix[j] * n] = Qp_[i + j * k];
1420 }
1421 }
1422 }
1423
1424 free(ix);
1425 free(x_);
1426 free(xp_);
1427 free(P_);
1428 free(Pp_);
1429 free(Qp_);
1430 free(H_);
1431 return info;
1432}
1433
1434/*============================================================================
1435 * DD residuals
1436 *===========================================================================*/
1437
1438/** Order in which ddres() is called within one ppp_rtk_pos() epoch */
1439enum ddres_order { FST = 1, SND, TRD };
1440
1441/**
1442 * @brief Test whether the two L6 channels currently carry the same facility.
1443 *
1444 * When both channels carry the same facility, either of them can supply any
1445 * satellite, so the reference satellite and each non-reference satellite may
1446 * pick their channel independently. With different facilities the two
1447 * correction sets are not interchangeable and the DD has to be kept within one
1448 * channel unless a common reference satellite exists.
1449 *
1450 * Mirrors upstream check_same_fac(), which compares the facility id of the two
1451 * CurrentCSSR entries directly.
1452 *
1453 * @param[in] clas CLAS context
1454 * @return 1 if both channels report the same facility, 0 otherwise
1455 */
1456static int check_same_fac(const clas_ctx_t* clas) { return clas->current[0].facility == clas->current[1].facility; }
1457
1458/**
1459 * @brief Check that a satellite yields a usable zero-difference signal.
1460 *
1461 * Validates the candidate non-reference satellite j against the reference
1462 * satellite i for system m and frequency f on one L6 channel's residual
1463 * vector, and reports the satellite/system/wavelength pair on success.
1464 *
1465 * @param[in] rtk RTK control struct
1466 * @param[in] obs Observation data
1467 * @param[in] nav Navigation data
1468 * @param[in] y ZD residuals of the channel under test
1469 * @param[in] nf Number of frequencies
1470 * @param[in] f Frequency index (f<nf: phase, f>=nf: code)
1471 * @param[in] i Reference satellite observation index
1472 * @param[in] j Candidate satellite observation index
1473 * @param[in] m System loop index
1474 * @param[in] sati Reference satellite number
1475 * @param[out] satj Candidate satellite number
1476 * @param[out] sysi Reference satellite system
1477 * @param[out] sysj Candidate satellite system
1478 * @param[out] lami Reference satellite wavelength (m)
1479 * @param[out] lamj Candidate satellite wavelength (m)
1480 * @return 1 if the pair is usable, 0 otherwise
1481 */
1482static int valid_zd_sig(rtk_t* rtk, const obsd_t* obs, const nav_t* nav, double* y, int nf, int f, int i, int j, int m,
1483 int sati, int* satj, int* sysi, int* sysj, double* lami, double* lamj) {
1484 prcopt_t* opt = &rtk->opt;
1485 int ff;
1486
1487 if (i == j) {
1488 return 0;
1489 }
1490 *satj = obs[j].sat;
1491 *sysi = rtk->ssat[sati - 1].sys;
1492 *sysj = rtk->ssat[*satj - 1].sys;
1493 if (!test_sys(*sysj, m, opt->qzsmodear, rtk->ssat[*satj - 1].code[f % nf])) {
1494 return 0;
1495 }
1496 if (!validobs(j, f, nf, y)) {
1497 return 0;
1498 }
1499
1500 ff = f % nf;
1501 *lami = sat_lambda(sati, ff, nav);
1502 *lamj = sat_lambda(*satj, ff, nav);
1503 if (*lami <= 0.0 || *lamj <= 0.0) {
1504 return 0;
1505 }
1506 return 1;
1507}
1508
1509/**
1510 * @brief Re-initialize the phase-bias state of one satellite from phase-code.
1511 *
1512 * Used when a satellite switches L6 channel: the accumulated bias belongs to
1513 * the previous channel's correction set and must not carry over. Mirrors the
1514 * second half of udbias_ppp(), restricted to a single satellite.
1515 *
1516 * @param[in,out] rtk RTK control struct
1517 * @param[in] nav Navigation data
1518 * @param[in] obs Observation data
1519 * @param[in] n Number of observations
1520 * @param[in] f Frequency index
1521 * @param[in] satj Satellite to reset
1522 */
1523static void set_init_pb(rtk_t* rtk, const nav_t* nav, const obsd_t* obs, int n, int f, int satj) {
1524 int i, j, isat;
1525 double offset, lami, cp, pr, com_bias;
1526 double* bias = zeros(n, 1);
1527
1528 /* estimate approximate phase-bias by phase - code */
1529 for (i = j = 0, offset = 0.0; i < n; i++) {
1530 cp = obs[i].L[f];
1531 pr = obs[i].P[f];
1532 lami = sat_lambda(obs[i].sat, f, nav);
1533 if (cp == 0.0 || pr == 0.0 || lami == 0.0) {
1534 continue;
1535 }
1536
1537 bias[i] = cp * lami - pr;
1538
1539 /* offset = sum of (bias - phase-bias) for all valid sats in meters */
1540 if (rtk->x[IB_RTK(obs[i].sat, f, &rtk->opt)] != 0.0) {
1541 offset += bias[i] - rtk->x[IB_RTK(obs[i].sat, f, &rtk->opt)] * lami;
1542 j++;
1543 }
1544 }
1545 com_bias = j > 0 ? offset / j : 0.0;
1546
1547 /* set initial state of the phase-bias of satj */
1548 for (i = 0; i < n; i++) {
1549 isat = obs[i].sat;
1550 if (isat != satj || bias[i] == 0.0) {
1551 continue;
1552 }
1553 lami = sat_lambda(isat, f, nav);
1554 if (lami == 0.0) {
1555 continue;
1556 }
1557 trace(NULL, 3, "set initial states of phase-bias, sat=%2d f=%2d\n", isat, f + 1);
1558 initx(rtk, (bias[i] - com_bias) / lami, SQR(rtk->opt.std[0]), IB_RTK(isat, f, &rtk->opt));
1559 rtk->ssat[isat - 1].lock[f] = -rtk->opt.minlock;
1560 rtk->ssat[isat - 1].pbreset[f] = 1;
1561 }
1562 free(bias);
1563}
1564
1565/**
1566 * @brief Compute DD residuals and measurement matrix, merging L6 channels.
1567 *
1568 * With opt->l6mrg set, residuals of both L6 channels are supplied and merged:
1569 * the reference satellite is preferentially one that both channels observe,
1570 * each non-reference satellite is assigned to a channel (priority channel for
1571 * l6mrg==1, load-balanced otherwise) with a fallback to the other channel, and
1572 * a satellite that changes channel has its phase bias re-initialized — only
1573 * when the channels carry different facilities: same-facility channels share
1574 * the correction basis, so a switch does not break phase-bias continuity.
1575 * With l6mrg clear only y[0]/pbslip[0] are read and the merge logic is
1576 * bypassed.
1577 *
1578 * @param[in,out] rtk RTK control struct
1579 * @param[in] nav Navigation data
1580 * @param[in] x State vector
1581 * @param[in,out] pbslip Per-channel phase bias slip corrections (NULL for postfit)
1582 * @param[in] P Covariance matrix (may be NULL)
1583 * @param[in] obs Observation data
1584 * @param[in] y Per-channel ZD residuals
1585 * @param[in] e Line-of-sight vectors (3*n)
1586 * @param[in] azel Azimuth/elevation (2*n)
1587 * @param[in] n Number of observations
1588 * @param[out] v DD residuals
1589 * @param[out] H Measurement matrix (may be NULL)
1590 * @param[out] R Measurement covariance
1591 * @param[out] vflg Residual flags
1592 * @param[in] niter Reference satellite search iteration
1593 * @param[in] pch Priority L6 channel
1594 * @param[in] order Call order within the epoch (FST/SND/TRD)
1595 * @return Number of DD residuals
1596 */
1597static int ddres(rtk_t* rtk, const nav_t* nav, double* x, double** pbslip, const double* P, const obsd_t* obs,
1598 double** y, double* e, double* azel, int n, double* v, double* H, double* R, int* vflg, int niter,
1599 int pch, int order) {
1600 prcopt_t* opt = &rtk->opt;
1601 double pos[3], lami, lamj, *Ri, *Rj, *Hi = NULL;
1602 double *tropu, *im, *dtdxu, didxi = 0.0, didxj = 0.0, fi, fj;
1603 /* nb is indexed once per (system group, freq, phase/code) DD block: the
1604 * system loop runs 6 groups, so size by 6 as upstream does (#307) */
1605 int i, j, k, l, m, f, nv = 0, nb[NFREQ * 6 * 2 + 1] = {0}, b = 0;
1606 int sati, satj, sysi, sysj, nf = NF_RTK(opt), flg = 0;
1607 int h, ch = -1, sch, refsatch, refchgflg, ch_sig_cnt[SSR_CH_NUM] = {0};
1608 int samefac = 0, nch = opt->l6mrg ? SSR_CH_NUM : 1;
1609 const clas_ctx_t* clas = (const clas_ctx_t*)nav->clas_ctx;
1610
1611 (void)P;
1612
1613 trace(NULL, 3, "ddres :niter=%2d nx=%d n=%d\n", niter, rtk->nx, n);
1614
1615 im = mat(n, 1);
1616 tropu = mat(n, 1);
1617 dtdxu = mat(n, 3);
1618 Ri = mat(n * nf * 2 * 2, 1);
1619 Rj = mat(n * nf * 2 * 2, 1);
1620 ecef2pos(x, pos);
1621
1622 if (opt->l6mrg) {
1623 memcpy(prtk_ctx.pre_use_ch, prtk_ctx.use_ch, sizeof(prtk_ctx.use_ch));
1624 for (f = 0; f < NFREQ; f++) {
1625 for (i = 0; i <= MAXSAT; i++) {
1626 prtk_ctx.use_ch[f][i] = -1;
1627 }
1628 }
1629 }
1630
1631 for (i = 0; i < MAXSAT; i++) {
1632 for (j = 0; j < NFREQ; j++) {
1633 rtk->ssat[i].resp[j] = rtk->ssat[i].resc[j] = 0.0;
1634 }
1635 }
1636
1637 /* compute ionospheric and tropospheric mapping factors */
1638 for (i = 0; i < n; i++) {
1639 if (opt->ionoopt >= IONOOPT_EST) {
1640 im[i] = ionmapf(pos, azel + i * 2);
1641 }
1642 if (opt->tropopt >= TROPOPT_EST) {
1643 tropu[i] = prectrop2(rtk->sol.time, pos, 0, azel + i * 2, opt, x, dtdxu + i * 3);
1644 }
1645 }
1646
1647 /* check if both channels carry the same facility */
1648 if (opt->l6mrg && clas) {
1649 samefac = check_same_fac(clas);
1650 }
1651
1652 /* system loop: 0=GPS/SBS, 1=GLO, 2=GAL, 3=BDS, 4=QZS, 5=GPS(L2C) */
1653 for (m = 0; m < 6; m++) {
1654 if (m == 1 && rtk->opt.glomodear == 0) {
1655 continue;
1656 }
1657 if (m == 4 && rtk->opt.qzsmodear != 1) {
1658 continue;
1659 }
1660
1661 /* count satellites of this system whose corrections are still fresh;
1662 * a single such satellite disables the freshness gate below so that a
1663 * lone satellite is not dropped outright (l6mrg only) */
1664 l = 0;
1665 if (opt->l6mrg) {
1666 for (j = 0; j < n; j++) {
1667 sati = obs[j].sat;
1668 if (!test_sys(rtk->ssat[sati - 1].sys, m, opt->qzsmodear, rtk->ssat[sati - 1].code[1])) {
1669 continue;
1670 }
1671 if (samefac) {
1672 /* same facility: either channel may supply the satellite */
1673 if (timediff(obs[0].time, nav->ssr_ch[0][sati - 1].t0[1]) >= 5.0 &&
1674 timediff(obs[0].time, nav->ssr_ch[1][sati - 1].t0[1]) >= 5.0) {
1675 continue;
1676 }
1677 } else {
1678 /* different facilities: both channels must be fresh */
1679 if (timediff(obs[0].time, nav->ssr_ch[0][sati - 1].t0[1]) >= 5.0) {
1680 continue;
1681 }
1682 if (timediff(obs[0].time, nav->ssr_ch[1][sati - 1].t0[1]) >= 5.0) {
1683 continue;
1684 }
1685 }
1686 l++;
1687 }
1688 }
1689
1690 for (f = opt->mode > PMODE_DGPS ? 0 : nf; f < nf * 2; f++) {
1691 /* search reference satellite with highest elevation. h==0 looks
1692 * for a satellite common to both channels; if none exists, h==1
1693 * accepts a satellite seen by one channel only. */
1694 for (h = 0, sch = -1, refchgflg = 0, refsatch = -1; h < ((opt->l6mrg && !samefac) ? SSR_CH_NUM : 1); h++) {
1695 for (i = -1, j = 0; j < n; j++) {
1696 flg = 0;
1697 sati = obs[j].sat;
1698 sysi = rtk->ssat[sati - 1].sys;
1699 if (!test_sys(sysi, m, rtk->opt.qzsmodear, rtk->ssat[sati - 1].code[f % nf])) {
1700 continue;
1701 }
1702 if (niter > 0) {
1703 for (k = 0; k < niter; k++) {
1704 if (prtk_ctx.refsat[NFREQ * 2 * MAXREF * m + NFREQ * 2 * k + f] == sati) {
1705 flg = 1;
1706 }
1707 }
1708 if (flg == 1) {
1709 continue;
1710 }
1711 }
1712 if (!opt->l6mrg) {
1713 if (!validobs(j, f, nf, y[0])) {
1714 continue;
1715 }
1716 } else if (samefac || h == 1) {
1717 /* refsat is taken from whichever single channel has it
1718 * with a fresh clock: in real time one channel can
1719 * carry the new 5 s clock before the other (#380) */
1720 for (ch = -1, k = 0; k < SSR_CH_NUM; k++) {
1721 if (!validobs(j, f, nf, y[k])) {
1722 continue;
1723 }
1724 if (l > 1 && timediff(obs[0].time, nav->ssr_ch[k][sati - 1].t0[1]) >= 5.0) {
1725 continue;
1726 }
1727 ch = k;
1728 break;
1729 }
1730 if (ch < 0) {
1731 continue;
1732 }
1733 } else {
1734 /* select a refsat common to both channels */
1735 if (!validobs(j, f, nf, y[0]) || !validobs(j, f, nf, y[1])) {
1736 continue;
1737 }
1738 if (l > 1 && timediff(obs[0].time, nav->ssr_ch[0][sati - 1].t0[1]) >= 5.0) {
1739 continue;
1740 }
1741 if (l > 1 && timediff(obs[0].time, nav->ssr_ch[1][sati - 1].t0[1]) >= 5.0) {
1742 continue;
1743 }
1744 }
1745 if (i < 0 || azel[1 + j * 2] >= azel[1 + i * 2]) {
1746 i = j;
1747 if (samefac) {
1748 refsatch = ch;
1749 }
1750 if (h == 1) {
1751 sch = ch;
1752 }
1753 }
1754 }
1755 if (i >= 0) {
1756 break;
1757 }
1758 }
1759 if (i < 0) {
1760 continue;
1761 }
1762
1763 sati = obs[i].sat;
1764 if (opt->l6mrg) {
1765 if (f < nf) {
1766 prtk_ctx.use_ch[f][sati] = 2; /* use_ch=2 means ref sat */
1767 if (prtk_ctx.refsat2[NFREQ * 2 * m + f] != sati && prtk_ctx.refsat2[NFREQ * 2 * m + f] != 0) {
1768 trace(NULL, 2, "refsat changed2 %s L%d sat=%2d -> %2d\n", time_str(rtk->sol.time, 0),
1769 f % nf + 1, prtk_ctx.refsat2[NFREQ * 2 * m + f], sati);
1770 refchgflg = 1;
1771 }
1772 }
1773 }
1774 prtk_ctx.refsat[NFREQ * 2 * MAXREF * m + NFREQ * 2 * niter + f] = sati;
1775 prtk_ctx.refsat2[NFREQ * 2 * m + f] = sati;
1776
1777 if (niter > 0) {
1778 trace(NULL, 2, "refsat changed %s L%d sat=%2d -> %2d\n", time_str(rtk->sol.time, 0), f % nf + 1,
1779 prtk_ctx.refsat[NFREQ * MAXREF * 2 * m + NFREQ * 2 * (niter - 1) + f],
1780 prtk_ctx.refsat[NFREQ * MAXREF * 2 * m + NFREQ * 2 * niter + f]);
1781 }
1782
1783 if (f < nf) {
1784 trace(NULL, 4, "refsat m:%d L%d %s: %d niter:%d\n", m, f + 1, time_str(rtk->sol.time, 0), sati, niter);
1785 }
1786
1787 /* apply reference satellite phase bias slip correction */
1788 if (f < nf && pbslip != NULL) {
1789 for (ch = 0; ch < nch; ch++) {
1790 x[IB_RTK(sati, f, opt)] += pbslip[ch][IB_RTK(sati, f, opt)];
1791 pbslip[ch][IB_RTK(sati, f, opt)] = 0.0;
1792 }
1793 }
1794
1795 /* if both channels share a facility, each satellite may pick its
1796 * own channel independently of the reference satellite's */
1797 if (opt->l6mrg && samefac) {
1798 sch = -1;
1799 }
1800
1801 /* form DD with each non-reference satellite */
1802 for (j = 0; j < n; j++) {
1803 if (!opt->l6mrg) {
1804 ch = 0;
1805 if (!valid_zd_sig(rtk, obs, nav, y[ch], nf, f, i, j, m, sati, &satj, &sysi, &sysj, &lami, &lamj)) {
1806 continue;
1807 }
1808 } else if (sch == -1) {
1809 if (opt->l6mrg == 1) {
1810 ch = pch; /* priority channel selection */
1811 } else {
1812 /* balanced channel selection */
1813 ch = ch_sig_cnt[0] <= ch_sig_cnt[1] ? 0 : 1;
1814 }
1815 if (!valid_zd_sig(rtk, obs, nav, y[ch], nf, f, i, j, m, sati, &satj, &sysi, &sysj, &lami, &lamj)) {
1816 ch = ch == 1 ? 0 : 1; /* fall back to the other channel */
1817 if (!valid_zd_sig(rtk, obs, nav, y[ch], nf, f, i, j, m, sati, &satj, &sysi, &sysj, &lami,
1818 &lamj)) {
1819 continue;
1820 }
1821 }
1822 } else {
1823 /* the reference satellite came from one channel only, so
1824 * the whole difference has to be formed on that channel */
1825 ch = sch;
1826 if (!valid_zd_sig(rtk, obs, nav, y[ch], nf, f, i, j, m, sati, &satj, &sysi, &sysj, &lami, &lamj)) {
1827 continue;
1828 }
1829 }
1830
1831 if (opt->l6mrg && f < nf) {
1832 prtk_ctx.use_ch[f][satj] = ch;
1833
1834 /* the phase bias of a satellite is tied to the corrections
1835 * of the channel it was computed on, so reset it when
1836 * (1) the reference satellite changed and this satellite
1837 * differs from the new reference's previous channel, or
1838 * (2) the satellite itself changed channel */
1839 if (order == FST && !samefac) {
1840 if ((prtk_ctx.pre_use_ch[f][sati] != ch && refchgflg) || (prtk_ctx.pre_use_ch[f][satj] != ch)) {
1841 set_init_pb(rtk, nav, obs, n, f, satj);
1842 }
1843 }
1844 }
1845
1846 /* apply non-ref satellite phase bias slip correction */
1847 if (f < nf && pbslip != NULL) {
1848 x[IB_RTK(satj, f, opt)] += pbslip[ch][IB_RTK(satj, f, opt)];
1849 pbslip[ch][IB_RTK(satj, f, opt)] = 0.0;
1850 }
1851
1852 if (H) {
1853 Hi = H + nv * rtk->nx;
1854 for (k = 0; k < rtk->nx; k++) {
1855 Hi[k] = 0.0;
1856 }
1857 }
1858
1859 /* DD residual; refsatch>=0 pins the reference satellite to the
1860 * channel it was found valid on (same-facility case) */
1861 double yi = refsatch < 0 ? y[ch][f + i * nf * 2] : y[refsatch][f + i * nf * 2];
1862 double yj = y[ch][f + j * nf * 2];
1863 double raw = yi - yj;
1864 double ion = 0.0, trop = 0.0, amb = 0.0;
1865 v[nv] = raw;
1866
1867 /* partial derivatives by rover position */
1868 if (H) {
1869 for (k = 0; k < 3; k++) {
1870 Hi[k] = -e[k + i * 3] + e[k + j * 3];
1871 }
1872 }
1873
1874 /* DD ionospheric delay term */
1875 if (opt->ionoopt == IONOOPT_EST || opt->ionoopt == IONOOPT_EST_ADPT) {
1876 fi = lami / lam_carr[0];
1877 fj = lamj / lam_carr[0];
1878 didxi = (f < nf ? -1.0 : 1.0) * fi * fi * im[i];
1879 didxj = (f < nf ? -1.0 : 1.0) * fj * fj * im[j];
1880 ion = didxi * x[II_RTK(sati, opt)] - didxj * x[II_RTK(satj, opt)];
1881 v[nv] -= ion;
1882 if (H) {
1883 Hi[II_RTK(sati, opt)] = didxi;
1884 Hi[II_RTK(satj, opt)] = -didxj;
1885 }
1886 }
1887
1888 /* DD tropospheric delay term */
1889 if (opt->tropopt == TROPOPT_EST || opt->tropopt == TROPOPT_ESTG) {
1890 trop = tropu[i] - tropu[j];
1891 v[nv] -= trop;
1892 for (k = 0; k < (opt->tropopt < TROPOPT_ESTG ? 1 : 3); k++) {
1893 if (!H) {
1894 continue;
1895 }
1896 Hi[IT_RTK(opt) + k] = (dtdxu[k + i * 3] - dtdxu[k + j * 3]);
1897 }
1898 }
1899
1900 /* DD phase-bias term */
1901 if (f < nf) {
1902 amb = lami * x[IB_RTK(sati, f, opt)] - lamj * x[IB_RTK(satj, f, opt)];
1903 v[nv] -= amb;
1904 if (H) {
1905 Hi[IB_RTK(sati, f, opt)] = lami;
1906 Hi[IB_RTK(satj, f, opt)] = -lamj;
1907 }
1908 }
1909
1910 /* measurement error variance */
1911 if (f < nf) { /* phase */
1912 rtk->ssat[satj - 1].resc[f] = v[nv];
1913 Ri[nv] = varerr(sati, sysi, azel[1 + i * 2], 0, opt);
1914 Rj[nv] = varerr(satj, sysj, azel[1 + j * 2], 0, opt);
1915 /* L2 scaling factor; MRTK sigcfg can place E5/L5 at slot 1. */
1916 if (f > 0 && (freq_idx2freq_num(sysi, f) == 2 || freq_idx2freq_num(sysj, f) == 2)) {
1917 Ri[nv] *= SQR(2.55 / 1.55);
1918 Rj[nv] *= SQR(2.55 / 1.55);
1919 }
1920 } else { /* code */
1921 rtk->ssat[satj - 1].resp[f - nf] = v[nv];
1922 Ri[nv] = varerr(sati, sysi, azel[1 + i * 2], 1, opt);
1923 Rj[nv] = varerr(satj, sysj, azel[1 + j * 2], 1, opt);
1924 }
1925
1926 /* set valid data flags */
1927 if (f < nf) {
1928 rtk->ssat[sati - 1].vsat[f] = 1;
1929 rtk->ssat[satj - 1].vsat[f] = 1;
1930 /* mark the satellite valid in the solution so the status
1931 * output (outsolstat $SAT, rtkoutstat $ION) reports it.
1932 * vs is otherwise only set by the SPP seed, which can omit
1933 * satellites that the PPP-RTK filter still uses (e.g. a
1934 * Galileo satellite usable on E1 only) -> they were missing
1935 * from the .stat. Mirrors upstream RTKLIB relpos ddres. */
1936 rtk->ssat[sati - 1].vs = 1;
1937 rtk->ssat[satj - 1].vs = 1;
1938 }
1939
1940 vflg[nv++] = (sati << 16) | (satj << 8) | ((f < nf ? 0 : 1) << 4) | (f % nf);
1941 nb[b]++;
1942
1943 if (ch == 0) {
1944 ch_sig_cnt[0]++;
1945 } else {
1946 ch_sig_cnt[1]++;
1947 }
1948 }
1949 b++;
1950 }
1951 } /* end of system loop */
1952
1953 /* measurement error covariance */
1954 ddcov(nb, b, Ri, Rj, nv, R);
1955
1956 free(Ri);
1957 free(Rj);
1958 free(im);
1959 free(tropu);
1960 free(dtdxu);
1961
1962 return nv;
1963}
1964
1965/*============================================================================
1966 * SD-to-DD transformation matrix
1967 *===========================================================================*/
1968
1969/**
1970 * @brief Build SD-to-DD transformation matrix D.
1971 * @param[in,out] rtk RTK control struct
1972 * @param[out] D Transformation matrix (nx x nx)
1973 * @param[out] indI Reference satellite bias indices
1974 * @param[out] indJ Non-reference satellite bias indices
1975 * @return Number of DD ambiguities (nb)
1976 */
1977static int ddmat(rtk_t* rtk, double* D, int* indI, int* indJ) {
1978 int i, j, k, m, f, nb = 0, nx = rtk->nx, na = rtk->na, code;
1979 int nf = NF_RTK(&rtk->opt);
1980
1981 trace(NULL, 3, "ddmat :\n");
1982
1983 for (i = 0; i < MAXSAT; i++) {
1984 for (j = 0; j < NFREQ; j++) {
1985 rtk->ssat[i].fix[j] = 0;
1986 }
1987 }
1988 for (i = 0; i < na; i++) {
1989 D[i + i * nx] = 1.0;
1990 }
1991
1992 for (m = 0; m < 6; m++) { /* 0:GPS/SBS, 1:GLO, 2:GAL, 3:BDS, 4:QZS, 5:GPS(L2C) */
1993 if (m == 1 && rtk->opt.glomodear == 0) {
1994 continue;
1995 }
1996 if (m == 4 && rtk->opt.qzsmodear != 1) {
1997 continue;
1998 }
1999
2000 for (f = 0, k = na; f < nf; f++, k += MAXSAT) {
2001 /* find reference satellite (highest elevation with lock) */
2002 for (i = k; i < k + MAXSAT; i++) {
2003 code = rtk->ssat[i - k].code[f];
2004 if (rtk->x[i] == 0.0 || !test_sys(rtk->ssat[i - k].sys, m, rtk->opt.qzsmodear, code) ||
2005 !rtk->ssat[i - k].vsat[f]) {
2006 continue;
2007 }
2008 if (rtk->ssat[i - k].lock[f] > 0 && !(rtk->ssat[i - k].slip[f] & 2) &&
2009 rtk->ssat[i - k].azel[1] >= rtk->opt.elmaskar) {
2010 rtk->ssat[i - k].fix[f] = 2; /* fix */
2011 break;
2012 } else {
2013 rtk->ssat[i - k].fix[f] = 1;
2014 }
2015 }
2016 /* form DD pairs */
2017 for (j = k; j < k + MAXSAT; j++) {
2018 code = rtk->ssat[j - k].code[f];
2019 if (i == j || rtk->x[j] == 0.0 || !test_sys(rtk->ssat[j - k].sys, m, rtk->opt.qzsmodear, code) ||
2020 !rtk->ssat[j - k].vsat[f]) {
2021 continue;
2022 }
2023 if (rtk->ssat[j - k].lock[f] > 0 && !(rtk->ssat[j - k].slip[f] & 2) &&
2024 rtk->ssat[j - k].azel[1] >= rtk->opt.elmaskar) {
2025 D[i + (na + nb) * nx] = 1.0;
2026 D[j + (na + nb) * nx] = -1.0;
2027 indI[nb] = i - na;
2028 indJ[nb] = j - na;
2029 nb++;
2030 rtk->ssat[j - k].fix[f] = 2; /* fix */
2031 } else {
2032 rtk->ssat[j - k].fix[f] = 1;
2033 }
2034 }
2035 }
2036 }
2037
2038 trace(NULL, 5, "D=\n");
2039 tracemat(NULL, 5, D, nx, na + nb, 2, 0);
2040
2041 return nb;
2042}
2043
2044/*============================================================================
2045 * Restore SD ambiguities from DD
2046 *===========================================================================*/
2047
2048/**
2049 * @brief Restore SD ambiguities from DD fixed values.
2050 * @param[in,out] rtk RTK control struct
2051 * @param[in] bias DD ambiguity biases
2052 * @param[in] nb Number of DD ambiguities
2053 * @param[out] xa Fixed state vector
2054 */
2055static void restamb(rtk_t* rtk, const double* bias, int nb, double* xa) {
2056 int i, n, m, f, index[MAXSAT], nv = 0, nf = NF_RTK(&rtk->opt);
2057
2058 trace(NULL, 3, "restamb :\n");
2059
2060 (void)nb;
2061
2062 /* init all fixed states to float state values */
2063 for (i = 0; i < rtk->nx; i++) {
2064 xa[i] = rtk->x[i];
2065 }
2066 /* overwrite non phase-bias states with fixed values */
2067 for (i = 0; i < rtk->na; i++) {
2068 xa[i] = rtk->xa[i];
2069 }
2070
2071 for (m = 0; m < 6; m++) {
2072 if (m == 1 && rtk->opt.glomodear == 0) {
2073 continue;
2074 }
2075 if (m == 4 && rtk->opt.qzsmodear != 1) {
2076 continue;
2077 }
2078
2079 for (f = 0; f < nf; f++) {
2080 for (n = i = 0; i < MAXSAT; i++) {
2081 if (!test_sys(rtk->ssat[i].sys, m, rtk->opt.qzsmodear, rtk->ssat[i].code[f]) ||
2082 rtk->ssat[i].fix[f] != 2) {
2083 continue;
2084 }
2085 index[n++] = IB_RTK(i + 1, f, &rtk->opt);
2086 rtk->ssat[i].lock[f]++; /* increment AR count */
2087 }
2088 if (n < 2) {
2089 continue;
2090 }
2091
2092 xa[index[0]] = rtk->x[index[0]];
2093
2094 for (i = 1; i < n; i++) {
2095 xa[index[i]] = xa[index[0]] - bias[nv++];
2096 }
2097 }
2098 }
2099}
2100
2101/*============================================================================
2102 * Hold integer ambiguity
2103 *===========================================================================*/
2104
2105/**
2106 * @brief Hold fixed ambiguities via pseudomeasurement constraint.
2107 * @param[in,out] rtk RTK control struct
2108 * @param[in] xa Fixed state vector
2109 */
2110static void holdamb(rtk_t* rtk, const double* xa) {
2111 double *v, *H, *R;
2112 int i, n, m, f, info, index[MAXSAT], nb = rtk->nx - rtk->na;
2113 int nv = 0, nf = NF_RTK(&rtk->opt);
2114
2115 trace(NULL, 3, "holdamb :\n");
2116
2117 v = mat(nb, 1);
2118 H = zeros(nb, rtk->nx);
2119
2120 for (m = 0; m < 6; m++) {
2121 if (m == 1 && rtk->opt.glomodear == 0) {
2122 continue;
2123 }
2124 if (m == 4 && rtk->opt.qzsmodear != 1) {
2125 continue;
2126 }
2127
2128 for (f = 0; f < nf; f++) {
2129 for (n = i = 0; i < MAXSAT; i++) {
2130 if (!test_sys(rtk->ssat[i].sys, m, rtk->opt.qzsmodear, rtk->ssat[i].code[f]) ||
2131 rtk->ssat[i].fix[f] != 2 || rtk->ssat[i].azel[1] < rtk->opt.elmaskhold || !rtk->ssat[i].vsat[f]) {
2132 continue;
2133 }
2134 index[n++] = IB_RTK(i + 1, f, &rtk->opt);
2135 rtk->ssat[i].fix[f] = 3; /* hold */
2136 }
2137 /* constraint to fixed ambiguity */
2138 for (i = 1; i < n; i++) {
2139 v[nv] = (xa[index[0]] - xa[index[i]]) - (rtk->x[index[0]] - rtk->x[index[i]]);
2140
2141 H[index[0] + nv * rtk->nx] = 1.0;
2142 H[index[i] + nv * rtk->nx] = -1.0;
2143 nv++;
2144 }
2145 }
2146 }
2147 if (nv > 0) {
2148 R = zeros(nv, nv);
2149 {
2150 double vh = rtk->opt.varholdamb > 0.0 ? rtk->opt.varholdamb : VARHOLDAMB;
2151 for (i = 0; i < nv; i++) {
2152 R[i + i * nv] = vh;
2153 }
2154 }
2155
2156 /* update states with constraints */
2157 if ((info = filter(rtk->x, rtk->P, H, v, R, rtk->nx, nv))) {
2158 trace(NULL, 2, "filter error (info=%d)\n", info);
2159 }
2160 free(R);
2161 }
2162 free(v);
2163 free(H);
2164}
2165
2166/*============================================================================
2167 * Sort satellites by elevation (ascending) for PAR
2168 *===========================================================================*/
2169
2170/**
2171 * @brief Sort satellites by elevation angle (ascending) for partial AR.
2172 * @param[in] rtk RTK control struct
2173 * @param[in] obs Observation data array
2174 * @param[out] isat Sorted satellite numbers
2175 * @param[out] el Sorted elevation angles (rad)
2176 * @param[in] n Number of observations
2177 */
2178static void elsort(rtk_t* rtk, const obsd_t* obs, int* isat, double* el, int n) {
2179 double el_;
2180 int i, j, k, m = 0, sati;
2181
2182 for (i = 0; i < n; i++) {
2183 sati = obs[i].sat;
2184 el_ = rtk->ssat[sati - 1].azel[1];
2185 if (m <= 0) {
2186 isat[m] = sati;
2187 el[m++] = el_;
2188 } else {
2189 for (j = 0; j < m; j++) {
2190 if (el_ < el[j]) {
2191 break;
2192 }
2193 }
2194 if (j >= n) {
2195 continue;
2196 }
2197 for (k = m; k > j; k--) {
2198 isat[k] = isat[k - 1];
2199 el[k] = el[k - 1];
2200 }
2201 isat[j] = sati;
2202 el[j] = el_;
2203 if (m < n) {
2204 m++;
2205 }
2206 }
2207 }
2208}
2209
2210/* compute single-difference PDOP */
2211static double sddop(const rtk_t* rtk, const double* H, const int* vflg, double elmask, int nv) {
2212 double *H_, Q[9], dop = 1000.0;
2213 int i, cnt, freq, type, sat;
2214
2215 H_ = mat(3, nv);
2216 for (i = cnt = 0; i < nv; i++) {
2217 sat = (vflg[i] >> 8) & 0xFF;
2218 type = (vflg[i] >> 4) & 0xF;
2219 freq = vflg[i] & 0xF;
2220 if (freq > 0 || type > 0 || rtk->ssat[sat - 1].vsat[0] == 0 || rtk->ssat[sat - 1].azel[1] <= elmask) {
2221 continue;
2222 }
2223 H_[cnt * 3] = H[i * rtk->nx];
2224 H_[cnt * 3 + 1] = H[i * rtk->nx + 1];
2225 H_[cnt * 3 + 2] = H[i * rtk->nx + 2];
2226 cnt++;
2227 }
2228 matmul("NT", 3, 3, cnt, 1.0, H_, H_, 0.0, Q);
2229 if (cnt > 0 && !matinv(Q, 3)) {
2230 dop = sqrt(Q[0] + Q[4] + Q[8]);
2231 }
2232 free(H_);
2233
2234 return dop;
2235}
2236
2237/* compute zero-difference PDOP */
2238static double zddop(const rtk_t* rtk, const double* azel, const obsd_t* obs, double elmask, int n) {
2239 double *azel_, dop[4];
2240 int i, j;
2241
2242 azel_ = mat(2, n);
2243 for (i = j = 0; i < n; i++) {
2244 if (rtk->ssat[obs[i].sat - 1].vsat[0] == 0) {
2245 continue;
2246 }
2247 azel_[j * 2] = azel[i * 2];
2248 azel_[j * 2 + 1] = azel[i * 2 + 1];
2249 j++;
2250 }
2251 dops(j, azel_, elmask, dop);
2252 free(azel_);
2253
2254 return dop[1];
2255}
2256
2257/* compute PDOP for AR gates */
2258static double calpdop(const rtk_t* rtk, const double* H, const double* azel, const int* vflg, const obsd_t* obs,
2259 double elmask, int n, int nv) {
2260 switch (rtk->opt.refdop) {
2261 case 0:
2262 return zddop(rtk, azel, obs, elmask, n);
2263 case 1:
2264 return sddop(rtk, H, vflg, elmask, nv);
2265 }
2266 return 0.0;
2267}
2268
2269/*============================================================================
2270 * LAMBDA integer ambiguity resolution
2271 *===========================================================================*/
2272
2273/**
2274 * @brief Resolve integer ambiguity by LAMBDA.
2275 * @param[in,out] rtk RTK control struct
2276 * @param[out] bias DD fixed ambiguity biases
2277 * @param[out] xa Fixed state vector
2278 * @return Number of resolved ambiguities (0 = failed)
2279 */
2280static int resamb_LAMBDA(rtk_t* rtk, double* bias, double* xa) {
2281 int i, j, nb, nc, info, nx = rtk->nx, na = rtk->na;
2282 int minamb, alpha = -1;
2283 int indI[MAXOBS * NFREQ], indJ[MAXOBS * NFREQ];
2284 double *D, *y, *b, *db, *Qb, *Qab, *QQ, s[2], thres = 0.0, ratio = 0.0;
2285 double* D_P;
2286
2287 trace(NULL, 3, "resamb_LAMBDA : nx=%d\n", nx);
2288
2289 rtk->sol.ratio = 0.0;
2290
2291 if (rtk->opt.mode <= PMODE_DGPS || rtk->opt.modear == ARMODE_OFF) {
2292 return 0;
2293 }
2294
2295 /* SD-to-DD transformation matrix */
2296 D = zeros(nx, nx);
2297 nb = ddmat(rtk, D, indI, indJ);
2298 minamb = rtk->opt.minamb > 0 ? rtk->opt.minamb : MINAMB_DEFAULT;
2299
2300 if (nb <= 1) {
2301 trace(NULL, 2, "no valid double-difference\n");
2302 free(D);
2303 return 0;
2304 }
2305 if (nb < minamb) {
2306 trace(NULL, 2, "less than minimum number of ambiguities nb=%2d minamb=%d\n", nb, minamb);
2307 free(D);
2308 return 0;
2309 }
2310
2311 y = mat(na + nb, 1);
2312 b = mat(nb, 2);
2313 db = mat(nb, 1);
2314 Qb = mat(nb, nb);
2315 Qab = mat(na, nb);
2316 QQ = mat(na, nb);
2317 nc = nx - na;
2318 D_P = mat(nb, nc);
2319
2320 for (i = 0; i < na; i++) {
2321 y[i] = rtk->x[i];
2322 for (j = 0; j < nb; j++) {
2323 Qab[i + j * na] = rtk->P[i + (na + indI[j]) * nx] - rtk->P[i + (na + indJ[j]) * nx];
2324 }
2325 }
2326 for (i = 0; i < nb; i++) {
2327 y[i + na] = rtk->x[na + indI[i]] - rtk->x[na + indJ[i]];
2328 for (j = 0; j < nc; j++) {
2329 D_P[i + j * nb] = rtk->P[na + indI[i] + (na + j) * nx] - rtk->P[na + indJ[i] + (na + j) * nx];
2330 }
2331 }
2332 for (i = 0; i < nb; i++) {
2333 for (j = 0; j < nb; j++) {
2334 Qb[i + j * nb] = D_P[i + indI[j] * nb] - D_P[i + indJ[j] * nb];
2335 }
2336 }
2337 free(D_P);
2338
2339 /* LAMBDA integer least-squares */
2340 if (!(info = lambda(nb, 2, y + na, Qb, b, s))) {
2341 trace(NULL, 5, "N(1)=\n");
2342 tracemat(NULL, 5, b, 1, nb, 10, 3);
2343 trace(NULL, 5, "N(2)=\n");
2344 tracemat(NULL, 5, b + nb, 1, nb, 10, 3);
2345
2346 rtk->sol.ratio = (s[0] > 0) ? (float)(s[1] / s[0]) : 0.0f;
2347 if (rtk->sol.ratio > 999.9) {
2348 rtk->sol.ratio = 999.9f;
2349 }
2350
2351 /* F-test threshold (significance level from config or default) */
2352 alpha = rtk->opt.alphaar > 0 ? rtk->opt.alphaar : ALPHAAR_DEFAULT;
2353 thres = (alpha < 6 && nb - 1 < 60) ? qf[alpha][nb - 1] : 1.5;
2354 ratio = (s[0] == 0.0) ? 0.0 : s[1] / s[0];
2355
2356 /* validation by ratio-test */
2357 if (s[0] <= 0.0 || ratio >= thres) {
2358 for (i = 0; i < na; i++) {
2359 rtk->xa[i] = rtk->x[i];
2360 for (j = 0; j < na; j++) {
2361 rtk->Pa[i + j * na] = rtk->P[i + j * nx];
2362 }
2363 }
2364 for (i = 0; i < nb; i++) {
2365 bias[i] = b[i];
2366 y[na + i] -= b[i];
2367 }
2368 if (!matinv(Qb, nb)) {
2369 /* transform float to fixed solution: xa=xa-Qab*Qb\‍(b0-b) */
2370 matmul("NN", nb, 1, nb, 1.0, Qb, y + na, 0.0, db);
2371 matmul("NN", na, 1, nb, -1.0, Qab, db, 1.0, rtk->xa);
2372
2373 /* covariance of fixed solution: Qa=Qa-Qab*Qb^-1*Qab' */
2374 matmul("NN", na, nb, nb, 1.0, Qab, Qb, 0.0, QQ);
2375 matmul("NT", na, na, nb, -1.0, QQ, Qab, 1.0, rtk->Pa);
2376
2377 trace(NULL, 3, "resamb : validation ok (nb=%d ratio=%.2f s=%.2f/%.2f)\n", nb, ratio, s[1], s[0]);
2378
2379 /* restore SD ambiguities */
2380 restamb(rtk, bias, nb, xa);
2381 } else {
2382 trace(NULL, 2, "lambda error nb=%d\n", nb);
2383 nb = 0;
2384 }
2385 } else {
2386 trace(NULL, 3,
2387 "ambiguity validation failed "
2388 "(nb=%d thres=%.2f ratio=%.2f s=%.2f/%.2f)\n",
2389 nb, thres, ratio, s[1], s[0]);
2390 nb = 0;
2391 }
2392 } else {
2393 trace(NULL, 2, "lambda error (info=%d)\n", info);
2394 nb = 0;
2395 }
2396
2397 free(D);
2398 free(y);
2399 free(b);
2400 free(db);
2401 free(Qb);
2402 free(Qab);
2403 free(QQ);
2404
2405 return nb; /* number of ambiguities */
2406}
2407
2408/*============================================================================
2409 * Public API: solution status output
2410 *===========================================================================*/
2411
2412/**
2413 * @brief Write PPP-RTK solution status to buffer.
2414 *
2415 * PPP-RTK uses a different state layout from standard PPP (no clock states,
2416 * iono at NP, trop at NP+NI). This function outputs state information using
2417 * the correct PPP-RTK indices (II_RTK, IT_RTK, IB_RTK).
2418 *
2419 * @param[in,out] rtk RTK control/result struct
2420 * @param[out] buff Output buffer
2421 * @return Number of bytes written
2422 */
2423extern int ppprtk_outstat(rtk_t* rtk, char* buff) {
2424 ssat_t* ssat;
2425 double tow, pos[3], vel[3], acc[3], *x;
2426 int i, j, week, nx = rtk->nx;
2427 char id[32], *p = buff;
2428
2429 if (!rtk->sol.stat) {
2430 return 0;
2431 }
2432
2433 tow = time2gpst(rtk->sol.time, &week);
2434 x = rtk->sol.stat == SOLQ_FIX ? rtk->xa : rtk->x;
2435
2436 /* receiver position */
2437 p += sprintf(p, "$POS,%d,%.3f,%d,%.4f,%.4f,%.4f,%.4f,%.4f,%.4f\n", week, tow, rtk->sol.stat, x[0], x[1], x[2],
2438 SQRT(rtk->P[0]), SQRT(rtk->P[1 + nx]), SQRT(rtk->P[2 + 2 * nx]));
2439
2440 /* receiver velocity and acceleration */
2441 if (rtk->opt.dynamics) {
2442 ecef2pos(rtk->sol.rr, pos);
2443 ecef2enu(pos, rtk->x + 3, vel);
2444 ecef2enu(pos, rtk->x + 6, acc);
2445 p += sprintf(p,
2446 "$VELACC,%d,%.3f,%d,%.4f,%.4f,%.4f,%.5f,%.5f,%.5f,%.4f,%.4f,"
2447 "%.4f,%.5f,%.5f,%.5f\n",
2448 week, tow, rtk->sol.stat, vel[0], vel[1], vel[2], acc[0], acc[1], acc[2], 0.0, 0.0, 0.0, 0.0, 0.0,
2449 0.0);
2450 }
2451 /* no clock states in PPP-RTK (absorbed by DD) */
2452
2453 /* tropospheric parameters */
2454 if (rtk->opt.tropopt == TROPOPT_EST || rtk->opt.tropopt == TROPOPT_ESTG) {
2455 i = IT_RTK(&rtk->opt);
2456 if (i < nx) {
2457 p += sprintf(p, "$TROP,%d,%.3f,%d,%d,%.4f,%.4f\n", week, tow, rtk->sol.stat, 1, x[i],
2458 SQRT(rtk->P[i + i * nx]));
2459 }
2460 }
2461 if (rtk->opt.tropopt == TROPOPT_ESTG) {
2462 i = IT_RTK(&rtk->opt);
2463 if (i + 2 < nx) {
2464 p += sprintf(p, "$TRPG,%d,%.3f,%d,%d,%.5f,%.5f,%.5f,%.5f\n", week, tow, rtk->sol.stat, 1, x[i + 1],
2465 x[i + 2], SQRT(rtk->P[(i + 1) + (i + 1) * nx]), SQRT(rtk->P[(i + 2) + (i + 2) * nx]));
2466 }
2467 }
2468 /* ionosphere parameters */
2469 if (rtk->opt.ionoopt == IONOOPT_EST || rtk->opt.ionoopt == IONOOPT_EST_ADPT) {
2470 for (i = 0; i < MAXSAT; i++) {
2471 ssat = rtk->ssat + i;
2472 if (!ssat->vs) {
2473 continue;
2474 }
2475 j = II_RTK(i + 1, &rtk->opt);
2476 if (j >= nx || rtk->x[j] == 0.0) {
2477 continue;
2478 }
2479 satno2id(i + 1, id);
2480 p += sprintf(p, "$ION,%d,%.3f,%d,%s,%.1f,%.1f,%.4f,%.4f\n", week, tow, rtk->sol.stat, id,
2481 ssat->azel[0] * R2D, ssat->azel[1] * R2D, x[j], SQRT(rtk->P[j + j * nx]));
2482 }
2483 }
2484 return (int)(p - buff);
2485}
2486
2487/*============================================================================
2488 * Public API: number of estimated states
2489 *===========================================================================*/
2490
2491/**
2492 * @brief Compute number of estimated states for PPP-RTK.
2493 * @param[in] opt Processing options
2494 * @return Total number of states
2495 */
2496extern int ppp_rtk_nx(const prcopt_t* opt) { return NX_RTK(opt); }
2497
2498extern int ppp_rtk_na(const prcopt_t* opt) { return NR_RTK(opt); }
2499
2500static void check_clas_facility(nav_t* nav, const clas_ctx_t* clas, const clas_grid_t* grid, const prcopt_t* opt) {
2501 static int savefacility[SSR_CH_NUM] = {-1, -1};
2502 int ch, nch = opt->l6mrg ? SSR_CH_NUM : 1;
2503
2504 (void)grid;
2505
2506 for (ch = 0; ch < nch; ch++) {
2507 int facility = clas->current[ch].facility;
2508 /* a channel without corrections carries no facility id: keep the last
2509 * known one so that a correction gap does not look like a facility change */
2510 if (!clas->current[ch].use) {
2511 continue;
2512 }
2513 if (savefacility[ch] != facility) {
2514 if (savefacility[ch] != -1) {
2515 nav->filreset = 1;
2516 }
2517 savefacility[ch] = facility;
2518 }
2519 }
2520}
2521
2522/*============================================================================
2523 * Public API: PPP-RTK positioning for one epoch
2524 *===========================================================================*/
2525
2526/**
2527 * @brief PPP-RTK positioning for one epoch (CLAS).
2528 *
2529 * Main entry point for the CLAS PPP-RTK positioning engine. Uses
2530 * double-differencing with CLAS grid corrections (STEC, ZWD) and iterative
2531 * Kalman filter with LAMBDA integer ambiguity resolution.
2532 *
2533 * Uses CLAS grid corrections with a short CSSR backup fallback matching claslib.
2534 *
2535 * @param[in,out] rtk RTK control/result struct
2536 * @param[in] obs Observation data for epoch
2537 * @param[in] n Number of observations
2538 * @param[in,out] nav Navigation data (includes SSR via nav->ssr_ch[][])
2539 */
2540extern void ppp_rtk_pos(rtk_t* rtk, const obsd_t* obs, int n, nav_t* nav) {
2541 static gtime_t regularly = {-1, 0.0};
2542 prcopt_t* opt = &rtk->opt;
2543 clas_ctx_t* clas;
2544 clas_grid_t* grid;
2545 double *v, *H, *R, *Pp, *bias, *xp, *xa, *azel, *e;
2546 double *y[SSR_CH_NUM] = {0}, *rs[SSR_CH_NUM] = {0}, *dts[SSR_CH_NUM] = {0}, *var[SSR_CH_NUM] = {0};
2547 double* pbslip[SSR_CH_NUM] = {0};
2548 double pos[3], pdop;
2549 int i, j, k, f, l, nv, nvtmp[SSR_CH_NUM], nvsum, info, nf = rtk->opt.nf, sati;
2550 int ch, pch = 0, nch = opt->l6mrg ? SSR_CH_NUM : 1;
2551 int stat = (rtk->opt.mode <= PMODE_DGPS) ? SOLQ_DGPS : SOLQ_FLOAT;
2552 int vflg[MAXOBS * NFREQ * 4 + 1], nb, svh[SSR_CH_NUM][MAXOBS];
2553 int nn, use_backup = 0;
2554 double cpc_[SSR_CH_NUM][MAXSAT * NFREQ];
2555 gtime_t pt0_[SSR_CH_NUM][MAXSAT];
2556
2557 trace(NULL, 2, "ppp_rtk_pos: time=%s nx=%d n=%d\n", time_str(obs[0].time, 0), rtk->nx, n);
2558 /* check CLAS context availability */
2559 if (!nav->clas_ctx) {
2560 trace(NULL, 2, "ppp_rtk_pos: no CLAS context\n");
2561 rtk->sol.stat = SOLQ_NONE;
2562 return;
2563 }
2564 clas = (clas_ctx_t*)nav->clas_ctx;
2565
2566 azel = zeros(2, n);
2567 e = mat(n, 3);
2568
2569 nv = n * nf * 2 * 3;
2570
2571 for (ch = 0; ch < nch; ch++) {
2572 y[ch] = zeros(nv, 1);
2573 rs[ch] = mat(6, n);
2574 dts[ch] = mat(2, n);
2575 var[ch] = mat(1, n);
2576 }
2577
2578 /* initialize satellite system info */
2579 for (i = 0; i < MAXSAT; i++) {
2580 rtk->ssat[i].sys = satsys(i + 1, NULL);
2581 for (j = 0; j < NFREQ; j++) {
2582 rtk->ssat[i].vsat[j] = 0;
2583 rtk->ssat[i].snr[j] = 0;
2584 }
2585 }
2586 for (i = 0; i < MAXSAT; i++) {
2587 rtk->ssat[i].fix[0] = 0;
2588 }
2589
2590 /* reset all states if time gap exceeds threshold */
2591 if (rtk->tt > (opt->maxobsloss_s > 0.0 ? opt->maxobsloss_s : MAXOBSLOSS_DEFAULT)) {
2592 for (i = 0; i < rtk->nx; i++) {
2593 rtk->x[i] = 0.0;
2594 }
2595 }
2596
2597 /* periodic filter reset: force a full state reset every opt->regularly
2598 * seconds (0 = disabled). Ported from upstream claslib ppp_rtk_pos()
2599 * (#264); mirrors the already-ported VRS path in relposvrs(). Returns
2600 * early to skip the epoch, so the buffers allocated above must be freed
2601 * here (unlike VRS, which allocates them after this check). */
2602 regularly = (regularly.time == -1 ? obs[0].time : regularly);
2603 if (opt->regularly != 0 && timediff(obs[0].time, regularly) >= (double)opt->regularly) {
2604 trace(NULL, 1, "ppp_rtk_pos: regularly reset filter, tow=%.1f\n", time2gpst(obs[0].time, NULL));
2605 for (i = 0; i < rtk->nx; i++) {
2606 rtk->x[i] = 0.0;
2607 }
2608 rtk->sol.stat = SOLQ_SINGLE;
2609 regularly = obs[0].time;
2610 nav->filreset = 0;
2611 prtk_ctx.float_count = 0;
2612 free(azel);
2613 free(e);
2614 for (ch = 0; ch < nch; ch++) {
2615 free(y[ch]);
2616 free(rs[ch]);
2617 free(dts[ch]);
2618 free(var[ch]);
2619 }
2620 return;
2621 }
2622
2623 /* reset filter if consecutive FLOAT epochs exceed floatcnt threshold */
2624 if (opt->floatcnt > 0 && prtk_ctx.float_count >= opt->floatcnt) {
2625 trace(NULL, 2, "ppp_rtk_pos: float persisted %d epochs, resetting filter\n", prtk_ctx.float_count);
2626 if (opt->dynamics) {
2627 /* kinematic: reset all states (position re-initialized from SPP) */
2628 for (i = 0; i < rtk->nx; i++) {
2629 rtk->x[i] = 0.0;
2630 }
2631 } else {
2632 /* static: keep position, reset iono/trop/bias */
2633 for (i = NP_RTK(opt); i < rtk->nx; i++) {
2634 rtk->x[i] = 0.0;
2635 }
2636 }
2637 prtk_ctx.float_count = 0;
2638 prtk_ctx.floatcnt_reset = 1;
2639 }
2640
2641 if (nav->filreset) {
2642 trace(NULL, 1, "ppp_rtk_pos: reset filter, tow=%.1f\n", time2gpst(obs[0].time, NULL));
2643 if (opt->dynamics) {
2644 for (i = 0; i < rtk->nx; i++) {
2645 rtk->x[i] = 0.0;
2646 }
2647 } else {
2648 for (i = NP_RTK(opt); i < rtk->nx; i++) {
2649 rtk->x[i] = 0.0;
2650 }
2651 }
2652 rtk->sol.stat = SOLQ_FLOAT;
2653 nav->filreset = 0;
2654 stat = SOLQ_FLOAT;
2655 }
2656
2657 /* temporal update of states */
2658 udstate_ppp(rtk, obs, n, nav);
2659
2660 /* get grid index */
2661 grid = &clas->grid[0]; /* all channels share one grid, as upstream does */
2662 ecef2pos(rtk->x, pos);
2663 if ((nn = clas_get_grid_index(clas, pos, grid, opt->gridsel, obs[0].time)) <= 0) {
2664 if (!clas_backup_valid(clas, obs[0].time, opt->l6mrg)) {
2665 trace(NULL, 2, "ppp_rtk_pos: no valid grid\n");
2666 free(azel);
2667 free(e);
2668 for (ch = 0; ch < nch; ch++) {
2669 free(y[ch]);
2670 free(rs[ch]);
2671 free(dts[ch]);
2672 free(var[ch]);
2673 }
2674 rtk->sol.stat = SOLQ_SINGLE;
2675 return;
2676 }
2677 clas_restore_backup(clas, obs[0].time, grid, opt->l6mrg);
2678 use_backup = 1;
2679 }
2680
2681 /* fetch corrections for all active channels, or restore a short-gap backup.
2682 * OR semantics as in claslib get_close_cssr(): a channel that cannot supply
2683 * corrections is skipped and the epoch fails only if every channel fails. */
2684 {
2685 int nch = opt->l6mrg ? SSR_CH_NUM : 1;
2686 int ch, nok = 0;
2687 int chok[SSR_CH_NUM] = {0};
2688 for (ch = 0; ch < nch && !use_backup; ch++) {
2689 /* re-fetch corrections closest to observation time,
2690 * falling back to L6 buffer time if obs lookup fails */
2691 if (grid->network > 0) {
2692 int bgrc = clas_bank_get_close(clas, obs[0].time, grid->network, ch, &clas->current[ch]);
2693 if (bgrc != 0) {
2694 /* fallback: try L6 buffer time */
2695 bgrc = clas_bank_get_close(clas, clas->l6buf[ch].time, grid->network, ch, &clas->current[ch]);
2696 }
2697 if (bgrc != 0) {
2698 trace(NULL, 3, "ppp_rtk_pos: bank lookup failed ch=%d\n", ch);
2699 continue;
2700 }
2701 clas_check_grid_status(clas, &clas->current[ch], ch);
2702 }
2703 chok[ch] = 1;
2704 nok++;
2705 }
2706 if (opt->l6mrg && nok == 1) {
2707 static int warned_single_ch = 0;
2708 if (!warned_single_ch) {
2709 warned_single_ch = 1;
2710 trace(NULL, 2,
2711 "ppp_rtk_pos: l6_merge running single-channel on ch=%d (only one L6 transmit pattern "
2712 "received so far)\n",
2713 chok[0] ? 0 : 1);
2714 }
2715 }
2716 if (!use_backup && nok == 0) {
2717 if (!clas_backup_valid(clas, obs[0].time, opt->l6mrg)) {
2718 free(azel);
2719 free(e);
2720 for (ch = 0; ch < nch; ch++) {
2721 free(y[ch]);
2722 free(rs[ch]);
2723 free(dts[ch]);
2724 free(var[ch]);
2725 }
2726 rtk->sol.stat = SOLQ_SINGLE;
2727 return;
2728 }
2729 clas_restore_backup(clas, obs[0].time, grid, opt->l6mrg);
2730 use_backup = 1;
2731 } else if (!use_backup) {
2732 clas_backup_current(clas, grid, opt->l6mrg);
2733 /* the snapshot of a channel that failed this epoch is stale: drop it so
2734 * that neither ssr_ch[ch] nor the clock merge below picks it up */
2735 for (ch = 0; ch < nch; ch++) {
2736 if (!chok[ch]) {
2737 clas_clear_current(clas, ch);
2738 }
2739 }
2740 }
2741 for (ch = 0; ch < nch; ch++) {
2742 clas_update_global(nav, &clas->current[ch], ch);
2743 clas_update_local(nav, &clas->current[ch], ch);
2744 }
2745 check_clas_facility(nav, clas, grid, opt);
2746 }
2747
2748 trace(NULL, 4, "x(0)=\n");
2749 tracemat(NULL, 4, rtk->x, 1, NR_RTK(opt), 13, 4);
2750
2751 /* satellite positions and clocks, computed separately per channel:
2752 * set_ssr_ch_idx() steers the SSR lookups inside satposs()/satcorr to the
2753 * channel being processed */
2754 for (i = 0; i < 3; i++) {
2755 prtk_ctx.osr_ctx.saved_rr[i] = rtk->x[i];
2756 }
2757 for (ch = 0; ch < nch; ch++) {
2758 set_ssr_ch_idx(ch);
2759 satposs(obs[0].time, obs, n, nav, rtk->opt.sateph, rs[ch], dts[ch], var[ch], svh[ch]);
2760 for (i = 0; i < MAXSAT; i++) {
2761 clas_osr_satcorr_update(obs[0].time, obs[0].time, i + 1, &rtk->opt, nav, ch, &prtk_ctx.osr_ctx);
2762 }
2763 }
2764
2765 xp = mat(rtk->nx, 1);
2766 Pp = zeros(rtk->nx, rtk->nx);
2767 xa = mat(rtk->nx, 1);
2768 for (ch = 0; ch < nch; ch++) {
2769 pbslip[ch] = zeros(rtk->nx, 1);
2770 }
2771
2772 matcpy(xp, rtk->x, rtk->nx, 1);
2773
2774 if ((rtk->opt.ionoopt == IONOOPT_EST_ADPT || rtk->opt.prnadpt) && prtk_ctx.Qp && prtk_ctx.Qp_nx == rtk->nx) {
2775 for (i = 0; i < rtk->nx * rtk->nx; i++) {
2776 prtk_ctx.Qp[i] = 0.0;
2777 }
2778 }
2779
2780 v = mat(nv, 1);
2781 R = mat(nv, nv);
2782 H = mat(rtk->nx, nv);
2783 bias = mat(rtk->nx, 1);
2784
2785 /* iterative measurement update */
2786 for (i = 0; i < rtk->opt.niter; i++) {
2787 for (k = 0; k < MAXREF; k++) {
2788 /* save persistent zdres state, and compute the prefit
2789 * zero-difference residuals separately on each channel */
2790 for (ch = 0, nvsum = 0; ch < nch; ch++) {
2791 for (j = 0; j < MAXSAT; j++) {
2792 pt0_[ch][j] = prtk_ctx.pt0[ch][j];
2793 for (f = 0; f < NFREQ; f++) {
2794 cpc_[ch][f * MAXSAT + j] = prtk_ctx.cpc[ch][f * MAXSAT + j];
2795 }
2796 }
2797 for (j = 0; j < 3; j++) {
2798 pbslip[ch][j] = xp[j];
2799 }
2800 load_ppp_osr_cpc(&prtk_ctx.osr_ctx, cpc_[ch], pt0_[ch]);
2801 nvtmp[ch] = clas_osr_zdres(obs, n, rs[ch], dts[ch], var[ch], svh[ch], nav, pbslip[ch], y[ch], e, azel,
2802 rtk, 1, &prtk_ctx.osr_ctx, grid, &clas->current[ch], rtk->ssat, &rtk->opt,
2803 &rtk->sol, NULL, ch);
2804 save_ppp_osr_cpc(&prtk_ctx.osr_ctx, cpc_[ch], pt0_[ch]);
2805 nvsum += nvtmp[ch];
2806 /* the channel contributing more measurements is preferred when
2807 * a satellite is available on both */
2808 if (ch != 0 && opt->l6mrg) {
2809 pch = nvtmp[0] >= nvtmp[1] ? 0 : 1;
2810 }
2811 }
2812
2813 if (!nvsum) {
2814 trace(NULL, 2, "rover initial position error\n");
2815 stat = SOLQ_NONE;
2816 break;
2817 }
2818 /* DD residuals and partial derivatives */
2819 if ((nv = ddres(rtk, nav, xp, pbslip, Pp, obs, y, e, azel, n, v, H, R, vflg, k, pch, FST)) <= 0) {
2820 trace(NULL, 2, "no double-differenced residual\n");
2821 stat = SOLQ_NONE;
2822 break;
2823 }
2824 /* Kalman filter measurement update */
2825 matcpy(Pp, rtk->P, rtk->nx, rtk->nx);
2826 matcpy(xp, rtk->x, rtk->nx, 1);
2827 if ((info = filter2(rtk, xp, Pp, H, v, R, rtk->nx, nv, vflg, 0)) != 0) {
2828 trace(NULL, 2, "ppp-rtk filter error %s info=%d\n", time_str(rtk->sol.time, 0), info);
2829 if (i == (rtk->opt.niter - 1) && k == (MAXREF - 1)) {
2830 stat = SOLQ_NONE;
2831 }
2832 continue;
2833 }
2834
2835 /* recompute zero-difference residuals (postfit) */
2836 for (ch = 0, nvsum = 0; ch < nch; ch++) {
2837 load_ppp_osr_cpc(&prtk_ctx.osr_ctx, cpc_[ch], pt0_[ch]);
2838 nvtmp[ch] = clas_osr_zdres(obs, n, rs[ch], dts[ch], var[ch], svh[ch], nav, xp, y[ch], e, azel, rtk, 0,
2839 &prtk_ctx.osr_ctx, grid, &clas->current[ch], rtk->ssat, &rtk->opt, &rtk->sol,
2840 NULL, ch);
2841 save_ppp_osr_cpc(&prtk_ctx.osr_ctx, cpc_[ch], pt0_[ch]);
2842 nvsum += nvtmp[ch];
2843 if (ch != 0 && opt->l6mrg) {
2844 pch = nvtmp[0] >= nvtmp[1] ? 0 : 1;
2845 }
2846 }
2847
2848 if (nvsum) {
2849 nv = ddres(rtk, nav, xp, NULL, Pp, obs, y, e, azel, n, v, NULL, R, vflg, k, pch, SND);
2850 /* validation of float solution */
2851 filter2(rtk, xp, Pp, H, v, R, rtk->nx, nv, vflg, 1);
2852 if (prtk_ctx.chisq < 100.0) {
2853 break;
2854 }
2855 }
2856 matcpy(xp, rtk->x, rtk->nx, 1);
2857 if (i == (rtk->opt.niter - 1) && k == (MAXREF - 1)) {
2858 stat = SOLQ_NONE;
2859 trace(NULL, 2, "measurement update failed\n");
2860 }
2861 }
2862 }
2863
2864 /* SSR age: the freshest correction over all active channels */
2865 rtk->sol.age = 1e4;
2866 for (ch = 0; ch < nch; ch++) {
2867 for (i = 0; i < n && i < MAXOBS; i++) {
2868 float age;
2869 sati = obs[i].sat;
2870 age = (float)timediff(obs[i].time, nav->ssr_ch[ch][sati - 1].t0[1]);
2871 if (rtk->sol.age > age) {
2872 rtk->sol.age = age;
2873 }
2874 }
2875 }
2876
2877 /* Upstream updates adaptive process noise from the epoch-local prefit Qp
2878 * before deciding whether the float state itself is accepted. Keep the same
2879 * ordering so rejection-heavy storm windows do not freeze the Q history. */
2880 if (rtk->opt.ionoopt == IONOOPT_EST_ADPT && prtk_ctx.Qp) {
2881 for (i = 0; i < n; i++) {
2882 j = II_RTK(obs[i].sat, &rtk->opt);
2883 if (j < prtk_ctx.Qp_nx) {
2884 rtk->Q[j + j * rtk->nx] =
2885 rtk->opt.forgetion * rtk->Q[j + j * rtk->nx] +
2886 (1.0 - rtk->opt.forgetion) * SQR(rtk->opt.afgainion) * prtk_ctx.Qp[j + j * prtk_ctx.Qp_nx];
2887 }
2888 }
2889 }
2890 if (rtk->opt.prnadpt && prtk_ctx.Qp) {
2891 int np = NP_RTK(&rtk->opt);
2892 for (i = 0; i < np; i++) {
2893 for (j = 0; j < np; j++) {
2894 if (i < prtk_ctx.Qp_nx && j < prtk_ctx.Qp_nx) {
2895 rtk->Q[i + j * rtk->nx] =
2896 rtk->opt.forgetpva * rtk->Q[i + j * rtk->nx] +
2897 (1.0 - rtk->opt.forgetpva) * SQR(rtk->opt.afgainpva) * prtk_ctx.Qp[i + j * prtk_ctx.Qp_nx];
2898 }
2899 }
2900 }
2901 }
2902
2903 if (stat != SOLQ_NONE) {
2904 /* update state and covariance matrix */
2905 matcpy(rtk->x, xp, rtk->nx, 1);
2906 matcpy(rtk->P, Pp, rtk->nx, rtk->nx);
2907
2908 /* update ambiguity control struct */
2909 rtk->sol.ns = 0;
2910 for (i = 0; i < n; i++) {
2911 sati = obs[i].sat;
2912 for (f = 0; f < nf; f++) {
2913 if (!rtk->ssat[sati - 1].vsat[f]) {
2914 continue;
2915 }
2916 rtk->ssat[sati - 1].lock[f]++;
2917 rtk->ssat[sati - 1].outc[f] = 0;
2918 if (f == 0) {
2919 rtk->sol.ns++;
2920 }
2921 }
2922 }
2923 trace(NULL, 5, "valid sat(L1=%d)\n", rtk->sol.ns);
2924
2925 /* lack of valid satellites */
2926 if (rtk->sol.ns < 4) {
2927 trace(NULL, 2, "lack of valid satellites.\n");
2928 stat = SOLQ_NONE;
2929 }
2930 }
2931
2932 /* keep symmetry of covariance matrix */
2933 for (i = 0; i < rtk->nx; i++) {
2934 for (j = 0; j < i; j++) {
2935 rtk->P[i + j * rtk->nx] = rtk->P[j + i * rtk->nx] =
2936 (rtk->P[i + j * rtk->nx] + rtk->P[j + i * rtk->nx]) * 0.5;
2937 }
2938 }
2939
2940 /* resolve integer ambiguity by LAMBDA */
2941 nb = 0;
2942 pdop = calpdop(rtk, H, azel, vflg, obs, rtk->opt.elmaskar, n, nv);
2943 if (stat != SOLQ_NONE && (rtk->opt.maxpdopar == 0.0 || pdop <= rtk->opt.maxpdopar)) {
2944 /* B2: position variance gate — skip AR if filter not yet converged */
2945 {
2946 double posvar = 0.0;
2947 for (i = 0; i < 3; i++) {
2948 posvar += rtk->P[i + i * rtk->nx];
2949 }
2950 posvar /= 3.0;
2951 if (rtk->opt.thresar[1] > 0.0 && posvar > rtk->opt.thresar[1]) {
2952 trace(NULL, 3, "ppprtk: skip AR (posvar=%.4f>thresar1=%.4f)\n", posvar, rtk->opt.thresar[1]);
2953 nb = 0;
2954 goto ppprtk_ar_done;
2955 }
2956 }
2957 nb = resamb_LAMBDA(rtk, bias, xa);
2958
2959 /* partial AR: exclude satellites to improve fix (posopt[6]=on) */
2960 if (nb <= 1 && rtk->opt.posopt[6]) {
2961 int isat[MAXOBS], exsat = 0, sati, nf2 = NF_RTK(&rtk->opt);
2962 double el[MAXOBS], ratio = 0.0;
2963 elsort(rtk, obs, isat, el, n);
2964 for (l = 0; l < rtk->opt.armaxdelsat; l++) {
2965 ratio = 0.0;
2966 for (i = 0; i < n; i++) {
2967 sati = isat[i];
2968 if (el[i] <= rtk->opt.elmaskar) {
2969 continue;
2970 }
2971 /* skip satellites with >1 frequency missing lock */
2972 for (f = j = 0; f < nf2; f++) {
2973 if (!rtk->ssat[sati - 1].vsat[f] || rtk->ssat[sati - 1].lock[f] <= 0) {
2974 j++;
2975 }
2976 }
2977 if (j > 1) {
2978 continue;
2979 }
2980 /* temporarily exclude this satellite */
2981 for (f = 0; f < nf2; f++) {
2982 rtk->ssat[sati - 1].vsat[f] = 0;
2983 }
2984 pdop = calpdop(rtk, H, azel, vflg, obs, rtk->opt.elmaskar, n, nv);
2985 if ((rtk->opt.maxpdopar == 0.0 || pdop <= rtk->opt.maxpdopar) &&
2986 (nb = resamb_LAMBDA(rtk, bias, xa)) > 1) {
2987 trace(NULL, 2, "PAR OK: excluded sat=%2d el=%.1f nb=%d ratio=%.2f\n", sati, el[i] * R2D, nb,
2988 rtk->sol.ratio);
2989 break;
2990 } else {
2991 /* restore satellite and track best ratio */
2992 for (f = 0; f < nf2; f++) {
2993 rtk->ssat[sati - 1].vsat[f] = 1;
2994 }
2995 if (ratio < rtk->sol.ratio) {
2996 exsat = sati;
2997 ratio = rtk->sol.ratio;
2998 }
2999 }
3000 }
3001 if (nb > 1) {
3002 break;
3003 }
3004 /* permanently exclude satellite with best ratio */
3005 if (exsat > 0) {
3006 for (f = 0; f < nf2; f++) {
3007 rtk->ssat[exsat - 1].vsat[f] = 0;
3008 }
3009 }
3010 }
3011 }
3012
3013 /* B3: arfilter — if AR still failing, exclude newly-locked sats and retry */
3014 if (nb <= 1 && rtk->opt.arfilter) {
3015 int nf2 = NF_RTK(&rtk->opt), rerun = 0, dly = 2, sat;
3016 for (sat = 0; sat < MAXSAT; sat++) {
3017 for (f = 0; f < nf2; f++) {
3018 if (rtk->ssat[sat].vsat[f] && rtk->ssat[sat].lock[f] == 0) {
3019 rtk->ssat[sat].lock[f] = -(rtk->opt.minlock + dly);
3020 dly += 2;
3021 rerun = 1;
3022 }
3023 }
3024 }
3025 if (rerun) {
3026 trace(NULL, 3, "ppprtk arfilter: rerun with newly-locked sats excluded\n");
3027 nb = resamb_LAMBDA(rtk, bias, xa);
3028 }
3029 }
3030 ppprtk_ar_done:;
3031
3032 if (nb > 1) {
3033 /* recompute zdres with fixed solution */
3034 for (ch = 0, nvsum = 0; ch < nch; ch++) {
3035 load_ppp_osr_cpc(&prtk_ctx.osr_ctx, cpc_[ch], pt0_[ch]);
3036 nvtmp[ch] = clas_osr_zdres(obs, n, rs[ch], dts[ch], var[ch], svh[ch], nav, xa, y[ch], e, azel, rtk, 0,
3037 &prtk_ctx.osr_ctx, grid, &clas->current[ch], rtk->ssat, &rtk->opt, &rtk->sol,
3038 NULL, ch);
3039 save_ppp_osr_cpc(&prtk_ctx.osr_ctx, cpc_[ch], pt0_[ch]);
3040 nvsum += nvtmp[ch];
3041 if (ch != 0 && opt->l6mrg) {
3042 pch = nvtmp[0] >= nvtmp[1] ? 0 : 1;
3043 }
3044 }
3045 if (nvsum) {
3046 /* post-fix DD residuals */
3047 nv = ddres(rtk, nav, xa, NULL, NULL, obs, y, e, azel, n, v, NULL, R, vflg, k, pch, TRD);
3048 /* validation of fixed solution */
3049 filter2(rtk, xp, Pp, H, v, R, rtk->nx, nv, vflg, 2);
3050
3051 {
3052 /* chi-square thresholds from rejionno4/5 or defaults */
3053 double thres_hold = rtk->opt.maxinno_ext[3] > 0.0 ? rtk->opt.maxinno_ext[3] : 0.5;
3054 double thres_fix = rtk->opt.maxinno_ext[4] > 0.0 ? rtk->opt.maxinno_ext[4] : 5.0;
3055
3056 if (prtk_ctx.chisq < thres_hold) {
3057 /* hold integer ambiguity */
3058 pdop = calpdop(rtk, H, azel, vflg, obs, rtk->opt.elmaskhold, n, nv);
3059 if (++rtk->nfix >= rtk->opt.minfix && rtk->opt.modear == ARMODE_FIXHOLD &&
3060 (rtk->opt.maxpdophold == 0.0 || pdop <= rtk->opt.maxpdophold)) {
3061 holdamb(rtk, xa);
3062 }
3063 }
3064 if (prtk_ctx.chisq < thres_fix) {
3065 stat = SOLQ_FIX;
3066 }
3067 }
3068 }
3069 }
3070 }
3071
3072 /* reset all states due to large difference between PPP-RTK and SPP positions */
3073 if (opt->maxdiffp > 0.0 && norm(rtk->x, 3) > 0.0 && norm(rtk->sol.rr, 3) > 0.0) {
3074 double dist2 =
3075 SQR(rtk->x[0] - rtk->sol.rr[0]) + SQR(rtk->x[1] - rtk->sol.rr[1]) + SQR(rtk->x[2] - rtk->sol.rr[2]);
3076 if (dist2 > SQR(opt->maxdiffp)) {
3077 ++prtk_ctx.cntdiffp;
3078 if (prtk_ctx.cntdiffp > opt->poserrcnt) {
3079 trace(NULL, 2, "reset all states due to large pos error dist=%.2f\n", sqrt(dist2));
3080 stat = SOLQ_NONE;
3081 for (i = 0; i < 3; i++) {
3082 rtk->x[i] = rtk->sol.rr[i];
3083 rtk->P[i + i * rtk->nx] = (double)rtk->sol.qr[i];
3084 }
3085 prtk_ctx.cntdiffp = 0;
3086 }
3087 } else {
3088 prtk_ctx.cntdiffp = 0;
3089 }
3090 }
3091
3092 /* save solution status */
3093 if (stat == SOLQ_FIX) {
3094 for (i = 0; i < 3; i++) {
3095 rtk->sol.rr[i] = rtk->xa[i];
3096 rtk->sol.qr[i] = (float)rtk->Pa[i + i * rtk->na];
3097 }
3098 rtk->sol.qr[3] = (float)rtk->Pa[1];
3099 rtk->sol.qr[4] = (float)rtk->Pa[1 + 2 * rtk->na];
3100 rtk->sol.qr[5] = (float)rtk->Pa[2];
3101 } else {
3102 for (i = 0; i < 3; i++) {
3103 rtk->sol.rr[i] = rtk->x[i];
3104 rtk->sol.qr[i] = (float)rtk->P[i + i * rtk->nx];
3105 }
3106 rtk->sol.qr[3] = (float)rtk->P[1];
3107 rtk->sol.qr[4] = (float)rtk->P[1 + 2 * rtk->nx];
3108 rtk->sol.qr[5] = (float)rtk->P[2];
3109 rtk->nfix = 0;
3110 }
3111
3112 /* save carrier phase observations */
3113 for (i = 0; i < n; i++) {
3114 for (j = 0; j < nf; j++) {
3115 if (obs[i].L[j] == 0.0) {
3116 continue;
3117 }
3118 rtk->ssat[obs[i].sat - 1].pt[obs[i].rcv - 1][j] = obs[i].time;
3119 rtk->ssat[obs[i].sat - 1].ph[obs[i].rcv - 1][j] = obs[i].L[j];
3120 }
3121 }
3122 /* save SNR */
3123 for (i = 0; i < n; i++) {
3124 for (j = 0; j < nf; j++) {
3125 rtk->ssat[obs[i].sat - 1].snr[j] = obs[i].SNR[j];
3126 }
3127 }
3128 /* update fix/slip status */
3129 for (i = 0; i < MAXSAT; i++) {
3130 for (j = 0; j < nf; j++) {
3131 if (rtk->ssat[i].fix[j] == 2 && stat != SOLQ_FIX) {
3132 rtk->ssat[i].fix[j] = 1;
3133 }
3134 if (rtk->ssat[i].slip[j] & 1) {
3135 rtk->ssat[i].slipc[j] = 1;
3136 }
3137 }
3138 }
3139
3140 if (stat != SOLQ_NONE) {
3141 rtk->sol.stat = stat;
3142 } else {
3143 rtk->sol.stat = SOLQ_SINGLE;
3144 }
3145
3146 /* update consecutive FLOAT counter for filter reset */
3147 if (stat == SOLQ_FIX) {
3148 prtk_ctx.float_count = 0;
3149 } else if (stat == SOLQ_FLOAT) {
3150 prtk_ctx.float_count++;
3151 }
3152
3153 /* reset states on failure */
3154 if (stat == SOLQ_NONE) {
3155 if (rtk->opt.dynamics) {
3156 for (i = 0; i < rtk->nx; i++) {
3157 rtk->x[i] = 0.0;
3158 }
3159 } else {
3160 for (i = NP_RTK(opt); i < rtk->nx; i++) {
3161 rtk->x[i] = 0.0;
3162 }
3163 }
3164 }
3165
3166 /* save persistent zdres state */
3167 for (ch = 0; ch < nch; ch++) {
3168 for (j = 0; j < MAXSAT; j++) {
3169 prtk_ctx.pt0[ch][j] = pt0_[ch][j];
3170 for (f = 0; f < NFREQ; f++) {
3171 prtk_ctx.cpc[ch][f * MAXSAT + j] = cpc_[ch][f * MAXSAT + j];
3172 }
3173 }
3174 }
3175
3176 /* compute DOPs (sol_t has no dop[] field in MRTKLIB, skip) */
3177
3178 free(azel);
3179 free(e);
3180 free(xp);
3181 free(Pp);
3182 free(xa);
3183 free(v);
3184 free(H);
3185 free(R);
3186 free(bias);
3187 for (ch = 0; ch < nch; ch++) {
3188 free(pbslip[ch]);
3189 free(y[ch]);
3190 free(rs[ch]);
3191 free(dts[ch]);
3192 free(var[ch]);
3193 }
3194}
3195
3196/*============================================================================
3197 * Public API: initialization
3198 *===========================================================================*/
3199
3200/**
3201 * @brief Initialize RTK control struct for PPP-RTK positioning.
3202 * @param[in,out] rtk RTK control struct
3203 * @param[in] opt Processing options
3204 */
3205extern void rtkinit_ppp_rtk(rtk_t* rtk, const prcopt_t* opt) {
3206 trace(NULL, 3, "rtkinit_ppp_rtk:\n");
3207 (void)rtk;
3208 (void)opt;
3209
3210 /* zero out file-scope persistent state */
3211 memset(&prtk_ctx, 0, sizeof(prtk_ctx));
3212 prtk_ctx.initialized = 1;
3213}
3214
3215/*============================================================================
3216 * Public API: cleanup
3217 *===========================================================================*/
3218
3219/**
3220 * @brief Free RTK control struct for PPP-RTK positioning.
3221 * @param[in,out] rtk RTK control struct
3222 */
3223extern void rtkfree_ppp_rtk(rtk_t* rtk) {
3224 trace(NULL, 3, "rtkfree_ppp_rtk:\n");
3225 (void)rtk;
3226
3227 /* reset file-scope persistent state */
3228 memset(&prtk_ctx, 0, sizeof(prtk_ctx));
3229}
tropmapf
double tropmapf(gtime_t time, const double pos[], const double azel[], double *mapfw)
Definition mrtk_atmos.c:256
ionmapf
double ionmapf(const double *pos, const double *azel)
Compute ionospheric delay mapping function by single layer model.
Definition mrtk_atmos.c:114
mrtk_atmos.h
MRTKLIB Atmosphere Module — Ionospheric and tropospheric models.
satno2id
void satno2id(int sat, char *id)
Definition mrtk_nav.c:817
code2obs
char * code2obs(uint8_t code)
Definition mrtk_obs.c:190
clas_check_grid_status
void clas_check_grid_status(clas_ctx_t *ctx, const clas_corr_t *corr, int ch)
Check and update grid status flags based on correction availability.
Definition mrtk_clas.c:1320
clas_restore_backup
void clas_restore_backup(clas_ctx_t *ctx, gtime_t time, clas_grid_t *grid, int l6mrg)
Definition mrtk_clas.c:1274
clas_bank_get_close
int clas_bank_get_close(const clas_ctx_t *ctx, gtime_t time, int network, int ch, clas_corr_t *corr)
Merge closest corrections from all banks into a clas_corr_t snapshot.
Definition mrtk_clas.c:1121
SQR
#define SQR(x)
Definition mrtk_clas.c:46
clas_update_global
void clas_update_global(nav_t *nav, const clas_corr_t *corr, int ch)
Apply global corrections (orbit/clock/bias) to nav_t.ssr[].
Definition mrtk_clas.c:1430
clas_backup_valid
int clas_backup_valid(const clas_ctx_t *ctx, gtime_t obstime, int l6mrg)
Save/restore the latest usable CSSR snapshot for short correction gaps.
Definition mrtk_clas.c:1233
clas_update_local
void clas_update_local(nav_t *nav, const clas_corr_t *corr, int ch)
Apply local corrections (trop/STEC) to nav_t.
Definition mrtk_clas.c:1493
clas_backup_current
void clas_backup_current(clas_ctx_t *ctx, const clas_grid_t *grid, int l6mrg)
Definition mrtk_clas.c:1250
clas_clear_current
void clas_clear_current(clas_ctx_t *ctx, int ch)
Drop the current correction snapshot of one channel.
Definition mrtk_clas.c:1266
mrtk_clas.h
CLAS (QZSS L6D) CSSR decoder and grid correction type definitions.
clas_get_grid_index
int clas_get_grid_index(clas_ctx_t *ctx, const double *pos, clas_grid_t *grid, int gridsel, gtime_t obstime)
Select surrounding grid points and compute interpolation weights.
Definition mrtk_clas_grid.c:459
clas_osr_selfreqpair
int clas_osr_selfreqpair(int sat, const prcopt_t *opt, const obsd_t *obs)
Select frequency pair for multi-frequency processing.
Definition mrtk_clas_osr.c:368
clas_osr_zdres
int clas_osr_zdres(const obsd_t *obs, int n, const double *rs, const double *dts, const double *vare, const int *svh, nav_t *nav, double *x, double *y, double *e, double *azel, rtk_t *rtk, int osrlog, clas_osr_ctx_t *osr_ctx, clas_grid_t *grid, clas_corr_t *corr, ssat_t *ssat, prcopt_t *opt, sol_t *sol, clas_osrd_t *osr, int ch)
Compute zero-differenced phase/code residuals with CLAS corrections.
Definition mrtk_clas_osr.c:1080
NP_RTK
#define NP_RTK(opt)
Definition mrtk_clas_osr.c:65
clas_osr_satcorr_update
int clas_osr_satcorr_update(gtime_t time, gtime_t teph, int sat, const prcopt_t *opt, const nav_t *nav, int ch, clas_osr_ctx_t *ctx)
Update per-satellite SIS correction state for IODE adjustment.
Definition mrtk_clas_osr.c:793
NR_RTK
#define NR_RTK(opt)
Definition mrtk_clas_osr.c:72
IB_RTK
#define IB_RTK(s, f, opt)
Definition mrtk_clas_osr.c:73
mrtk_const.h
MRTKLIB Constants — Physical, navigation, and system constants.
EFACT_IRN
#define EFACT_IRN
Definition mrtk_const.h:178
EFACT_CMP
#define EFACT_CMP
Definition mrtk_const.h:175
EFACT_GAL
#define EFACT_GAL
Definition mrtk_const.h:169
EFACT_QZS
#define EFACT_QZS
Definition mrtk_const.h:172
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
SQRT
#define SQRT(x)
Definition mrtk_coords.c:543
ecef2pos
void ecef2pos(const double *r, double *pos)
Transform ECEF position to geodetic position.
Definition mrtk_coords.c:102
ecef2enu
void ecef2enu(const double *pos, const double *r, double *e)
Transform ECEF vector to local tangential coordinate.
Definition mrtk_coords.c:162
mrtk_coords.h
MRTKLIB Coordinates Module — Coordinate transformations and geometry.
satsys
int satsys(int sat, int *prn)
Definition mrtk_nav.c:117
set_ssr_ch_idx
void set_ssr_ch_idx(int ch)
Set SSR channel index for satpos_ssr() (multi-L6E).
Definition mrtk_eph.c:761
SYS_GPS
#define SYS_GPS
Definition mrtk_eph.c:65
SYS_SBS
#define SYS_SBS
Definition mrtk_eph.c:66
SYS_QZS
#define SYS_QZS
Definition mrtk_eph.c:69
satposs
void satposs(gtime_t teph, const obsd_t *obs, int n, const nav_t *nav, int ephopt, double *rs, double *dts, double *var, int *svh)
Compute satellite positions, velocities and clocks.
Definition mrtk_eph.c:940
SYS_IRN
#define SYS_IRN
Definition mrtk_eph.c:71
SYS_CMP
#define SYS_CMP
Definition mrtk_eph.c:70
SYS_GAL
#define SYS_GAL
Definition mrtk_eph.c:68
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
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
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.
freq_num2freq
double freq_num2freq(int sys, int freq_num, int fcn)
Convert frequency number to frequency in Hz.
Definition mrtk_nav.c:680
freq_idx2freq_num
int freq_idx2freq_num(int sys, int freq_idx)
Convert frequency index to frequency number using obsdef table.
Definition mrtk_nav.c:649
mrtk_nav.h
MRTKLIB Navigation Module — Navigation data types and satellite/nav utility functions.
FREQ5
#define FREQ5
Definition mrtk_obs.c:69
FREQ2
#define FREQ2
Definition mrtk_obs.c:68
FREQ8
#define FREQ8
Definition mrtk_obs.c:72
FREQ6
#define FREQ6
Definition mrtk_obs.c:70
R2D
#define R2D
Definition mrtk_obs.c:51
FREQ7
#define FREQ7
Definition mrtk_obs.c:71
FREQ1
#define FREQ1
Definition mrtk_obs.c:67
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
PMODE_PPP_FIXED
#define PMODE_PPP_FIXED
Definition mrtk_opt.h:53
IONOOPT_EST_ADPT
#define IONOOPT_EST_ADPT
Definition mrtk_opt.h:100
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
VAR_VEL
#define VAR_VEL
Definition mrtk_ppp.c:195
VAR_GRA
#define VAR_GRA
Definition mrtk_ppp.c:199
VAR_POS
#define VAR_POS
Definition mrtk_ppp.c:194
VAR_ACC
#define VAR_ACC
Definition mrtk_ppp.c:196
floatcnt_reset
int floatcnt_reset
Definition mrtk_ppp_rtk.c:194
NX_RTK
#define NX_RTK(opt)
Definition mrtk_ppp_rtk.c:104
ddres
static int ddres(rtk_t *rtk, const nav_t *nav, double *x, double **pbslip, const double *P, const obsd_t *obs, double **y, double *e, double *azel, int n, double *v, double *H, double *R, int *vflg, int niter, int pch, int order)
Compute DD residuals and measurement matrix, merging L6 channels.
Definition mrtk_ppp_rtk.c:1597
VARHOLDAMB
#define VARHOLDAMB
Definition mrtk_ppp_rtk.c:76
detslp_ll
static void detslp_ll(rtk_t *rtk, const obsd_t *obs, int n)
Detect cycle slip by LLI flag.
Definition mrtk_ppp_rtk.c:751
ppp_rtk_pos
void ppp_rtk_pos(rtk_t *rtk, const obsd_t *obs, int n, nav_t *nav)
PPP-RTK positioning for one epoch (CLAS).
Definition mrtk_ppp_rtk.c:2540
resamb_LAMBDA
static int resamb_LAMBDA(rtk_t *rtk, double *bias, double *xa)
Resolve integer ambiguity by LAMBDA.
Definition mrtk_ppp_rtk.c:2280
pt0
gtime_t pt0[SSR_CH_NUM][MAXSAT]
Definition mrtk_ppp_rtk.c:168
rtkinit_ppp_rtk
void rtkinit_ppp_rtk(rtk_t *rtk, const prcopt_t *opt)
Initialize RTK control struct for PPP-RTK positioning.
Definition mrtk_ppp_rtk.c:3205
gfmeas_L1L2
static double gfmeas_L1L2(const obsd_t *obs, const nav_t *nav)
L1/L2 geometry-free phase measurement.
Definition mrtk_ppp_rtk.c:772
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
check_clas_facility
static void check_clas_facility(nav_t *nav, const clas_ctx_t *clas, const clas_grid_t *grid, const prcopt_t *opt)
Definition mrtk_ppp_rtk.c:2500
holdamb
static void holdamb(rtk_t *rtk, const double *xa)
Hold fixed ambiguities via pseudomeasurement constraint.
Definition mrtk_ppp_rtk.c:2110
BASELEN
#define BASELEN
Definition mrtk_ppp_rtk.c:71
prtk_ctx
static struct @371154060365014213017353377321306024134051277272 prtk_ctx
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 i, int j, int m, int sati, int *satj, int *sysi, int *sysj, double *lami, double *lamj)
Check that a satellite yields a usable zero-difference signal.
Definition mrtk_ppp_rtk.c:1482
initx
static void initx(rtk_t *rtk, double xi, double var, int i)
Initialize a state and its covariance.
Definition mrtk_ppp_rtk.c:384
varerr
static double varerr(int sat, int sys, double el, int type, const prcopt_t *opt)
Compute measurement error variance.
Definition mrtk_ppp_rtk.c:334
ppp_rtk_na
int ppp_rtk_na(const prcopt_t *opt)
Compute number of non-ambiguity states for PPP-RTK.
Definition mrtk_ppp_rtk.c:2498
sddop
static double sddop(const rtk_t *rtk, const double *H, const int *vflg, double elmask, int nv)
Definition mrtk_ppp_rtk.c:2211
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)
Inner KF measurement update with residual test.
Definition mrtk_ppp_rtk.c:1268
NI_RTK
#define NI_RTK(opt)
Definition mrtk_ppp_rtk.c:100
refsat2
int refsat2[5 *2 *6+1]
Definition mrtk_ppp_rtk.c:163
udpos_ppp
static void udpos_ppp(rtk_t *rtk, const obsd_t *obs, int n, double tt)
Temporal update of position states.
Definition mrtk_ppp_rtk.c:403
save_ppp_osr_cpc
static void save_ppp_osr_cpc(const clas_osr_ctx_t *ctx, double *cpc, gtime_t *pt0)
Definition mrtk_ppp_rtk.c:212
prectrop2
static double prectrop2(gtime_t time, const double *pos, int r, const double *azel, const prcopt_t *opt, const double *x, double *dtdx)
Compute precise tropospheric delay for KF.
Definition mrtk_ppp_rtk.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)
Validate residuals and reject outliers.
Definition mrtk_ppp_rtk.c:1123
chisq
double chisq
Definition mrtk_ppp_rtk.c:181
obs_code_for_freq
static uint8_t obs_code_for_freq(int sat, int f)
Return default observation code for frequency index and system.
Definition mrtk_ppp_rtk.c:235
gfmeas_L1L5
static double gfmeas_L1L5(const obsd_t *obs, const nav_t *nav)
L1/L5 geometry-free phase measurement.
Definition mrtk_ppp_rtk.c:784
filter2
static int filter2(rtk_t *rtk, double *x, double *P, const double *H, const double *v, const double *R, const int n, const int m, const int *vflg, const int flg)
KF measurement update with state compaction.
Definition mrtk_ppp_rtk.c:1370
ppprtk_outstat
int ppprtk_outstat(rtk_t *rtk, char *buff)
Write PPP-RTK solution status to buffer.
Definition mrtk_ppp_rtk.c:2423
udtrop
static void udtrop(rtk_t *rtk, double tt)
Temporal update of tropospheric parameters.
Definition mrtk_ppp_rtk.c:643
II_RTK
#define II_RTK(s, opt)
Definition mrtk_ppp_rtk.c:107
cntdiffp
int cntdiffp
Definition mrtk_ppp_rtk.c:191
LAM_CARR_L5
#define LAM_CARR_L5
Definition mrtk_ppp_rtk.c:84
NP_RTK
#define NP_RTK(opt)
Definition mrtk_ppp_rtk.c:99
check_same_fac
static int check_same_fac(const clas_ctx_t *clas)
Test whether the two L6 channels currently carry the same facility.
Definition mrtk_ppp_rtk.c:1456
MINAMB_DEFAULT
#define MINAMB_DEFAULT
Definition mrtk_ppp_rtk.c:78
ddmat
static int ddmat(rtk_t *rtk, double *D, int *indI, int *indJ)
Build SD-to-DD transformation matrix D.
Definition mrtk_ppp_rtk.c:1977
float_count
int float_count
Definition mrtk_ppp_rtk.c:188
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
NR_RTK
#define NR_RTK(opt)
Definition mrtk_ppp_rtk.c:102
IT_RTK
#define IT_RTK(opt)
Definition mrtk_ppp_rtk.c:108
ddcov
static void ddcov(const int *nb, int n, const double *Ri, const double *Rj, int nv, double *R)
Build DD covariance from SD variances.
Definition mrtk_ppp_rtk.c:1062
calpdop
static double calpdop(const rtk_t *rtk, const double *H, const double *azel, const int *vflg, const obsd_t *obs, double elmask, int n, int nv)
Definition mrtk_ppp_rtk.c:2258
SQR
#define SQR(x)
Definition mrtk_ppp_rtk.c:58
cpc
double cpc[SSR_CH_NUM][MAXSAT *5]
Definition mrtk_ppp_rtk.c:167
restamb
static void restamb(rtk_t *rtk, const double *bias, int nb, double *xa)
Restore SD ambiguities from DD fixed values.
Definition mrtk_ppp_rtk.c:2055
load_ppp_osr_cpc
static void load_ppp_osr_cpc(clas_osr_ctx_t *ctx, const double *cpc, const gtime_t *pt0)
Definition mrtk_ppp_rtk.c:199
zddop
static double zddop(const rtk_t *rtk, const double *azel, const obsd_t *obs, double elmask, int n)
Definition mrtk_ppp_rtk.c:2238
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
MAXOBSLOSS_DEFAULT
#define MAXOBSLOSS_DEFAULT
Definition mrtk_ppp_rtk.c:79
elsort
static void elsort(rtk_t *rtk, const obsd_t *obs, int *isat, double *el, int n)
Sort satellites by elevation angle (ascending) for partial AR.
Definition mrtk_ppp_rtk.c:2178
detslp_gf
static void detslp_gf(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
Detect cycle slip by geometry-free phase jump.
Definition mrtk_ppp_rtk.c:804
chisqr
static const double chisqr[100]
Definition mrtk_ppp_rtk.c:114
udion
static void udion(rtk_t *rtk, double tt, double bl, const obsd_t *obs, int ns)
Temporal update of ionospheric parameters.
Definition mrtk_ppp_rtk.c:679
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
initialized
int initialized
Definition mrtk_ppp_rtk.c:196
use_ch
int use_ch[5][MAXSAT+1]
Definition mrtk_ppp_rtk.c:174
NF_RTK
#define NF_RTK(opt)
Definition mrtk_ppp_rtk.c:98
udbias_ppp
static void udbias_ppp(rtk_t *rtk, double tt, const obsd_t *obs, int n, const nav_t *nav)
Temporal update of phase biases.
Definition mrtk_ppp_rtk.c:856
osr_ctx
clas_osr_ctx_t osr_ctx
Definition mrtk_ppp_rtk.c:178
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
ppp_rtk_nx
int ppp_rtk_nx(const prcopt_t *opt)
Compute number of estimated states for PPP-RTK.
Definition mrtk_ppp_rtk.c:2496
udstate_ppp
static void udstate_ppp(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
Temporal update of all PPP-RTK states.
Definition mrtk_ppp_rtk.c:971
pre_use_ch
int pre_use_ch[5][MAXSAT+1]
Definition mrtk_ppp_rtk.c:175
ALPHAAR_DEFAULT
#define ALPHAAR_DEFAULT
Definition mrtk_ppp_rtk.c:77
sat_lambda
static double sat_lambda(int sat, int f, const nav_t *nav)
Compute satellite-specific wavelength for frequency index.
Definition mrtk_ppp_rtk.c:274
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)
Re-initialize the phase-bias state of one satellite from phase-code.
Definition mrtk_ppp_rtk.c:1523
IB_RTK
#define IB_RTK(s, f, opt)
Definition mrtk_ppp_rtk.c:109
qf
static const double qf[6][60]
Definition mrtk_ppp_rtk.c:125
rtkfree_ppp_rtk
void rtkfree_ppp_rtk(rtk_t *rtk)
Free RTK control struct for PPP-RTK positioning.
Definition mrtk_ppp_rtk.c:3223
Qp_nx
int Qp_nx
Definition mrtk_ppp_rtk.c:185
mrtk_ppp_rtk.h
PPP-RTK positioning engine declarations.
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
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.
clas_corr_t::facility
int facility
Definition mrtk_clas.h:368
clas_corr_t::use
int use
Definition mrtk_clas.h:371
clas_ctx_t
Definition mrtk_clas.h:580
clas_ctx_t::grid
clas_grid_t grid[2]
Definition mrtk_clas.h:592
clas_ctx_t::current
clas_corr_t current[2]
Definition mrtk_clas.h:588
clas_ctx_t::l6buf
clas_l6buf_t l6buf[2]
Definition mrtk_clas.h:608
clas_grid_t
Definition mrtk_clas.h:462
clas_grid_t::network
int network
Definition mrtk_clas.h:467
clas_l6buf_t::time
gtime_t time
Definition mrtk_clas.h:565
clas_osr_ctx_t
Definition mrtk_clas.h:922
clas_osr_ctx_t::cpctmp
double cpctmp[MAXSAT *(5+5)]
Definition mrtk_clas.h:925
clas_osr_ctx_t::pt0tmp
gtime_t pt0tmp[MAXSAT *(5+5)]
Definition mrtk_clas.h:926
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::glo_fcn
int glo_fcn[32]
Definition mrtk_nav.h:426
nav_s::clas_ctx
void * clas_ctx
Definition mrtk_nav.h:439
nav_s::ssr_ch
ssr_t ssr_ch[SSR_CH_NUM][MAXSAT]
Definition mrtk_nav.h:434
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
prcopt_t
Definition mrtk_opt.h:171
prcopt_t::stats_errtrop
double stats_errtrop
Definition mrtk_opt.h:329
prcopt_t::prnadpt
int prnadpt
Definition mrtk_opt.h:261
prcopt_t::niter
int niter
Definition mrtk_opt.h:191
prcopt_t::minfix
int minfix
Definition mrtk_opt.h:185
prcopt_t::qzsmodear
int qzsmodear
Definition mrtk_opt.h:254
prcopt_t::nf
int nf
Definition mrtk_opt.h:174
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::ionoopt
int ionoopt
Definition mrtk_opt.h:187
prcopt_t::armaxdelsat
int armaxdelsat
Definition mrtk_opt.h:256
prcopt_t::stats_prnpositv
double stats_prnpositv
Definition mrtk_opt.h:272
prcopt_t::stats_erriono
double stats_erriono
Definition mrtk_opt.h:328
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::floatcnt
int floatcnt
Definition mrtk_opt.h:257
prcopt_t::eratio
double eratio[5]
Definition mrtk_opt.h:200
prcopt_t::tropopt
int tropopt
Definition mrtk_opt.h:188
prcopt_t::arfilter
int arfilter
Definition mrtk_opt.h:292
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::stats_tconstiono
double stats_tconstiono
Definition mrtk_opt.h:273
prcopt_t::thresar
double thresar[8]
Definition mrtk_opt.h:211
prcopt_t::maxout
int maxout
Definition mrtk_opt.h:183
prcopt_t::maxobsloss_s
double maxobsloss_s
Definition mrtk_opt.h:265
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::minamb
int minamb
Definition mrtk_opt.h:255
prcopt_t::gridsel
int gridsel
Definition mrtk_opt.h:252
prcopt_t::varholdamb
double varholdamb
Definition mrtk_opt.h:263
prcopt_t::forgetion
double forgetion
Definition mrtk_opt.h:266
prcopt_t::ru
double ru[3]
Definition mrtk_opt.h:220
prcopt_t::std
double std[3]
Definition mrtk_opt.h:208
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::maxinno
double maxinno
Definition mrtk_opt.h:217
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::na
int na
Definition mrtk_rtkpos.h:52
rtk_t::P
double * P
Definition mrtk_rtkpos.h:54
rtk_t::xa
double * xa
Definition mrtk_rtkpos.h:56
rtk_t::nfix
int nfix
Definition mrtk_rtkpos.h:57
rtk_t::ssat
ssat_t ssat[MAXSAT]
Definition mrtk_rtkpos.h:60
rtk_t::sol
sol_t sol
Definition mrtk_rtkpos.h:50
rtk_t::Q
double * Q
Definition mrtk_rtkpos.h:55
sol_t::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::age
float age
Definition mrtk_sol.h:56
sol_t::ns
uint8_t ns
Definition mrtk_sol.h:55
sol_t::time
gtime_t time
Definition mrtk_sol.h:45
ssat_t
Definition mrtk_sol.h:98
ssat_t::pt
gtime_t pt[2][5]
Definition mrtk_sol.h:116
ssat_t::vsat
uint8_t vsat[5]
Definition mrtk_sol.h:104
ssat_t::slipc
uint32_t slipc[5]
Definition mrtk_sol.h:111
ssat_t::vs
uint8_t vs
Definition mrtk_sol.h:100
ssat_t::ph
double ph[2][5]
Definition mrtk_sol.h:117
ssat_t::outc
uint32_t outc[5]
Definition mrtk_sol.h:110
ssat_t::rejc
uint32_t rejc[5]
Definition mrtk_sol.h:112
ssat_t::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
ssr_t::t0
gtime_t t0[6]
Definition mrtk_nav.h:176
free
#define free(x)
Definition toml.c:53