diff --git a/app/qtapp/appcmn_qt/navi_post_opt.cpp b/app/qtapp/appcmn_qt/navi_post_opt.cpp
index ca9c89be6..d2324e866 100644
--- a/app/qtapp/appcmn_qt/navi_post_opt.cpp
+++ b/app/qtapp/appcmn_qt/navi_post_opt.cpp
@@ -77,6 +77,7 @@ OptDialog::OptDialog(QWidget *parent, int opts)
ui->cBIonosphereOption->setItemData(IONOOPT_EST, tr("Estimate ionospheric parameter STEC (slant total electron content)"), Qt::ToolTipRole);
ui->cBIonosphereOption->setItemData(IONOOPT_TEC, tr("Use IONEX TEC grid data"), Qt::ToolTipRole);
ui->cBIonosphereOption->setItemData(IONOOPT_QZS, tr("Apply broadcast ionosphere model provided by QZSS"), Qt::ToolTipRole);
+ ui->cBIonosphereOption->setItemData(IONOOPT_VTEC, tr("Apply SSR IONEX VTEC"), Qt::ToolTipRole);
ui->cBTroposphereOption->setItemData(TROPOPT_OFF, tr("Not apply troposphere correction"), Qt::ToolTipRole);
ui->cBTroposphereOption->setItemData(TROPOPT_SAAS, tr("Apply Saastamoinen model"), Qt::ToolTipRole);
diff --git a/app/qtapp/appcmn_qt/navi_post_opt.ui b/app/qtapp/appcmn_qt/navi_post_opt.ui
index a69e2df57..04111e209 100644
--- a/app/qtapp/appcmn_qt/navi_post_opt.ui
+++ b/app/qtapp/appcmn_qt/navi_post_opt.ui
@@ -112,6 +112,11 @@
QZSS Broadcast
+ -
+
+ SSR VTEC
+
+
-
diff --git a/app/qtapp/rtknavi_qt/mondlg.cpp b/app/qtapp/rtknavi_qt/mondlg.cpp
index 158d53547..b97fbe866 100644
--- a/app/qtapp/rtknavi_qt/mondlg.cpp
+++ b/app/qtapp/rtknavi_qt/mondlg.cpp
@@ -421,7 +421,7 @@ void MonitorDialog::showRtk()
unsigned int nmsg[3][10] = {{0}};
char tstr[40], id[8], s1[40] = "-", s2[40] = "-", s3[40] = "-";
char file[1024] = "";
- const QString ionoopt[] = {tr("OFF"), tr("Broadcast"), tr("SBAS"), tr("Dual-Frequency"), tr("Estimate STEC"), tr("IONEX TEC"), tr("QZSS LEX")};
+ const QString ionoopt[] = {tr("OFF"), tr("Broadcast"), tr("SBAS"), tr("Dual-Frequency"), tr("Estimate STEC"), tr("IONEX TEC"), tr("QZSS LEX"), tr("SSR VTEC")};
const QString tropopt[] = {tr("OFF"), tr("Saastamoinen"), tr("SBAS"), tr("Estimate ZTD"), tr("Estimate ZTD+Grad")};
const QString ephopt [] = {tr("Broadcast"), tr("Precise"), tr("Broadcast+SBAS"), tr("Broadcast+SSR APC"), tr("Broadcast+SSR CoM"), tr("QZSS LEX")};
@@ -530,7 +530,7 @@ void MonitorDialog::showRtk()
ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1, %2").arg(rtk->opt.dynamics ? tr("ON") : tr("OFF"), rtk->opt.tidecorr ? tr("ON") : tr("OFF")));
ui->tWConsole->item(row, 0)->setText(tr("Ionosphere/Troposphere Model"));
- ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1, %2").arg(ionoopt[rtk->opt.ionoopt < 0 || rtk->opt.ionoopt > 6 ? 0 : rtk->opt.ionoopt],
+ ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1, %2").arg(ionoopt[rtk->opt.ionoopt < 0 || rtk->opt.ionoopt > 7 ? 0 : rtk->opt.ionoopt],
tropopt[rtk->opt.tropopt < 0 || rtk->opt.tropopt > 4 ? 0 : rtk->opt.tropopt]));
ui->tWConsole->item(row, 0)->setText(tr("Satellite Ephemeris"));
diff --git a/app/winapp/rtknavi/mondlg.cpp b/app/winapp/rtknavi/mondlg.cpp
index 51897b931..814f6097a 100644
--- a/app/winapp/rtknavi/mondlg.cpp
+++ b/app/winapp/rtknavi/mondlg.cpp
@@ -372,7 +372,7 @@ void __fastcall TMonitorDialog::ShowRtk(void)
int cputime,nb[3]={0},nmsg[3][10]={{0}},ne;
char tstr[40],*ant,id[8],s1[40]="-",s2[40]="-",s3[40]="-";
char file[1024]="";
- const char *ionoopt[]={"OFF","Broadcast","SBAS","Dual-Frequency","Estimate STEC","IONEX TEC","QZSS LEX",""};
+ const char *ionoopt[]={"OFF","Broadcast","SBAS","Dual-Frequency","Estimate STEC","IONEX TEC","QZSS LEX","SSR VTEC",""};
const char *tropopt[]={"OFF","Saastamoinen","SBAS","Estimate ZTD","Estimate ZTD+Grad",""};
const char *ephopt []={"Broadcast","Precise","Broadcast+SBAS","Broadcat+SSR APC","Broadcast+SSR CoM","QZSS LEX",""};
diff --git a/app/winapp/rtknavi/naviopt.dfm b/app/winapp/rtknavi/naviopt.dfm
index 0be0683e1..1a5fc2f3d 100644
--- a/app/winapp/rtknavi/naviopt.dfm
+++ b/app/winapp/rtknavi/naviopt.dfm
@@ -238,7 +238,8 @@ object OptDialog: TOptDialog
'Iono-Free LC'
'Estimate STEC'
'IONEX TEC'
- 'QZSS Broadcast')
+ 'QZSS Broadcast'
+ 'SSR VTEC')
end
object TropOpt: TComboBox
Left = 248
diff --git a/app/winapp/rtkpost/postopt.dfm b/app/winapp/rtkpost/postopt.dfm
index 14ba9f8aa..b6442234b 100644
--- a/app/winapp/rtkpost/postopt.dfm
+++ b/app/winapp/rtkpost/postopt.dfm
@@ -257,7 +257,8 @@ object OptDialog: TOptDialog
'Iono-Free LC'
'Estimate STEC'
'IONEX TEC'
- 'QZSS Broadcast')
+ 'QZSS Broadcast'
+ 'SSR VTEC')
end
object TropOpt: TComboBox
Left = 248
diff --git a/src/ionex.c b/src/ionex.c
index c5555e139..f439e0c6c 100644
--- a/src/ionex.c
+++ b/src/ionex.c
@@ -21,7 +21,7 @@
#define SQR(x) ((x)*(x))
#define MIN(x,y) ((x)<=(y)?(x):(y))
#define VAR_NOTEC SQR(30.0) /* variance of no tec */
-#define VAR_SSR_VTEC SQR(10.0); /* variance of SSR VTEC */
+#define VAR_SSR_VTEC SQR(3.0) // Variance of SSR VTEC in TEUC^2.
#define MIN_EL 0.0 /* min elevation angle (rad) */
#define MIN_HGT -1000.0 /* min user height (m) */
@@ -500,9 +500,6 @@ extern int ionvtec(gtime_t time, const nav_t *nav, const double *pos,
const double *azel, double freq, double *delay, double *var)
{
const double re=6371.0; /* fallback single-layer height (km) */
- double pppos[2],fs,vtec,P[17][17],cosl[17],sinl[17];
- double sinp,cosp,x;
- int i,j,k,nmax,mmax;
char tstr[40];
trace(4,"ionvtec: time=%s pos=%.3f %.3f azel=%.3f %.3f\n",
@@ -516,52 +513,64 @@ extern int ionvtec(gtime_t time, const nav_t *nav, const double *pos,
}
if (pos[2]vtec.nlay;i++) {
- /* ionospheric pierce point position */
- ionppp(pos,azel,re,nav->vtec.hgt[i],pppos);
+ double tow = time2gpst(time, NULL), tod = fmod(tow, 86400);
+ double vtec = 0.0, vartec = 0.0;
+ for (int i = 0; i < nav->vtec.nlay; i++) {
+ // Ionospheric pierce point position.
+ double pppos[3];
+ double fs = ionppp(pos, azel, re, nav->vtec.hgt[i], pppos);
- nmax=nav->vtec.nmax[i];
- mmax=nav->vtec.mmax[i];
- sinp=sin(pppos[0]);
- cosp=cos(pppos[0]);
+ // Mean sun fixed longitude pierce point phase shifted to the local
+ // time of maximum TEC, 14:00, and modulo 2 * PI.
+ double pplon = pppos[1] + (tod - 50400) * PI / 43200;
+ pplon -= floor(pplon / (2 * PI)) * 2 * PI;
- /* fully normalized associated Legendre polynomials */
- P[0][0]=1.0;
- for (j=1;j<=nmax;j++) {
- /* diagonal term */
- P[j][j]=sqrt((2*j+1)/(2.0*j))*cosp*P[j-1][j-1];
- /* one above diagonal */
- P[j][j-1]=sqrt(2*j+1)*sinp*P[j-1][j-1];
- }
- for (j=2;j<=nmax;j++) {
- for (k=0;k<=j-2;k++) {
- P[j][k]=sqrt((2*j+1)*(2*j-1)/((double)(j-k)*(j+k)))*sinp*P[j-1][k]
- -sqrt((2*j+1)*(j+k-1)*(j-k-1)/((double)(j-k)*(j+k)*(2*j-3)))*P[j-2][k];
+ int nmax = nav->vtec.nmax[i];
+ int mmax = nav->vtec.mmax[i];
+
+ // Fully normalized associated Legendre polynomials.
+ double sinp = sin(pppos[0]), cosp = cos(pppos[0]);
+ double P[17][17];
+ P[0][0] = 1.0;
+ P[1][0] = sqrt(3) * sinp;
+ P[1][1] = sqrt(3) * cosp;
+ // Diagonal term.
+ for (int n = 2; n <= nmax; n++) {
+ for (int m = 0; m < n; m++) {
+ double anm = sqrt((2 * n + 1) * (2 * n - 1) / ((double)(n - m) * (n + m)));
+ double bnm = sqrt((2 * n + 1) * (n + m - 1) * (n - m - 1) /
+ ((double)(n - m) * (n + m) * (2 * n - 3)));
+ P[n][m] = anm * sinp * P[n - 1][m] - bnm * P[n - 2][m];
}
+ P[n][n] = cosp * sqrt((2 * n + 1) / (2.0 * n)) * P[n - 1][n - 1];
}
- /* cosine and sine of longitude multiples */
- cosl[0]=1.0; sinl[0]=0.0;
- for (j=1;j<=mmax;j++) {
- cosl[j]=cos(j*pppos[1]);
- sinl[j]=sin(j*pppos[1]);
+
+ // Cosine and sine of longitude multiples.
+ double cosl[17], sinl[17];
+ cosl[0] = 1.0; sinl[0] = 0.0;
+ for (int m = 1; m <= mmax; m++) {
+ cosl[m] = cos(m * pplon);
+ sinl[m] = sin(m * pplon);
}
- /* evaluate spherical harmonic expansion (TECU) */
- x=0.0;
- for (j=0;j<=nmax;j++) {
- x+=nav->vtec.cosC[i][j][0]*P[j][0];
- for (k=1;k<=MIN(j,mmax);k++) {
- x+=nav->vtec.cosC[i][j][k]*P[j][k]*cosl[k]+
- nav->vtec.sinC[i][j][k]*P[j][k]*sinl[k];
+ // Evaluate spherical harmonic expansion (TECU).
+ double x = 0.0;
+ for (int n = 0; n <= nmax; n++) {
+ x += nav->vtec.cosC[i][n][0] * P[n][0];
+ for (int m = 1; m <= MIN(n, mmax); m++) {
+ x += nav->vtec.cosC[i][n][m] * P[n][m] * cosl[m] +
+ nav->vtec.sinC[i][n][m] * P[n][m] * sinl[m];
}
}
- vtec+=x;
+ // Mapping function.
+ vtec += fs * x;
+ vartec += SQR(fs * nav->vtec.qi);
+ vartec += SQR(fs) * VAR_SSR_VTEC; // Additional fixed variance for now.
}
- /* convert TECU to metres on L1, then scale to signal frequency */
- fs=ionmapf(pos,azel); /* mapping function */
- *delay=40.3E16/SQR(FREQL1)*vtec*SQR(FREQL1/freq)*fs;
- *var=VAR_SSR_VTEC; /* use fixed variance for now */
- trace(4,"ionvtec: pppos=%.3f %.3f vtec=%.1f delay=%.3f var=%.2f, freq=%.1f\n",
- pppos[0]*R2D,pppos[1]*R2D,vtec,*delay,*var,freq);
+ // Guard against negative vtec.
+ if (vtec < 0) vtec = 0;
+ // Convert TECU to metres at the signal frequency.
+ *delay = vtec * 40.3E16 / SQR(freq);
+ *var = vartec * SQR(40.3E16) / SQR(SQR(freq));
+ trace(4,"ionvtec: vtec=%.1f delay=%.3f var=%.2f, freq=%.1f\n", vtec, *delay, *var, freq);
return 1;
}
diff --git a/src/options.c b/src/options.c
index 2204e07f3..c30a622ba 100644
--- a/src/options.c
+++ b/src/options.c
@@ -45,7 +45,7 @@ static char snrmask_[NFREQ][1024];
#define MODOPT "0:single,1:dgps,2:kinematic,3:static,4:static-start,5:movingbase,6:fixed,7:ppp-kine,8:ppp-static,9:ppp-fixed"
#define FRQOPT "1:l1,2:l1+l2,3:l1+l2+l5,4:l1+l2+l5+l6"
#define TYPOPT "0:forward,1:backward,2:combined,3:combined-nophasereset"
-#define IONOPT "0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc"
+#define IONOPT "0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc,7:ssr-vtec"
#define TRPOPT "0:off,1:saas,2:sbas,3:est-ztd,4:est-ztdgrad"
#define EPHOPT "0:brdc,1:precise,2:brdc+sbas,3:brdc+ssrapc,4:brdc+ssrcom"
#define NAVOPT "1:gps+2:sbas+4:glo+8:gal+16:qzs+32:bds+64:navic"
diff --git a/src/pntpos.c b/src/pntpos.c
index 5f6095cb8..c5ded7b61 100644
--- a/src/pntpos.c
+++ b/src/pntpos.c
@@ -221,6 +221,11 @@ extern int ionocorr(gtime_t time, const nav_t *nav, int sat, const double *pos,
if (iontec(time,nav,pos,azel,1,ion,var)) return 1;
err=1;
}
+ // IONEX VTEC model.
+ if (ionoopt == IONOOPT_VTEC) {
+ if (ionvtec(time, nav, pos, azel, FREQL1, ion, var)) return 1;
+ err=1;
+ }
/* QZSS broadcast ionosphere model */
if (ionoopt==IONOOPT_QZS&&norm(nav->ion_qzs,8)>0.0) {
*ion=ionmodel(time,nav->ion_qzs,pos,azel);
diff --git a/src/ppp.c b/src/ppp.c
index a4a056ec3..dd822d9a9 100644
--- a/src/ppp.c
+++ b/src/ppp.c
@@ -713,7 +713,7 @@ static void udiono_ppp(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
int sys=satsys(sat,NULL);
double P0_corr=obs[i].P[0];
double Pf_corr=obs[i].P[f2];
- if (&rtk->opt.sateph==EPHOPT_SSRAPC||&rtk->opt.sateph==EPHOPT_SSRCOM) {
+ if (rtk->opt.sateph == EPHOPT_SSRAPC || rtk->opt.sateph == EPHOPT_SSRCOM) {
/* apply SSR correction */
P0_corr-=nav->ssr[obs->sat-1].cbias[obs[i].code[0]-1];
Pf_corr-=nav->ssr[obs->sat-1].cbias[obs[i].code[f2]-1];
diff --git a/src/rtcm3.c b/src/rtcm3.c
index 8de486eef..eb145bdfa 100644
--- a/src/rtcm3.c
+++ b/src/rtcm3.c
@@ -2059,8 +2059,8 @@ static int decode_ssr7(rtcm_t *rtcm, int sys, int subtype)
/* decode SSR 8: VTEC ionosphere --------------------------------------------- */
static int decode_ssr8(rtcm_t *rtcm, int subtype)
{
- double udint,hgt,cosC,sinC,tow;
- int i=0,j,k,l,type,sync,iod,udi,nlay,nmax,mmax,np,offp,qi;
+ double udint,hgt,cosC,sinC;
+ int i=0,type,sync,iod,udi,nlay,nmax,mmax,qi;
char *msg,tstr[40];
type=getbitu(rtcm->buff,24,12);
@@ -2086,7 +2086,7 @@ static int decode_ssr8(rtcm_t *rtcm, int subtype)
return -1;
}
- for (j=0;jlen*8;j++) {
+ for (int j=0;jlen*8;j++) {
hgt =getbitu(rtcm->buff,i,8)*10; i+=8; /* height (km) DF472 */
nmax =getbitu(rtcm->buff,i, 4)+1; i+= 4; /* degree DF473 */
mmax =getbitu(rtcm->buff,i, 4)+1; i+= 4; /* order DF474 */
@@ -2099,17 +2099,17 @@ static int decode_ssr8(rtcm_t *rtcm, int subtype)
memset(rtcm->nav.vtec.cosC, 0, sizeof(rtcm->nav.vtec.cosC));
memset(rtcm->nav.vtec.sinC, 0, sizeof(rtcm->nav.vtec.sinC));
/* cosine coefficients: for order o=0..mmax, degree n=o..nmax */
- for (k=0;k<=mmax;k++) {
- for (l=k;l<=nmax;l++) {
+ for (int m=0;m<=mmax;m++) {
+ for (int n=m;n<=nmax&&i+16<=rtcm->len*8;n++) {
cosC=getbits(rtcm->buff,i,16); i+=16;
- rtcm->nav.vtec.cosC[j][l][k]=cosC/200.0;
+ rtcm->nav.vtec.cosC[j][n][m]=cosC/200.0;
}
}
/* sine coefficients: for order o=1..mmax, degree n=o..nmax */
- for (k=1;k<=mmax;k++) {
- for (l=k;l<=nmax;l++) {
+ for (int m=1;m<=mmax;m++) {
+ for (int n=m;n<=nmax&&i+16<=rtcm->len*8;n++) {
sinC=getbits(rtcm->buff,i,16); i+=16;
- rtcm->nav.vtec.sinC[j][l][k]=sinC/200.0;
+ rtcm->nav.vtec.sinC[j][n][m]=sinC/200.0;
}
}
}
@@ -2717,6 +2717,7 @@ static int decode_type4076(rtcm_t *rtcm)
case 125: return decode_ssr3(rtcm,SYS_SBS,subtype);
case 126: return decode_ssr7(rtcm,SYS_SBS,subtype);
case 127: return decode_ssr5(rtcm,SYS_SBS,subtype);
+ case 201: return decode_ssr8(rtcm, subtype);
}
trace(2,"rtcm3 4076: unsupported message subtype=%d\n",subtype);
return 0;
diff --git a/src/rtklib.h b/src/rtklib.h
index b11da21c4..743b46db2 100644
--- a/src/rtklib.h
+++ b/src/rtklib.h
@@ -431,6 +431,7 @@ extern "C" {
#define IONOOPT_EST 4 /* ionosphere option: estimation */
#define IONOOPT_TEC 5 /* ionosphere option: IONEX TEC model */
#define IONOOPT_QZS 6 /* ionosphere option: QZSS broadcast model */
+#define IONOOPT_VTEC 7 /* ionosphere option: IONEX VTEC model */
#define TROPOPT_OFF 0 /* troposphere option: correction off */
#define TROPOPT_SAAS 1 /* troposphere option: Saastamoinen model */