Corrections in processing of Galileo Nav Messages

app/convbin/convbin.c

@@ -37,6 +37,8 @@ static const char rcsid[]="$Id: convbin.c,v 1.1 2008/07/17 22:13:04 ttaka Exp $"
 #define PRGNAME   "CONVBIN"
 #define TRACEFILE "convbin.trace"
 
+int galmessagetype;
+
 /* help text -----------------------------------------------------------------*/
 static const char *help[]={
 "",
@@ -165,6 +167,8 @@ static int convbin(int format, rnxopt_t *opt, const char *ifile, char **file,
     char *extnav=opt->rnxver<=2.99||opt->navsys==SYS_GPS?"N":"P";
     char *extlog=format==STRFMT_LEXR?"lex":"sbs";
 
+    galmessagetype = GALMESS_INAV;
+
     def=!file[0]&&!file[1]&&!file[2]&&!file[3]&&!file[4]&&!file[5]&&!file[6];
 
     for (i=0;i<7;i++) ofile[i]=ofile_[i];

app/rnx2rtkp/rnx2rtkp.c

@@ -22,6 +22,8 @@ static const char rcsid[]="$Id: rnx2rtkp.c,v 1.1 2008/07/17 21:55:16 ttaka Exp $
 #define PROGNAME    "rnx2rtkp"          /* program name */
 #define MAXFILE     8                   /* max number of input files */
 
+int galmessagetype;
+
 /* help text -----------------------------------------------------------------*/
 static const char *help[]={
 "",
@@ -97,6 +99,8 @@ int main(int argc, char **argv)
     int i,j,n,ret;
     char *infile[MAXFILE],*outfile="";
 
+    galmessagetype = GALMESS_INAV;
+
     prcopt.mode  =PMODE_KINEMA;
     prcopt.navsys=SYS_GPS|SYS_GLO;
     prcopt.refpos=1;
@@ -162,6 +166,8 @@ int main(int argc, char **argv)
         showmsg("error : no input file");
         return -2;
     }
+
+    if (prcopt.ionoopt==IONOOPT_IFLC) galmessagetype=GALMESS_FNAV;
     ret=postpos(ts,te,tint,0.0,&prcopt,&solopt,&filopt,infile,n,outfile,"","");
 
     if (!ret) fprintf(stderr,"%40s\r","");

app/rtkrcv/rtkrcv.c

@@ -138,6 +138,8 @@ static const char *pathopts[]={         /* path options help */
 #define SOLOPT  "0:llh,1:xyz,2:enu,3:nmea"
 #define MSGOPT  "0:all,1:rover,2:base,3:corr"
 
+int galmessagetype;
+
 static opt_t rcvopts[]={
     {"console-passwd",  2,  (void *)passwd,              ""     },
     {"console-timetype",3,  (void *)&timetype,           TIMOPT },
@@ -1369,6 +1371,8 @@ int main(int argc, char **argv)
     int i,start=0,port=0,outstat=0,trace=0;
     char *dev="",file[MAXSTR]="";
 
+    galmessagetype = GALMESS_INAV;
+
     for (i=1;i<argc;i++) {
         if      (!strcmp(argv[i],"-s")) start=1;
         else if (!strcmp(argv[i],"-p")&&i+1<argc) port=atoi(argv[++i]);

src/ephemeris.c

@@ -96,6 +96,15 @@ static double var_uraeph(int ura)
     };
     return ura<0||15<ura?SQR(6144.0):SQR(ura_value[ura]);
 }
+/* variance by sisa ephemeris (ref [7] Issue 1.3 5.1.11) --------------------------*/
+static double var_sisaeph(int ura)
+{
+    if (ura<50) return SQR(ura*0.01);
+    if (ura<75) return SQR(0.50+(ura-50)*0.02);
+    if (ura<100) return SQR(1.0+(ura-75)*0.04);
+    if (ura<=125) return SQR(2.0+(ura-100)*0.16);
+    return -1.0;
+}
 /* variance by ura ssr (ref [4]) ---------------------------------------------*/
 static double var_urassr(int ura)
 {
@@ -246,7 +255,11 @@ extern void eph2pos(gtime_t time, const eph_t *eph, double *rs, double *dts,
     *dts-=2.0*sqrt(mu*eph->A)*eph->e*sinE/SQR(CLIGHT);
 
     /* position and clock error variance */
-    *var=var_uraeph(eph->sva);
+    if (sys==SYS_GAL){
+    	*var=var_sisaeph(eph->sva);
+    } else {
+    	*var=var_uraeph(eph->sva);
+    }
 }
 /* glonass orbit differential equations --------------------------------------*/
 static void deq(const double *x, double *xdot, const double *acc)
@@ -382,25 +395,38 @@ extern void seph2pos(gtime_t time, const seph_t *seph, double *rs, double *dts,
     *var=var_uraeph(seph->sva);
 }
 /* select ephememeris --------------------------------------------------------*/
-static eph_t *seleph(gtime_t time, int sat, int iode, const nav_t *nav)
+extern eph_t *seleph(gtime_t time, int sat, int iode, const nav_t *nav)
 {
     double t,tmax,tmin;
-    int i,j=-1;
-    
+    int i,j=-1,sys;
+
     trace(4,"seleph  : time=%s sat=%2d iode=%d\n",time_str(time,3),sat,iode);
-
-    switch (satsys(sat,NULL)) {
+
+    sys=satsys(sat,NULL);
+    switch (sys) {
         case SYS_QZS: tmax=MAXDTOE_QZS+1.0; break;
-        case SYS_GAL: tmax=MAXDTOE_GAL+1.0; break;
+        case SYS_GAL: tmax=MAXDTOE_GAL; break;
         case SYS_CMP: tmax=MAXDTOE_CMP+1.0; break;
         default: tmax=MAXDTOE+1.0; break;
     }
     tmin=tmax+1.0;
-    
+
     for (i=0;i<nav->n;i++) {
         if (nav->eph[i].sat!=sat) continue;
         if (iode>=0&&nav->eph[i].iode!=iode) continue;
-        if ((t=fabs(timediff(nav->eph[i].toe,time)))>tmax) continue;
+        if (sys==SYS_GAL){
+        	if (galmessagetype==GALMESS_FNAV){
+        		/* Check if F/NAV ephemeris */
+        		if (!(nav->eph[i].code & 0x02)) continue;
+        	} else {
+        		/* Check if I/NAV ephemeris */
+        		if (!(nav->eph[i].code & 0x05)) continue;
+        	}
+        	if ((t=timediff(time,nav->eph[i].toe))>tmax) continue;
+        	if (t<0) continue;
+        } else {
+        	if ((t=fabs(timediff(nav->eph[i].toe,time)))>tmax) continue;
+        }
         if (iode>=0) return nav->eph+i;
         if (t<=tmin) {j=i; tmin=t;} /* toe closest to time */
     }
@@ -499,7 +525,6 @@ static int ephpos(gtime_t time, gtime_t teph, int sat, const nav_t *nav,
 
     if (sys==SYS_GPS||sys==SYS_GAL||sys==SYS_QZS||sys==SYS_CMP) {
         if (!(eph=seleph(teph,sat,iode,nav))) return 0;
-
         eph2pos(time,eph,rs,dts,var);
         time=timeadd(time,tt);
         eph2pos(time,eph,rst,dtst,var);

src/pntpos.c

@@ -203,7 +203,9 @@ static int rescode(int iter, const obsd_t *obs, int n, const double *rs,
                    double *resp, int *ns)
 {
     double r,dion,dtrp,vmeas,vion,vtrp,rr[3],pos[3],dtr,e[3],P,lam_L1;
-    int i,j,nv=0,sys,mask[4]={0};
+    int i,j,nv=0,sys,mask[4]={0},sva=-1;
+
+    eph_t *eph;
 
     trace(3,"resprng : n=%d\n",n);
 
@@ -230,8 +232,13 @@ static int rescode(int iter, const obsd_t *obs, int n, const double *rs,
         /* psudorange with code bias correction */
         if ((P=prange(obs+i,nav,azel+i*2,iter,opt,&vmeas))==0.0) continue;
 
+        /* if Gal satellite, check SISA is not NAPA*/
+        if (sys==SYS_GAL){
+        	if (!(eph=seleph(obs[i].time,obs[i].sat,-1,nav))) continue;
+        	sva = eph->sva;
+        }
         /* excluded satellite? */
-        if (satexclude(obs[i].sat,svh[i],opt)) continue;
+        if (satexclude(obs[i].sat,svh[i],sva,opt)) continue;
 
         /* ionospheric corrections */
         if (!ionocorr(obs[i].time,nav,obs[i].sat,pos,azel+i*2,

src/ppp.c

@@ -884,7 +884,9 @@ static int res_ppp(int iter, const obsd_t *obs, int n, const double *rs,
     prcopt_t *opt=&rtk->opt;
     double r,rr[3],disp[3],pos[3],e[3],meas[2],dtdx[3],dantr[NFREQ]={0};
     double dants[NFREQ]={0},var[MAXOBS*2],dtrp=0.0,vart=0.0,varm[2]={0};
-    int i,j,k,sat,sys,nv=0,nx=rtk->nx,brk,tideopt;
+    int i,j,k,sat,sys,nv=0,nx=rtk->nx,brk,tideopt,sva=-1;
+
+    eph_t *eph;
 
     trace(3,"res_ppp : n=%d nx=%d\n",n,nx);
 
@@ -910,8 +912,14 @@ static int res_ppp(int iter, const obsd_t *obs, int n, const double *rs,
         if ((r=geodist(rs+i*6,rr,e))<=0.0||
             satazel(pos,e,azel+i*2)<opt->elmin) continue;
 
+        /* if Gal satellite, check SISA is not NAPA*/
+        if (sys==SYS_GAL){
+        	if (!(eph=seleph(obs[i].time,obs[i].sat,-1,nav))) continue;
+        	sva = eph->sva;
+        }
+
         /* excluded satellite? */
-        if (satexclude(obs[i].sat,svh[i],opt)) continue;
+        if (satexclude(obs[i].sat,svh[i],sva,opt)) continue;
 
         /* tropospheric delay correction */
         if (opt->tropopt==TROPOPT_SAAS) {

src/rinex.c

@@ -184,6 +184,36 @@ static int uraindex(double value)
     for (i=0;i<15;i++) if (ura_eph[i]>=value) break;
     return i;
 }
+/* SISA index to SISA value (m) */
+static double sisavalue(int sva)
+{
+	if (sva<50) return sva*0.01;
+	if (sva<75) return 0.50+(sva-50)*0.02;
+	if (sva<100) return 1.0+(sva-75)*0.04;
+	if (sva<=125) return 2.0+(sva-100)*0.16;
+	return -1.0;
+}
+/* SISA value (m) to SISA index */
+static int sisaindex(double value)
+{
+    /* 0.0 m is an impossible value due to the statistical meaning of the SISA
+     * and is identified as NAPA due to some receivers setting the SISA field to
+     * 0.0 m when the SISA is set to NAPA (NAPA is undefined in the RINEX format
+     * definition). */
+    if (value<0.01) return 255;
+    /* index     SISA value
+     * 0..49      0m to 0.49m with 1cm resolution
+     * 50..74     0.50cm to 0.98m with 2cm resolution
+     * 75..99     1.0m to 1.96m with 4cm resolution
+     * 100..125   2.0m to 6.0m with 16cm resolution
+     * 126..254   spare
+     * 255        NAPA */
+    if (value<0.5) return (int)(value/0.01);
+    if (value<1.0) return (int)(50+(value-0.5)/0.02);
+    if (value<2.0) return (int)(75+(value-1.0)/0.04);
+    if (value<=6.0) return (int)(100+(value-2.0)/0.16);
+    return 255;
+}
 /* initialize station parameter ----------------------------------------------*/
 static void init_sta(sta_t *sta)
 {
@@ -1066,7 +1096,7 @@ static int decode_eph(double ver, int sat, gtime_t toc, const double *data,
                                       /* bit   4-5: E5a HS */
                                       /* bit     6: E5b DVS */
                                       /* bit   7-8: E5b HS */
-        eph->sva =uraindex(data[23]); /* ura (m->index) */
+        eph->sva =sisaindex(data[23]); /* ura (m->index) */
 
         eph->tgd[0]=   data[25];      /* BGD E5a/E1 */
         eph->tgd[1]=   data[26];      /* BGD E5b/E1 */
@@ -2270,7 +2300,11 @@ extern int outrnxnavb(FILE *fp, const rnxopt_t *opt, const eph_t *eph)
     outnavf(fp,eph->flag   );
     fprintf(fp,"\n%s",sep  );
 
-    outnavf(fp,uravalue(eph->sva));
+    if (sys==SYS_GAL){
+    	outnavf(fp,sisavalue(eph->sva));
+    } else {
+    	outnavf(fp,uravalue(eph->sva));
+    }
     outnavf(fp,eph->svh    );
     outnavf(fp,eph->tgd[0] ); /* GPS/QZS:TGD, GAL:BGD E5a/E1, BDS: TGD1 B1/B3 */
     if (sys==SYS_GAL||sys==SYS_CMP) {

src/rtkcmn.c

@@ -468,10 +468,11 @@ extern void satno2id(int sat, char *id)
 * test excluded satellite
 * args   : int    sat       I   satellite number
 *          int    svh       I   sv health flag
+*          int    sva       I   SV URA/SISA
 *          prcopt_t *opt    I   processing options (NULL: not used)
 * return : status (1:excluded,0:not excluded)
 *-----------------------------------------------------------------------------*/
-extern int satexclude(int sat, int svh, const prcopt_t *opt)
+extern int satexclude(int sat, int svh, const int sva, const prcopt_t *opt)
 {
     int sys=satsys(sat,NULL);
 
@@ -483,6 +484,11 @@ extern int satexclude(int sat, int svh, const prcopt_t *opt)
         if (!(sys&opt->navsys)) return 1; /* unselected sat sys */
     }
     if (sys==SYS_QZS) svh&=0xFE; /* mask QZSS LEX health */
+    /* Exclude Galileo satellites with NAPA or invalid SISA */
+    if ((sys==SYS_GAL) && (sva>125)){
+    	trace(3,"Galileo satellite with NAPA or invalid SISA: sat=%3d",sat);
+    	return 1;
+    }
     if (svh) {
         trace(3,"unhealthy satellite: sat=%3d svh=%02X\n",sat,svh);
         return 1;
@@ -2393,7 +2399,7 @@ static int cmpeph(const void *p1, const void *p2)
 static void uniqeph(nav_t *nav)
 {
     eph_t *nav_eph;
-    int i,j;
+    int i,j,sys;
 
     trace(3,"uniqeph: n=%d\n",nav->n);
 
@@ -2402,8 +2408,11 @@ static void uniqeph(nav_t *nav)
     qsort(nav->eph,nav->n,sizeof(eph_t),cmpeph);
 
     for (i=1,j=0;i<nav->n;i++) {
+    	sys = satsys(nav->eph[i].sat,NULL);
         if (nav->eph[i].sat!=nav->eph[j].sat||
-            nav->eph[i].iode!=nav->eph[j].iode) {
+            nav->eph[i].iode!=nav->eph[j].iode||
+            ((sys==SYS_GAL)&&(nav->eph[i].code!=nav->eph[j].code))||
+            (nav->eph[i].toe.sec!=nav->eph[j].toe.sec)) {
             nav->eph[++j]=nav->eph[i];
         }
     }

src/rtklib.h

@@ -137,7 +137,7 @@ extern "C" {
 #endif
 #ifdef ENAGAL
 #define MINPRNGAL   1                   /* min satellite PRN number of Galileo */
-#define MAXPRNGAL   30                  /* max satellite PRN number of Galileo */
+#define MAXPRNGAL   36                  /* max satellite PRN number of Galileo */
 #define NSATGAL    (MAXPRNGAL-MINPRNGAL+1) /* number of Galileo satellites */
 #define NSYSGAL     1
 #else
@@ -199,7 +199,7 @@ extern "C" {
 #define DTTOL       0.005               /* tolerance of time difference (s) */
 #define MAXDTOE     7200.0              /* max time difference to GPS Toe (s) */
 #define MAXDTOE_QZS 7200.0              /* max time difference to QZSS Toe (s) */
-#define MAXDTOE_GAL 10800.0             /* max time difference to Galileo Toe (s) */
+#define MAXDTOE_GAL 14400.0             /* max time difference to Galileo Toe (s) */
 #define MAXDTOE_CMP 21600.0             /* max time difference to BeiDou Toe (s) */
 #define MAXDTOE_GLO 1800.0              /* max time difference to GLONASS Toe (s) */
 #define MAXDTOE_SBS 360.0               /* max time difference to SBAS Toe (s) */
@@ -338,6 +338,9 @@ extern "C" {
 #define IONOOPT_LEX 7                   /* ionosphere option: QZSS LEX ionospehre */
 #define IONOOPT_STEC 8                  /* ionosphere option: SLANT TEC model */
 
+#define GALMESS_INAV 0                  /* Galileo message type: FNAV */
+#define GALMESS_FNAV 1                  /* Galileo message type: INAV */
+
 #define TROPOPT_OFF 0                   /* troposphere option: correction off */
 #define TROPOPT_SAAS 1                  /* troposphere option: Saastamoinen model */
 #define TROPOPT_SBAS 2                  /* troposphere option: SBAS model */
@@ -1260,6 +1263,7 @@ extern const sbsigpband_t igpband1[][8]; /* SBAS IGP band 0-8 */
 extern const sbsigpband_t igpband2[][5]; /* SBAS IGP band 9-10 */
 extern const char *formatstrs[];        /* stream format strings */
 extern opt_t sysopts[];                 /* system options table */
+extern int galmessagetype;              /* GAL message type */
 
 /* satellites, systems, codes functions --------------------------------------*/
 extern int  satno   (int sys, int prn);
@@ -1268,7 +1272,7 @@ extern int  satid2no(const char *id);
 extern void satno2id(int sat, char *id);
 extern unsigned char obs2code(const char *obs, int *freq);
 extern char *code2obs(unsigned char code, int *freq);
-extern int  satexclude(int sat, int svh, const prcopt_t *opt);
+extern int  satexclude(int sat, int svh, const int sva, const prcopt_t *opt);
 extern int  testsnr(int base, int freq, double el, double snr,
                     const snrmask_t *mask);
 extern void setcodepri(int sys, int freq, const char *pri);
@@ -1491,6 +1495,8 @@ extern int tle_name_read(const char *file, tle_t *tle);
 extern int tle_pos(gtime_t time, const char *name, const char *satno,
                    const char *desig, const tle_t *tle, const erp_t *erp,
                    double *rs);
+extern eph_t *seleph(const gtime_t time, const int sat, const int iode,
+                     const nav_t *nav);
 
 /* receiver raw data functions -----------------------------------------------*/
 extern unsigned int getbitu(const unsigned char *buff, int pos, int len);

src/rtkpos.c

@@ -851,7 +851,9 @@ static int zdres(int base, const obsd_t *obs, int n, const double *rs,
 {
     double r,rr_[3],pos[3],dant[NFREQ]={0},disp[3];
     double zhd,zazel[]={0.0,90.0*D2R};
-    int i,nf=NF(opt);
+    int i,nf=NF(opt),sys,sva=-1;
+
+    eph_t *eph;
 
     trace(3,"zdres   : n=%d\n",n);
 
@@ -874,8 +876,14 @@ static int zdres(int base, const obsd_t *obs, int n, const double *rs,
         if ((r=geodist(rs+i*6,rr_,e+i*3))<=0.0) continue;
         if (satazel(pos,e+i*3,azel+i*2)<opt->elmin) continue;
 
+        /* if Gal satellite, check SISA is not NAPA*/
+        if (!(sys=satsys(obs[i].sat,NULL))) continue;
+        if (sys==SYS_GAL){
+        	if (!(eph=seleph(obs[i].time,obs[i].sat,-1,nav))) continue;
+        	sva = eph->sva;
+        }
         /* excluded satellite? */
-        if (satexclude(obs[i].sat,svh[i],opt)) continue;
+        if (satexclude(obs[i].sat,svh[i],sva,opt)) continue;
 
         /* satellite clock-bias */
         r+=-CLIGHT*dts[i*2];