Ssr vtec misc

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);

app/qtapp/appcmn_qt/navi_post_opt.ui

@@ -112,6 +112,11 @@
            <string>QZSS Broadcast</string>
           </property>
          </item>
+         <item>
+          <property name="text">
+           <string>SSR VTEC</string>
+          </property>
+         </item>
         </widget>
        </item>
        <item row="4" column="2">

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")};
 
@@ -536,7 +536,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"));

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",""};
 

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

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

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) */
 
@@ -499,10 +499,7 @@ extern int iontec(gtime_t time, const nav_t *nav, const double *pos,
 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;
+    const double re = 6370.0; // Spherical earths radius (km).
     char tstr[40];
 
     trace(4,"ionvtec: time=%s pos=%.3f %.3f azel=%.3f %.3f\n",
@@ -516,52 +513,70 @@ extern int ionvtec(gtime_t time, const nav_t *nav, const double *pos,
     }
     if (pos[2]<MIN_HGT||azel[1]<=0.0) return 0;
 
-    vtec=0.0;
-    for (i=0;i<nav->vtec.nlay;i++) {
-        /* ionospheric pierce point position */
-        ionppp(pos,azel,re,nav->vtec.hgt[i],pppos);
+    // Convert the station position to the spherical earth model.
+    double r[3]; // ecef
+    pos2ecef(pos, r);
+    double spos[3], r2 = SQR(r[0]) + SQR(r[1]);
+    spos[0] = r2 > 1E-12 ? atan(r[2] / sqrt(r2)) : (r[2] > 0.0 ? PI / 2.0 : -PI / 2.0);
+    spos[1] = r2 > 1E-12 ? atan2(r[1], r[0]) : 0.0;
+    spos[2] = sqrt(r2 + SQR(r[2])) - re * 1000;
 
-        nmax=nav->vtec.nmax[i];
-        mmax=nav->vtec.mmax[i];
-        sinp=sin(pppos[0]);
-        cosp=cos(pppos[0]);
+    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(spos, azel, re, nav->vtec.hgt[i], pppos);
 
-       /* 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];
+        // Mean sun fixed longitude pierce point phase shifted to the local
+        // time of maximum TEC, 14:00, and modulo 2 * PI.
+        double pplon = fmod(pppos[1] + (tod - 50400) * PI / 43200, 2 * PI);
+        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;
+        for (int n = 2; n <= nmax; n++) {
+            for (int m = 0; m < n - 1; 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 - 1] = sqrt(2 * n + 1) * sinp * P[n - 1][n - 1];
+            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] * cosl[m] +
+                      nav->vtec.sinC[i][n][m] * sinl[m]) * P[n][m];
             }
         }
-        vtec+=x;
+        // Guard against a negative vtec, per layer.
+        if (x < 0) x = 0;
+        // 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);
+    // 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;
 }

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"

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);

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;j<nlay&&i+16<=rtcm->len*8;j++) {
+    for (int j=0;j<nlay&&i+16<=rtcm->len*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;

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 */