MRTKLIB 0.7.14
Modernized GNSS Positioning Library
Loading...
Searching...
No Matches
mrtk_ppp.c
Go to the documentation of this file.
1/*------------------------------------------------------------------------------
2 * mrtk_ppp.c : precise point positioning functions
3 *
4 * Copyright (C) 2026 H.SHIONO (MRTKLIB Project)
5 * Copyright (C) 2023-2025 Cabinet Office, Japan
6 * Copyright (C) 2024-2025 Lighthouse Technology & Consulting Co. Ltd.
7 * Copyright (C) 2023-2025 Japan Aerospace Exploration Agency
8 * Copyright (C) 2023-2025 TOSHIBA ELECTRONIC TECHNOLOGIES CORPORATION
9 * Copyright (C) 2015- Mitsubishi Electric Corp.
10 * Copyright (C) 2014 Geospatial Information Authority of Japan
11 * Copyright (C) 2014 T.SUZUKI
12 * Copyright (C) 2007-2023 T.TAKASU
13 *
14 * SPDX-License-Identifier: BSD-2-Clause
15 *
16 * options : -DIERS_MODEL use IERS tide model
17 * -DOUTSTAT_AMB output ambiguity parameters to solution status
18 *
19 * references :
20 * [1] D.D.McCarthy, IERS Technical Note 21, IERS Conventions 1996, July 1996
21 * [2] D.D.McCarthy and G.Petit, IERS Technical Note 32, IERS Conventions
22 * 2003, November 2003
23 * [3] D.A.Vallado, Fundamentals of Astrodynamics and Applications 2nd ed,
24 * Space Technology Library, 2004
25 * [4] J.Kouba, A Guide to using International GNSS Service (IGS) products,
26 * May 2009
27 * [5] RTCM Paper, April 12, 2010, Proposed SSR Messages for SV Orbit Clock,
28 * Code Biases, URA
29 * [6] MacMillan et al., Atmospheric gradients and the VLBI terrestrial and
30 * celestial reference frames, Geophys. Res. Let., 1997
31 * [7] G.Petit and B.Luzum (eds), IERS Technical Note No. 36, IERS
32 * Conventions (2010), 2010
33 * [8] J.Kouba, A simplified yaw-attitude model for eclipsing GPS satellites,
34 * GPS Solutions, 13:1-12, 2009
35 * [9] F.Dilssner, GPS IIF-1 satellite antenna phase center and attitude
36 * modeling, InsideGNSS, September, 2010
37 * [10] F.Dilssner, The GLONASS-M satellite yaw-attitude model, Advances in
38 * Space Research, 2010
39 * [11] IGS MGEX (http://igs.org/mgex)
40 * [12] CAO IS-QZSS-MDC-004, TBD, 2025
41 *
42 * version : $Revision:$ $Date:$
43 * history : 2010/07/20 1.0 new
44 * added api:
45 * tidedisp()
46 * 2010/12/11 1.1 enable exclusion of eclipsing satellite
47 * 2012/02/01 1.2 add gps-glonass h/w bias correction
48 * move windupcorr() to rtkcmn.c
49 * 2013/03/11 1.3 add otl and pole tides corrections
50 * involve iers model with -DIERS_MODEL
51 * change initial variances
52 * suppress acos domain error
53 * 2013/09/01 1.4 pole tide model by iers 2010
54 * add mode of ionosphere model off
55 * 2014/05/23 1.5 add output of trop gradient in solution status
56 * 2014/10/13 1.6 fix bug on P0(a[3]) computation in tide_oload()
57 * fix bug on m2 computation in tide_pole()
58 * 2015/03/19 1.7 fix bug on ionosphere correction for GLO and BDS
59 * 2015/05/10 1.8 add function to detect slip by MW-LC jump
60 * fix ppp solutin problem with large clock variance
61 * 2015/06/08 1.9 add precise satellite yaw-models
62 * cope with day-boundary problem of satellite clock
63 * 2015/07/31 1.10 fix bug on nan-solution without glonass nav-data
64 * pppoutsolsat() -> pppoutstat()
65 * 2015/11/13 1.11 add L5-receiver-dcb estimation
66 * merge post-residual validation by rnx2rtkp_test
67 * support support option opt->pppopt=-GAP_RESION=nnnn
68 * 2016/01/22 1.12 delete support for yaw-model bug
69 * add support for ura of ephemeris
70 * 2018/10/10 1.13 support api change of satexclude()
71 * 2020/11/30 1.14 use sat2freq() to get carrier frequency
72 * use E1-E5b for Galileo iono-free LC
73 * 2023/02/01 1.15 branch from ver.2.4.3b34 for MADOCALIB
74 * 2023/03/08 1.16 change for suppress -Wstringop-overflow warnings.
75 * 2024/01/10 1.17 support MADOCA-PPP ionospheric corrections, add
76 * udbsys_ppp(), NM, IM
77 * change udstate_ppp(), ppp_res()
78 * delete L5-receiver-dcb estimation, uddcb_ppp(), ND, ID
79 * fix bug on L1 iono computation in udbias_ppp()
80 * considered signal selection for PCO and PCV, add
81 * signal_sel_pco(), signal_sel_pcv()
82 * 2024/03/15 1.18 change the value of $ION from vertical to slant.
83 * 2024/05/24 1.19 estimate receiver clock of QZS.
84 * delete processing regarding MADOCA-PPP ionospheric
85 * corrections system bias.
86 * delete udbsys_ppp().
87 * 2024/06/06 1.20 coonsider URA ratio in the calculation of the
88 * variance of phase and code residuals.
89 * 2024/07/23 1.21 support cycle slip detection by ssr phase bias
90 * discontinuity indicator.
91 * apply phase bias correction only when performing
92 * ambiguity resolution.
93 * change the sign of the phase bias correction to
94 * correspond to IS-QZSS-MDC (ref.[12]).
95 * 2024/09/27 1.22 delete NM,IM
96 * 2024/11/13 1.23 apply phase bias correction even when ambiguity
97 * resolution is not performed.
98 * 2025/03/18 1.24 support triple/quad-frequency ppp.
99 * support bds-2 and bds-3.
100 *-----------------------------------------------------------------------------*/
101#include "mrtklib/mrtk_ppp.h"
102
103#include <math.h>
104#include <stdio.h>
105#include <stdlib.h>
106#include <string.h>
107
108#include "mrtklib/mrtk_antenna.h"
109#include "mrtklib/mrtk_astro.h"
110#include "mrtklib/mrtk_atmos.h"
111#include "mrtklib/mrtk_bias_sinex.h"
112#include "mrtklib/mrtk_const.h"
113#include "mrtklib/mrtk_context.h"
114#include "mrtklib/mrtk_coords.h"
115#include "mrtklib/mrtk_has.h"
116#include "mrtklib/mrtk_ionex.h"
117#include "mrtklib/mrtk_madoca.h"
118#include "mrtklib/mrtk_mat.h"
119#include "mrtklib/mrtk_peph.h"
120#include "mrtklib/mrtk_ppp_ar.h"
121#include "mrtklib/mrtk_ppp_rtk.h"
122#include "mrtklib/mrtk_sbas.h"
123#include "mrtklib/mrtk_tides.h"
124#include "mrtklib/mrtk_trace.h"
125
126/* forward declarations */
127extern int satsys(int sat, int* prn);
128extern int satsys_bd2(int sat, int* prn);
129extern void satno2id(int sat, char* id);
130extern int satid2no(const char* id);
131extern double sat2freq(int sat, uint8_t code, const nav_t* nav);
132extern int satexclude(int sat, double var, int svh, const prcopt_t* opt);
133extern char* time_str(gtime_t t, int n);
134extern void time2str(gtime_t t, char* s, int n);
135extern double time2gpst(gtime_t t, int* week);
136extern void sunmoonpos(gtime_t tut, const double* erpv, double* rsun, double* rmoon, double* gmst);
137extern double geodist(const double* rs, const double* rr, double* e);
138extern double satazel(const double* pos, const double* e, double* azel);
139extern void satantoff(gtime_t time, const double* rs, int sat, const nav_t* nav, double* dant);
140extern void ecef2pos(const double* r, double* pos);
141extern void ecef2enu(const double* pos, const double* r, double* e);
142extern void xyz2enu(const double* pos, double* E);
143extern void covecef(const double* pos, const double* Q, double* P);
144extern void covenu(const double* pos, const double* P, double* Q);
145extern double dot(const double* a, const double* b, int n);
146extern double norm(const double* a, int n);
147extern int normv3(const double* a, double* b);
148extern void cross3(const double* a, const double* b, double* c);
149extern void matcpy(double* A, const double* B, int n, int m);
150extern void matmul(const char* tr, int n, int k, int m, double alpha, const double* A, const double* B, double beta,
151 double* C);
152extern double* mat(int n, int m);
153extern int* imat(int n, int m);
154extern double* zeros(int n, int m);
155extern double* eye(int n);
156extern int filter(double* x, double* P, const double* H, const double* v, const double* R, int n, int m);
157extern double sbstropcorr(gtime_t time, const double* pos, const double* azel, double* var);
158extern double ionmodel(gtime_t t, const double* ion, const double* pos, const double* azel);
159extern int iontec(gtime_t time, const nav_t* nav, const double* pos, const double* azel, int opt, double* delay,
160 double* var);
161extern int sbsioncorr(gtime_t time, const nav_t* nav, const double* pos, const double* azel, double* delay,
162 double* var);
163extern double tropmodel(gtime_t time, const double* pos, const double* azel, double humi);
164extern double tropmapf(gtime_t time, const double* pos, const double* azel, double* mapfw);
165extern int testsnr(int base, int freq, double el, double snr, const snrmask_t* mask);
166extern void antmodel(const pcv_t* pcv, const double* del, const double* azel, int opt, double* dant);
167extern void antmodel_s(const pcv_t* pcv, double nadir, double* dant);
168extern void satposs(gtime_t teph, const obsd_t* obs, int n, const nav_t* nav, int ephopt, double* rs, double* dts,
169 double* var, int* svh);
170extern gtime_t gpst2utc(gtime_t t);
171extern void tidedisp(gtime_t tutc, const double* rr, int opt, const erp_t* erp, const double* odisp, double* dr);
172extern unsigned int tickget(void);
173extern char* code2obs(uint8_t code);
174extern int mcssr_sel_biascode(const int sys, const int code);
175/* freq_idx2freq_num, freq_num2freq, freq_idx2ant_idx declared in mrtk_nav.h */
176extern int const_iono_corr(rtk_t* rtk, const obsd_t* obs, const nav_t* nav, const int n, const double* azel,
177 const int* exc, const double* rr, const double* x, double* v, double* H, double* var);
178
179#define SQR(x) ((x) * (x))
180#define SQRT(x) ((x) <= 0.0 || (x) != (x) ? 0.0 : sqrt(x))
181#define MAX(x, y) ((x) > (y) ? (x) : (y))
182#define MIN(x, y) ((x) < (y) ? (x) : (y))
183#define ROUND(x) (int)floor((x) + 0.5)
184
185#define PPP_NFREQ 4 /* number of freq. for PPP */
186#define SOLQ_FIX_WL 99 /* solution status: wide-lane fixed */
187
188#define MAX_ITER 8 /* max number of iterations */
189#define MIN_NSAT_SOL 4 /* min satellite number for solution */
190#define THRES_REJECT 4.0 /* reject threshold of posfit-res (sigma) */
191
192#define THRES_MW_JUMP 10.0
193
194#define VAR_POS SQR(60.0) /* init variance receiver position (m^2) */
195#define VAR_VEL SQR(10.0) /* init variance of receiver vel ((m/s)^2) */
196#define VAR_ACC SQR(10.0) /* init variance of receiver acc ((m/ss)^2) */
197#define VAR_CLK SQR(60.0) /* init variance receiver clock (m^2) */
198#define VAR_ZTD SQR(0.12) /* init variance ztd (m^2) */
199#define VAR_GRA SQR(0.01) /* init variance gradient (m^2) */
200#define VAR_IFB SQR(30.0) /* init variance ifb (m^2) */
201#define VAR_BIAS SQR(60.0) /* init variance phase-bias (m^2) */
202#define VAR_IONO SQR(60.0) /* init variance iono-delay */
203#define VAR_GLO_IFB SQR(0.6) /* variance of glonass ifb */
204
205#define ERR_SAAS 0.3 /* saastamoinen model error std (m) */
206#define ERR_BRDCI 0.5 /* broadcast iono model error factor */
207#define REL_HUMI 0.7 /* relative humidity for saastamoinen model */
208#define GAP_RESION 120 /* default gap to reset ionos parameters (ep) */
209
210/* number and index of states */
211#define NF(opt) ((opt)->ionoopt == IONOOPT_IFLC ? 1 : (opt)->nf)
212#define NP(opt) ((opt)->dynamics ? 9 : 3)
213#define NC(opt) (NSYS + 1) /* BDS3 and BDS2 */
214static inline int NT(const prcopt_t* opt) {
215 if (opt->tropopt < TROPOPT_EST) return 0;
216 if (opt->tropopt == TROPOPT_EST) return 1;
217 return 3;
218}
219#define NI(opt) ((opt)->ionoopt == IONOOPT_EST ? MAXSAT : 0)
220#define N3(opt) (NF(opt) >= 3 ? NC(opt) : 0) /* 3rd-frequency receivre bias */
221#define N4(opt) (NF(opt) >= 4 ? NC(opt) : 0) /* 4th-frequency receivre bias */
222#define NR(opt) (NP(opt) + NC(opt) + NT(opt) + NI(opt) + N3(opt) + N4(opt))
223#define NB(opt) (NF(opt) * MAXSAT)
224#define NX(opt) (NR(opt) + NB(opt))
225#define IC(s, opt) (NP(opt) + (s))
226#define IT(opt) (NP(opt) + NC(opt))
227#define II(s, opt) (NP(opt) + NC(opt) + NT(opt) + (s) - 1)
228#define I3(s, opt) (NP(opt) + NC(opt) + NT(opt) + NI(opt) + (s))
229#define I4(s, opt) (NP(opt) + NC(opt) + NT(opt) + NI(opt) + N3(opt) + (s))
230#define IB(s, f, opt) (NR(opt) + MAXSAT * (f) + (s) - 1)
231
232/* standard deviation of state -----------------------------------------------*/
233static double STD(rtk_t* rtk, int i) {
234 if (rtk->sol.stat == SOLQ_FIX) {
235 return SQRT(rtk->Pa[i + i * rtk->nx]);
236 }
237 return SQRT(rtk->P[i + i * rtk->nx]);
238}
239/* write solution status for PPP ---------------------------------------------*/
240extern int pppoutstat(rtk_t* rtk, char* buff) {
241 ssat_t* ssat;
242 double tow, pos[3], vel[3], acc[3], *x;
243 int i, j, week;
244 char id[32], *p = buff;
245
246 if (!rtk->sol.stat) {
247 return 0;
248 }
249
250 /* PPP-RTK uses different state layout (no clocks, different indices) */
251 if (rtk->opt.mode == PMODE_PPP_RTK) {
252 return ppprtk_outstat(rtk, buff);
253 }
254
255 trace(NULL, 3, "pppoutstat:\n");
256
257 tow = time2gpst(rtk->sol.time, &week);
258
259 x = rtk->sol.stat == SOLQ_FIX ? rtk->xa : rtk->x;
260
261 /* receiver position */
262 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],
263 STD(rtk, 0), STD(rtk, 1), STD(rtk, 2));
264
265 /* receiver velocity and acceleration */
266 if (rtk->opt.dynamics) {
267 ecef2pos(rtk->sol.rr, pos);
268 ecef2enu(pos, rtk->x + 3, vel);
269 ecef2enu(pos, rtk->x + 6, acc);
270 p += sprintf(p,
271 "$VELACC,%d,%.3f,%d,%.4f,%.4f,%.4f,%.5f,%.5f,%.5f,%.4f,%.4f,"
272 "%.4f,%.5f,%.5f,%.5f\n",
273 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,
274 0.0);
275 }
276 /* receiver clocks */
277 i = IC(0, &rtk->opt);
278 p += sprintf(p, "$CLK,%d,%.3f,%d,%d,%.3f,%.3f,%.3f,%.3f\n", week, tow, rtk->sol.stat, 1, x[i] * 1E9 / CLIGHT,
279 x[i + 1] * 1E9 / CLIGHT, STD(rtk, i) * 1E9 / CLIGHT, STD(rtk, i + 1) * 1E9 / CLIGHT);
280
281 /* tropospheric parameters */
282 if (rtk->opt.tropopt == TROPOPT_EST || rtk->opt.tropopt == TROPOPT_ESTG) {
283 i = IT(&rtk->opt);
284 p += sprintf(p, "$TROP,%d,%.3f,%d,%d,%.4f,%.4f\n", week, tow, rtk->sol.stat, 1, x[i], STD(rtk, i));
285 }
286 if (rtk->opt.tropopt == TROPOPT_ESTG) {
287 i = IT(&rtk->opt);
288 p += sprintf(p, "$TRPG,%d,%.3f,%d,%d,%.5f,%.5f,%.5f,%.5f\n", week, tow, rtk->sol.stat, 1, x[i + 1], x[i + 2],
289 STD(rtk, i + 1), STD(rtk, i + 2));
290 }
291 /* ionosphere parameters */
292 if (rtk->opt.ionoopt == IONOOPT_EST) {
293 ecef2pos(rtk->sol.rr, pos);
294 for (i = 0; i < MAXSAT; i++) {
295 ssat = rtk->ssat + i;
296 if (!ssat->vs) {
297 continue;
298 }
299 j = II(i + 1, &rtk->opt);
300 if (rtk->x[j] == 0.0) {
301 continue;
302 }
303 satno2id(i + 1, id);
304 p += sprintf(p, "$ION,%d,%.3f,%d,%s,%.1f,%.1f,%.4f,%.4f\n", week, tow, rtk->sol.stat, id,
305 rtk->ssat[i].azel[0] * R2D, rtk->ssat[i].azel[1] * R2D, x[j], STD(rtk, j));
306 }
307 }
308#ifdef OUTSTAT_AMB
309 /* ambiguity parameters */
310 for (i = 0; i < MAXSAT; i++)
311 for (j = 0; j < NF(&rtk->opt); j++) {
312 k = IB(i + 1, j, &rtk->opt);
313 if (rtk->x[k] == 0.0) continue;
314 satno2id(i + 1, id);
315 p +=
316 sprintf(p, "$AMB,%d,%.3f,%d,%s,%d,%.4f,%.4f\n", week, tow, rtk->sol.stat, id, j + 1, x[k], STD(rtk, k));
317 }
318#endif
319 return (int)(p - buff);
320}
321/* exclude meas of eclipsing satellite (block IIA) ---------------------------*/
322static void testeclipse(const obsd_t* obs, int n, const nav_t* nav, double* rs) {
323 double rsun[3], esun[3], r, ang, erpv[5] = {0}, cosa;
324 int i, j;
325 const char* type;
326
327 trace(NULL, 3, "testeclipse:\n");
328
329 /* unit vector of sun direction (ecef) */
330 sunmoonpos(gpst2utc(obs[0].time), erpv, rsun, NULL, NULL);
331 normv3(rsun, esun);
332
333 for (i = 0; i < n; i++) {
334 type = nav->pcvs[obs[i].sat - 1].type;
335
336 if ((r = norm(rs + i * 6, 3)) <= 0.0) {
337 continue;
338 }
339
340 /* only block IIA */
341 if (*type && !strstr(type, "BLOCK IIA")) {
342 continue;
343 }
344
345 /* sun-earth-satellite angle */
346 cosa = dot(rs + i * 6, esun, 3) / r;
347 cosa = cosa < -1.0 ? -1.0 : (cosa > 1.0 ? 1.0 : cosa);
348 ang = acos(cosa);
349
350 /* test eclipse */
351 if (ang < PI / 2.0 || r * sin(ang) > RE_WGS84) {
352 continue;
353 }
354
355 trace(NULL, 3, "eclipsing sat excluded %s sat=%2d\n", time_str(obs[0].time, 0), obs[i].sat);
356
357 for (j = 0; j < 3; j++) {
358 rs[j + i * 6] = 0.0;
359 }
360 }
361}
362/* nominal yaw-angle ---------------------------------------------------------*/
363static double yaw_nominal(double beta, double mu) {
364 if (fabs(beta) < 1E-12 && fabs(mu) < 1E-12) {
365 return PI;
366 }
367 return atan2(-tan(beta), sin(mu)) + PI;
368}
369/* yaw-angle of satellite ----------------------------------------------------*/
370extern int yaw_angle(int sat, const char* type, int opt, double beta, double mu, double* yaw) {
371 *yaw = yaw_nominal(beta, mu);
372 return 1;
373}
374/* satellite attitude model --------------------------------------------------*/
375static int sat_yaw(gtime_t time, int sat, const char* type, int opt, const double* rs, double* exs, double* eys) {
376 double rsun[3], ri[6], es[3], esun[3], n[3], p[3], en[3], ep[3], ex[3], E, beta, mu;
377 double yaw, cosy, siny, erpv[5] = {0};
378 int i;
379
380 sunmoonpos(gpst2utc(time), erpv, rsun, NULL, NULL);
381
382 /* beta and orbit angle */
383 matcpy(ri, rs, 6, 1);
384 ri[3] -= OMGE * ri[1];
385 ri[4] += OMGE * ri[0];
386 cross3(ri, ri + 3, n);
387 cross3(rsun, n, p);
388 if (!normv3(rs, es) || !normv3(rsun, esun) || !normv3(n, en) || !normv3(p, ep)) {
389 return 0;
390 }
391 beta = PI / 2.0 - acos(dot(esun, en, 3));
392 E = acos(dot(es, ep, 3));
393 mu = PI / 2.0 + (dot(es, esun, 3) <= 0 ? -E : E);
394 if (mu < -PI / 2.0) {
395 mu += 2.0 * PI;
396 } else if (mu >= PI / 2.0) {
397 mu -= 2.0 * PI;
398 }
399
400 /* yaw-angle of satellite */
401 if (!yaw_angle(sat, type, opt, beta, mu, &yaw)) {
402 return 0;
403 }
404
405 /* satellite fixed x,y-vector */
406 cross3(en, es, ex);
407 cosy = cos(yaw);
408 siny = sin(yaw);
409 for (i = 0; i < 3; i++) {
410 exs[i] = -siny * en[i] + cosy * ex[i];
411 eys[i] = -cosy * en[i] - siny * ex[i];
412 }
413 return 1;
414}
415/* phase windup model --------------------------------------------------------*/
416static int model_phw(gtime_t time, int sat, const char* type, int opt, const double* rs, const double* rr,
417 double* phw) {
418 double exs[3], eys[3], ek[3], exr[3], eyr[3], eks[3], ekr[3], E[9];
419 double dr[3], ds[3], drs[3], r[3], pos[3], cosp, ph;
420 int i;
421
422 if (opt <= 0) {
423 return 1; /* no phase windup */
424 }
425
426 /* satellite yaw attitude model */
427 if (!sat_yaw(time, sat, type, opt, rs, exs, eys)) {
428 return 0;
429 }
430
431 /* unit vector satellite to receiver */
432 for (i = 0; i < 3; i++) {
433 r[i] = rr[i] - rs[i];
434 }
435 if (!normv3(r, ek)) {
436 return 0;
437 }
438
439 /* unit vectors of receiver antenna */
440 ecef2pos(rr, pos);
441 xyz2enu(pos, E);
442 exr[0] = E[1];
443 exr[1] = E[4];
444 exr[2] = E[7]; /* x = north */
445 eyr[0] = -E[0];
446 eyr[1] = -E[3];
447 eyr[2] = -E[6]; /* y = west */
448
449 /* phase windup effect */
450 cross3(ek, eys, eks);
451 cross3(ek, eyr, ekr);
452 for (i = 0; i < 3; i++) {
453 ds[i] = exs[i] - ek[i] * dot(ek, exs, 3) - eks[i];
454 dr[i] = exr[i] - ek[i] * dot(ek, exr, 3) + ekr[i];
455 }
456 cosp = dot(ds, dr, 3) / norm(ds, 3) / norm(dr, 3);
457 if (cosp < -1.0) {
458 cosp = -1.0;
459 } else if (cosp > 1.0) {
460 cosp = 1.0;
461 }
462 ph = acos(cosp) / 2.0 / PI;
463 cross3(ds, dr, drs);
464 if (dot(ek, drs, 3) < 0.0) {
465 ph = -ph;
466 }
467
468 *phw = ph + floor(*phw - ph + 0.5); /* in cycle */
469 return 1;
470}
471/* measurement error variance ------------------------------------------------*/
472static double varerr(int sat, int sys, double el, int idx, int type, const prcopt_t* opt) {
473 double fact = 1.0, sinel = sin(el);
474
475 if (type == 1) {
476 fact *= opt->eratio[idx == 0 ? 0 : 1];
477 }
478 fact *= sys == SYS_GLO
479 ? EFACT_GLO
480 : (sys == SYS_GAL
481 ? EFACT_GAL
482 : (sys == SYS_QZS
483 ? EFACT_QZS
484 : (sys == SYS_CMP
485 ? EFACT_CMP
486 : (sys == SYS_IRN ? EFACT_IRN : (sys == SYS_SBS ? EFACT_SBS : EFACT_GPS)))));
487 if (opt->ionoopt == IONOOPT_IFLC) {
488 fact *= 3.0;
489 }
490 return SQR(fact * opt->err[1]) + SQR(fact * opt->err[2] / sinel);
491}
492/* initialize state and covariance -------------------------------------------*/
493static void initx(rtk_t* rtk, double xi, double var, int i) {
494 int j;
495 rtk->x[i] = xi;
496 for (j = 0; j < rtk->nx; j++) {
497 rtk->P[i + j * rtk->nx] = rtk->P[j + i * rtk->nx] = i == j ? var : 0.0;
498 }
499}
500/* geometry-free phase measurement -------------------------------------------*/
501static double gfmeas(const obsd_t* obs, const nav_t* nav, int f) {
502 double freq1, freq2;
503
504 freq1 = sat2freq(obs->sat, obs->code[0], nav);
505 freq2 = sat2freq(obs->sat, obs->code[f], nav);
506 if (freq1 == 0.0 || freq2 == 0.0 || obs->L[0] == 0.0 || obs->L[f] == 0.0) {
507 return 0.0;
508 }
509 return (obs->L[0] / freq1 - obs->L[f] / freq2) * CLIGHT;
510}
511/* Melbourne-Wubbena linear combination --------------------------------------*/
512static double mwmeas(const obsd_t* obs, const nav_t* nav, int f) {
513 double freq1, freq2;
514
515 freq1 = sat2freq(obs->sat, obs->code[0], nav);
516 freq2 = sat2freq(obs->sat, obs->code[f], nav);
517
518 if (freq1 == 0.0 || freq2 == 0.0 || obs->L[0] == 0.0 || obs->L[f] == 0.0 || obs->P[0] == 0.0 || obs->P[f] == 0.0) {
519 return 0.0;
520 }
521 return (obs->L[0] - obs->L[f]) * CLIGHT / (freq1 - freq2) -
522 (freq1 * obs->P[0] + freq2 * obs->P[f]) / (freq1 + freq2);
523}
524/* notify the user once per system when IGS-PPP auto-selects a non-nominal 2nd band
525 * for the iono-free combination (#135). trace() is a no-op in release builds, so the
526 * primary channel is stderr. Per-process state is adequate for the post path. */
527static void ppp_note_autosig(int sys, uint8_t code) {
528 static int reported = 0; /* bitmask of SYS_* already notified */
529 const char* name;
530
531 if (reported & sys) {
532 return;
533 }
534 reported |= sys;
535
536 switch (sys) {
537 case SYS_GPS:
538 name = "GPS";
539 break;
540 case SYS_GLO:
541 name = "GLO";
542 break;
543 case SYS_GAL:
544 name = "GAL";
545 break;
546 case SYS_QZS:
547 name = "QZS";
548 break;
549 case SYS_CMP:
550 name = "BDS";
551 break;
552 case SYS_IRN:
553 name = "IRN";
554 break;
555 default:
556 name = "???";
557 break;
558 }
559 fprintf(stderr, "note: auto signal: %s iono-free 2nd freq = L%s (nominal band absent in data)\n", name,
560 code2obs(code));
561 trace(NULL, 1, "corr_meas: auto IFLC 2nd band sys=%s code=L%s\n", name, code2obs(code));
562}
563/* antenna corrected measurements --------------------------------------------*/
564static int corr_meas(const obsd_t* obs, const nav_t* nav, const double* azel, const prcopt_t* opt, const double* dantr,
565 const double* dants, double phw, double* L, double* P, double* Lc, double* Pc, int tl) {
566 char satid[8], *tstr = time_str(obs->time, 3);
567 double freq[NFREQ] = {0}, C1, C2, cb = 0.0, pb = 0.0;
568 int i, i2 = 1, nfread, sys = satsys(obs->sat, NULL), ssrcode = CODE_NONE, ant_idx;
569 int ssr_bias = 0;
570 int corr = opt->correction; /* correction source (CORR_???); selects the bias path below */
571
572 if (strstr(opt->pppopt, "-SSR_BIAS")) {
573 ssr_bias = 1;
574 }
575
576 satno2id(obs->sat, satid);
577 /* #135: for IGS-product PPP, prepare all main frequency slots (not just nf) so a
578 * receiver that delivers its 2nd band in slot 2 (e.g. GAL E5b / BDS B2I, with
579 * slot-1 E5a/B3I absent) still yields a usable second observation for the
580 * iono-free pair below. Other correction sources keep the nf-bounded loop. */
581 nfread = (corr == CORR_IGS && NFREQ > opt->nf) ? NFREQ : opt->nf;
582 for (i = 0; i < nfread && i < NFREQ; i++) {
583 L[i] = P[i] = cb = pb = 0.0;
584 freq[i] = sat2freq(obs->sat, obs->code[i], nav);
585
586 if (freq[i] == 0.0 || obs->L[i] == 0.0 || obs->P[i] == 0.0) {
587 continue;
588 }
589 if (testsnr(0, 0, azel[1], obs->SNR[i] * SNR_UNIT, &opt->snrmask)) {
590 continue;
591 }
592
593 /* antenna index */
594 ant_idx = freq_idx2ant_idx(satsys_bd2(obs->sat, NULL), i);
595
596 /* antenna phase center and phase windup correction */
597 L[i] = obs->L[i] * CLIGHT / freq[i] - dants[ant_idx] - dantr[ant_idx] - phw * CLIGHT / freq[i];
598 P[i] = obs->P[i] - dants[ant_idx] - dantr[ant_idx];
599
600 if (corr == CORR_IGS || corr == CORR_IGS_RTS) {
601 /* RTKLIB-style float PPP measurement model, shared by IGS precise
602 * products (files, correction=igs) and IGS-RTS RTCM-SSR streams
603 * (correction=igs-rts, #138). The satellite orbit/clock already come
604 * from satpos_ssr (EPHOPT_SSRAPC) for igs-rts; here only the receiver
605 * measurement model differs from MADOCA. Apply optional P1-C1 / P2-C2
606 * DCB from nav->cbias (GPS/GLO) if present; never discard a measurement
607 * for a missing satellite bias. The satellite phase bias is absorbed by
608 * the float phase-ambiguity state. The per-signal SSR code/phase bias
609 * (MADOCA convention) is deliberately NOT applied here: RTCM-SSR uses a
610 * different bias convention and applying it degrades the solution
611 * (verified against upstream RTKLIB rnx2rtkp). */
612 if (corr == CORR_IGS && opt->ionoopt == IONOOPT_EST) {
613 /* #142: uncombined IGS PPP-AR. Apply the per-signal absolute
614 * Observable-Specific Bias (OSB) from the Bias-SINEX product
615 * (nav->osb, populated by udsatcb()/udsatpb()): satellite code
616 * bias to the pseudorange and satellite phase bias to the carrier
617 * (both m, keyed by obs->code[]). Removing the satellite phase
618 * bias makes the float ambiguity integer-recoverable so
619 * ppp_amb_ILS can fix WL/NL. Gated on uncombined mode so the
620 * float iono-free path (correction=igs, ionosphere=dual-freq)
621 * stays bit-identical (#142 before/after guard). Sign: nav->osb
622 * stores the OSB already negated (udsatcb/udsatpb apply *=-1), so
623 * adding it here yields the formal observed-minus-OSB correction
624 * (same consumer convention as the OSR engine). */
625 int c = obs->code[i] - 1;
626 if (c >= 0) {
627 if (nav->osb.vspb[obs->sat - 1][c]) {
628 L[i] += nav->osb.spb[obs->sat - 1][c];
629 }
630 if (nav->osb.vscb[obs->sat - 1][c]) {
631 P[i] += nav->osb.scb[obs->sat - 1][c];
632 } else if (opt->unbias) {
633 P[i] = 0.0;
634 trace(NULL, tl > 0 ? 3 : 4,
635 "corr_meas: %s satellite code bias does not exist. %s obscode=C%s\n", tstr, satid,
636 code2obs(obs->code[i]));
637 continue;
638 }
639 /* Apply a receiver-specific OSB when available, otherwise
640 * its system-wide OSB. With unbias enabled, a usable code
641 * observation must have both satellite and receiver OSBs.
642 * This gates P[i] only -- L[i] (phase) already carries the
643 * satellite phase bias above and must not be dropped here. */
644 {
645 int sysno = getsysno(obs->sat);
646
647 if (nav->osb.vrsatcb[obs->sat - 1][c]) {
648 P[i] += nav->osb.rsatcb[obs->sat - 1][c];
649 } else if (sysno >= 0 && nav->osb.vrsyscb[sysno][c]) {
650 P[i] += nav->osb.rsyscb[sysno][c];
651 } else if (opt->unbias) {
652 P[i] = 0.0;
653 trace(NULL, tl > 0 ? 3 : 4,
654 "corr_meas: %s receiver code bias does not exist. %s obscode=C%s\n", tstr, satid,
655 code2obs(obs->code[i]));
656 continue;
657 }
658 }
659 } else if (opt->unbias) {
660 P[i] = 0.0;
661 trace(NULL, tl > 0 ? 3 : 4, "corr_meas: %s invalid observation code. %s\n", tstr, satid);
662 }
663 } else if (sys == SYS_GPS || sys == SYS_GLO) {
664 if (obs->code[i] == CODE_L1C) {
665 P[i] += nav->cbias[obs->sat - 1][1];
666 }
667 if (obs->code[i] == CODE_L2C) {
668 P[i] += nav->cbias[obs->sat - 1][2];
669 }
670 }
671 } else if (corr == CORR_NONE) {
672 /* no satellite bias correction (SPP/DGPS/RTK route through pntpos/relpos,
673 * not here; branch kept for completeness) */
674 } else {
675 /* SSR path (qzs-madoca / gal-has / bds-b2b; also CORR_AUTO fallback).
676 * Satellite cbias/pbias come from nav->ssr_ch; a measurement with no
677 * valid SSR bias is dropped (required for PPP-AR, MADOCALIB design).
678 * Note: igs-rts does NOT take this path -- it uses the RTKLIB-style
679 * measurement model above (RTCM-SSR has a different bias convention). */
680 ssrcode =
681 (corr == CORR_GAL_HAS) ? has_sel_biascode(sys, obs->code[i]) : mcssr_sel_biascode(sys, obs->code[i]);
682
683 /* for backward compatible to IS-QZSS-MDC-003 (MADOCA-specific) */
684 if (corr != CORR_GAL_HAS && sys == SYS_GPS && ssrcode == CODE_L5Q) {
685 if (nav->ssr_ch[0][obs->sat - 1].pbias[ssrcode - 1] == 0.0 &&
686 nav->ssr_ch[0][obs->sat - 1].cbias[ssrcode - 1] == 0.0) {
687 ssrcode = CODE_L5X;
688 }
689 }
690
691 if (ssrcode != CODE_NONE) {
692 if (nav->ssr_ch[0][obs->sat - 1].pbias[ssrcode - 1] != 0.0) {
693 pb = nav->ssr_ch[0][obs->sat - 1].pbias[ssrcode - 1];
694 }
695 if (nav->ssr_ch[0][obs->sat - 1].cbias[ssrcode - 1] != 0.0) {
696 cb = nav->ssr_ch[0][obs->sat - 1].cbias[ssrcode - 1];
697 }
698 }
699 if (cb != 0.0) {
700 P[i] += cb;
701 trace(NULL, tl > 0 ? 4 : 4, "corr_meas: %s cbias %s obscode=C%s ssrcode=C%s cbias=%7.3f\n", tstr, satid,
702 code2obs(obs->code[i]), code2obs(ssrcode), cb);
703 } else {
704 /* ssrcode==CODE_NONE has no SSR bias mapping for this signal: drop
705 * the measurement (and short-circuit to avoid vcbias[-1] OOB). */
706 if (ssrcode == CODE_NONE || !nav->ssr_ch[0][obs->sat - 1].vcbias[ssrcode - 1]) {
707 P[i] = 0.0;
708 trace(NULL, tl > 0 ? 3 : 4,
709 "corr_meas: %s cbias does not exist. %s obscode=C%s ssrcode=C%s cbias=%7.3f\n", tstr, satid,
710 code2obs(obs->code[i]), code2obs(ssrcode), cb);
711 }
712 }
713 if (cb <= (SSR_INVALID_CBIAS + 0.005)) {
714 P[i] = 0.0;
715 trace(NULL, tl > 0 ? 2 : 4, "corr_meas: %s invalid cbias %s obscode=C%s ssrcode=C%s cbias=%7.3f\n",
716 tstr, satid, code2obs(obs->code[i]), code2obs(ssrcode), cb);
717 }
718
719 if (pb != 0.0) {
720 L[i] += pb;
721 trace(NULL, tl > 0 ? 4 : 4, "corr_meas: %s pbias %s obscode=L%s ssrcode=L%s pbias=%7.3f\n", tstr, satid,
722 code2obs(obs->code[i]), code2obs(ssrcode), pb);
723 } else if (sys != SYS_GLO && corr != CORR_GAL_HAS) {
724 /* HAS Initial Service may broadcast no phase bias; for float PPP the
725 * carrier is kept and the float ambiguity absorbs the satellite phase
726 * bias (phase bias is optional for HAS, unlike the required code bias). */
727 if (ssrcode == CODE_NONE || !nav->ssr_ch[0][obs->sat - 1].vpbias[ssrcode - 1]) {
728 L[i] = 0.0;
729 trace(NULL, tl > 0 ? 3 : 4,
730 "corr_meas: %s pbias does not exist. %s obscode=C%s ssrcode=C%s pbias=%7.3f\n", tstr, satid,
731 code2obs(obs->code[i]), code2obs(ssrcode), pb);
732 }
733 }
734 if (pb <= (SSR_INVALID_PBIAS + 0.00005)) {
735 L[i] = 0.0;
736 trace(NULL, tl > 0 ? 2 : 4, "corr_meas: %s invalid pbias %s obscode=L%s ssrcode=L%s pbias=%7.3f\n",
737 tstr, satid, code2obs(obs->code[i]), code2obs(ssrcode), pb);
738 }
739 }
740 }
741
742 /* iono-free LC: pair = slot 0 + the 2nd frequency. #135: for IGS-product PPP,
743 * if the nominal slot-1 band is absent (F9P-class GAL E5b / BDS B2I), fall back
744 * to the first populated higher slot. No-op for other cases (i2 stays 1). */
745 if (corr == CORR_IGS && (freq[1] == 0.0 || P[1] == 0.0 || L[1] == 0.0)) {
746 int j;
747 for (j = 2; j < NFREQ; j++) {
748 if (freq[j] != 0.0 && P[j] != 0.0 && L[j] != 0.0) {
749 i2 = j;
750 ppp_note_autosig(sys, obs->code[i2]); /* notify user of the fallback */
751 break;
752 }
753 }
754 }
755 *Lc = *Pc = 0.0;
756 if (freq[0] == 0.0 || freq[i2] == 0.0) {
757 return 0;
758 }
759 C1 = SQR(freq[0]) / (SQR(freq[0]) - SQR(freq[i2]));
760 C2 = -SQR(freq[i2]) / (SQR(freq[0]) - SQR(freq[i2]));
761
762 if (P[0] == 0.0 || P[i2] == 0.0 || L[0] == 0.0 || L[i2] == 0.0) {
763 trace(NULL, tl > 0 ? 2 : 4, "corr_meas: no sufficient obs data %s %s \n", tstr, satid);
764 return 0;
765 }
766 *Lc = C1 * L[0] + C2 * L[i2];
767 *Pc = C1 * P[0] + C2 * P[i2];
768
769 return 1;
770}
771/* detect cycle slip by LLI --------------------------------------------------*/
772/* detect cycle slip by Doppler/phase rate comparison (demo5, single-receiver) */
773static void detslp_dop(rtk_t* rtk, const obsd_t* obs, int n, const nav_t* nav) {
774 (void)nav;
775 int i, f, sat, ndop = 0, nf = rtk->opt.nf < NFREQ ? rtk->opt.nf : NFREQ;
776 double dph, dpt, tt, mean_dop = 0.0;
777 double dopdif[MAXOBS][NFREQ] = {{0}}, tts[MAXOBS][NFREQ] = {{0}};
778
779 if (rtk->opt.thresdop <= 0.0) {
780 return;
781 }
782
783 for (i = 0; i < n && i < MAXOBS; i++) {
784 sat = obs[i].sat;
785 for (f = 0; f < nf; f++) {
786 if (obs[i].L[f] == 0.0 || obs[i].D[f] == 0.0 || rtk->ssat[sat - 1].ph[0][f] == 0.0) {
787 continue;
788 }
789 tt = timediff(obs[i].time, rtk->ssat[sat - 1].pt[0][f]);
790 if (fabs(tt) < DTTOL) {
791 continue;
792 }
793 dph = (obs[i].L[f] - rtk->ssat[sat - 1].ph[0][f]) / tt;
794 dpt = -obs[i].D[f];
795 dopdif[i][f] = dph - dpt;
796 tts[i][f] = tt;
797 if (fabs(dopdif[i][f]) < 3.0 * rtk->opt.thresdop) {
798 mean_dop += dopdif[i][f];
799 ndop++;
800 }
801 }
802 }
803 if (ndop == 0) {
804 return;
805 }
806 mean_dop /= ndop;
807
808 for (i = 0; i < n && i < MAXOBS; i++) {
809 sat = obs[i].sat;
810 for (f = 0; f < nf; f++) {
811 if (dopdif[i][f] == 0.0) {
812 continue;
813 }
814 if (fabs(dopdif[i][f] - mean_dop) > rtk->opt.thresdop) {
815 rtk->ssat[sat - 1].slip[f] |= 1;
816 trace(NULL, 3, "ppp detslp_dop: sat=%2d F=%d dL=%.3f off=%.3f tt=%.2f\n", sat, f + 1,
817 dopdif[i][f] - mean_dop, mean_dop, tts[i][f]);
818 }
819 }
820 }
821}
822/* detect cycle slip by observation code change (demo5, single-receiver) */
823static void detslp_code(rtk_t* rtk, const obsd_t* obs, int n) {
824 int i, f, sat, nf = rtk->opt.nf < NFREQ ? rtk->opt.nf : NFREQ;
825 for (i = 0; i < n && i < MAXOBS; i++) {
826 sat = obs[i].sat;
827 for (f = 0; f < nf; f++) {
828 uint8_t code = obs[i].code[f];
829 if (code == CODE_NONE) {
830 continue;
831 }
832 uint8_t ccode = rtk->ssat[sat - 1].codeprev[f][0];
833 if (code != ccode) {
834 rtk->ssat[sat - 1].codeprev[f][0] = code;
835 if (ccode != CODE_NONE) {
836 rtk->ssat[sat - 1].slip[f] |= 1;
837 trace(NULL, 3, "ppp detslp_code: sat=%2d F=%d %s->%s\n", sat, f + 1, code2obs(ccode),
838 code2obs(code));
839 }
840 }
841 }
842 }
843}
844static void detslp_ll(rtk_t* rtk, const obsd_t* obs, int n) {
845 int i, j;
846 char satid[4];
847 trace(NULL, 3, "detslp_ll: n=%d\n", n);
848
849 /* Note, When using gcc, use constants to suppress
850 -Wstringop-overflow warnings. */
851 for (i = 0; i < n && i < MAXOBS; i++) {
852 for (j = 0; j < rtk->opt.nf && j < NFREQ; j++) {
853 if (obs[i].L[j] == 0.0 || !(obs[i].LLI[j] & 3)) {
854 continue;
855 }
856 satno2id(obs[i].sat, satid);
857 trace(NULL, 3, "detslp_ll: slip detected %s sat=%2d f=%d\n", satid, obs[i].sat, j + 1);
858
859 rtk->ssat[obs[i].sat - 1].slip[j] = 1;
860 }
861 }
862}
863/* detect cycle slip by geometry free phase jump -----------------------------*/
864static void detslp_gf(rtk_t* rtk, const obsd_t* obs, int n, const nav_t* nav) {
865 double g0, g1;
866 int i, f;
867 char satid[4];
868 trace(NULL, 3, "detslp_gf: n=%d\n", n);
869
870 for (i = 0; i < n && i < MAXOBS; i++) {
871 for (f = 1; f < rtk->opt.nf && f < NFREQ; f++) {
872 if ((g1 = gfmeas(obs + i, nav, f)) == 0.0) {
873 continue;
874 }
875
876 g0 = rtk->ssat[obs[i].sat - 1].gf[f];
877 rtk->ssat[obs[i].sat - 1].gf[f] = g1;
878
879 trace(NULL, 4, "detslip_gf: sat=%2d f=%d gf0=%8.3f gf1=%8.3f\n", obs[i].sat, f + 1, g0, g1);
880
881 if (g0 != 0.0 && fabs(g1 - g0) > rtk->opt.thresslip) {
882 satno2id(obs[i].sat, satid);
883 trace(NULL, 2, "detslip_gf: slip detected %s sat=%2d f=%d gf=%8.3f->%8.3f\n", satid, obs[i].sat, f + 1,
884 g0, g1);
885
886 rtk->ssat[obs[i].sat - 1].slip[0] |= 1;
887 rtk->ssat[obs[i].sat - 1].slip[f] |= 1;
888 }
889 }
890 }
891}
892/* detect slip by Melbourne-Wubbena linear combination jump ------------------*/
893static void detslp_mw(rtk_t* rtk, const obsd_t* obs, int n, const nav_t* nav) {
894 double w0, w1;
895 int i, f;
896 char satid[4];
897 trace(NULL, 3, "detslp_mw: n=%d\n", n);
898
899 for (i = 0; i < n && i < MAXOBS; i++) {
900 for (f = 1; f < rtk->opt.nf && f < NFREQ; f++) {
901 if ((w1 = mwmeas(obs + i, nav, f)) == 0.0) {
902 continue;
903 }
904
905 w0 = rtk->ssat[obs[i].sat - 1].mw[f];
906 rtk->ssat[obs[i].sat - 1].mw[f] = w1;
907
908 trace(NULL, 4, "detslip_mw: sat=%2d f=%d mw0=%8.3f mw1=%8.3f\n", obs[i].sat, f + 1, w0, w1);
909
910 if (w0 != 0.0 && fabs(w1 - w0) > THRES_MW_JUMP) {
911 satno2id(obs[i].sat, satid);
912 trace(NULL, 2, "detslip_mw: slip detected %s sat=%2d f=%d mw=%8.3f->%8.3f\n", satid, obs[i].sat, f + 1,
913 w0, w1);
914
915 rtk->ssat[obs[i].sat - 1].slip[0] |= 1;
916 rtk->ssat[obs[i].sat - 1].slip[f] |= 1;
917 }
918 }
919 }
920}
921/* detect cycle slip by ssr phase bias discontinuity indicator ---------------*/
922static void detslp_ssr(rtk_t* rtk, const obsd_t* obs, int n, const nav_t* nav) {
923 int discont0, discont1, ssrcode;
924 int sys, i, j;
925 char satid[4];
926 trace(NULL, 3, "detslp_ssr: n=%d\n", n);
927
928 if (rtk->opt.sateph != EPHOPT_SSRAPC && rtk->opt.sateph != EPHOPT_SSRCOM) {
929 return;
930 }
931
932 for (i = 0; i < n && i < MAXOBS; i++) {
933 if (!nav->ssr_ch[0][obs[i].sat - 1].t0[5].time) {
934 continue;
935 }
936
937 sys = satsys(obs[i].sat, NULL);
938
939 /* detect slip by discontinuity indicator */
940 for (j = 0; j < rtk->opt.nf && j < NFREQ; j++) {
941 ssrcode = mcssr_sel_biascode(sys, obs[i].code[j]);
942 if (ssrcode == CODE_NONE) {
943 continue;
944 }
945 if (nav->ssr_ch[0][obs[i].sat - 1].pbias[ssrcode - 1] == 0.0) {
946 continue;
947 }
948
949 discont0 = rtk->ssat[obs[i].sat - 1].discont[j];
950 discont1 = nav->ssr_ch[0][obs[i].sat - 1].discnt[ssrcode - 1];
951 rtk->ssat[obs[i].sat - 1].discont[j] = discont1;
952
953 if (discont0 != discont1) {
954 satno2id(obs[i].sat, satid);
955 trace(NULL, 2,
956 "detslp_ssr: discont detected sys=%2d %s sat=%2d f=%d code=%2d discont0=%d discont1=%d\n", sys,
957 satid, obs[i].sat, j + 1, ssrcode, discont0, discont1);
958
959 rtk->ssat[obs[i].sat - 1].slip[j] |= 1;
960 }
961 }
962 }
963}
964/* temporal update of position -----------------------------------------------*/
965static void udpos_ppp(rtk_t* rtk) {
966 double *F, *P, *FP, *x, *xp, pos[3], Q[9] = {0}, Qv[9];
967 int i, j, *ix, nx;
968
969 trace(NULL, 3, "udpos_ppp:\n");
970
971 /* fixed mode */
972 if (rtk->opt.mode == PMODE_PPP_FIXED) {
973 for (i = 0; i < 3; i++) {
974 initx(rtk, rtk->opt.ru[i], 1E-8, i);
975 }
976 return;
977 }
978 /* initialize position for first epoch */
979 if (norm(rtk->x, 3) <= 0.0) {
980 for (i = 0; i < 3; i++) {
981 initx(rtk, rtk->sol.rr[i], VAR_POS, i);
982 }
983 if (rtk->opt.dynamics) {
984 for (i = 3; i < 6; i++) {
985 initx(rtk, rtk->sol.rr[i], VAR_VEL, i);
986 }
987 for (i = 6; i < 9; i++) {
988 initx(rtk, 1E-6, VAR_ACC, i);
989 }
990 }
991 }
992 /* static ppp mode */
993 if (rtk->opt.mode == PMODE_PPP_STATIC) {
994 for (i = 0; i < 3; i++) {
995 rtk->P[i * (1 + rtk->nx)] += SQR(rtk->opt.prn[5]) * fabs(rtk->tt);
996 }
997 return;
998 }
999 /* kinmatic mode without dynamics */
1000 if (!rtk->opt.dynamics) {
1001 for (i = 0; i < 3; i++) {
1002 initx(rtk, rtk->sol.rr[i], VAR_POS, i);
1003 }
1004 return;
1005 }
1006 /* generate valid state index */
1007 ix = imat(rtk->nx, 1);
1008 for (i = nx = 0; i < rtk->nx; i++) {
1009 if (rtk->x[i] != 0.0 && rtk->P[i + i * rtk->nx] > 0.0) {
1010 ix[nx++] = i;
1011 }
1012 }
1013 if (nx < 9) {
1014 free(ix);
1015 return;
1016 }
1017 /* state transition of position/velocity/acceleration */
1018 F = eye(nx);
1019 P = mat(nx, nx);
1020 FP = mat(nx, nx);
1021 x = mat(nx, 1);
1022 xp = mat(nx, 1);
1023
1024 for (i = 0; i < 6; i++) {
1025 F[i + (i + 3) * nx] = rtk->tt;
1026 }
1027 for (i = 0; i < 3; i++) {
1028 F[i + (i + 6) * nx] = SQR(rtk->tt) / 2.0;
1029 }
1030 for (i = 0; i < nx; i++) {
1031 x[i] = rtk->x[ix[i]];
1032 for (j = 0; j < nx; j++) {
1033 P[i + j * nx] = rtk->P[ix[i] + ix[j] * rtk->nx];
1034 }
1035 }
1036 /* x=F*x, P=F*P*F+Q */
1037 matmul("NN", nx, 1, nx, 1.0, F, x, 0.0, xp);
1038 matmul("NN", nx, nx, nx, 1.0, F, P, 0.0, FP);
1039 matmul("NT", nx, nx, nx, 1.0, FP, F, 0.0, P);
1040
1041 for (i = 0; i < nx; i++) {
1042 rtk->x[ix[i]] = xp[i];
1043 for (j = 0; j < nx; j++) {
1044 rtk->P[ix[i] + ix[j] * rtk->nx] = P[i + j * nx];
1045 }
1046 }
1047 /* process noise added to only acceleration */
1048 Q[0] = Q[4] = SQR(rtk->opt.prn[3]) * fabs(rtk->tt);
1049 Q[8] = SQR(rtk->opt.prn[4]) * fabs(rtk->tt);
1050 ecef2pos(rtk->x, pos);
1051 covecef(pos, Q, Qv);
1052 for (i = 0; i < 3; i++) {
1053 for (j = 0; j < 3; j++) {
1054 rtk->P[i + 6 + (j + 6) * rtk->nx] += Qv[i + j * 3];
1055 }
1056 }
1057 free(ix);
1058 free(F);
1059 free(P);
1060 free(FP);
1061 free(x);
1062 free(xp);
1063}
1064/* temporal update of clock --------------------------------------------------*/
1065static void udclk_ppp(rtk_t* rtk) {
1066 double dtr;
1067 int i;
1068
1069 trace(NULL, 3, "udclk_ppp:\n");
1070
1071 /* initialize every epoch for clock (white noise) */
1072 for (i = 0; i < NC(rtk->opt); i++) {
1073 if (rtk->opt.sateph == EPHOPT_PREC) {
1074 /* time of prec ephemeris is based gpst */
1075 /* negelect receiver inter-system bias */
1076 dtr = rtk->sol.dtr[0];
1077 } else {
1078 dtr = i == 0 ? rtk->sol.dtr[0] : rtk->sol.dtr[0] + rtk->sol.dtr[i];
1079 }
1080 initx(rtk, CLIGHT * dtr, VAR_CLK, IC(i, &rtk->opt));
1081 }
1082}
1083/* temporal update of tropospheric parameters --------------------------------*/
1084static void udtrop_ppp(rtk_t* rtk) {
1085 double pos[3], azel[] = {0.0, PI / 2.0}, ztd, var;
1086 int i = IT(&rtk->opt), j;
1087
1088 trace(NULL, 3, "udtrop_ppp:\n");
1089
1090 if (rtk->x[i] == 0.0) {
1091 ecef2pos(rtk->sol.rr, pos);
1092 ztd = sbstropcorr(rtk->sol.time, pos, azel, &var);
1093 initx(rtk, ztd, VAR_ZTD, i);
1094
1095 if (rtk->opt.tropopt >= TROPOPT_ESTG) {
1096 for (j = i + 1; j < i + 3; j++) {
1097 initx(rtk, 1E-6, VAR_GRA, j);
1098 }
1099 }
1100 } else {
1101 rtk->P[i + i * rtk->nx] += SQR(rtk->opt.prn[2]) * fabs(rtk->tt);
1102 trace(NULL, 2, "udtrop_ppp: rtk->opt.prn[2]=%.3e\n", rtk->opt.prn[2]);
1103 if (rtk->opt.tropopt >= TROPOPT_ESTG) {
1104 for (j = i + 1; j < i + 3; j++) {
1105 rtk->P[j + j * rtk->nx] += SQR(rtk->opt.prn[2] * 0.1) * fabs(rtk->tt);
1106 }
1107 }
1108 }
1109}
1110/* temporal update of ionospheric parameters ---------------------------------*/
1111static void udiono_ppp(rtk_t* rtk, const obsd_t* obs, int n, const nav_t* nav) {
1112 double freq1, freq2, ionp, ionc, sinel, pos[3], *azel, el;
1113 char *p, satid[4];
1114 int i, j, gap_resion = GAP_RESION;
1115 double L[NFREQ], P[NFREQ], Lc, Pc;
1116 double dantr[NFREQPCV] = {0}, dants[NFREQPCV] = {0}, dion;
1117 uint8_t* slip;
1118 trace(NULL, 3, "udiono_ppp:\n");
1119
1120 if ((p = strstr(rtk->opt.pppopt, "-GAP_RESION="))) {
1121 sscanf(p, "-GAP_RESION=%d", &gap_resion);
1122 }
1123 for (i = 0; i < MAXSAT; i++) {
1124 j = II(i + 1, &rtk->opt);
1125 if (rtk->x[j] != 0.0 && (int)rtk->ssat[i].outc[0] > gap_resion) {
1126 rtk->x[j] = 0.0;
1127
1128 if ((int)rtk->ssat[i].outc[0] == gap_resion + 1 && rtk->ssat[i].azel[1] >= rtk->opt.elmin) {
1129 satno2id(i + 1, satid);
1130 el = rtk->ssat[i].azel[1] * R2D;
1131 trace(NULL, 3, "udiono_ppp(%4d): %s reset ionospheric estimates. %s outc=%d el=%2.0f\n", __LINE__,
1132 time_str(obs[0].time, 0), satid, (int)rtk->ssat[i].outc[0], el);
1133 }
1134 }
1135 }
1136 for (i = 0; i < n; i++) {
1137 j = II(obs[i].sat, &rtk->opt);
1138 corr_meas(obs + i, nav, rtk->ssat[obs[i].sat - 1].azel, &rtk->opt, dantr, dants, 0.0, L, P, &Lc, &Pc, 0);
1139 freq1 = sat2freq(obs[i].sat, obs[i].code[0], nav);
1140 freq2 = sat2freq(obs[i].sat, obs[i].code[1], nav);
1141 slip = rtk->ssat[obs[i].sat - 1].slip;
1142
1143 if (rtk->x[j] == 0.0) {
1144 if (P[0] == 0.0 || P[1] == 0.0 || freq1 == 0.0 || freq2 == 0.0) {
1145 continue;
1146 }
1147 ionp = (P[0] - P[1]) / (SQR(FREQ1 / freq1) - SQR(FREQ1 / freq2));
1148 ecef2pos(rtk->sol.rr, pos);
1149 azel = rtk->ssat[obs[i].sat - 1].azel;
1150 initx(rtk, ionp, VAR_IONO, j);
1151
1152 if (L[0] == 0.0 || L[1] == 0.0 || slip[0] || slip[1]) {
1153 rtk->ssat[obs[i].sat - 1].ionc = 0.0;
1154 } else {
1155 rtk->ssat[obs[i].sat - 1].ionc = -1.0 * (L[0] - L[1]) / (SQR(FREQ1 / freq1) - SQR(FREQ1 / freq2));
1156 }
1157 if (azel[1] >= rtk->opt.elmin) {
1158 satno2id(obs[i].sat, satid);
1159 trace(NULL, 3, "udiono_ppp(%4d): %s initialize ionospheric estimates. %s el=%2.0f x=%8.4f\n", __LINE__,
1160 time_str(obs[i].time, 0), satid, azel[1] * R2D, rtk->x[j]);
1161 }
1162 } else {
1163 sinel = sin(MAX(rtk->ssat[obs[i].sat - 1].azel[1], 5.0 * D2R));
1164 rtk->P[j + j * rtk->nx] += SQR(rtk->opt.prn[1] / sinel) * fabs(rtk->tt);
1165
1166 if (L[0] == 0.0 || L[1] == 0.0 || freq1 == 0.0 || freq2 == 0.0 || slip[0] || slip[1]) {
1167 rtk->ssat[obs[i].sat - 1].ionc = 0.0;
1168 continue;
1169 }
1170 ionc = -1.0 * (L[0] - L[1]) / (SQR(FREQ1 / freq1) - SQR(FREQ1 / freq2));
1171 dion = 0.0;
1172 if (rtk->ssat[obs[i].sat - 1].ionc != 0.0) {
1173 dion = ionc - rtk->ssat[obs[i].sat - 1].ionc;
1174 }
1175 rtk->ssat[obs[i].sat - 1].ionc = ionc;
1176 rtk->x[j] += dion;
1177 }
1178 }
1179}
1180/* temporal update of 3rd/4th-freqency-receiver-bias parameters --------------*/
1181static void udifb_ppp(rtk_t* rtk) {
1182 int i, j;
1183
1184 trace(NULL, 3, "udifb_ppp:\n");
1185
1186 /* 3rd-freqency-receiver-bias */
1187 for (i = 0; i < N3(&rtk->opt); i++) {
1188 j = I3(i, &rtk->opt);
1189
1190 if (rtk->x[j] == 0.0) {
1191 initx(rtk, 1e-6, VAR_IFB, j);
1192 } else {
1193 /* add process noise */
1194 rtk->P[j + j * rtk->nx] += SQR(rtk->opt.prn[6]) * fabs(rtk->tt);
1195 }
1196 }
1197
1198 /* 4th-freqency-receiver-bias */
1199 for (i = 0; i < N4(&rtk->opt); i++) {
1200 j = I4(i, &rtk->opt);
1201
1202 if (rtk->x[j] == 0.0) {
1203 initx(rtk, 1e-6, VAR_IFB, j);
1204 } else {
1205 /* add process noise */
1206 rtk->P[j + j * rtk->nx] += SQR(rtk->opt.prn[6]) * fabs(rtk->tt);
1207 }
1208 }
1209}
1210/* temporal update of phase biases -------------------------------------------*/
1211static void udbias_ppp(rtk_t* rtk, const obsd_t* obs, int n, const nav_t* nav) {
1212 double L[NFREQ], P[NFREQ], Lc, Pc, bias[MAXOBS], offset = 0.0, pos[3] = {0};
1213 double freq1, freq2, freq3, ion, dantr[NFREQPCV] = {0}, dants[NFREQPCV] = {0};
1214 int i, j, k, f, sat, slip[MAXOBS] = {0}, clk_jump = 0;
1215
1216 trace(NULL, 3, "udbias : n=%d\n", n);
1217
1218 /* handle day-boundary clock jump */
1219 if (rtk->opt.clockjump) {
1220 clk_jump = ROUND(time2gpst(obs[0].time, NULL) * 10) % 864000 == 0;
1221 }
1222 for (i = 0; i < MAXSAT; i++) {
1223 for (j = 0; j < rtk->opt.nf && j < NFREQ; j++) {
1224 rtk->ssat[i].slip[j] = 0;
1225 }
1226 }
1227 /* detect cycle slip by LLI */
1228 detslp_ll(rtk, obs, n);
1229
1230 /* detect cycle slip by geometry-free phase jump */
1231 detslp_gf(rtk, obs, n, nav);
1232
1233 /* detect slip by Melbourne-Wubbena linear combination jump */
1234 detslp_mw(rtk, obs, n, nav);
1235
1236 /* detect cycle slip by ssr */
1237 detslp_ssr(rtk, obs, n, nav);
1238
1239 /* detect cycle slip by Doppler rate (demo5) */
1240 detslp_dop(rtk, obs, n, nav);
1241
1242 /* detect cycle slip by observation code change (demo5) */
1243 detslp_code(rtk, obs, n);
1244
1245 ecef2pos(rtk->sol.rr, pos);
1246
1247 for (f = 0; f < NF(&rtk->opt); f++) {
1248 /* reset phase-bias if expire obs outage counter */
1249 for (i = 0; i < MAXSAT; i++) {
1250 if (++rtk->ssat[i].outc[f] > (uint32_t)rtk->opt.maxout || rtk->opt.modear == ARMODE_INST || clk_jump) {
1251 initx(rtk, 0.0, 0.0, IB(i + 1, f, &rtk->opt));
1252 }
1253 }
1254 for (i = k = 0; i < n && i < MAXOBS; i++) {
1255 sat = obs[i].sat;
1256 j = IB(sat, f, &rtk->opt);
1257 corr_meas(obs + i, nav, rtk->ssat[sat - 1].azel, &rtk->opt, dantr, dants, 0.0, L, P, &Lc, &Pc, 0);
1258
1259 bias[i] = 0.0;
1260
1261 if (rtk->opt.ionoopt == IONOOPT_IFLC) {
1262 bias[i] = Lc - Pc;
1263 slip[i] = rtk->ssat[sat - 1].slip[0] || rtk->ssat[sat - 1].slip[1];
1264 } else if (L[f] != 0.0 && P[f] != 0.0) {
1265 freq1 = sat2freq(sat, obs[i].code[0], nav);
1266 freq2 = sat2freq(sat, obs[i].code[1], nav);
1267 freq3 = sat2freq(sat, obs[i].code[f], nav);
1268 slip[i] = rtk->ssat[sat - 1].slip[f];
1269 if (obs[i].P[0] == 0.0 || obs[i].P[1] == 0.0 || obs[i].P[f] == 0.0 || freq1 == 0.0 || freq2 == 0.0 ||
1270 freq3 == 0.0) {
1271 continue;
1272 }
1273 ion = (obs[i].P[0] - obs[i].P[1]) / (1.0 - SQR(freq1 / freq2));
1274 bias[i] = L[f] - P[f] + 2.0 * ion * SQR(freq1 / freq3);
1275 }
1276 if (rtk->x[j] == 0.0 || slip[i] || bias[i] == 0.0) {
1277 continue;
1278 }
1279
1280 offset += bias[i] - rtk->x[j];
1281 k++;
1282 }
1283 /* correct phase-code jump to ensure phase-code coherency */
1284 if (k >= 2 && fabs(offset / k) > 0.0005 * CLIGHT) {
1285 for (i = 0; i < MAXSAT; i++) {
1286 j = IB(i + 1, f, &rtk->opt);
1287 if (rtk->x[j] != 0.0) {
1288 rtk->x[j] += offset / k;
1289 }
1290 }
1291 trace(NULL, 2, "phase-code jump corrected: %s n=%2d dt=%12.9fs\n", time_str(rtk->sol.time, 0), k,
1292 offset / k / CLIGHT);
1293 }
1294 for (i = 0; i < n && i < MAXOBS; i++) {
1295 sat = obs[i].sat;
1296 j = IB(sat, f, &rtk->opt);
1297
1298 rtk->P[j + j * rtk->nx] += SQR(rtk->opt.prn[0]) * fabs(rtk->tt);
1299
1300 if (bias[i] == 0.0 || (rtk->x[j] != 0.0 && !slip[i])) {
1301 continue;
1302 }
1303
1304 /* reinitialize phase-bias if detecting cycle slip */
1305 initx(rtk, bias[i], VAR_BIAS, IB(sat, f, &rtk->opt));
1306
1307 /* reset fix flags */
1308 for (k = 0; k < MAXSAT; k++) {
1309 rtk->ambc[sat - 1].flags[k] = 0;
1310 }
1311
1312 trace(NULL, 5, "udbias_ppp: sat=%2d bias=%.3f\n", sat, bias[i]);
1313 }
1314 }
1315}
1316/* temporal update of states --------------------------------------------------*/
1317static void udstate_ppp(rtk_t* rtk, const obsd_t* obs, int n, const nav_t* nav) {
1318 trace(NULL, 3, "udstate_ppp: n=%d\n", n);
1319
1320 /* temporal update of position */
1321 udpos_ppp(rtk);
1322
1323 /* temporal update of clock */
1324 udclk_ppp(rtk);
1325
1326 /* temporal update of tropospheric parameters */
1327 if (rtk->opt.tropopt == TROPOPT_EST || rtk->opt.tropopt == TROPOPT_ESTG) {
1328 udtrop_ppp(rtk);
1329 }
1330 /* temporal update of 3rd/4th-freqency-receiver-bias parameters */
1331 if (NF(&rtk->opt) >= 3) {
1332 udifb_ppp(rtk);
1333 }
1334 /* temporal update of phase-bias */
1335 udbias_ppp(rtk, obs, n, nav);
1336
1337 /* temporal update of ionospheric parameters */
1338 if (rtk->opt.ionoopt == IONOOPT_EST) {
1339 udiono_ppp(rtk, obs, n, nav);
1340 }
1341}
1342/* satellite antenna phase center variation ----------------------------------*/
1343static void satantpcv(const double* rs, const double* rr, const pcv_t* pcv, double* dant) {
1344 double ru[3], rz[3], eu[3], ez[3], nadir, cosa;
1345 int i;
1346
1347 return; /* disable for madoca-ppp */
1348
1349 for (i = 0; i < 3; i++) {
1350 ru[i] = rr[i] - rs[i];
1351 rz[i] = -rs[i];
1352 }
1353 if (!normv3(ru, eu) || !normv3(rz, ez)) {
1354 return;
1355 }
1356
1357 cosa = dot(eu, ez, 3);
1358 cosa = cosa < -1.0 ? -1.0 : (cosa > 1.0 ? 1.0 : cosa);
1359 nadir = acos(cosa);
1360
1361 antmodel_s(pcv, nadir, dant);
1362}
1363/* precise tropospheric model ------------------------------------------------*/
1364static double trop_model_prec(gtime_t time, const double* pos, const double* azel, const double* x, double* dtdx,
1365 double* var) {
1366 const double zazel[] = {0.0, PI / 2.0};
1367 double zhd, m_h, m_w, cotz, grad_n, grad_e;
1368
1369 /* zenith hydrostatic delay */
1370 zhd = tropmodel(time, pos, zazel, 0.0);
1371
1372 /* mapping function */
1373 m_h = tropmapf(time, pos, azel, &m_w);
1374
1375 if (azel[1] > 0.0) {
1376 /* m_w=m_0+m_0*cot(el)*(Gn*cos(az)+Ge*sin(az)): ref [6] */
1377 cotz = 1.0 / tan(azel[1]);
1378 grad_n = m_w * cotz * cos(azel[0]);
1379 grad_e = m_w * cotz * sin(azel[0]);
1380 m_w += grad_n * x[1] + grad_e * x[2];
1381 dtdx[1] = grad_n * (x[0] - zhd);
1382 dtdx[2] = grad_e * (x[0] - zhd);
1383 }
1384 dtdx[0] = m_w;
1385 *var = SQR(0.01);
1386 return m_h * zhd + m_w * (x[0] - zhd);
1387}
1388/* tropospheric model ---------------------------------------------------------*/
1389static int model_trop(gtime_t time, const double* pos, const double* azel, const prcopt_t* opt, const double* x,
1390 double* dtdx, const nav_t* nav, double* dtrp, double* var) {
1391 double trp[3] = {0};
1392
1393 if (opt->tropopt == TROPOPT_SAAS) {
1394 *dtrp = tropmodel(time, pos, azel, REL_HUMI);
1395 *var = SQR(ERR_SAAS);
1396 return 1;
1397 }
1398 if (opt->tropopt == TROPOPT_SBAS) {
1399 *dtrp = sbstropcorr(time, pos, azel, var);
1400 return 1;
1401 }
1402 if (opt->tropopt == TROPOPT_EST || opt->tropopt == TROPOPT_ESTG) {
1403 matcpy(trp, x + IT(opt), opt->tropopt == TROPOPT_EST ? 1 : 3, 1);
1404 *dtrp = trop_model_prec(time, pos, azel, trp, dtdx, var);
1405 return 1;
1406 }
1407 return 0;
1408}
1409/* ionospheric model ---------------------------------------------------------*/
1410static int model_iono(gtime_t time, const double* pos, const double* azel, const prcopt_t* opt, int sat,
1411 const double* x, const nav_t* nav, double* dion, double* var) {
1412 if (opt->ionoopt == IONOOPT_SBAS) {
1413 return sbsioncorr(time, nav, pos, azel, dion, var);
1414 }
1415 if (opt->ionoopt == IONOOPT_TEC) {
1416 return iontec(time, nav, pos, azel, 1, dion, var);
1417 }
1418 if (opt->ionoopt == IONOOPT_BRDC) {
1419 *dion = ionmodel(time, nav->ion_gps, pos, azel);
1420 *var = SQR(*dion * ERR_BRDCI);
1421 return 1;
1422 }
1423 if (opt->ionoopt == IONOOPT_EST) {
1424 *dion = x[II(sat, opt)];
1425 *var = 0.0;
1426 return 1;
1427 }
1428 if (opt->ionoopt == IONOOPT_IFLC) {
1429 *dion = *var = 0.0;
1430 return 1;
1431 }
1432 return 0;
1433}
1434/* phase and code residuals --------------------------------------------------*/
1435static int ppp_res(int post, const obsd_t* obs, int n, const double* rs, const double* dts, const double* var_rs,
1436 const int* svh, const double* dr, int* exc, const nav_t* nav, const double* x, const double* cov,
1437 rtk_t* rtk, double* v, double* H, double* R, double* azel, double* resrms) {
1438 prcopt_t* opt = &rtk->opt;
1439 double y, r, cdtr, bias, C = 0.0, rr[3], pos[3], e[3], dtdx[3] = {0.0}, L[NFREQ], P[NFREQ], Lc, Pc;
1440 double var[MAXOBS * 2 * NFREQ], dtrp = 0.0, dion = 0.0, vart = 0.0, vari = 0.0, freq, ifb, ifb_std, bias_std,
1441 dion_std, dtrp_std;
1442 double dantr[NFREQPCV] = {0}, dants[NFREQPCV] = {0};
1443 double ve[MAXOBS * 2 * NFREQ] = {0}, vmax = 0, resc_rms = 0, resp_rms = 0, resc_ave = 0, resp_ave = 0;
1444 char str[32], satstr[6];
1445 int ne = 0, obsi[MAXOBS * 2 * NFREQ] = {0}, frqi[MAXOBS * 2 * NFREQ], maxobs, maxfrq, rej;
1446 int i, j, k, sat, sys, nv = 0, nx = rtk->nx, stat = 1, nresc = 0, nresp = 0;
1447
1448 time2str(obs[0].time, str, 2);
1449
1450 for (i = 0; i < MAXSAT; i++) {
1451 for (j = 0; j < opt->nf && j < NFREQ; j++) {
1452 rtk->ssat[i].vsat[j] = 0;
1453 }
1454 }
1455
1456 for (i = 0; i < 3; i++) {
1457 rr[i] = x[i] + dr[i];
1458 }
1459 ecef2pos(rr, pos);
1460
1461 for (i = 0; i < n && i < MAXOBS; i++) {
1462 sat = obs[i].sat;
1463 satno2id(sat, satstr);
1464
1465 if ((r = geodist(rs + i * 6, rr, e)) <= 0.0 || satazel(pos, e, azel + i * 2) < opt->elmin) {
1466 exc[i] = 1;
1467 continue;
1468 }
1469 if (!(sys = satsys(sat, NULL)) || !rtk->ssat[sat - 1].vs || satexclude(obs[i].sat, var_rs[i], svh[i], opt) ||
1470 exc[i]) {
1471 exc[i] = 1;
1472 continue;
1473 }
1474 /* tropospheric and ionospheric model */
1475 if (!model_trop(obs[i].time, pos, azel + i * 2, opt, x, dtdx, nav, &dtrp, &vart) ||
1476 !model_iono(obs[i].time, pos, azel + i * 2, opt, sat, x, nav, &dion, &vari)) {
1477 continue;
1478 }
1479 /* receiver antenna model */
1480 if (opt->posopt[0]) {
1481 satantpcv(rs + i * 6, rr, nav->pcvs + sat - 1, dants);
1482 }
1483 antmodel(opt->pcvr, opt->antdel[0], azel + i * 2, opt->posopt[1], dantr);
1484
1485 /* phase windup model */
1486 if (!model_phw(rtk->sol.time, sat, nav->pcvs[sat - 1].type, opt->posopt[2] ? 2 : 0, rs + i * 6, rr,
1487 &rtk->ssat[sat - 1].phw)) {
1488 continue;
1489 }
1490 /* corrected phase and code measurements */
1491 if (!corr_meas(obs + i, nav, azel + i * 2, &rtk->opt, dantr, dants, rtk->ssat[sat - 1].phw, L, P, &Lc, &Pc,
1492 post)) {
1493 continue;
1494 }
1495 trace(NULL, post > 0 ? 3 : 4, "ppp_res: %s %s sat=%3d nf=%d\n", str, satstr, sat, NF(opt));
1496
1497 /* stack phase and code residuals {L1,P1,L2,P2,...} */
1498 for (j = 0; j < 2 * NF(opt); j++) {
1499 bias = ifb = 0.0;
1500 bias_std = ifb_std = dion_std = dtrp_std = 0.0;
1501
1502 if (opt->ionoopt == IONOOPT_IFLC) {
1503 if ((y = j % 2 == 0 ? Lc : Pc) == 0.0) {
1504 continue;
1505 }
1506 } else {
1507 if ((y = j % 2 == 0 ? L[j / 2] : P[j / 2]) == 0.0) {
1508 continue;
1509 }
1510
1511 if ((freq = sat2freq(sat, obs[i].code[j / 2], nav)) == 0.0) {
1512 continue;
1513 }
1514 C = SQR(FREQ1 / freq) * (j % 2 == 0 ? -1.0 : 1.0);
1515 }
1516 for (k = 0; k < nx; k++) {
1517 H[k + nx * nv] = k < 3 ? -e[k] : 0.0;
1518 }
1519
1520 /* receiver clock */
1521 switch (satsys_bd2(sat, NULL)) {
1522 case SYS_GLO:
1523 k = 1;
1524 break;
1525 case SYS_GAL:
1526 k = 2;
1527 break;
1528 case SYS_CMP:
1529 k = 3;
1530 break; /* BDS-3 */
1531 case SYS_IRN:
1532 k = 4;
1533 break;
1534 case SYS_QZS:
1535 k = 5;
1536 break;
1537 case SYS_BD2:
1538 k = 6;
1539 break; /* BDS-2 */
1540 default:
1541 k = 0;
1542 break;
1543 }
1544 cdtr = x[IC(k, opt)];
1545 H[IC(k, opt) + nx * nv] = 1.0;
1546
1547 if (j / 2 == 2 && j % 2 == 1) { /* 3rd-frequency-receiver-bias */
1548 ifb = x[I3(k, opt)];
1549 ifb_std = SQRT(cov[I3(k, opt) + I3(k, opt) * rtk->nx]);
1550 H[I3(k, opt) + nx * nv] = 1.0;
1551 }
1552 if (j / 2 == 3 && j % 2 == 1) { /* 4th-frequency-receiver-bias */
1553 ifb = x[I4(k, opt)];
1554 ifb_std = SQRT(cov[I4(k, opt) + I4(k, opt) * rtk->nx]);
1555 H[I4(k, opt) + nx * nv] = 1.0;
1556 }
1557
1558 if (opt->tropopt == TROPOPT_EST || opt->tropopt == TROPOPT_ESTG) {
1559 dtrp_std = SQRT(cov[IT(opt) + IT(opt) * rtk->nx]);
1560 for (k = 0; k < (opt->tropopt >= TROPOPT_ESTG ? 3 : 1); k++) {
1561 H[IT(opt) + k + nx * nv] = dtdx[k];
1562 }
1563 }
1564 if (opt->ionoopt == IONOOPT_EST) {
1565 if (rtk->x[II(sat, opt)] == 0.0) {
1566 continue;
1567 }
1568 dion_std = SQRT(cov[II(sat, opt) + II(sat, opt) * rtk->nx]);
1569 H[II(sat, opt) + nx * nv] = C;
1570 }
1571 if (j % 2 == 0) { /* phase bias */
1572 if ((bias = x[IB(sat, j / 2, opt)]) == 0.0) {
1573 continue;
1574 }
1575 bias_std = SQRT(cov[IB(obs[i].sat, j / 2, opt) + IB(obs[i].sat, j / 2, opt) * rtk->nx]);
1576 ;
1577 H[IB(sat, j / 2, opt) + nx * nv] = 1.0;
1578 }
1579 /* residual */
1580 v[nv] = y - (r + cdtr - CLIGHT * dts[i * 2] + dtrp + C * dion + bias + ifb);
1581
1582 if (j % 2 == 0) {
1583 rtk->ssat[sat - 1].resc[j / 2] = v[nv];
1584 } else {
1585 rtk->ssat[sat - 1].resp[j / 2] = v[nv];
1586 }
1587
1588 /* variance */
1589 var[nv] = varerr(obs[i].sat, sys, azel[1 + i * 2], j / 2, j % 2, opt) + vart + SQR(C) * vari +
1590 opt->uraratio * var_rs[i];
1591 if (sys == SYS_GLO && j % 2 == 1) {
1592 var[nv] += VAR_GLO_IFB;
1593 }
1594
1595 trace(NULL, post > 0 ? 3 : 4,
1596 "ppp_res: post=%d %s %s sat=%3d %s%s res=%9.4f sig=%9.4f el=%4.1f cdtr=%7.4f "
1597 "C=%7.4f dion=%7.4f (%5.2f) dtrp=%7.4f (%5.2f) bias=%9.4f (%5.2f) ifb=%9.4f (%5.2f) fix=%d\n",
1598 post, str, satstr, sat, j % 2 ? "C" : "L", code2obs(obs[i].code[j / 2]), v[nv], sqrt(var[nv]),
1599 azel[1 + i * 2] * R2D, cdtr, C, dion, dion_std, dtrp, dtrp_std, bias, bias_std, ifb, ifb_std,
1600 rtk->ssat[sat - 1].fix[j / 2]);
1601
1602 /* reject satellite by pre-fit residuals */
1603 if (!post && opt->maxinno > 0.0 && fabs(v[nv]) > opt->maxinno) {
1604 trace(NULL, 2, "outlier (%d) rejected %s %s sat=%2d %s%s res=%9.4f el=%4.1f\n", post, str, satstr, sat,
1605 j % 2 ? "C" : "L", code2obs(obs[i].code[j / 2]), v[nv], azel[1 + i * 2] * R2D);
1606 exc[i] = 1;
1607 rtk->ssat[sat - 1].rejc[j % 2]++;
1608 continue;
1609 }
1610 /* record large post-fit residuals */
1611 if (post && fabs(v[nv]) > sqrt(var[nv]) * THRES_REJECT) {
1612 obsi[ne] = i;
1613 frqi[ne] = j;
1614 ve[ne] = v[nv];
1615 ne++;
1616 }
1617 if (j % 2 == 0) {
1618 rtk->ssat[sat - 1].vsat[j / 2] = 1;
1619 }
1620 if (j % 2 == 0) {
1621 resc_rms += (v[nv] * v[nv]);
1622 resc_ave += v[nv];
1623 nresc++;
1624 } else {
1625 resp_rms += (v[nv] * v[nv]);
1626 resp_ave += v[nv];
1627 nresp++;
1628 }
1629 nv++;
1630 }
1631 }
1632 resp_rms = nresp ? sqrt(resp_rms / nresp) : 0.0;
1633 resp_ave = nresp ? resp_ave / nresp : 0.0;
1634 resc_rms = nresc ? sqrt(resc_rms / nresc) : 0.0;
1635 resc_ave = nresc ? resc_ave / nresc : 0.0;
1636 resrms[0] = resp_rms;
1637 resrms[1] = resc_rms;
1638 trace(NULL, 2,
1639 "ppp_res: post=%d, %s, nresp=, %2d, resp_rms=, %7.3f, resp_ave=, %7.3f, nresc=, %2d, resc_rms=, %7.3f, "
1640 "resc_ave=, %7.3f\n",
1641 post, str, nresp, resp_rms, resp_ave, nresc, resc_rms, resc_ave);
1642
1643 /* reject satellite with large and max post-fit residual */
1644 if (post && ne > 0) {
1645 vmax = ve[0];
1646 maxobs = obsi[0];
1647 maxfrq = frqi[0];
1648 rej = 0;
1649 for (j = 1; j < ne; j++) {
1650 if (fabs(vmax) >= fabs(ve[j])) {
1651 continue;
1652 }
1653 vmax = ve[j];
1654 maxobs = obsi[j];
1655 maxfrq = frqi[j];
1656 rej = j;
1657 }
1658 sat = obs[maxobs].sat;
1659 satno2id(sat, satstr);
1660 trace(NULL, 2, "outlier (%d) rejected %s %s sat=%2d %s%d res=%9.4f el=%4.1f\n", post, str, satstr, sat,
1661 maxfrq % 2 ? "P" : "L", maxfrq / 2 + 1, vmax, azel[1 + maxobs * 2] * R2D);
1662 exc[maxobs] = 1;
1663 rtk->ssat[sat - 1].rejc[maxfrq % 2]++;
1664 stat = 0;
1665 ve[rej] = 0;
1666 }
1667
1668 /* constraint to ionospheric correction */
1669 if (!post && opt->ionocorr) {
1670 nv += const_iono_corr(rtk, obs, nav, n, azel, exc, rr, x, v + nv, H + nv * rtk->nx, var + nv);
1671 }
1672
1673 for (i = 0; i < nv; i++) {
1674 for (j = 0; j < nv; j++) {
1675 R[i + j * nv] = i == j ? var[i] : 0.0;
1676 }
1677 }
1678 return post ? stat : nv;
1679}
1680/* number of estimated states ------------------------------------------------*/
1681extern int pppnx(const prcopt_t* opt) { return NX(opt); }
1682/* update solution status ----------------------------------------------------*/
1683static void update_stat(rtk_t* rtk, const obsd_t* obs, int n, int stat) {
1684 const prcopt_t* opt = &rtk->opt;
1685 int i, j;
1686
1687 /* test # of valid satellites */
1688 rtk->sol.ns = 0;
1689 for (i = 0; i < n && i < MAXOBS; i++) {
1690 for (j = 0; j < opt->nf && j < NFREQ; j++) {
1691 if (!rtk->ssat[obs[i].sat - 1].vsat[j]) {
1692 continue;
1693 }
1694 rtk->ssat[obs[i].sat - 1].lock[j]++;
1695 rtk->ssat[obs[i].sat - 1].outc[j] = 0;
1696 if (j == 0) {
1697 rtk->sol.ns++;
1698 }
1699 }
1700 }
1701
1702 if (stat == SOLQ_FIX || stat == SOLQ_FIX_WL) {
1703 for (i = 0; i < 3; i++) {
1704 rtk->sol.rr[i] = rtk->xa[i];
1705 rtk->sol.qr[i] = (float)rtk->Pa[i + i * rtk->na];
1706 }
1707 rtk->sol.qr[3] = (float)rtk->Pa[1];
1708 rtk->sol.qr[4] = (float)rtk->Pa[1 + 2 * rtk->na];
1709 rtk->sol.qr[5] = (float)rtk->Pa[2];
1710 } else {
1711 for (i = 0; i < 3; i++) {
1712 rtk->sol.rr[i] = rtk->x[i];
1713 rtk->sol.qr[i] = (float)rtk->P[i + i * rtk->nx];
1714 }
1715 rtk->sol.qr[3] = (float)rtk->P[1];
1716 rtk->sol.qr[4] = (float)rtk->P[2 + rtk->nx];
1717 rtk->sol.qr[5] = (float)rtk->P[2];
1718 }
1719 if (stat == SOLQ_FIX_WL) {
1720 stat = SOLQ_PPP;
1721 }
1722 rtk->sol.stat = rtk->sol.ns < MIN_NSAT_SOL ? SOLQ_NONE : stat;
1723
1724 rtk->sol.dtr[0] = rtk->x[IC(0, opt)] / CLIGHT; /* [s] */
1725 for (i = 1; i < 7; i++) {
1726 rtk->sol.dtr[i] = (rtk->x[IC(i, opt)] - rtk->x[IC(0, opt)]) / CLIGHT;
1727 }
1728
1729 for (i = 0; i < n && i < MAXOBS; i++) {
1730 for (j = 0; j < opt->nf && j < NFREQ; j++) {
1731 rtk->ssat[obs[i].sat - 1].snr[j] = obs[i].SNR[j];
1732 }
1733 }
1734 /* update previous phase/time for Doppler slip detection */
1735 for (i = 0; i < n && i < MAXOBS; i++) {
1736 for (j = 0; j < opt->nf && j < NFREQ; j++) {
1737 if (obs[i].L[j] != 0.0) {
1738 rtk->ssat[obs[i].sat - 1].pt[0][j] = obs[i].time;
1739 rtk->ssat[obs[i].sat - 1].ph[0][j] = obs[i].L[j];
1740 }
1741 }
1742 }
1743 for (i = 0; i < MAXSAT; i++) {
1744 for (j = 0; j < opt->nf && j < NFREQ; j++) {
1745 if (rtk->ssat[i].slip[j] & 3) {
1746 rtk->ssat[i].slipc[j]++;
1747 }
1748 }
1749 }
1750}
1751/* test hold ambiguity -------------------------------------------------------*/
1752static int test_hold_amb(rtk_t* rtk) {
1753 int i, j, stat = 0;
1754
1755 /* no fix-and-hold mode */
1756 if (rtk->opt.modear != ARMODE_FIXHOLD) {
1757 return 0;
1758 }
1759
1760 /* reset # of continuous fixed if new ambiguity introduced */
1761 for (i = 0; i < MAXSAT; i++) {
1762 if (rtk->ssat[i].fix[0] != 2) {
1763 continue;
1764 }
1765 for (j = 0; j < MAXSAT; j++) {
1766 if (rtk->ssat[j].fix[0] != 2) {
1767 continue;
1768 }
1769 if (!rtk->ambc[j].flags[i] || !rtk->ambc[i].flags[j]) {
1770 stat = 1;
1771 }
1772 rtk->ambc[j].flags[i] = rtk->ambc[i].flags[j] = 1;
1773 }
1774 }
1775 if (stat) {
1776 rtk->nfix = 0;
1777 return 0;
1778 }
1779 /* test # of continuous fixed */
1780 return ++rtk->nfix >= ROUND(rtk->opt.minfix / fabs(rtk->tt));
1781}
1782/* precise point positioning -------------------------------------------------*/
1783extern void pppos(mrtk_ctx_t* ctx, rtk_t* rtk, const obsd_t* obs, int n, const nav_t* nav) {
1784 const prcopt_t* opt = &rtk->opt;
1785 double *rs, *dts, *var, *v, *H, *R, *azel, *xp, *Pp, dr[3] = {0}, std[3] = {0.0};
1786 double resrms[2] = {0}, resrms_ar[2] = {0};
1787 char str[32];
1788 int i, j, nv, info, svh[MAXOBS], exc[MAXOBS] = {0}, stat = SOLQ_SINGLE, ret;
1789
1790 time2str(obs[0].time, str, 2);
1791 trace(NULL, 3, "pppos : time=%s nx=%d n=%d\n", str, rtk->nx, n);
1792
1793 rs = mat(6, n);
1794 dts = mat(2, n);
1795 var = mat(1, n);
1796 azel = zeros(2, n);
1797
1798 for (i = 0; i < MAXSAT; i++) {
1799 for (j = 0; j < opt->nf && j < NFREQ; j++) {
1800 rtk->ssat[i].fix[j] = 0;
1801 }
1802 }
1803
1804 /* temporal update of ekf states */
1805 udstate_ppp(rtk, obs, n, nav);
1806
1807 /* satellite positions and clocks */
1808 satposs(obs[0].time, obs, n, nav, rtk->opt.sateph, rs, dts, var, svh);
1809
1810 /* exclude measurements of eclipsing satellite (block IIA) */
1811 if (rtk->opt.posopt[3]) {
1812 testeclipse(obs, n, nav, rs);
1813 }
1814 /* earth tides correction */
1815 if (opt->tidecorr) {
1816 tidedisp(gpst2utc(obs[0].time), rtk->x, opt->tidecorr == 1 ? 1 : 7, &nav->erp, opt->odisp[0], dr);
1817 }
1818 nv = n * rtk->opt.nf * 2 + MAXSAT + 3;
1819 xp = mat(rtk->nx, 1);
1820 Pp = zeros(rtk->nx, rtk->nx);
1821 v = mat(nv, 1);
1822 H = mat(rtk->nx, nv);
1823 R = mat(nv, nv);
1824
1825 for (i = 0; i < MAX_ITER; i++) {
1826 matcpy(xp, rtk->x, rtk->nx, 1);
1827 matcpy(Pp, rtk->P, rtk->nx, rtk->nx);
1828
1829 /* prefit residuals */
1830 if (!(nv = ppp_res(0, obs, n, rs, dts, var, svh, dr, exc, nav, xp, Pp, rtk, v, H, R, azel, resrms))) {
1831 trace(NULL, 2, "%s ppp (%d) no valid obs data\n", str, i + 1);
1832 break;
1833 }
1834 /* measurement update of ekf states */
1835 if ((info = filter(xp, Pp, H, v, R, rtk->nx, nv))) {
1836 trace(NULL, 2, "%s ppp (%d) filter error info=%d\n", str, i + 1, info);
1837 break;
1838 }
1839 /* postfit residuals */
1840 if (ppp_res(i + 1, obs, n, rs, dts, var, svh, dr, exc, nav, xp, Pp, rtk, v, H, R, azel, resrms)) {
1841 matcpy(rtk->x, xp, rtk->nx, 1);
1842 matcpy(rtk->P, Pp, rtk->nx, rtk->nx);
1843 stat = SOLQ_PPP;
1844 break;
1845 }
1846 }
1847 if (i >= MAX_ITER) {
1848 trace(NULL, 2, "%s ppp (%d) iteration overflows\n", str, i);
1849 }
1850 if (stat == SOLQ_PPP) {
1851 /* ambiguity resolution in ppp */
1852 if ((ret = ppp_ar(ctx, rtk, obs, n, exc, nav, azel, xp, Pp)) &&
1853 ppp_res(9, obs, n, rs, dts, var, svh, dr, exc, nav, xp, Pp, rtk, v, H, R, azel, resrms_ar)) {
1854 matcpy(rtk->xa, xp, rtk->nx, 1);
1855 matcpy(rtk->Pa, Pp, rtk->nx, rtk->nx);
1856
1857 for (i = 0; i < 3; i++) {
1858 std[i] = sqrt(Pp[i + i * rtk->nx]);
1859 }
1860
1861 if (ret > 1) {
1862 stat = SOLQ_FIX; /* narrow-lane fixed */
1863 } else {
1864 stat = SOLQ_FIX_WL; /* only wide-lane fixed */
1865 }
1866
1867 trace(NULL, 2,
1868 "pppos(%4d): time=%s stat=%d std=%.3f resp=%7.3f -> %7.3f (%6.3f) resc=%7.3f -> %7.3f (%6.3f)\n",
1869 __LINE__, str, stat, norm(std, 3), resrms[0], resrms_ar[0], resrms_ar[0] / resrms[0], resrms[1],
1870 resrms_ar[1], resrms_ar[1] / resrms[1]);
1871 } else {
1872 rtk->nfix = 0;
1873 }
1874 /* update solution status */
1875 update_stat(rtk, obs, n, stat);
1876
1877 /* hold fixed ambiguities */
1878 if (stat == SOLQ_FIX && test_hold_amb(rtk)) {
1879 matcpy(rtk->x, xp, rtk->nx, 1);
1880 matcpy(rtk->P, Pp, rtk->nx, rtk->nx);
1881 trace(NULL, 3, "pppos(%4d): %s hold ambiguity nfix=%4d / %4d\n", __LINE__, str, rtk->nfix, rtk->opt.minfix);
1882 }
1883 }
1884 free(rs);
1885 free(dts);
1886 free(var);
1887 free(azel);
1888 free(xp);
1889 free(Pp);
1890 free(v);
1891 free(H);
1892 free(R);
1893}
mrtk_antenna.h
MRTKLIB Antenna Module — Antenna phase center model and parameter I/O.
mrtk_astro.h
MRTKLIB Astronomy Module — Sun and moon position computation.
ERR_SAAS
#define ERR_SAAS
Definition mrtk_atmos.c:294
REL_HUMI
#define REL_HUMI
Definition mrtk_atmos.c:296
ERR_BRDCI
#define ERR_BRDCI
Definition mrtk_atmos.c:295
mrtk_atmos.h
MRTKLIB Atmosphere Module — Ionospheric and tropospheric models.
getsysno
int getsysno(int sat)
Get system code number based on satellite number.
Definition mrtk_bias_sinex.c:47
code2obs
char * code2obs(uint8_t code)
Definition mrtk_obs.c:190
mrtk_bias_sinex.h
MRTKLIB Bias SINEX Module — BIAS-SINEX data types and I/O functions.
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
SSR_INVALID_PBIAS
#define SSR_INVALID_PBIAS
Definition mrtk_const.h:427
EFACT_GAL
#define EFACT_GAL
Definition mrtk_const.h:169
EFACT_QZS
#define EFACT_QZS
Definition mrtk_const.h:172
SSR_INVALID_CBIAS
#define SSR_INVALID_CBIAS
Definition mrtk_const.h:426
mrtk_context.h
MRTKLIB Context — runtime context for processing, trace, and error handling.
D2R
#define D2R
Definition mrtk_coords.c:36
CLIGHT
#define CLIGHT
Definition mrtk_coords.c:41
PI
#define PI
Definition mrtk_coords.c:35
MAXSAT
#define MAXSAT
Definition mrtk_coords.c:45
RE_WGS84
#define RE_WGS84
Definition mrtk_coords.c:38
OMGE
#define OMGE
Definition mrtk_coords.c:40
mrtk_coords.h
MRTKLIB Coordinates Module — Coordinate transformations and geometry.
EPHOPT_PREC
#define EPHOPT_PREC
Definition mrtk_eph.c:84
SYS_GPS
#define SYS_GPS
Definition mrtk_eph.c:65
SYS_SBS
#define SYS_SBS
Definition mrtk_eph.c:66
EPHOPT_SSRAPC
#define EPHOPT_SSRAPC
Definition mrtk_eph.c:86
SYS_QZS
#define SYS_QZS
Definition mrtk_eph.c:69
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
EPHOPT_SSRCOM
#define EPHOPT_SSRCOM
Definition mrtk_eph.c:87
nav_t
struct nav_s nav_t
Definition mrtk_eph.h:40
SNR_UNIT
#define SNR_UNIT
Definition mrtk_foundation.h:54
CODE_L5Q
#define CODE_L5Q
Definition mrtk_foundation.h:210
MAXOBS
#define MAXOBS
Definition mrtk_foundation.h:151
CODE_L5X
#define CODE_L5X
Definition mrtk_foundation.h:211
CODE_L1C
#define CODE_L1C
Definition mrtk_foundation.h:186
CODE_L2C
#define CODE_L2C
Definition mrtk_foundation.h:199
DTTOL
#define DTTOL
Definition mrtk_foundation.h:158
CODE_NONE
#define CODE_NONE
Definition mrtk_foundation.h:185
NFREQPCV
#define NFREQPCV
Definition mrtk_foundation.h:48
has_sel_biascode
int has_sel_biascode(int sys, uint8_t code)
Select the HAS bias slot obs code for an observation code.
Definition mrtk_has.c:902
mrtk_has.h
MRTKLIB Galileo HAS Module — E6-B C/NAV HAS page collection and Message Type 1 (MT1) decode into ssr_...
mrtk_ionex.h
MRTKLIB IONEX Module — IONEX ionospheric TEC map I/O and interpolation.
mrtk_madoca.h
MRTKLIB MADOCA-PPP Module — QZSS L6E CSSR message decode functions.
mrtk_mat.h
MRTKLIB Matrix Module — Matrix/vector operations and estimation.
freq_idx2ant_idx
int freq_idx2ant_idx(int sys, int freq_idx)
Convert frequency index to antenna PCV index.
Definition mrtk_nav.c:743
satexclude
int satexclude(int sat, double var, int svh, const struct prcopt_t *opt)
Test excluded satellite.
Definition mrtk_nav.c:1169
R2D
#define R2D
Definition mrtk_obs.c:51
SYS_BD2
#define SYS_BD2
Definition mrtk_obs.c:62
FREQ1
#define FREQ1
Definition mrtk_obs.c:67
PMODE_PPP_RTK
#define PMODE_PPP_RTK
Definition mrtk_opt.h:54
SOLQ_SINGLE
#define SOLQ_SINGLE
Definition mrtk_opt.h:75
IONOOPT_EST
#define IONOOPT_EST
Definition mrtk_opt.h:96
TROPOPT_SBAS
#define TROPOPT_SBAS
Definition mrtk_opt.h:104
TROPOPT_SAAS
#define TROPOPT_SAAS
Definition mrtk_opt.h:103
CORR_IGS_RTS
#define CORR_IGS_RTS
Definition mrtk_opt.h:118
ARMODE_FIXHOLD
#define ARMODE_FIXHOLD
Definition mrtk_opt.h:131
CORR_GAL_HAS
#define CORR_GAL_HAS
Definition mrtk_opt.h:120
TROPOPT_EST
#define TROPOPT_EST
Definition mrtk_opt.h:105
IONOOPT_SBAS
#define IONOOPT_SBAS
Definition mrtk_opt.h:94
PMODE_PPP_STATIC
#define PMODE_PPP_STATIC
Definition mrtk_opt.h:52
PMODE_PPP_FIXED
#define PMODE_PPP_FIXED
Definition mrtk_opt.h:53
CORR_NONE
#define CORR_NONE
Definition mrtk_opt.h:116
TROPOPT_ESTG
#define TROPOPT_ESTG
Definition mrtk_opt.h:106
SOLQ_NONE
#define SOLQ_NONE
Definition mrtk_opt.h:70
SOLQ_PPP
#define SOLQ_PPP
Definition mrtk_opt.h:76
IONOOPT_BRDC
#define IONOOPT_BRDC
Definition mrtk_opt.h:93
CORR_IGS
#define CORR_IGS
Definition mrtk_opt.h:117
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
IONOOPT_TEC
#define IONOOPT_TEC
Definition mrtk_opt.h:97
mrtk_peph.h
MRTKLIB Precise Ephemeris Module — Earth rotation parameter types and ERP I/O functions.
satno2id
void satno2id(int sat, char *id)
Definition mrtk_nav.c:817
geodist
double geodist(const double *rs, const double *rr, double *e)
Definition mrtk_coords.c:493
detslp_ll
static void detslp_ll(rtk_t *rtk, const obsd_t *obs, int n)
Definition mrtk_ppp.c:844
udifb_ppp
static void udifb_ppp(rtk_t *rtk)
Definition mrtk_ppp.c:1181
satsys
int satsys(int sat, int *prn)
Definition mrtk_nav.c:117
NX
#define NX(opt)
Definition mrtk_ppp.c:224
SOLQ_FIX_WL
#define SOLQ_FIX_WL
Definition mrtk_ppp.c:186
mat
double * mat(int n, int m)
Definition mrtk_mat.c:76
THRES_REJECT
#define THRES_REJECT
Definition mrtk_ppp.c:190
VAR_CLK
#define VAR_CLK
Definition mrtk_ppp.c:197
VAR_IONO
#define VAR_IONO
Definition mrtk_ppp.c:202
detslp_dop
static void detslp_dop(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
Definition mrtk_ppp.c:773
const_iono_corr
int const_iono_corr(rtk_t *rtk, const obsd_t *obs, const nav_t *nav, const int n, const double *azel, const int *exc, const double *rr, const double *x, double *v, double *H, double *var)
Definition mrtk_ppp_iono.c:259
initx
static void initx(rtk_t *rtk, double xi, double var, int i)
Definition mrtk_ppp.c:493
IC
#define IC(s, opt)
Definition mrtk_ppp.c:225
zeros
double * zeros(int n, int m)
Definition mrtk_mat.c:108
VAR_VEL
#define VAR_VEL
Definition mrtk_ppp.c:195
sunmoonpos
void sunmoonpos(gtime_t tut, const double *erpv, double *rsun, double *rmoon, double *gmst)
Definition mrtk_astro.c:135
satazel
double satazel(const double *pos, const double *e, double *azel)
Definition mrtk_coords.c:517
matcpy
void matcpy(double *A, const double *B, int n, int m)
Definition mrtk_mat.c:205
tidedisp
void tidedisp(gtime_t tutc, const double *rr, int opt, const erp_t *erp, const double *odisp, double *dr)
Definition mrtk_tides.c:261
tickget
unsigned int tickget(void)
Definition mrtk_time.c:566
test_hold_amb
static int test_hold_amb(rtk_t *rtk)
Definition mrtk_ppp.c:1752
cross3
void cross3(const double *a, const double *b, double *c)
Definition mrtk_mat.c:172
xyz2enu
void xyz2enu(const double *pos, double *E)
Definition mrtk_coords.c:142
covecef
void covecef(const double *pos, const double *Q, double *P)
Definition mrtk_coords.c:202
tropmodel
double tropmodel(gtime_t time, const double *pos, const double *azel, double humi)
Definition mrtk_atmos.c:162
NC
#define NC(opt)
Definition mrtk_ppp.c:213
pppoutstat
int pppoutstat(rtk_t *rtk, char *buff)
Write PPP solution status to buffer.
Definition mrtk_ppp.c:240
tropmapf
double tropmapf(gtime_t time, const double *pos, const double *azel, double *mapfw)
mcssr_sel_biascode
int mcssr_sel_biascode(const int sys, const int code)
Definition mrtk_madoca.c:988
iontec
int iontec(gtime_t time, const nav_t *nav, const double *pos, const double *azel, int opt, double *delay, double *var)
Definition mrtk_ionex.c:481
udclk_ppp
static void udclk_ppp(rtk_t *rtk)
Definition mrtk_ppp.c:1065
testsnr
int testsnr(int base, int freq, double el, double snr, const snrmask_t *mask)
Definition mrtk_obs.c:132
sbsioncorr
int sbsioncorr(gtime_t time, const nav_t *nav, const double *pos, const double *azel, double *delay, double *var)
Definition mrtk_sbas.c:823
GAP_RESION
#define GAP_RESION
Definition mrtk_ppp.c:208
antmodel
void antmodel(const pcv_t *pcv, const double *del, const double *azel, int opt, double *dant)
Definition mrtk_antenna.c:481
testeclipse
static void testeclipse(const obsd_t *obs, int n, const nav_t *nav, double *rs)
Definition mrtk_ppp.c:322
imat
int * imat(int n, int m)
Definition mrtk_mat.c:92
time2str
void time2str(gtime_t t, char *s, int n)
Definition mrtk_time.c:507
gpst2utc
gtime_t gpst2utc(gtime_t t)
Definition mrtk_time.c:422
NF
#define NF(opt)
Definition mrtk_ppp.c:211
covenu
void covenu(const double *pos, const double *P, double *Q)
Definition mrtk_coords.c:188
IT
#define IT(opt)
Definition mrtk_ppp.c:226
update_stat
static void update_stat(rtk_t *rtk, const obsd_t *obs, int n, int stat)
Definition mrtk_ppp.c:1683
filter
int filter(double *x, double *P, const double *H, const double *v, const double *R, int n, int m)
Definition mrtk_mat.c:524
yaw_angle
int yaw_angle(int sat, const char *type, int opt, double beta, double mu, double *yaw)
Definition mrtk_ppp.c:370
udbias_ppp
static void udbias_ppp(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
Definition mrtk_ppp.c:1211
STD
static double STD(rtk_t *rtk, int i)
Definition mrtk_ppp.c:233
ppp_note_autosig
static void ppp_note_autosig(int sys, uint8_t code)
Definition mrtk_ppp.c:527
pppos
void pppos(mrtk_ctx_t *ctx, rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
Precise point positioning.
Definition mrtk_ppp.c:1783
time2gpst
double time2gpst(gtime_t t, int *week)
Definition mrtk_time.c:181
satid2no
int satid2no(const char *id)
Definition mrtk_nav.c:753
I4
#define I4(s, opt)
Definition mrtk_ppp.c:229
pppnx
int pppnx(const prcopt_t *opt)
Number of estimated states for PPP.
Definition mrtk_ppp.c:1681
N3
#define N3(opt)
Definition mrtk_ppp.c:220
VAR_GRA
#define VAR_GRA
Definition mrtk_ppp.c:199
N4
#define N4(opt)
Definition mrtk_ppp.c:221
ionmodel
double ionmodel(gtime_t t, const double *ion, const double *pos, const double *azel)
Definition mrtk_atmos.c:63
VAR_IFB
#define VAR_IFB
Definition mrtk_ppp.c:200
SQR
#define SQR(x)
Definition mrtk_ppp.c:179
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)
Definition mrtk_eph.c:940
sat_yaw
static int sat_yaw(gtime_t time, int sat, const char *type, int opt, const double *rs, double *exs, double *eys)
Definition mrtk_ppp.c:375
MAX
#define MAX(x, y)
Definition mrtk_ppp.c:181
sbstropcorr
double sbstropcorr(gtime_t time, const double *pos, const double *azel, double *var)
Definition mrtk_sbas.c:930
model_iono
static int model_iono(gtime_t time, const double *pos, const double *azel, const prcopt_t *opt, int sat, const double *x, const nav_t *nav, double *dion, double *var)
Definition mrtk_ppp.c:1410
code2obs
char * code2obs(uint8_t code)
Definition mrtk_obs.c:190
model_phw
static int model_phw(gtime_t time, int sat, const char *type, int opt, const double *rs, const double *rr, double *phw)
Definition mrtk_ppp.c:416
norm
double norm(const double *a, int n)
Definition mrtk_mat.c:165
satsys_bd2
int satsys_bd2(int sat, int *prn)
Definition mrtk_nav.c:159
udpos_ppp
static void udpos_ppp(rtk_t *rtk)
Definition mrtk_ppp.c:965
time_str
char * time_str(gtime_t t, int n)
Definition mrtk_time.c:530
normv3
int normv3(const double *a, double *b)
Definition mrtk_mat.c:183
VAR_GLO_IFB
#define VAR_GLO_IFB
Definition mrtk_ppp.c:203
detslp_gf
static void detslp_gf(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
Definition mrtk_ppp.c:864
MAX_ITER
#define MAX_ITER
Definition mrtk_ppp.c:188
ROUND
#define ROUND(x)
Definition mrtk_ppp.c:183
IB
#define IB(s, f, opt)
Definition mrtk_ppp.c:230
udtrop_ppp
static void udtrop_ppp(rtk_t *rtk)
Definition mrtk_ppp.c:1084
yaw_nominal
static double yaw_nominal(double beta, double mu)
Definition mrtk_ppp.c:363
I3
#define I3(s, opt)
Definition mrtk_ppp.c:228
satexclude
int satexclude(int sat, double var, int svh, const prcopt_t *opt)
mwmeas
static double mwmeas(const obsd_t *obs, const nav_t *nav, int f)
Definition mrtk_ppp.c:512
model_trop
static int model_trop(gtime_t time, const double *pos, const double *azel, const prcopt_t *opt, const double *x, double *dtdx, const nav_t *nav, double *dtrp, double *var)
Definition mrtk_ppp.c:1389
detslp_ssr
static void detslp_ssr(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
Definition mrtk_ppp.c:922
matmul
void matmul(const char *tr, int n, int k, int m, double alpha, const double *A, const double *B, double beta, double *C)
Definition mrtk_mat.c:310
dot
double dot(const double *a, const double *b, int n)
Definition mrtk_mat.c:151
eye
double * eye(int n)
Definition mrtk_mat.c:129
varerr
static double varerr(int sat, int sys, double el, int idx, int type, const prcopt_t *opt)
Definition mrtk_ppp.c:472
VAR_POS
#define VAR_POS
Definition mrtk_ppp.c:194
VAR_ZTD
#define VAR_ZTD
Definition mrtk_ppp.c:198
SQRT
#define SQRT(x)
Definition mrtk_ppp.c:180
gfmeas
static double gfmeas(const obsd_t *obs, const nav_t *nav, int f)
Definition mrtk_ppp.c:501
udiono_ppp
static void udiono_ppp(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
Definition mrtk_ppp.c:1111
corr_meas
static int corr_meas(const obsd_t *obs, const nav_t *nav, const double *azel, const prcopt_t *opt, const double *dantr, const double *dants, double phw, double *L, double *P, double *Lc, double *Pc, int tl)
Definition mrtk_ppp.c:564
ecef2pos
void ecef2pos(const double *r, double *pos)
Definition mrtk_coords.c:102
MIN_NSAT_SOL
#define MIN_NSAT_SOL
Definition mrtk_ppp.c:189
THRES_MW_JUMP
#define THRES_MW_JUMP
Definition mrtk_ppp.c:192
detslp_code
static void detslp_code(rtk_t *rtk, const obsd_t *obs, int n)
Definition mrtk_ppp.c:823
VAR_ACC
#define VAR_ACC
Definition mrtk_ppp.c:196
udstate_ppp
static void udstate_ppp(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
Definition mrtk_ppp.c:1317
VAR_BIAS
#define VAR_BIAS
Definition mrtk_ppp.c:201
ppp_res
static int ppp_res(int post, const obsd_t *obs, int n, const double *rs, const double *dts, const double *var_rs, const int *svh, const double *dr, int *exc, const nav_t *nav, const double *x, const double *cov, rtk_t *rtk, double *v, double *H, double *R, double *azel, double *resrms)
Definition mrtk_ppp.c:1435
antmodel_s
void antmodel_s(const pcv_t *pcv, double nadir, double *dant)
Definition mrtk_antenna.c:507
trop_model_prec
static double trop_model_prec(gtime_t time, const double *pos, const double *azel, const double *x, double *dtdx, double *var)
Definition mrtk_ppp.c:1364
satantoff
void satantoff(gtime_t time, const double *rs, int sat, const nav_t *nav, double *dant)
Definition mrtk_peph.c:785
satantpcv
static void satantpcv(const double *rs, const double *rr, const pcv_t *pcv, double *dant)
Definition mrtk_ppp.c:1343
sat2freq
double sat2freq(int sat, uint8_t code, const nav_t *nav)
Definition mrtk_obs.c:244
ecef2enu
void ecef2enu(const double *pos, const double *r, double *e)
Definition mrtk_coords.c:162
detslp_mw
static void detslp_mw(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
Definition mrtk_ppp.c:893
II
#define II(s, opt)
Definition mrtk_ppp.c:227
mrtk_ppp.h
MRTKLIB Precise Point Positioning (PPP) Module.
ppp_ar
int ppp_ar(mrtk_ctx_t *ctx, rtk_t *rtk, const obsd_t *obs, int n, int *exc, const nav_t *nav, const double *azel, double *x, double *P)
Ambiguity resolution in PPP.
Definition mrtk_ppp_ar.c:1003
mrtk_ppp_ar.h
MRTKLIB PPP Ambiguity Resolution Module.
NT
#define NT(opt)
Definition mrtk_ppp_iono.c:63
ppprtk_outstat
int ppprtk_outstat(rtk_t *rtk, char *buff)
Write PPP-RTK solution status to buffer.
Definition mrtk_ppp_rtk.c:2423
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
mrtk_sbas.h
MRTKLIB SBAS Module — SBAS message decoding and corrections.
mrtk_tides.h
MRTKLIB Tides Module — Tidal displacement corrections.
timediff
double timediff(gtime_t t1, gtime_t t2)
Compute time difference (t1 - t2).
Definition mrtk_time.c:279
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.
ambc_t::flags
char flags[MAXSAT]
Definition mrtk_rtkpos.h:46
erp_t
ERP (earth rotation parameter) type.
Definition mrtk_peph.h:59
gtime_t
GNSS time representation.
Definition mrtk_time.h:49
gtime_t::time
time_t time
Definition mrtk_time.h:50
mrtk_ctx_t
Runtime context for MRTKLIB processing.
Definition mrtk_context.h:73
nav_s::erp
erp_t erp
Definition mrtk_nav.h:413
nav_s::cbias
double cbias[MAXSAT][3]
Definition mrtk_nav.h:427
nav_s::pcvs
pcv_t pcvs[MAXSAT]
Definition mrtk_nav.h:430
nav_s::ssr_ch
ssr_t ssr_ch[SSR_CH_NUM][MAXSAT]
Definition mrtk_nav.h:434
nav_s::ion_gps
double ion_gps[8]
Definition mrtk_nav.h:421
nav_s::osb
osb_t osb
Definition mrtk_nav.h:441
obsd_t
Observation data record.
Definition mrtk_obs.h:44
obsd_t::D
float D[5+5]
Definition mrtk_obs.h:52
obsd_t::time
gtime_t time
Definition mrtk_obs.h:45
obsd_t::sat
uint8_t sat
Definition mrtk_obs.h:46
obsd_t::L
double L[5+5]
Definition mrtk_obs.h:50
obsd_t::P
double P[5+5]
Definition mrtk_obs.h:51
obsd_t::SNR
uint16_t SNR[5+5]
Definition mrtk_obs.h:47
obsd_t::code
uint8_t code[5+5]
Definition mrtk_obs.h:49
osb_t::rsyscb
double rsyscb[MAXBSNXSYS][MAXCODE]
Definition mrtk_nav.h:203
osb_t::scb
double scb[MAXSAT][MAXCODE]
Definition mrtk_nav.h:201
osb_t::vrsatcb
int vrsatcb[MAXSAT][MAXCODE]
Definition mrtk_nav.h:208
osb_t::vspb
int vspb[MAXSAT][MAXCODE]
Definition mrtk_nav.h:206
osb_t::vscb
int vscb[MAXSAT][MAXCODE]
Definition mrtk_nav.h:205
osb_t::vrsyscb
int vrsyscb[MAXBSNXSYS][MAXCODE]
Definition mrtk_nav.h:207
osb_t::spb
double spb[MAXSAT][MAXCODE]
Definition mrtk_nav.h:202
osb_t::rsatcb
double rsatcb[MAXSAT][MAXCODE]
Definition mrtk_nav.h:204
pcv_t
Definition mrtk_nav.h:45
pcv_t::type
char type[MAXANT]
Definition mrtk_nav.h:47
prcopt_t
Definition mrtk_opt.h:171
prcopt_t::clockjump
int clockjump
Definition mrtk_opt.h:324
prcopt_t::ionocorr
int ionocorr
Definition mrtk_opt.h:248
prcopt_t::thresdop
double thresdop
Definition mrtk_opt.h:215
prcopt_t::snrmask
snrmask_t snrmask
Definition mrtk_opt.h:177
prcopt_t::minfix
int minfix
Definition mrtk_opt.h:185
prcopt_t::elmin
double elmin
Definition mrtk_opt.h:176
prcopt_t::nf
int nf
Definition mrtk_opt.h:174
prcopt_t::ionoopt
int ionoopt
Definition mrtk_opt.h:187
prcopt_t::pcvr
pcv_t pcvr[2]
Definition mrtk_opt.h:224
prcopt_t::thresslip
double thresslip
Definition mrtk_opt.h:214
prcopt_t::prn
double prn[7]
Definition mrtk_opt.h:209
prcopt_t::dynamics
int dynamics
Definition mrtk_opt.h:189
prcopt_t::err
double err[8]
Definition mrtk_opt.h:201
prcopt_t::eratio
double eratio[5]
Definition mrtk_opt.h:200
prcopt_t::tropopt
int tropopt
Definition mrtk_opt.h:188
prcopt_t::odisp
double odisp[2][6 *11]
Definition mrtk_opt.h:239
prcopt_t::unbias
int unbias
Definition mrtk_opt.h:231
prcopt_t::pppopt
char pppopt[256]
Definition mrtk_opt.h:243
prcopt_t::posopt
int posopt[13]
Definition mrtk_opt.h:237
prcopt_t::uraratio
double uraratio
Definition mrtk_opt.h:249
prcopt_t::maxout
int maxout
Definition mrtk_opt.h:183
prcopt_t::modear
int modear
Definition mrtk_opt.h:179
prcopt_t::mode
int mode
Definition mrtk_opt.h:172
prcopt_t::tidecorr
int tidecorr
Definition mrtk_opt.h:190
prcopt_t::antdel
double antdel[2][3]
Definition mrtk_opt.h:223
prcopt_t::ru
double ru[3]
Definition mrtk_opt.h:220
prcopt_t::correction
int correction
Definition mrtk_opt.h:300
prcopt_t::sateph
int sateph
Definition mrtk_opt.h:178
prcopt_t::maxinno
double maxinno
Definition mrtk_opt.h:217
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::ambc
ambc_t ambc[MAXSAT]
Definition mrtk_rtkpos.h:59
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
snrmask_t
SNR mask type.
Definition mrtk_obs.h:67
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::dtr
double dtr[NSYS+1]
Definition mrtk_sol.h:52
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::mw
double mw[5]
Definition mrtk_sol.h:114
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::discont
int discont[5]
Definition mrtk_sol.h:122
ssat_t::vs
uint8_t vs
Definition mrtk_sol.h:100
ssat_t::ph
double ph[2][5]
Definition mrtk_sol.h:117
ssat_t::outc
uint32_t outc[5]
Definition mrtk_sol.h:110
ssat_t::rejc
uint32_t rejc[5]
Definition mrtk_sol.h:112
ssat_t::codeprev
uint8_t codeprev[5][2]
Definition mrtk_sol.h:126
ssat_t::azel
double azel[2]
Definition mrtk_sol.h:101
ssat_t::phw
double phw
Definition mrtk_sol.h:115
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::ionc
double ionc
Definition mrtk_sol.h:123
ssat_t::fix
uint8_t fix[5]
Definition mrtk_sol.h:106
ssat_t::resc
double resc[5]
Definition mrtk_sol.h:103
ssr_t::t0
gtime_t t0[6]
Definition mrtk_nav.h:176
ssr_t::vpbias
int vpbias[MAXCODE]
Definition mrtk_nav.h:191
ssr_t::vcbias
int vcbias[MAXCODE]
Definition mrtk_nav.h:190
ssr_t::pbias
double pbias[MAXCODE]
Definition mrtk_nav.h:188
ssr_t::cbias
float cbias[MAXCODE]
Definition mrtk_nav.h:187
ssr_t::discnt
int discnt[MAXCODE]
Definition mrtk_nav.h:192
free
#define free(x)
Definition toml.c:53