From 70e8c3737b49900cd32cc3b65f43ed34996c7874 Mon Sep 17 00:00:00 2001 From: Miquel Garcia Date: Thu, 31 Jul 2025 09:24:55 +0200 Subject: [PATCH 1/6] chore(rtcm): remove trailing spaces --- src/rtcm3.c | 294 +++++++++++++++--------------- src/rtcm3e.c | 400 ++++++++++++++++++++--------------------- util/testeph/dumpssr.c | 30 ++-- 3 files changed, 362 insertions(+), 362 deletions(-) diff --git a/src/rtcm3.c b/src/rtcm3.c index dd71d1eec..62dddc576 100644 --- a/src/rtcm3.c +++ b/src/rtcm3.c @@ -166,7 +166,7 @@ static void adjweek(rtcm_t *rtcm, double tow) { double tow_p; int week; - + /* if no time, get cpu time */ if (rtcm->time.time==0) rtcm->time=utc2gpst(timeget()); tow_p=time2gpst(rtcm->time,&week); @@ -188,7 +188,7 @@ static void adjday_glot(rtcm_t *rtcm, double tod) gtime_t time; double tow,tod_p; int week; - + if (rtcm->time.time==0) rtcm->time=utc2gpst(timeget()); time=timeadd(gpst2utc(rtcm->time),10800.0); /* glonass time */ tow=time2gpst(time,&week); @@ -223,12 +223,12 @@ static double snratio(double snr) static int obsindex(obs_t *obs, gtime_t time, int sat) { int i,j; - + for (i=0;in;i++) { if (obs->data[i].sat==sat) return i; /* field already exists */ } if (i>=MAXOBS) return -1; /* overflow */ - + /* add new field */ obs->data[i].time=time; obs->data[i].sat=sat; @@ -245,7 +245,7 @@ static int test_staid(rtcm_t *rtcm, int staid) { char *p; int type,id; - + /* test station id option */ if ((p=strstr(rtcm->opt,"-STA="))&&sscanf(p,"-STA=%d",&id)==1) { if (staid!=id) return 0; @@ -257,7 +257,7 @@ static int test_staid(rtcm_t *rtcm, int staid) else if (staid!=rtcm->staid) { type=getbitu(rtcm->buff,24,12); trace(2,"rtcm3 %d staid invalid id=%d %d\n",type,staid,rtcm->staid); - + /* reset station id if station id error */ rtcm->staid=0; return 0; @@ -270,11 +270,11 @@ static int decode_head1001(rtcm_t *rtcm, int *sync) double tow; char *msg,tstr[40]; int i=24,staid,nsat,type; - + if (rtcm->obsflag) rtcm->obs.n = rtcm->obsflag = 0; type=getbitu(rtcm->buff,i,12); i+=12; - + if (i+52<=rtcm->len*8) { staid=getbitu(rtcm->buff,i,12); i+=12; tow =getbitu(rtcm->buff,i,30)*0.001; i+=30; @@ -287,12 +287,12 @@ static int decode_head1001(rtcm_t *rtcm, int *sync) } /* test station ID */ if (!test_staid(rtcm,staid)) return -1; - + adjweek(rtcm,tow); - + time2str(rtcm->time,tstr,2); trace(4,"decode_head1001: time=%s nsat=%d sync=%d\n",tstr,nsat,*sync); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," staid=%4d %s nsat=%2d sync=%d",staid,tstr,nsat,*sync); @@ -312,9 +312,9 @@ static int decode_type1002(rtcm_t *rtcm) { double pr1,cnr1,tt,cp1,freq=FREQL1; int i=24+64,j,index,nsat,sync,prn,code,sat,ppr1,lock1,amb,sys; - + if ((nsat=decode_head1001(rtcm,&sync))<0) return -1; - + for (j=0;jlen*8;j++) { prn =getbitu(rtcm->buff,i, 6); i+= 6; code =getbitu(rtcm->buff,i, 1); i+= 1; @@ -338,7 +338,7 @@ static int decode_type1002(rtcm_t *rtcm) if ((index=obsindex(&rtcm->obs,rtcm->time,sat))<0) continue; pr1=pr1*0.02+amb*PRUNIT_GPS; rtcm->obs.data[index].P[0]=pr1; - + if (ppr1!=(int)0xFFF80000) { cp1=adjcp(rtcm,sat,0,ppr1*0.0005*freq/CLIGHT); rtcm->obs.data[index].L[0]=pr1*freq/CLIGHT+cp1; @@ -364,9 +364,9 @@ static int decode_type1004(rtcm_t *rtcm) double pr1,cnr1,cnr2,tt,cp1,cp2,freq[2]={FREQL1,FREQL2}; int i=24+64,j,index,nsat,sync,prn,sat,code1,code2,pr21,ppr1,ppr2; int lock1,lock2,amb,sys; - + if ((nsat=decode_head1001(rtcm,&sync))<0) return -1; - + for (j=0;jlen*8;j++) { prn =getbitu(rtcm->buff,i, 6); i+= 6; code1=getbitu(rtcm->buff,i, 1); i+= 1; @@ -395,7 +395,7 @@ static int decode_type1004(rtcm_t *rtcm) if ((index=obsindex(&rtcm->obs,rtcm->time,sat))<0) continue; pr1=pr1*0.02+amb*PRUNIT_GPS; rtcm->obs.data[index].P[0]=pr1; - + if (ppr1!=(int)0xFFF80000) { cp1=adjcp(rtcm,sat,0,ppr1*0.0005*freq[0]/CLIGHT); rtcm->obs.data[index].L[0]=pr1*freq[0]/CLIGHT+cp1; @@ -403,7 +403,7 @@ static int decode_type1004(rtcm_t *rtcm) rtcm->obs.data[index].LLI[0]=lossoflock(rtcm,sat,0,lock1); rtcm->obs.data[index].SNR[0]=snratio(cnr1*0.25); rtcm->obs.data[index].code[0]=code1?CODE_L1P:CODE_L1C; - + if (pr21!=(int)0xFFFFE000) { rtcm->obs.data[index].P[1]=pr1+pr21*0.02; } @@ -429,7 +429,7 @@ static int decode_type1005(rtcm_t *rtcm) double rr[3],re[3],pos[3]; char *msg; int i=24+12,j,staid,itrf; - + if (i+140==rtcm->len*8) { staid=getbitu(rtcm->buff,i,12); i+=12; itrf =getbitu(rtcm->buff,i, 6); i+= 6+4; @@ -450,7 +450,7 @@ static int decode_type1005(rtcm_t *rtcm) } /* test station id */ if (!test_staid(rtcm,staid)) return -1; - + sprintf(rtcm->sta.name,"%04d",staid); rtcm->sta.deltype=0; /* xyz */ for (j=0;j<3;j++) { @@ -467,7 +467,7 @@ static int decode_type1006(rtcm_t *rtcm) double rr[3],re[3],pos[3],anth; char *msg; int i=24+12,j,staid,itrf; - + if (i+156<=rtcm->len*8) { staid=getbitu(rtcm->buff,i,12); i+=12; itrf =getbitu(rtcm->buff,i, 6); i+= 6+4; @@ -489,7 +489,7 @@ static int decode_type1006(rtcm_t *rtcm) } /* test station id */ if (!test_staid(rtcm,staid)) return -1; - + sprintf(rtcm->sta.name,"%04d",staid); rtcm->sta.deltype=1; /* xyz */ for (j=0;j<3;j++) { @@ -506,9 +506,9 @@ static int decode_type1007(rtcm_t *rtcm) char des[32]=""; char *msg; int i=24+12,j,staid,n,setup; - + n=getbitu(rtcm->buff,i+12,8); - + if (i+28+8*n<=rtcm->len*8) { staid=getbitu(rtcm->buff,i,12); i+=12+8; for (j=0;jsta.name,"%04d",staid); strncpy(rtcm->sta.antdes,des,n); rtcm->sta.antdes[n]='\0'; rtcm->sta.antsetup=setup; @@ -539,10 +539,10 @@ static int decode_type1008(rtcm_t *rtcm) char des[32]="",sno[32]=""; char *msg; int i=24+12,j,staid,n,m,setup; - + n=getbitu(rtcm->buff,i+12,8); m=getbitu(rtcm->buff,i+28+8*n,8); - + if (i+36+8*(n+m)<=rtcm->len*8) { staid=getbitu(rtcm->buff,i,12); i+=12+8; for (j=0;jsta.name,"%04d",staid); strncpy(rtcm->sta.antdes,des,n); rtcm->sta.antdes[n]='\0'; rtcm->sta.antsetup=setup; @@ -576,11 +576,11 @@ static int decode_head1009(rtcm_t *rtcm, int *sync) double tod; char *msg,tstr[40]; int i=24,staid,nsat,type; - + if (rtcm->obsflag) rtcm->obs.n = rtcm->obsflag = 0; type=getbitu(rtcm->buff,i,12); i+=12; - + if (i+49<=rtcm->len*8) { staid=getbitu(rtcm->buff,i,12); i+=12; tod =getbitu(rtcm->buff,i,27)*0.001; i+=27; /* sec in a day */ @@ -593,12 +593,12 @@ static int decode_head1009(rtcm_t *rtcm, int *sync) } /* test station ID */ if (!test_staid(rtcm,staid)) return -1; - + adjday_glot(rtcm,tod); - + time2str(rtcm->time,tstr,2); trace(4,"decode_head1009: time=%s nsat=%d sync=%d\n",tstr,nsat,*sync); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," staid=%4d %s nsat=%2d sync=%d",staid,tstr,nsat,*sync); @@ -618,9 +618,9 @@ static int decode_type1010(rtcm_t *rtcm) { double pr1,cnr1,tt,cp1,freq1; int i=24+61,j,index,nsat,sync,prn,sat,code,fcn,ppr1,lock1,amb,sys=SYS_GLO; - + if ((nsat=decode_head1009(rtcm,&sync))<0) return -1; - + for (j=0;jlen*8;j++) { prn =getbitu(rtcm->buff,i, 6); i+= 6; code =getbitu(rtcm->buff,i, 1); i+= 1; @@ -642,7 +642,7 @@ static int decode_type1010(rtcm_t *rtcm) if ((index=obsindex(&rtcm->obs,rtcm->time,sat))<0) continue; pr1=pr1*0.02+amb*PRUNIT_GLO; rtcm->obs.data[index].P[0]=pr1; - + if (ppr1!=(int)0xFFF80000) { freq1=code2freq(SYS_GLO,CODE_L1C,fcn-7); cp1=adjcp(rtcm,sat,0,ppr1*0.0005*freq1/CLIGHT); @@ -668,9 +668,9 @@ static int decode_type1012(rtcm_t *rtcm) double pr1,cnr1,cnr2,tt,cp1,cp2,freq1,freq2; int i=24+61,j,index,nsat,sync,prn,sat,fcn,code1,code2,pr21,ppr1,ppr2; int lock1,lock2,amb,sys=SYS_GLO; - + if ((nsat=decode_head1009(rtcm,&sync))<0) return -1; - + for (j=0;jlen*8;j++) { prn =getbitu(rtcm->buff,i, 6); i+= 6; code1=getbitu(rtcm->buff,i, 1); i+= 1; @@ -697,7 +697,7 @@ static int decode_type1012(rtcm_t *rtcm) if ((index=obsindex(&rtcm->obs,rtcm->time,sat))<0) continue; pr1=pr1*0.02+amb*PRUNIT_GLO; rtcm->obs.data[index].P[0]=pr1; - + if (ppr1!=(int)0xFFF80000) { freq1=code2freq(SYS_GLO,CODE_L1C,fcn-7); cp1=adjcp(rtcm,sat,0,ppr1*0.0005*freq1/CLIGHT); @@ -706,7 +706,7 @@ static int decode_type1012(rtcm_t *rtcm) rtcm->obs.data[index].LLI[0]=lossoflock(rtcm,sat,0,lock1); rtcm->obs.data[index].SNR[0]=snratio(cnr1*0.25); rtcm->obs.data[index].code[0]=code1?CODE_L1P:CODE_L1C; - + if (pr21!=(int)0xFFFFE000) { rtcm->obs.data[index].P[1]=pr1+pr21*0.02; } @@ -734,7 +734,7 @@ static int decode_type1019(rtcm_t *rtcm) double toc,sqrtA,tt; char *msg; int i=24+12,prn,sat,week,sys=SYS_GPS; - + if (i+476<=rtcm->len*8) { prn =getbitu(rtcm->buff,i, 6); i+= 6; week =getbitu(rtcm->buff,i,10); i+=10; @@ -775,7 +775,7 @@ static int decode_type1019(rtcm_t *rtcm) sys=SYS_SBS; prn+=80; } trace(4,"decode_type1019: prn=%d iode=%d toe=%.0f\n",prn,eph.iode,eph.toes); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," prn=%2d iode=%3d iodc=%3d week=%d toe=%6.0f toc=%6.0f svh=%02X", @@ -851,7 +851,7 @@ static int decode_type1020(rtcm_t *rtcm) return -1; } trace(4,"decode_type1020: prn=%d tk=%02.0f:%02.0f:%02.0f\n",prn,tk_h,tk_m,tk_s); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," prn=%2d tk=%02.0f:%02.0f:%02.0f frq=%2d Bn=%d tb=%d", @@ -872,7 +872,7 @@ static int decode_type1020(rtcm_t *rtcm) if (toetod+43200.0) toe-=86400.0; geph.toe=utc2gpst(gpst2time(week,tow+toe)); /* utc->gpst */ - + if (!strstr(rtcm->opt,"-EPHALL")) { if (fabs(timediff(geph.toe,rtcm->nav.geph[prn-1].toe))<1.0&& geph.svh==rtcm->nav.geph[prn-1].svh) return 0; /* unchanged */ @@ -929,7 +929,7 @@ static int decode_type1029(rtcm_t *rtcm) { char *msg; int i=24+12,j,staid,mjd,tod,nchar,cunit; - + if (i+60<=rtcm->len*8) { staid=getbitu(rtcm->buff,i,12); i+=12; mjd =getbitu(rtcm->buff,i,16); i+=16; @@ -944,12 +944,12 @@ static int decode_type1029(rtcm_t *rtcm) if (i+nchar*8>rtcm->len*8) { trace(2,"rtcm3 1029 length error: len=%d nchar=%d\n",rtcm->len,nchar); return -1; - } + } for (j=0;jmsg[j]=getbitu(rtcm->buff,i,8); i+=8; } rtcm->msg[j]='\0'; - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," staid=%4d text=%s",staid,rtcm->msg); @@ -980,13 +980,13 @@ static int decode_type1033(rtcm_t *rtcm) char des[32]="",sno[32]="",rec[32]="",ver[32]="",rsn[32]=""; char *msg; int i=24+12,j,staid,n,m,n1,n2,n3,setup; - + n =getbitu(rtcm->buff,i+12,8); m =getbitu(rtcm->buff,i+28+8*n,8); n1=getbitu(rtcm->buff,i+36+8*(n+m),8); n2=getbitu(rtcm->buff,i+44+8*(n+m+n1),8); n3=getbitu(rtcm->buff,i+52+8*(n+m+n1+n2),8); - + if (i+60+8*(n+m+n1+n2+n3)<=rtcm->len*8) { staid=getbitu(rtcm->buff,i,12); i+=12+8; for (j=0;jsta.name,"%04d",staid); strncpy(rtcm->sta.antdes, des,n ); rtcm->sta.antdes [n] ='\0'; rtcm->sta.antsetup=setup; @@ -1027,7 +1027,7 @@ static int decode_type1033(rtcm_t *rtcm) strncpy(rtcm->sta.rectype,rec,n1); rtcm->sta.rectype[n1]='\0'; strncpy(rtcm->sta.recver, ver,n2); rtcm->sta.recver [n2]='\0'; strncpy(rtcm->sta.recsno, rsn,n3); rtcm->sta.recsno [n3]='\0'; - + trace(3,"rtcm3 1033: ant=%s:%s rec=%s:%s:%s\n",des,sno,rec,ver,rsn); return 5; } @@ -1068,7 +1068,7 @@ static int decode_type1041(rtcm_t *rtcm) double toc,sqrtA,tt; char *msg; int i=24+12,prn,sat,week,sys=SYS_IRN; - + if (i+482-12<=rtcm->len*8) { prn =getbitu(rtcm->buff,i, 6); i+= 6; week =getbitu(rtcm->buff,i,10); i+=10; @@ -1102,7 +1102,7 @@ static int decode_type1041(rtcm_t *rtcm) return -1; } trace(4,"decode_type1041: prn=%d iode=%d toe=%.0f\n",prn,eph.iode,eph.toes); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," prn=%2d iode=%3d week=%d toe=%6.0f toc=%6.0f svh=%02X", @@ -1138,7 +1138,7 @@ static int decode_type1044(rtcm_t *rtcm) double toc,sqrtA,tt; char *msg; int i=24+12,prn,sat,week,sys=SYS_QZS; - + if (i+473<=rtcm->len*8) { prn =getbitu(rtcm->buff,i, 4)+192; i+= 4; toc =getbitu(rtcm->buff,i,16)*16.0; i+=16; @@ -1175,7 +1175,7 @@ static int decode_type1044(rtcm_t *rtcm) return -1; } trace(4,"decode_type1044: prn=%d iode=%d toe=%.0f\n",prn,eph.iode,eph.toes); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," prn=%3d iode=%3d iodc=%3d week=%d toe=%6.0f toc=%6.0f svh=%02X", @@ -1212,7 +1212,7 @@ static int decode_type1045(rtcm_t *rtcm) double toc,sqrtA,tt; char *msg; int i=24+12,prn,sat,week,e5a_hs,e5a_dvs,rsv,sys=SYS_GAL; - + if (strstr(rtcm->opt,"-GALINAV")) return 0; if (i+484<=rtcm->len*8) { @@ -1250,7 +1250,7 @@ static int decode_type1045(rtcm_t *rtcm) return -1; } trace(4,"decode_type1045: prn=%d iode=%d toe=%.0f\n",prn,eph.iode,eph.toes); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," prn=%2d iode=%3d week=%d toe=%6.0f toc=%6.0f hs=%d dvs=%d", @@ -1291,7 +1291,7 @@ static int decode_type1046(rtcm_t *rtcm) double toc,sqrtA,tt; char *msg; int i=24+12,prn,sat,week,e5b_hs,e5b_dvs,e1_hs,e1_dvs,sys=SYS_GAL; - + if (strstr(rtcm->opt,"-GALFNAV")) return 0; if (i+492<=rtcm->len*8) { @@ -1331,7 +1331,7 @@ static int decode_type1046(rtcm_t *rtcm) return -1; } trace(4,"decode_type1046: prn=%d iode=%d toe=%.0f\n",prn,eph.iode,eph.toes); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," prn=%2d iode=%3d week=%d toe=%6.0f toc=%6.0f hs=%d %d dvs=%d %d", @@ -1372,7 +1372,7 @@ static int decode_type1042(rtcm_t *rtcm) double toc,sqrtA,tt; char *msg; int i=24+12,prn,sat,week,sys=SYS_CMP; - + if (i+499<=rtcm->len*8) { prn =getbitu(rtcm->buff,i, 6); i+= 6; week =getbitu(rtcm->buff,i,13); i+=13; @@ -1408,7 +1408,7 @@ static int decode_type1042(rtcm_t *rtcm) return -1; } trace(4,"decode_type1042: prn=%d iode=%d toe=%.0f\n",prn,eph.iode,eph.toes); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," prn=%2d iode=%3d iodc=%3d week=%d toe=%6.0f toc=%6.0f svh=%02X", @@ -1443,9 +1443,9 @@ static int decode_ssr_epoch(rtcm_t *rtcm, int sys, int subtype) { double tod,tow; int i=24+12; - + if (subtype==0) { /* RTCM SSR */ - + if (sys==SYS_GLO) { tod=getbitu(rtcm->buff,i,17); i+=17; adjday_glot(rtcm,tod); @@ -1468,7 +1468,7 @@ static int decode_ssr1_head(rtcm_t *rtcm, int sys, int subtype, int *sync, { char *msg,tstr[40]; int i=24+12,nsat,udi,provid=0,solid=0,ns; - + if (subtype==0) { /* RTCM SSR */ ns=(sys==SYS_QZS)?4:6; if (i+((sys==SYS_GLO)?53:50+ns)>rtcm->len*8) return -1; @@ -1491,11 +1491,11 @@ static int decode_ssr1_head(rtcm_t *rtcm, int sys, int subtype, int *sync, } nsat =getbitu(rtcm->buff,i,ns); i+=ns; *udint=ssrudint[udi]; - + time2str(rtcm->time,tstr,2); trace(4,"decode_ssr1_head: time=%s sys=%d subtype=%d nsat=%d sync=%d iod=%d" " provid=%d solid=%d\n",tstr,sys,subtype,nsat,*sync,*iod,provid,solid); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," %s nsat=%2d iod=%2d udi=%2d sync=%d",tstr,nsat,*iod,udi, @@ -1510,7 +1510,7 @@ static int decode_ssr2_head(rtcm_t *rtcm, int sys, int subtype, int *sync, { char *msg,tstr[40]; int i=24+12,nsat,udi,provid=0,solid=0,ns; - + if (subtype==0) { /* RTCM SSR */ ns=(sys==SYS_QZS)?4:6; if (i+((sys==SYS_GLO)?52:49+ns)>rtcm->len*8) return -1; @@ -1527,11 +1527,11 @@ static int decode_ssr2_head(rtcm_t *rtcm, int sys, int subtype, int *sync, solid =getbitu(rtcm->buff,i, 4); i+= 4; /* solution ID */ nsat =getbitu(rtcm->buff,i,ns); i+=ns; *udint=ssrudint[udi]; - + time2str(rtcm->time,tstr,2); trace(4,"decode_ssr2_head: time=%s sys=%d subtype=%d nsat=%d sync=%d iod=%d" " provid=%d solid=%d\n",tstr,sys,subtype,nsat,*sync,*iod,provid,solid); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," %s nsat=%2d iod=%2d udi=%2d sync=%d",tstr,nsat,*iod,udi, @@ -1545,9 +1545,9 @@ static int decode_ssr1(rtcm_t *rtcm, int sys, int subtype) { double udint,deph[3],ddeph[3]; int i,j,k,type,sync,iod,nsat,prn,sat,iode,iodcrc=0,refd=0,np,ni,nj,offp; - + type=getbitu(rtcm->buff,24,12); - + if ((nsat=decode_ssr1_head(rtcm,sys,subtype,&sync,&iod,&udint,&refd,&i))<0) { trace(2,"rtcm3 %d length error: len=%d\n",type,rtcm->len); return -1; @@ -1576,7 +1576,7 @@ static int decode_ssr1(rtcm_t *rtcm, int sys, int subtype) ddeph[0]=getbits(rtcm->buff,i,21)*1E-6; i+=21; ddeph[1]=getbits(rtcm->buff,i,19)*4E-6; i+=19; ddeph[2]=getbits(rtcm->buff,i,19)*4E-6; i+=19; - + if (!(sat=satno(sys,prn))) { trace(2,"rtcm3 %d satellite number error: prn=%d\n",type,prn); continue; @@ -1587,7 +1587,7 @@ static int decode_ssr1(rtcm_t *rtcm, int sys, int subtype) rtcm->ssr[sat-1].iode=iode; /* SBAS/BDS: toe/t0 modulo */ rtcm->ssr[sat-1].iodcrc=iodcrc; /* SBAS/BDS: IOD CRC */ rtcm->ssr[sat-1].refd=refd; - + for (k=0;k<3;k++) { rtcm->ssr[sat-1].deph [k]=deph [k]; rtcm->ssr[sat-1].ddeph[k]=ddeph[k]; @@ -1601,9 +1601,9 @@ static int decode_ssr2(rtcm_t *rtcm, int sys, int subtype) { double udint,dclk[3]; int i,j,k,type,sync,iod,nsat,prn,sat,np,offp; - + type=getbitu(rtcm->buff,24,12); - + if ((nsat=decode_ssr2_head(rtcm,sys,subtype,&sync,&iod,&udint,&i))<0) { trace(2,"rtcm3 %d length error: len=%d\n",type,rtcm->len); return -1; @@ -1627,7 +1627,7 @@ static int decode_ssr2(rtcm_t *rtcm, int sys, int subtype) dclk[0]=getbits(rtcm->buff,i,22)*1E-4; i+=22; dclk[1]=getbits(rtcm->buff,i,21)*1E-6; i+=21; dclk[2]=getbits(rtcm->buff,i,27)*2E-8; i+=27; - + if (!(sat=satno(sys,prn))) { trace(2,"rtcm3 %d satellite number error: prn=%d\n",type,prn); continue; @@ -1635,7 +1635,7 @@ static int decode_ssr2(rtcm_t *rtcm, int sys, int subtype) rtcm->ssr[sat-1].t0 [1]=rtcm->time; rtcm->ssr[sat-1].udi[1]=udint; rtcm->ssr[sat-1].iod[1]=iod; - + for (k=0;k<3;k++) { rtcm->ssr[sat-1].dclk[k]=dclk[k]; } @@ -1649,9 +1649,9 @@ static int decode_ssr3(rtcm_t *rtcm, int sys, int subtype) const uint8_t *sigs; double udint,bias,cbias[MAXCODE]; int i,j,k,type,mode,sync,iod,nsat,prn,sat,nbias,np,offp; - + type=getbitu(rtcm->buff,24,12); - + if ((nsat=decode_ssr2_head(rtcm,sys,subtype,&sync,&iod,&udint,&i))<0) { trace(2,"rtcm3 %d length error: len=%d\n",type,rtcm->len); return -1; @@ -1673,7 +1673,7 @@ static int decode_ssr3(rtcm_t *rtcm, int sys, int subtype) for (j=0;jlen*8;j++) { prn =getbitu(rtcm->buff,i,np)+offp; i+=np; nbias=getbitu(rtcm->buff,i, 5); i+= 5; - + for (k=0;klen*8;k++) { mode=getbitu(rtcm->buff,i, 5); i+= 5; @@ -1692,7 +1692,7 @@ static int decode_ssr3(rtcm_t *rtcm, int sys, int subtype) rtcm->ssr[sat-1].t0 [4]=rtcm->time; rtcm->ssr[sat-1].udi[4]=udint; rtcm->ssr[sat-1].iod[4]=iod; - + for (k=0;kssr[sat-1].cbias[k]=(float)cbias[k]; } @@ -1705,9 +1705,9 @@ static int decode_ssr4(rtcm_t *rtcm, int sys, int subtype) { double udint,deph[3],ddeph[3],dclk[3]; int i,j,k,type,nsat,sync,iod,prn,sat,iode,iodcrc=0,refd=0,np,ni,nj,offp; - + type=getbitu(rtcm->buff,24,12); - + if ((nsat=decode_ssr1_head(rtcm,sys,subtype,&sync,&iod,&udint,&refd,&i))<0) { trace(2,"rtcm3 %d length error: len=%d\n",type,rtcm->len); return -1; @@ -1736,11 +1736,11 @@ static int decode_ssr4(rtcm_t *rtcm, int sys, int subtype) ddeph[0]=getbits(rtcm->buff,i,21)*1E-6; i+=21; ddeph[1]=getbits(rtcm->buff,i,19)*4E-6; i+=19; ddeph[2]=getbits(rtcm->buff,i,19)*4E-6; i+=19; - + dclk [0]=getbits(rtcm->buff,i,22)*1E-4; i+=22; dclk [1]=getbits(rtcm->buff,i,21)*1E-6; i+=21; dclk [2]=getbits(rtcm->buff,i,27)*2E-8; i+=27; - + if (!(sat=satno(sys,prn))) { trace(2,"rtcm3 %d satellite number error: prn=%d\n",type,prn); continue; @@ -1751,7 +1751,7 @@ static int decode_ssr4(rtcm_t *rtcm, int sys, int subtype) rtcm->ssr[sat-1].iode=iode; rtcm->ssr[sat-1].iodcrc=iodcrc; rtcm->ssr[sat-1].refd=refd; - + for (k=0;k<3;k++) { rtcm->ssr[sat-1].deph [k]=deph [k]; rtcm->ssr[sat-1].ddeph[k]=ddeph[k]; @@ -1766,9 +1766,9 @@ static int decode_ssr5(rtcm_t *rtcm, int sys, int subtype) { double udint; int i,j,type,nsat,sync,iod,prn,sat,ura,np,offp; - + type=getbitu(rtcm->buff,24,12); - + if ((nsat=decode_ssr2_head(rtcm,sys,subtype,&sync,&iod,&udint,&i))<0) { trace(2,"rtcm3 %d length error: len=%d\n",type,rtcm->len); return -1; @@ -1790,7 +1790,7 @@ static int decode_ssr5(rtcm_t *rtcm, int sys, int subtype) for (j=0;jlen*8;j++) { prn=getbitu(rtcm->buff,i,np)+offp; i+=np; ura=getbitu(rtcm->buff,i, 6); i+= 6; - + if (!(sat=satno(sys,prn))) { trace(2,"rtcm3 %d satellite number error: prn=%d\n",type,prn); continue; @@ -1808,9 +1808,9 @@ static int decode_ssr6(rtcm_t *rtcm, int sys, int subtype) { double udint,hrclk; int i,j,type,nsat,sync,iod,prn,sat,np,offp; - + type=getbitu(rtcm->buff,24,12); - + if ((nsat=decode_ssr2_head(rtcm,sys,subtype,&sync,&iod,&udint,&i))<0) { trace(2,"rtcm3 %d length error: len=%d\n",type,rtcm->len); return -1; @@ -1832,7 +1832,7 @@ static int decode_ssr6(rtcm_t *rtcm, int sys, int subtype) for (j=0;jlen*8;j++) { prn =getbitu(rtcm->buff,i,np)+offp; i+=np; hrclk=getbits(rtcm->buff,i,22)*1E-4; i+=22; - + if (!(sat=satno(sys,prn))) { trace(2,"rtcm3 %d satellite number error: prn=%d\n",type,prn); continue; @@ -1852,7 +1852,7 @@ static int decode_ssr7_head(rtcm_t *rtcm, int sys, int subtype, int *sync, { char *msg,tstr[40]; int i=24+12,nsat,udi,provid=0,solid=0,ns; - + if (subtype==0) { /* RTCM SSR */ ns=(sys==SYS_QZS)?4:6; if (i+((sys==SYS_GLO)?54:51+ns)>rtcm->len*8) return -1; @@ -1871,11 +1871,11 @@ static int decode_ssr7_head(rtcm_t *rtcm, int sys, int subtype, int *sync, *mw =getbitu(rtcm->buff,i, 1); i+= 1; /* MW consistency indicator */ nsat =getbitu(rtcm->buff,i,ns); i+=ns; *udint=ssrudint[udi]; - + time2str(rtcm->time,tstr,2); trace(4,"decode_ssr7_head: time=%s sys=%d subtype=%d nsat=%d sync=%d iod=%d" " provid=%d solid=%d\n",tstr,sys,subtype,nsat,*sync,*iod,provid,solid); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," %s nsat=%2d iod=%2d udi=%2d sync=%d",tstr,nsat,*iod,udi, @@ -1891,9 +1891,9 @@ static int decode_ssr7(rtcm_t *rtcm, int sys, int subtype) double udint,bias,std=0.0,pbias[MAXCODE],stdpb[MAXCODE]; int i,j,k,type,mode,sync,iod,nsat,prn,sat,nbias,np,mw,offp,sii,swl; int dispe,sdc,yaw_ang,yaw_rate; - + type=getbitu(rtcm->buff,24,12); - + if ((nsat=decode_ssr7_head(rtcm,sys,subtype,&sync,&iod,&udint,&dispe,&mw, &i))<0) { trace(2,"rtcm3 %d length error: len=%d\n",type,rtcm->len); @@ -1917,7 +1917,7 @@ static int decode_ssr7(rtcm_t *rtcm, int sys, int subtype) nbias =getbitu(rtcm->buff,i, 5); i+= 5; yaw_ang =getbitu(rtcm->buff,i, 9); i+= 9; yaw_rate=getbits(rtcm->buff,i, 8); i+= 8; - + for (k=0;klen*8;k++) { mode=getbitu(rtcm->buff,i, 5); i+= 5; @@ -1945,7 +1945,7 @@ static int decode_ssr7(rtcm_t *rtcm, int sys, int subtype) rtcm->ssr[sat-1].iod[5]=iod; rtcm->ssr[sat-1].yaw_ang =yaw_ang / 256.0*180.0; /* (deg) */ rtcm->ssr[sat-1].yaw_rate=yaw_rate/8192.0*180.0; /* (deg/s) */ - + for (k=0;kssr[sat-1].pbias[k]=pbias[k]; rtcm->ssr[sat-1].stdpb[k]=(float)stdpb[k]; @@ -1958,18 +1958,18 @@ static void sigindex(int sys, const uint8_t *code, int n, const char *opt, int *idx) { int i,nex,pri,pri_h[8]={0},index[8]={0},ex[32]={0}; - + /* test code priority */ for (i=0;i=NFREQ) { /* save as extended signal if idx >= NFREQ */ ex[i]=1; continue; } /* code priority */ pri=getcodepri(sys,code[i],opt); - + /* select highest priority signal */ if (pri>pri_h[idx[i]]) { if (index[idx[i]]) ex[index[idx[i]]-1]=1; @@ -2002,9 +2002,9 @@ static void save_msm_obs(rtcm_t *rtcm, int sys, msm_h_t *h, const double *r, uint8_t code[32]; char *msm_type="",*q=NULL; int i,j,k,type,prn,sat,fcn,index=0,idx[32]; - + type=getbitu(rtcm->buff,24,12); - + switch (sys) { case SYS_GPS: msm_type=q=rtcm->msmtype[0]; break; case SYS_GLO: msm_type=q=rtcm->msmtype[1]; break; @@ -2029,27 +2029,27 @@ static void save_msm_obs(rtcm_t *rtcm, int sys, msm_h_t *h, const double *r, /* signal to rinex obs type */ code[i]=obs2code(sig[i]); idx[i]=code2idx(sys,code[i]); - + if (code[i]!=CODE_NONE) { if (q) q+=sprintf(q,"L%s%s",sig[i],insig-1?",":""); } else { if (q) q+=sprintf(q,"(%d)%s",h->sigs[i],insig-1?",":""); - + trace(2,"rtcm3 %d: unknown signal id=%2d\n",type,h->sigs[i]); } } trace(3,"rtcm3 %d: signals=%s\n",type,msm_type); - + /* get signal index */ sigindex(sys,code,h->nsig,rtcm->opt,idx); - + for (i=j=0;insat;i++) { - + prn=h->sats[i]; if (sys==SYS_QZS) prn+=MINPRNQZS-1; else if (sys==SYS_SBS) prn+=MINPRNSBS-1; - + if ((sat=satno(sys,prn))) { tt=timediff(rtcm->obs.data[0].time,rtcm->time); if (rtcm->obsflag||fabs(tt)>1E-9) { @@ -2078,10 +2078,10 @@ static void save_msm_obs(rtcm_t *rtcm, int sys, msm_h_t *h, const double *r, } for (k=0;knsig;k++) { if (!h->cellmask[k+i*h->nsig]) continue; - + if (sat&&index>=0&&idx[k]>=0) { freq=fcn<-7?0.0:code2freq(sys,code[k],fcn); - + /* pseudorange (m) */ if (r[i]!=0.0&&pr[j]>-1E12) { rtcm->obs.data[index].P[idx[k]]=r[i]+pr[j]; @@ -2112,15 +2112,15 @@ static int decode_msm_head(rtcm_t *rtcm, int sys, int *sync, int *iod, double tow,tod; char *msg,tstr[40]; int i=24,j,dow,mask,staid,type,ncell=0; - + if (rtcm->obsflag) rtcm->obs.n = rtcm->obsflag = 0; type=getbitu(rtcm->buff,i,12); i+=12; - + *h=h0; if (i+157<=rtcm->len*8) { staid =getbitu(rtcm->buff,i,12); i+=12; - + if (sys==SYS_GLO) { dow =getbitu(rtcm->buff,i, 3); i+= 3; tod =getbitu(rtcm->buff,i,27)*0.001; i+=27; @@ -2157,7 +2157,7 @@ static int decode_msm_head(rtcm_t *rtcm, int sys, int *sync, int *iod, } /* test station id */ if (!test_staid(rtcm,staid)) return -1; - + if (h->nsat*h->nsig>64) { trace(2,"rtcm3 %d number of sats and sigs error: nsat=%d nsig=%d\n", type,h->nsat,h->nsig); @@ -2173,11 +2173,11 @@ static int decode_msm_head(rtcm_t *rtcm, int sys, int *sync, int *iod, if (h->cellmask[j]) ncell++; } *hsize=i; - + time2str(rtcm->time,tstr,2); trace(4,"decode_head_msm: time=%s sys=%d staid=%d nsat=%d nsig=%d sync=%d iod=%d ncell=%d\n", tstr,sys,staid,h->nsat,h->nsig,*sync,*iod,ncell); - + if (rtcm->outtype) { msg=rtcm->msgtype+strlen(rtcm->msgtype); sprintf(msg," staid=%4d %s nsat=%2d nsig=%2d iod=%2d ncell=%2d sync=%d", @@ -2200,12 +2200,12 @@ static int decode_msm4(rtcm_t *rtcm, int sys) msm_h_t h={0}; double r[64],pr[64],cp[64],cnr[64]; int i,j,type,sync,iod,ncell,rng,rng_m,prv,cpv,lock[64],half[64]; - + type=getbitu(rtcm->buff,24,12); - + /* decode msm header */ if ((ncell=decode_msm_head(rtcm,sys,&sync,&iod,&h,&i))<0) return -1; - + if (i+h.nsat*18+ncell*48>rtcm->len*8) { trace(2,"rtcm3 %d length error: nsat=%d ncell=%d len=%d\n",type,h.nsat, ncell,rtcm->len); @@ -2214,7 +2214,7 @@ static int decode_msm4(rtcm_t *rtcm, int sys) } for (j=0;jbuff,i, 8); i+= 8; @@ -2244,7 +2244,7 @@ static int decode_msm4(rtcm_t *rtcm, int sys) } /* save obs data in msm message */ save_msm_obs(rtcm,sys,&h,r,pr,cp,NULL,NULL,cnr,lock,NULL,half); - + rtcm->obsflag=!sync; return sync?0:1; } @@ -2255,12 +2255,12 @@ static int decode_msm5(rtcm_t *rtcm, int sys) double r[64],rr[64],pr[64],cp[64],rrf[64],cnr[64]; int i,j,type,sync,iod,ncell,rng,rng_m,rate,prv,cpv,rrv,lock[64]; int ex[64],half[64]; - + type=getbitu(rtcm->buff,24,12); - + /* decode msm header */ if ((ncell=decode_msm_head(rtcm,sys,&sync,&iod,&h,&i))<0) return -1; - + if (i+h.nsat*36+ncell*63>rtcm->len*8) { trace(2,"rtcm3 %d length error: nsat=%d ncell=%d len=%d\n",type,h.nsat, ncell,rtcm->len); @@ -2271,7 +2271,7 @@ static int decode_msm5(rtcm_t *rtcm, int sys) r[j]=rr[j]=0.0; ex[j]=15; } for (j=0;jbuff,i, 8); i+= 8; @@ -2312,7 +2312,7 @@ static int decode_msm5(rtcm_t *rtcm, int sys) } /* save obs data in msm message */ save_msm_obs(rtcm,sys,&h,r,pr,cp,rr,rrf,cnr,lock,ex,half); - + rtcm->obsflag=!sync; return sync?0:1; } @@ -2322,12 +2322,12 @@ static int decode_msm6(rtcm_t *rtcm, int sys) msm_h_t h={0}; double r[64],pr[64],cp[64],cnr[64]; int i,j,type,sync,iod,ncell,rng,rng_m,prv,cpv,lock[64],half[64]; - + type=getbitu(rtcm->buff,24,12); - + /* decode msm header */ if ((ncell=decode_msm_head(rtcm,sys,&sync,&iod,&h,&i))<0) return -1; - + if (i+h.nsat*18+ncell*65>rtcm->len*8) { trace(2,"rtcm3 %d length error: nsat=%d ncell=%d len=%d\n",type,h.nsat, ncell,rtcm->len); @@ -2336,7 +2336,7 @@ static int decode_msm6(rtcm_t *rtcm, int sys) } for (j=0;jbuff,i, 8); i+= 8; @@ -2366,7 +2366,7 @@ static int decode_msm6(rtcm_t *rtcm, int sys) } /* save obs data in msm message */ save_msm_obs(rtcm,sys,&h,r,pr,cp,NULL,NULL,cnr,lock,NULL,half); - + rtcm->obsflag=!sync; return sync?0:1; } @@ -2377,12 +2377,12 @@ static int decode_msm7(rtcm_t *rtcm, int sys) double r[64],rr[64],pr[64],cp[64],rrf[64],cnr[64]; int i,j,type,sync,iod,ncell,rng,rng_m,rate,prv,cpv,rrv,lock[64]; int ex[64],half[64]; - + type=getbitu(rtcm->buff,24,12); - + /* decode msm header */ if ((ncell=decode_msm_head(rtcm,sys,&sync,&iod,&h,&i))<0) return -1; - + if (i+h.nsat*36+ncell*80>rtcm->len*8) { trace(2,"rtcm3 %d length error: nsat=%d ncell=%d len=%d\n",type,h.nsat, ncell,rtcm->len); @@ -2393,7 +2393,7 @@ static int decode_msm7(rtcm_t *rtcm, int sys) r[j]=rr[j]=0.0; ex[j]=15; } for (j=0;jbuff,i, 8); i+= 8; @@ -2440,7 +2440,7 @@ static int decode_msm7(rtcm_t *rtcm, int sys) } /* save obs data in msm message */ save_msm_obs(rtcm,sys,&h,r,pr,cp,rr,rrf,cnr,lock,ex,half); - + rtcm->obsflag=!sync; return sync?0:1; } @@ -2448,7 +2448,7 @@ static int decode_msm7(rtcm_t *rtcm, int sys) static int decode_type1230(rtcm_t *rtcm) { int i=24+12,j,staid,align,mask,bias; - + if (i+20>=rtcm->len*8) { trace(2,"rtcm3 1230: length error len=%d\n",rtcm->len); return -1; @@ -2456,14 +2456,14 @@ static int decode_type1230(rtcm_t *rtcm) staid=getbitu(rtcm->buff,i,12); i+=12; align=getbitu(rtcm->buff,i, 1); i+= 1+3; mask =getbitu(rtcm->buff,i, 4); i+= 4; - + if (rtcm->outtype) { sprintf(rtcm->msgtype+strlen(rtcm->msgtype), " staid=%4d align=%d mask=0x%X",staid,align,mask); } /* test station ID */ if (!test_staid(rtcm,staid)) return -1; - + rtcm->sta.glo_cp_align=align; for (j=0;j<4;j++) { rtcm->sta.glo_cp_bias[j]=0.0; @@ -2481,9 +2481,9 @@ static int decode_type1230(rtcm_t *rtcm) static int decode_type4073(rtcm_t *rtcm) { int i=24+12,subtype; - + subtype=getbitu(rtcm->buff,i,4); i+=4; - + if (rtcm->outtype) { sprintf(rtcm->msgtype+strlen(rtcm->msgtype)," subtype=%d",subtype); } @@ -2494,14 +2494,14 @@ static int decode_type4073(rtcm_t *rtcm) static int decode_type4076(rtcm_t *rtcm) { int i=24+12,ver,subtype; - + if (i+3+8>=rtcm->len*8) { trace(2,"rtcm3 4076: length error len=%d\n",rtcm->len); return -1; } ver =getbitu(rtcm->buff,i,3); i+=3; subtype=getbitu(rtcm->buff,i,8); i+=8; - + if (rtcm->outtype) { sprintf(rtcm->msgtype+strlen(rtcm->msgtype)," ver=%d subtype=%3d",ver, subtype); @@ -2558,9 +2558,9 @@ extern int decode_rtcm3(rtcm_t *rtcm) { double tow; int ret=0,type=getbitu(rtcm->buff,24,12),week; - + trace(3,"decode_rtcm3: len=%3d type=%d\n",rtcm->len,type); - + if (rtcm->outtype) { sprintf(rtcm->msgtype,"RTCM %4d (%4d):",type,rtcm->len); } diff --git a/src/rtcm3e.c b/src/rtcm3e.c index 923bd1e97..4841351ea 100644 --- a/src/rtcm3e.c +++ b/src/rtcm3e.c @@ -164,7 +164,7 @@ static int to_msm_lock(double lock) static int to_msm_lock_ex(double lock) { int lock_ms = (int)(lock * 1000.0); - + if (lock<0.0 ) return 0; if (lock<0.064 ) return lock_ms; if (lock<0.128 ) return (lock_ms+64 )/2; @@ -247,14 +247,14 @@ static void gen_obs_gps(rtcm_t *rtcm, const obsd_t *data, int *code1, int *pr1, { double lam1,lam2,pr1c=0.0,ppr; int lt1,lt2; - + lam1=CLIGHT/FREQL1; lam2=CLIGHT/FREQL2; *pr1=*amb=0; if (ppr1) *ppr1=0xFFF80000; /* invalid values */ if (pr21) *pr21=0xFFFFE000; if (ppr2) *ppr2=0xFFF80000; - + /* L1 peudorange */ if (data->P[0]!=0.0&&data->code[0]) { *amb=(int)floor(data->P[0]/PRUNIT_GPS); @@ -278,7 +278,7 @@ static void gen_obs_gps(rtcm_t *rtcm, const obsd_t *data, int *code1, int *pr1, } lt1=locktime(data->time,rtcm->lltime[data->sat-1] ,data->LLI[0]); lt2=locktime(data->time,rtcm->lltime[data->sat-1]+1,data->LLI[1]); - + if (lock1) *lock1=to_lock(lt1); if (lock2) *lock2=to_lock(lt2); if (cnr1 ) *cnr1=ROUND(data->SNR[0]/0.25); @@ -293,7 +293,7 @@ static void gen_obs_glo(rtcm_t *rtcm, const obsd_t *data, int fcn, int *code1, { double lam1=0.0,lam2=0.0,pr1c=0.0,ppr; int lt1,lt2; - + if (fcn>=0) { /* fcn+7 */ lam1=CLIGHT/(FREQ1_GLO+DFRQ1_GLO*(fcn-7)); lam2=CLIGHT/(FREQ2_GLO+DFRQ2_GLO*(fcn-7)); @@ -302,7 +302,7 @@ static void gen_obs_glo(rtcm_t *rtcm, const obsd_t *data, int fcn, int *code1, if (ppr1) *ppr1=0xFFF80000; /* invalid values */ if (pr21) *pr21=0xFFFFE000; if (ppr2) *ppr2=0xFFF80000; - + /* L1 pseudorange */ if (data->P[0]!=0.0) { *amb=(int)floor(data->P[0]/PRUNIT_GLO); @@ -327,7 +327,7 @@ static void gen_obs_glo(rtcm_t *rtcm, const obsd_t *data, int fcn, int *code1, } lt1=locktime(data->time,rtcm->lltime[data->sat-1] ,data->LLI[0]); lt2=locktime(data->time,rtcm->lltime[data->sat-1]+1,data->LLI[1]); - + if (lock1) *lock1=to_lock(lt1); if (lock2) *lock2=to_lock(lt2); if (cnr1 ) *cnr1=ROUND(data->SNR[0]/0.25); @@ -340,12 +340,12 @@ static int encode_head(int type, rtcm_t *rtcm, int sys, int sync, int nsat) { double tow; int i=24,week,epoch; - + trace(4,"encode_head: type=%d sync=%d sys=%d nsat=%d\n",type,sync,sys,nsat); - + setbitu(rtcm->buff,i,12,type ); i+=12; /* message no */ setbitu(rtcm->buff,i,12,rtcm->staid); i+=12; /* ref station id */ - + if (sys==SYS_GLO) { tow=time2gpst(timeadd(gpst2utc(rtcm->time),10800.0),&week); epoch=ROUND(fmod(tow,86400.0)/0.001); @@ -366,9 +366,9 @@ static int encode_head(int type, rtcm_t *rtcm, int sys, int sync, int nsat) static int encode_type1001(rtcm_t *rtcm, int sync) { int code1,pr1,ppr1,lock1,amb; - + trace(3,"encode_type1001: sync=%d\n",sync); - + int nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat,NULL); @@ -377,19 +377,19 @@ static int encode_type1001(rtcm_t *rtcm, int sync) } /* encode header */ int i = encode_head(1001,rtcm,SYS_GPS,sync,nsat); - + nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat,&prn); if (!(sys&(SYS_GPS|SYS_SBS))) continue; nsat++; - + if (sys==SYS_SBS) prn-=80; /* 40-58: sbas 120-138 */ - + /* generate obs field data gps */ gen_obs_gps(rtcm,rtcm->obs.data+j,&code1,&pr1,&ppr1,&lock1,&amb,NULL, NULL,NULL,NULL,NULL,NULL); - + setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i, 1,code1); i+= 1; setbitu(rtcm->buff,i,24,pr1 ); i+=24; @@ -403,9 +403,9 @@ static int encode_type1001(rtcm_t *rtcm, int sync) static int encode_type1002(rtcm_t *rtcm, int sync) { int code1,pr1,ppr1,lock1,amb,cnr1; - + trace(3,"encode_type1002: sync=%d\n",sync); - + int nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat,NULL); @@ -414,19 +414,19 @@ static int encode_type1002(rtcm_t *rtcm, int sync) } /* encode header */ int i = encode_head(1002,rtcm,SYS_GPS,sync,nsat); - + nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat,&prn); if (!(sys&(SYS_GPS|SYS_SBS))) continue; nsat++; - + if (sys==SYS_SBS) prn-=80; /* 40-58: sbas 120-138 */ - + /* generate obs field data gps */ gen_obs_gps(rtcm,rtcm->obs.data+j,&code1,&pr1,&ppr1,&lock1,&amb,&cnr1, NULL,NULL,NULL,NULL,NULL); - + setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i, 1,code1); i+= 1; setbitu(rtcm->buff,i,24,pr1 ); i+=24; @@ -442,9 +442,9 @@ static int encode_type1002(rtcm_t *rtcm, int sync) static int encode_type1003(rtcm_t *rtcm, int sync) { int code1,pr1,ppr1,lock1,amb,code2,pr21,ppr2,lock2; - + trace(3,"encode_type1003: sync=%d\n",sync); - + int nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat,NULL); @@ -453,19 +453,19 @@ static int encode_type1003(rtcm_t *rtcm, int sync) } /* encode header */ int i = encode_head(1003,rtcm,SYS_GPS,sync,nsat); - + nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat,&prn); if (!(sys&(SYS_GPS|SYS_SBS))) continue; nsat++; - + if (sys==SYS_SBS) prn-=80; /* 40-58: sbas 120-138 */ - + /* generate obs field data gps */ gen_obs_gps(rtcm,rtcm->obs.data+j,&code1,&pr1,&ppr1,&lock1,&amb, NULL,&code2,&pr21,&ppr2,&lock2,NULL); - + setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i, 1,code1); i+= 1; setbitu(rtcm->buff,i,24,pr1 ); i+=24; @@ -483,9 +483,9 @@ static int encode_type1003(rtcm_t *rtcm, int sync) static int encode_type1004(rtcm_t *rtcm, int sync) { int code1,pr1,ppr1,lock1,amb,cnr1,code2,pr21,ppr2,lock2,cnr2; - + trace(3,"encode_type1004: sync=%d\n",sync); - + int nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat,NULL); @@ -494,19 +494,19 @@ static int encode_type1004(rtcm_t *rtcm, int sync) } /* encode header */ int i = encode_head(1004,rtcm,SYS_GPS,sync,nsat); - + nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat,&prn); if (!(sys&(SYS_GPS|SYS_SBS))) continue; nsat++; - + if (sys==SYS_SBS) prn-=80; /* 40-58: sbas 120-138 */ - + /* generate obs field data gps */ gen_obs_gps(rtcm,rtcm->obs.data+j,&code1,&pr1,&ppr1,&lock1,&amb, &cnr1,&code2,&pr21,&ppr2,&lock2,&cnr2); - + setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i, 1,code1); i+= 1; setbitu(rtcm->buff,i,24,pr1 ); i+=24; @@ -528,9 +528,9 @@ static int encode_type1005(rtcm_t *rtcm, int sync) { double *p=rtcm->sta.pos; int i=24; - + trace(3,"encode_type1005: sync=%d\n",sync); - + setbitu(rtcm->buff,i,12,1005 ); i+=12; /* message no */ setbitu(rtcm->buff,i,12,rtcm->staid); i+=12; /* ref station id */ setbitu(rtcm->buff,i, 6,0 ); i+= 6; /* itrf realization year */ @@ -552,9 +552,9 @@ static int encode_type1006(rtcm_t *rtcm, int sync) { double *p=rtcm->sta.pos; int i=24,hgt=0; - + trace(3,"encode_type1006: sync=%d\n",sync); - + if (0.0<=rtcm->sta.hgt&&rtcm->sta.hgt<=6.5535) { hgt=ROUND(rtcm->sta.hgt/0.0001); } @@ -583,12 +583,12 @@ static int encode_type1007(rtcm_t *rtcm, int sync) { int i=24,j,antsetup=rtcm->sta.antsetup; int n=MIN((int)strlen(rtcm->sta.antdes),31); - + trace(3,"encode_type1007: sync=%d\n",sync); - + setbitu(rtcm->buff,i,12,1007 ); i+=12; /* message no */ setbitu(rtcm->buff,i,12,rtcm->staid); i+=12; /* ref station id */ - + /* antenna descriptor */ setbitu(rtcm->buff,i,8,n); i+=8; for (j=0;jsta.antsetup; int n=MIN((int)strlen(rtcm->sta.antdes),31); int m=MIN((int)strlen(rtcm->sta.antsno),31); - + trace(3,"encode_type1008: sync=%d\n",sync); - + setbitu(rtcm->buff,i,12,1008 ); i+=12; /* message no */ setbitu(rtcm->buff,i,12,rtcm->staid); i+=12; /* ref station id */ - + /* antenna descriptor */ setbitu(rtcm->buff,i,8,n); i+=8; for (j=0;jbuff,i,8,rtcm->sta.antdes[j]); i+=8; } setbitu(rtcm->buff,i,8,antsetup); i+=8; /* antenna setup id */ - + /* antenna serial number */ setbitu(rtcm->buff,i,8,m); i+=8; for (j=0;jobs.n&&nsatobs.data[j].sat, prn; @@ -645,12 +645,12 @@ static int encode_type1009(rtcm_t *rtcm, int sync) int fcn = fcn_glo(sat,rtcm); if (fcn < 0) continue; /* fcn+7 */ nsat++; - + /* generate obs field data glonass */ int code1,pr1,ppr1,lock1,amb; gen_obs_glo(rtcm,rtcm->obs.data+j,fcn,&code1,&pr1,&ppr1,&lock1,&amb, NULL,NULL,NULL,NULL,NULL,NULL); - + setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i, 1,code1); i+= 1; setbitu(rtcm->buff,i, 5,fcn ); i+= 5; /* fcn+7 */ @@ -665,7 +665,7 @@ static int encode_type1009(rtcm_t *rtcm, int sync) static int encode_type1010(rtcm_t *rtcm, int sync) { trace(3,"encode_type1010: sync=%d\n",sync); - + int nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat; @@ -675,7 +675,7 @@ static int encode_type1010(rtcm_t *rtcm, int sync) } /* encode header */ int i = encode_head(1010,rtcm,SYS_GLO,sync,nsat); - + nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat, prn; @@ -683,12 +683,12 @@ static int encode_type1010(rtcm_t *rtcm, int sync) int fcn = fcn_glo(sat,rtcm); if (fcn < 0) continue; /* fcn+7 */ nsat++; - + /* generate obs field data glonass */ int code1,pr1,ppr1,lock1,amb,cnr1; gen_obs_glo(rtcm,rtcm->obs.data+j,fcn,&code1,&pr1,&ppr1,&lock1,&amb, &cnr1,NULL,NULL,NULL,NULL,NULL); - + setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i, 1,code1); i+= 1; setbitu(rtcm->buff,i, 5,fcn ); i+= 5; /* fcn+7 */ @@ -705,7 +705,7 @@ static int encode_type1010(rtcm_t *rtcm, int sync) static int encode_type1011(rtcm_t *rtcm, int sync) { trace(3,"encode_type1011: sync=%d\n",sync); - + int nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat; @@ -723,12 +723,12 @@ static int encode_type1011(rtcm_t *rtcm, int sync) int fcn = fcn_glo(sat,rtcm); if (fcn < 0) continue; /* fcn+7 */ nsat++; - + /* generate obs field data glonass */ int code1,pr1,ppr1,lock1,amb,code2,pr21,ppr2,lock2; gen_obs_glo(rtcm,rtcm->obs.data+j,fcn,&code1,&pr1,&ppr1,&lock1,&amb, NULL,&code2,&pr21,&ppr2,&lock2,NULL); - + setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i, 1,code1); i+= 1; setbitu(rtcm->buff,i, 5,fcn ); i+= 5; /* fcn+7 */ @@ -747,7 +747,7 @@ static int encode_type1011(rtcm_t *rtcm, int sync) static int encode_type1012(rtcm_t *rtcm, int sync) { trace(3,"encode_type1012: sync=%d\n",sync); - + int nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat; @@ -757,7 +757,7 @@ static int encode_type1012(rtcm_t *rtcm, int sync) } /* encode header */ int i = encode_head(1012,rtcm,SYS_GLO,sync,nsat); - + nsat = 0; for (int j=0;jobs.n&&nsatobs.data[j].sat, prn; @@ -765,12 +765,12 @@ static int encode_type1012(rtcm_t *rtcm, int sync) int fcn = fcn_glo(sat,rtcm); if (fcn < 0) continue; /* fcn+7 */ nsat++; - + /* generate obs field data glonass */ int code1,pr1,ppr1,lock1,amb,cnr1,code2,pr21,ppr2,lock2,cnr2; gen_obs_glo(rtcm,rtcm->obs.data+j,fcn,&code1,&pr1,&ppr1,&lock1,&amb, &cnr1,&code2,&pr21,&ppr2,&lock2,&cnr2); - + setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i, 1,code1); i+= 1; setbitu(rtcm->buff,i, 5,fcn ); i+= 5; /* fcn+7 */ @@ -795,9 +795,9 @@ static int encode_type1019(rtcm_t *rtcm, int sync) uint32_t sqrtA,e; int i=24,prn,week,toe,toc,i0,OMG0,omg,M0,deln,idot,OMGd,crs,crc; int cus,cuc,cis,cic,af0,af1,af2,tgd; - + trace(3,"encode_type1019: sync=%d\n",sync); - + if (satsys(rtcm->ephsat,&prn)!=SYS_GPS) return 0; eph=rtcm->nav.eph+rtcm->ephsat-1; if (eph->sat!=rtcm->ephsat) return 0; @@ -823,7 +823,7 @@ static int encode_type1019(rtcm_t *rtcm, int sync) af1 =ROUND(eph->f1 /P2_43); af2 =ROUND(eph->f2 /P2_55); tgd =ROUND(eph->tgd[0]/P2_31); - + setbitu(rtcm->buff,i,12,1019 ); i+=12; setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i,10,week ); i+=10; @@ -866,31 +866,31 @@ static int encode_type1020(rtcm_t *rtcm, int sync) double ep[6]; int i=24,j,prn,tk_h,tk_m,tk_s,tb,pos[3],vel[3],acc[3],gamn,taun,dtaun; int fcn,NT; - + trace(3,"encode_type1020: sync=%d\n",sync); - + if (satsys(rtcm->ephsat,&prn)!=SYS_GLO) return 0; geph=rtcm->nav.geph+prn-1; if (geph->sat!=rtcm->ephsat) return 0; fcn=geph->frq+7; - + /* time of frame within day (utc(su) + 3 hr) */ time=timeadd(gpst2utc(geph->tof),10800.0); time2epoch(time,ep); tk_h=(int)ep[3]; tk_m=(int)ep[4]; tk_s=ROUND(ep[5]/30.0); - + /* # of days since jan 1 in leap year */ ep[0]=floor(ep[0]/4.0)*4.0; ep[1]=ep[2]=1.0; ep[3]=ep[4]=ep[5]=0.0; NT=(int)floor(timediff(time,epoch2time(ep))/86400.+1.0); - + /* index of time interval within day (utc(su) + 3 hr) */ time=timeadd(gpst2utc(geph->toe),10800.0); time2epoch(time,ep); tb=ROUND((ep[3]*3600.0+ep[4]*60.0+ep[5])/900.0); - + for (j=0;j<3;j++) { pos[j]=ROUND(geph->pos[j]/P2_11/1E3); vel[j]=ROUND(geph->vel[j]/P2_20/1E3); @@ -909,7 +909,7 @@ static int encode_type1020(rtcm_t *rtcm, int sync) int P2 = (geph->flags >> 4) & 1; int P1 = (geph->flags >> 2) & 3; int P = geph->flags & 3; - + setbitu(rtcm->buff,i,12,1020 ); i+=12; setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i, 5,fcn ); i+= 5; @@ -961,18 +961,18 @@ static int encode_type1033(rtcm_t *rtcm, int sync) int I=MIN((int)strlen(rtcm->sta.rectype),31); int J=MIN((int)strlen(rtcm->sta.recver ),31); int K=MIN((int)strlen(rtcm->sta.recsno ),31); - + trace(3,"encode_type1033: sync=%d\n",sync); - + setbitu(rtcm->buff,i,12,1033 ); i+=12; setbitu(rtcm->buff,i,12,rtcm->staid); i+=12; - + setbitu(rtcm->buff,i,8,n); i+= 8; for (j=0;jbuff,i,8,rtcm->sta.antdes[j]); i+=8; } setbitu(rtcm->buff,i,8,antsetup); i+= 8; - + setbitu(rtcm->buff,i,8,m); i+= 8; for (j=0;jbuff,i,8,rtcm->sta.antsno[j]); i+=8; @@ -999,9 +999,9 @@ static int encode_type1041(rtcm_t *rtcm, int sync) uint32_t sqrtA,e; int i=24,prn,week,toe,toc,i0,OMG0,omg,M0,deln,idot,OMGd,crs,crc; int cus,cuc,cis,cic,af0,af1,af2,tgd; - + trace(3,"encode_type1041: sync=%d\n",sync); - + if (satsys(rtcm->ephsat,&prn)!=SYS_IRN) return 0; eph=rtcm->nav.eph+rtcm->ephsat-1; if (eph->sat!=rtcm->ephsat) return 0; @@ -1027,7 +1027,7 @@ static int encode_type1041(rtcm_t *rtcm, int sync) af1 =ROUND(eph->f1 /P2_43); af2 =ROUND(eph->f2 /P2_55); tgd =ROUND(eph->tgd[0]/P2_31); - + setbitu(rtcm->buff,i,12,1041 ); i+=12; setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i,10,week ); i+=10; @@ -1065,9 +1065,9 @@ static int encode_type1044(rtcm_t *rtcm, int sync) uint32_t sqrtA,e; int i=24,prn,week,toe,toc,i0,OMG0,omg,M0,deln,idot,OMGd,crs,crc; int cus,cuc,cis,cic,af0,af1,af2,tgd; - + trace(3,"encode_type1044: sync=%d\n",sync); - + if (satsys(rtcm->ephsat,&prn)!=SYS_QZS) return 0; eph=rtcm->nav.eph+rtcm->ephsat-1; if (eph->sat!=rtcm->ephsat) return 0; @@ -1093,7 +1093,7 @@ static int encode_type1044(rtcm_t *rtcm, int sync) af1 =ROUND(eph->f1 /P2_43); af2 =ROUND(eph->f2 /P2_55); tgd =ROUND(eph->tgd[0]/P2_31); - + setbitu(rtcm->buff,i,12,1044 ); i+=12; setbitu(rtcm->buff,i, 4,prn-192 ); i+= 4; setbitu(rtcm->buff,i,16,toc ); i+=16; @@ -1134,9 +1134,9 @@ static int encode_type1045(rtcm_t *rtcm, int sync) uint32_t sqrtA,e; int i=24,prn,week,toe,toc,i0,OMG0,omg,M0,deln,idot,OMGd,crs,crc; int cus,cuc,cis,cic,af0,af1,af2,bgd1,bgd2,oshs,osdvs; - + trace(3,"encode_type1045: sync=%d\n",sync); - + if (satsys(rtcm->ephsat,&prn)!=SYS_GAL) return 0; eph=rtcm->nav.eph+rtcm->ephsat-1+MAXSAT; /* F/NAV */ if (eph->sat!=rtcm->ephsat) return 0; @@ -1204,9 +1204,9 @@ static int encode_type1046(rtcm_t *rtcm, int sync) uint32_t sqrtA,e; int i=24,prn,week,toe,toc,i0,OMG0,omg,M0,deln,idot,OMGd,crs,crc; int cus,cuc,cis,cic,af0,af1,af2,bgd1,bgd2,oshs1,osdvs1,oshs2,osdvs2; - + trace(3,"encode_type1046: sync=%d\n",sync); - + if (satsys(rtcm->ephsat,&prn)!=SYS_GAL) return 0; eph=rtcm->nav.eph+rtcm->ephsat-1; /* I/NAV */ if (eph->sat!=rtcm->ephsat) return 0; @@ -1278,9 +1278,9 @@ static int encode_type1042(rtcm_t *rtcm, int sync) uint32_t sqrtA,e; int i=24,prn,week,toe,toc,i0,OMG0,omg,M0,deln,idot,OMGd,crs,crc; int cus,cuc,cis,cic,af0,af1,af2,tgd1,tgd2; - + trace(3,"encode_type1042: sync=%d\n",sync); - + if (satsys(rtcm->ephsat,&prn)!=SYS_CMP) return 0; eph=rtcm->nav.eph+rtcm->ephsat-1; if (eph->sat!=rtcm->ephsat) return 0; @@ -1307,7 +1307,7 @@ static int encode_type1042(rtcm_t *rtcm, int sync) af2 =ROUND(eph->f2 /P2_66); tgd1 =ROUND(eph->tgd[0]/1E-10); tgd2 =ROUND(eph->tgd[1]/1E-10); - + setbitu(rtcm->buff,i,12,1042 ); i+=12; setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i,13,week ); i+=13; @@ -1347,9 +1347,9 @@ static int encode_type63(rtcm_t *rtcm, int sync) uint32_t sqrtA,e; int i=24,prn,week,toe,toc,i0,OMG0,omg,M0,deln,idot,OMGd,crs,crc; int cus,cuc,cis,cic,af0,af1,af2,tgd1,tgd2; - + trace(3,"encode_type63: sync=%d\n",sync); - + if (satsys(rtcm->ephsat,&prn)!=SYS_CMP) return 0; eph=rtcm->nav.eph+rtcm->ephsat-1; if (eph->sat!=rtcm->ephsat) return 0; @@ -1376,7 +1376,7 @@ static int encode_type63(rtcm_t *rtcm, int sync) af2 =ROUND(eph->f2 /P2_66); tgd1 =ROUND(eph->tgd[0]/1E-10); tgd2 =ROUND(eph->tgd[1]/1E-10); - + setbitu(rtcm->buff,i,12,63 ); i+=12; setbitu(rtcm->buff,i, 6,prn ); i+= 6; setbitu(rtcm->buff,i,13,week ); i+=13; @@ -1416,10 +1416,10 @@ static int encode_ssr_head(int type, rtcm_t *rtcm, int sys, int subtype, { double tow; int i=24,msgno,epoch,week,udi,ns; - + trace(4,"encode_ssr_head: type=%d sys=%d subtype=%d nsat=%d sync=%d iod=%d " "udint=%.0f\n",type,sys,subtype,nsat,sync,iod,udint); - + if (subtype==0) { /* RTCM SSR */ ns=(sys==SYS_QZS)?4:6; switch (sys) { @@ -1435,7 +1435,7 @@ static int encode_ssr_head(int type, rtcm_t *rtcm, int sys, int subtype, return 0; } setbitu(rtcm->buff,i,12,msgno); i+=12; /* message type */ - + if (sys==SYS_GLO) { tow=time2gpst(timeadd(gpst2utc(rtcm->time),10800.0),&week); epoch=ROUND(tow)%86400; @@ -1509,9 +1509,9 @@ static int encode_ssr1(rtcm_t *rtcm, int sys, int subtype, int sync) { double udint=0.0; int i,j,iod=0,nsat,prn,iode,iodcrc,refd=0,np,ni,nj,offp,deph[3],ddeph[3]; - + trace(3,"encode_ssr1: sys=%d subtype=%d sync=%d\n",sys,subtype,sync); - + switch (sys) { case SYS_GPS: np=6; ni= 8; nj= 0; offp= 0; break; case SYS_GLO: np=5; ni= 8; nj= 0; offp= 0; break; @@ -1536,13 +1536,13 @@ static int encode_ssr1(rtcm_t *rtcm, int sys, int subtype, int sync) } /* encode SSR header */ i=encode_ssr_head(1,rtcm,sys,subtype,nsat,sync,iod,udint,refd,0,0); - + for (j=0;jssr[j].update) continue; - + iode=rtcm->ssr[j].iode; /* SBAS/BDS: toe/t0 modulo */ iodcrc=rtcm->ssr[j].iodcrc; /* SBAS/BDS: IOD CRC */ - + if (subtype>0) { /* IGS SSR */ iode&=0xFF; } @@ -1552,7 +1552,7 @@ static int encode_ssr1(rtcm_t *rtcm, int sys, int subtype, int sync) ddeph[0]=ROUND(rtcm->ssr[j].ddeph[0]/1E-6); ddeph[1]=ROUND(rtcm->ssr[j].ddeph[1]/4E-6); ddeph[2]=ROUND(rtcm->ssr[j].ddeph[2]/4E-6); - + setbitu(rtcm->buff,i,np,prn-offp); i+=np; /* satellite ID */ setbitu(rtcm->buff,i,ni,iode ); i+=ni; /* IODE */ setbitu(rtcm->buff,i,nj,iodcrc ); i+=nj; /* IODCRC */ @@ -1571,9 +1571,9 @@ static int encode_ssr2(rtcm_t *rtcm, int sys, int subtype, int sync) { double udint=0.0; int i,j,iod=0,nsat,prn,np,offp,dclk[3]; - + trace(3,"encode_ssr2: sys=%d subtype=%d sync=%d\n",sys,subtype,sync); - + switch (sys) { case SYS_GPS: np=6; offp= 0; break; case SYS_GLO: np=5; offp= 0; break; @@ -1597,14 +1597,14 @@ static int encode_ssr2(rtcm_t *rtcm, int sys, int subtype, int sync) } /* encode SSR header */ i=encode_ssr_head(2,rtcm,sys,subtype,nsat,sync,iod,udint,0,0,0); - + for (j=0;jssr[j].update) continue; - + dclk[0]=ROUND(rtcm->ssr[j].dclk[0]/1E-4); dclk[1]=ROUND(rtcm->ssr[j].dclk[1]/1E-6); dclk[2]=ROUND(rtcm->ssr[j].dclk[2]/2E-8); - + setbitu(rtcm->buff,i,np,prn-offp); i+=np; /* satellite ID */ setbits(rtcm->buff,i,22,dclk[0] ); i+=22; /* delta clock C0 */ setbits(rtcm->buff,i,21,dclk[1] ); i+=21; /* delta clock C1 */ @@ -1620,9 +1620,9 @@ static int encode_ssr3(rtcm_t *rtcm, int sys, int subtype, int sync) double udint=0.0; int i,j,k,iod=0,nsat,prn,nbias,np,offp; int code[MAXCODE],bias[MAXCODE]; - + trace(3,"encode_ssr3: sys=%d subtype=%d sync=%d\n",sys,subtype,sync); - + switch (sys) { case SYS_GPS: np=6; offp= 0; codes=codes_gps; break; case SYS_GLO: np=5; offp= 0; codes=codes_glo; break; @@ -1646,10 +1646,10 @@ static int encode_ssr3(rtcm_t *rtcm, int sys, int subtype, int sync) } /* encode SSR header */ i=encode_ssr_head(3,rtcm,sys,subtype,nsat,sync,iod,udint,0,0,0); - + for (j=0;jssr[j].update) continue; - + for (k=nbias=0;k<32;k++) { if (!codes[k]||rtcm->ssr[j].cbias[codes[k]-1]==0.0) continue; code[nbias]=k; @@ -1657,7 +1657,7 @@ static int encode_ssr3(rtcm_t *rtcm, int sys, int subtype, int sync) } setbitu(rtcm->buff,i,np,prn-offp); i+=np; /* satellite ID */ setbitu(rtcm->buff,i, 5,nbias); i+= 5; /* number of code biases */ - + for (k=0;kbuff,i, 5,code[k]); i+= 5; /* signal indicator */ setbits(rtcm->buff,i,14,bias[k]); i+=14; /* code bias */ @@ -1672,9 +1672,9 @@ static int encode_ssr4(rtcm_t *rtcm, int sys, int subtype, int sync) double udint=0.0; int i,j,iod=0,nsat,prn,iode,iodcrc,refd=0,np,ni,nj,offp; int deph[3],ddeph[3],dclk[3]; - + trace(3,"encode_ssr4: sys=%d subtype=%d sync=%d\n",sys,subtype,sync); - + switch (sys) { case SYS_GPS: np=6; ni= 8; nj= 0; offp= 0; break; case SYS_GLO: np=5; ni= 8; nj= 0; offp= 0; break; @@ -1699,13 +1699,13 @@ static int encode_ssr4(rtcm_t *rtcm, int sys, int subtype, int sync) } /* encode SSR header */ i=encode_ssr_head(4,rtcm,sys,subtype,nsat,sync,iod,udint,refd,0,0); - + for (j=0;jssr[j].update) continue; - + iode=rtcm->ssr[j].iode; iodcrc=rtcm->ssr[j].iodcrc; - + if (subtype>0) { /* IGS SSR */ iode&=0xFF; } @@ -1718,7 +1718,7 @@ static int encode_ssr4(rtcm_t *rtcm, int sys, int subtype, int sync) dclk [0]=ROUND(rtcm->ssr[j].dclk [0]/1E-4); dclk [1]=ROUND(rtcm->ssr[j].dclk [1]/1E-6); dclk [2]=ROUND(rtcm->ssr[j].dclk [2]/2E-8); - + setbitu(rtcm->buff,i,np,prn-offp); i+=np; /* satellite ID */ setbitu(rtcm->buff,i,ni,iode ); i+=ni; /* IODE */ setbitu(rtcm->buff,i,nj,iodcrc ); i+=nj; /* IODCRC */ @@ -1740,9 +1740,9 @@ static int encode_ssr5(rtcm_t *rtcm, int sys, int subtype, int sync) { double udint=0.0; int i,j,nsat,iod=0,prn,ura,np,offp; - + trace(3,"encode_ssr5: sys=%d subtype=%d sync=%d\n",sys,subtype,sync); - + switch (sys) { case SYS_GPS: np=6; offp= 0; break; case SYS_GLO: np=5; offp= 0; break; @@ -1766,10 +1766,10 @@ static int encode_ssr5(rtcm_t *rtcm, int sys, int subtype, int sync) } /* encode ssr header */ i=encode_ssr_head(5,rtcm,sys,subtype,nsat,sync,iod,udint,0,0,0); - + for (j=0;jssr[j].update) continue; - + ura=rtcm->ssr[j].ura; setbitu(rtcm->buff,i,np,prn-offp); i+=np; /* satellite id */ setbitu(rtcm->buff,i, 6,ura ); i+= 6; /* ssr ura */ @@ -1782,9 +1782,9 @@ static int encode_ssr6(rtcm_t *rtcm, int sys, int subtype, int sync) { double udint=0.0; int i,j,nsat,iod=0,prn,hrclk,np,offp; - + trace(3,"encode_ssr6: sys=%d subtype=%d sync=%d\n",sys,subtype,sync); - + switch (sys) { case SYS_GPS: np=6; offp= 0; break; case SYS_GLO: np=5; offp= 0; break; @@ -1808,12 +1808,12 @@ static int encode_ssr6(rtcm_t *rtcm, int sys, int subtype, int sync) } /* encode SSR header */ i=encode_ssr_head(6,rtcm,sys,subtype,nsat,sync,iod,udint,0,0,0); - + for (j=0;jssr[j].update) continue; - + hrclk=ROUND(rtcm->ssr[j].hrclk/1E-4); - + setbitu(rtcm->buff,i,np,prn-offp); i+=np; /* satellite ID */ setbits(rtcm->buff,i,22,hrclk ); i+=22; /* high rate clock corr */ } @@ -1827,9 +1827,9 @@ static int encode_ssr7(rtcm_t *rtcm, int sys, int subtype, int sync) double udint=0.0; int i,j,k,iod=0,nsat,prn,nbias,np,offp; int code[MAXCODE],pbias[MAXCODE],stdpb[MAXCODE],yaw_ang,yaw_rate; - + trace(3,"encode_ssr7: sys=%d subtype=%d sync=%d\n",sys,subtype,sync); - + switch (sys) { case SYS_GPS: np=6; offp= 0; codes=codes_gps; break; case SYS_GLO: np=5; offp= 0; codes=codes_glo; break; @@ -1853,10 +1853,10 @@ static int encode_ssr7(rtcm_t *rtcm, int sys, int subtype, int sync) } /* encode SSR header */ i=encode_ssr_head(7,rtcm,sys,subtype,nsat,sync,iod,udint,0,0,0); - + for (j=0;jssr[j].update) continue; - + for (k=nbias=0;k<32;k++) { if (!codes[k]||rtcm->ssr[j].pbias[codes[k]-1]==0.0) continue; code[nbias]=k; @@ -1869,7 +1869,7 @@ static int encode_ssr7(rtcm_t *rtcm, int sys, int subtype, int sync) setbitu(rtcm->buff,i, 5,nbias); i+= 5; /* number of code biases */ setbitu(rtcm->buff,i, 9,yaw_ang); i+= 9; /* yaw angle */ setbits(rtcm->buff,i, 8,yaw_rate); i+= 8; /* yaw rate */ - + for (k=0;kbuff,i, 5,code[k] ); i+= 5; /* signal indicator */ setbitu(rtcm->buff,i, 1,0 ); i+= 1; /* integer-indicator */ @@ -1888,12 +1888,12 @@ static int encode_ssr7(rtcm_t *rtcm, int sys, int subtype, int sync) static int to_satid(int sys, int sat) { int prn; - + if (satsys(sat,&prn)!=sys) return 0; - + if (sys==SYS_QZS) prn-=MINPRNQZS-1; else if (sys==SYS_SBS) prn-=MINPRNSBS-1; - + return prn; } /* observation code to MSM signal ID -----------------------------------------*/ @@ -1902,7 +1902,7 @@ static int to_sigid(int sys, uint8_t code) const char **msm_sig; char *sig; int i; - + /* signal conversion for undefined signal by rtcm */ if (sys==SYS_GPS) { if (code==CODE_L1Y) code=CODE_L1P; @@ -1914,7 +1914,7 @@ static int to_sigid(int sys, uint8_t code) else if (code==CODE_L2N) code=CODE_L2P; } if (!*(sig=code2obs(code))) return 0; - + switch (sys) { case SYS_GPS: msm_sig=msm_sig_gps; break; case SYS_GLO: msm_sig=msm_sig_glo; break; @@ -1936,16 +1936,16 @@ static void gen_msm_index(const rtcm_t *rtcm, int sys, int *nsat, int *nsig, uint8_t *cell_ind) { int i,j,sat,sig,cell; - + *nsat=*nsig=*ncell=0; - + /* generate satellite and signal index */ for (i=0;iobs.n;i++) { if (!(sat=to_satid(sys,rtcm->obs.data[i].sat))) continue; - + for (j=0;jobs.data[i].code[j]))) continue; - + sat_ind[sat-1]=sig_ind[sig-1]=1; } } @@ -1958,10 +1958,10 @@ static void gen_msm_index(const rtcm_t *rtcm, int sys, int *nsat, int *nsig, /* generate cell index */ for (i=0;iobs.n;i++) { if (!(sat=to_satid(sys,rtcm->obs.data[i].sat))) continue; - + for (j=0;jobs.data[i].code[j]))) continue; - + cell=sig_ind[sig-1]-1+(sat_ind[sat-1]-1)*(*nsig); cell_ind[cell]=1; } @@ -1977,20 +1977,20 @@ static void gen_msm_sat(rtcm_t *rtcm, int sys, int nsat, const uint8_t *sat_ind, obsd_t *data; double freq; int i,j,k,sat,sig,fcn; - + for (i=0;i<64;i++) rrng[i]=rrate[i]=0.0; - + for (i=0;iobs.n;i++) { data=rtcm->obs.data+i; fcn=fcn_glo(data->sat,rtcm); /* fcn+7 */ - + if (!(sat=to_satid(sys,data->sat))) continue; - + for (j=0;jcode[j]))) continue; k=sat_ind[sat-1]-1; freq=code2freq(sys,data->code[j],fcn-7); - + /* rough range (ms) and rough phase-range-rate (m/s) */ if (rrng[k]==0.0&&data->P[j]!=0.0) { rrng[k]=ROUND( data->P[j]/RANGE_MS/P2_10)*RANGE_MS*P2_10; @@ -2013,7 +2013,7 @@ static void gen_msm_sig(rtcm_t *rtcm, int sys, int nsat, int nsig, int ncell, obsd_t *data; double freq,lambda,psrng_s,phrng_s,rate_s,lt; int i,j,k,sat,sig,fcn,cell,LLI; - + for (i=0;iobs.n;i++) { data=rtcm->obs.data+i; fcn=fcn_glo(data->sat,rtcm); /* fcn+7 */ - + if (!(sat=to_satid(sys,data->sat))) continue; - + for (j=0;jcode[j]))) continue; - + k=sat_ind[sat-1]-1; if ((cell=cell_ind[sig_ind[sig-1]-1+k*nsig])>=64) continue; - + freq=code2freq(sys,data->code[j],fcn-7); lambda=freq==0.0?0.0:CLIGHT/freq; psrng_s=data->P[j]==0.0?0.0:data->P[j]-rrng[k]; phrng_s=data->L[j]==0.0||lambda<=0.0?0.0: data->L[j]*lambda-rrng [k]; rate_s =data->D[j]==0.0||lambda<=0.0?0.0:-data->D[j]*lambda-rrate[k]; - + /* subtract phase - psudorange integer cycle offset */ LLI=data->LLI[j]; if ((LLI&1)||fabs(phrng_s-rtcm->cp[data->sat-1][j])>1171.0) { @@ -2044,9 +2044,9 @@ static void gen_msm_sig(rtcm_t *rtcm, int sys, int nsat, int nsig, int ncell, LLI|=1; } phrng_s-=rtcm->cp[data->sat-1][j]; - + lt=locktime_d(data->time,rtcm->lltime[data->sat-1]+j,LLI); - + if (psrng&&psrng_s!=0.0) psrng[cell-1]=psrng_s; if (phrng&&phrng_s!=0.0) phrng[cell-1]=phrng_s; if (rate &&rate_s !=0.0) rate [cell-1]=rate_s; @@ -2067,7 +2067,7 @@ static int encode_msm_head(int type, rtcm_t *rtcm, int sys, int sync, int *nsat, uint8_t sat_ind[64]={0},sig_ind[32]={0},cell_ind[32*64]={0}; uint32_t dow,epoch; int i=24,j,nsig=0; - + switch (sys) { case SYS_GPS: type+=1070; break; case SYS_GLO: type+=1080; break; @@ -2080,7 +2080,7 @@ static int encode_msm_head(int type, rtcm_t *rtcm, int sys, int sync, int *nsat, } /* generate msm satellite, signal and cell index */ gen_msm_index(rtcm,sys,nsat,&nsig,ncell,sat_ind,sig_ind,cell_ind); - + if (sys==SYS_GLO) { /* GLONASS time (dow + tod-ms) */ tow=time2gpst(timeadd(gpst2utc(rtcm->time),10800.0),NULL); @@ -2106,7 +2106,7 @@ static int encode_msm_head(int type, rtcm_t *rtcm, int sys, int sync, int *nsat, setbitu(rtcm->buff,i, 2,0 ); i+= 2; /* external clock indicator */ setbitu(rtcm->buff,i, 1,0 ); i+= 1; /* smoothing indicator */ setbitu(rtcm->buff,i, 3,0 ); i+= 3; /* smoothing interval */ - + /* satellite mask */ for (j=0;j<64;j++) { setbitu(rtcm->buff,i,1,sat_ind[j]?1:0); i+=1; @@ -2121,11 +2121,11 @@ static int encode_msm_head(int type, rtcm_t *rtcm, int sys, int sync, int *nsat, } /* generate msm satellite data fields */ gen_msm_sat(rtcm,sys,*nsat,sat_ind,rrng,rrate,info); - + /* generate msm signal data fields */ gen_msm_sig(rtcm,sys,*nsat,nsig,*ncell,sat_ind,sig_ind,cell_ind,rrng,rrate, psrng,phrng,rate,lock,half,cnr); - + return i; } /* encode rough range integer ms ---------------------------------------------*/ @@ -2134,7 +2134,7 @@ static int encode_msm_int_rrng(rtcm_t *rtcm, int i, const double *rrng, { uint32_t int_ms; int j; - + for (j=0;jRANGE_MS*255.0) { mod_ms=0; @@ -2174,7 +2174,7 @@ static int encode_msm_mod_rrng(rtcm_t *rtcm, int i, const double *rrng, static int encode_msm_info(rtcm_t *rtcm, int i, const uint8_t *info, int nsat) { int j; - + for (j=0;jbuff,i,4,info[j]); i+=4; } @@ -2184,7 +2184,7 @@ static int encode_msm_info(rtcm_t *rtcm, int i, const uint8_t *info, int nsat) static int encode_msm_rrate(rtcm_t *rtcm, int i, const double *rrate, int nsat) { int j,rrate_val; - + for (j=0;j8191.0) { char tstr[40]; @@ -2203,7 +2203,7 @@ static int encode_msm_rrate(rtcm_t *rtcm, int i, const double *rrate, int nsat) static int encode_msm_psrng(rtcm_t *rtcm, int i, const double *psrng, int ncell) { int j,psrng_val; - + for (j=0;jbuff,i,4,lock_val); i+=4; @@ -2305,7 +2305,7 @@ static int encode_msm_lock_ex(rtcm_t *rtcm, int i, const double *lock, int ncell) { int j,lock_val; - + for (j=0;jbuff,i,10,lock_val); i+=10; @@ -2317,7 +2317,7 @@ static int encode_msm_half_amb(rtcm_t *rtcm, int i, const uint8_t *half, int ncell) { int j; - + for (j=0;jbuff,i,1,half[j]); i+=1; } @@ -2327,7 +2327,7 @@ static int encode_msm_half_amb(rtcm_t *rtcm, int i, const uint8_t *half, static int encode_msm_cnr(rtcm_t *rtcm, int i, const float *cnr, int ncell) { int j,cnr_val; - + for (j=0;jbuff,i,6,cnr_val); i+=6; @@ -2338,7 +2338,7 @@ static int encode_msm_cnr(rtcm_t *rtcm, int i, const float *cnr, int ncell) static int encode_msm_cnr_ex(rtcm_t *rtcm, int i, const float *cnr, int ncell) { int j,cnr_val; - + for (j=0;jbuff,i,10,cnr_val); i+=10; @@ -2349,7 +2349,7 @@ static int encode_msm_cnr_ex(rtcm_t *rtcm, int i, const float *cnr, int ncell) static int encode_msm_rate(rtcm_t *rtcm, int i, const double *rate, int ncell) { int j,rate_val; - + for (j=0;jnbit=i; return 1; } @@ -2395,9 +2395,9 @@ static int encode_msm2(rtcm_t *rtcm, int sys, int sync) double rrng[64],rrate[64],phrng[64],lock[64]; uint8_t half[64]; int i,nsat,ncell; - + trace(3,"encode_msm2: sys=%d sync=%d\n",sys,sync); - + /* encode msm header */ if (!(i=encode_msm_head(2,rtcm,sys,sync,&nsat,&ncell,rrng,rrate,NULL,NULL, phrng,NULL,lock,half,NULL))) { @@ -2405,12 +2405,12 @@ static int encode_msm2(rtcm_t *rtcm, int sys, int sync) } /* encode msm satellite data */ i=encode_msm_mod_rrng(rtcm,i,rrng ,nsat ); /* rough range modulo 1 ms */ - + /* encode msm signal data */ i=encode_msm_phrng (rtcm,i,phrng,ncell); /* fine phase-range */ i=encode_msm_lock (rtcm,i,lock ,ncell); /* lock-time indicator */ i=encode_msm_half_amb(rtcm,i,half ,ncell); /* half-cycle-amb indicator */ - + rtcm->nbit=i; return 1; } @@ -2420,9 +2420,9 @@ static int encode_msm3(rtcm_t *rtcm, int sys, int sync) double rrng[64],rrate[64],psrng[64],phrng[64],lock[64]; uint8_t half[64]; int i,nsat,ncell; - + trace(3,"encode_msm3: sys=%d sync=%d\n",sys,sync); - + /* encode msm header */ if (!(i=encode_msm_head(3,rtcm,sys,sync,&nsat,&ncell,rrng,rrate,NULL,psrng, phrng,NULL,lock,half,NULL))) { @@ -2430,13 +2430,13 @@ static int encode_msm3(rtcm_t *rtcm, int sys, int sync) } /* encode msm satellite data */ i=encode_msm_mod_rrng(rtcm,i,rrng ,nsat ); /* rough range modulo 1 ms */ - + /* encode msm signal data */ i=encode_msm_psrng (rtcm,i,psrng,ncell); /* fine pseudorange */ i=encode_msm_phrng (rtcm,i,phrng,ncell); /* fine phase-range */ i=encode_msm_lock (rtcm,i,lock ,ncell); /* lock-time indicator */ i=encode_msm_half_amb(rtcm,i,half ,ncell); /* half-cycle-amb indicator */ - + rtcm->nbit=i; return 1; } @@ -2447,9 +2447,9 @@ static int encode_msm4(rtcm_t *rtcm, int sys, int sync) float cnr[64]; uint8_t half[64]; int i,nsat,ncell; - + trace(3,"encode_msm4: sys=%d sync=%d\n",sys,sync); - + /* encode msm header */ if (!(i=encode_msm_head(4,rtcm,sys,sync,&nsat,&ncell,rrng,rrate,NULL,psrng, phrng,NULL,lock,half,cnr))) { @@ -2458,7 +2458,7 @@ static int encode_msm4(rtcm_t *rtcm, int sys, int sync) /* encode msm satellite data */ i=encode_msm_int_rrng(rtcm,i,rrng ,nsat ); /* rough range integer ms */ i=encode_msm_mod_rrng(rtcm,i,rrng ,nsat ); /* rough range modulo 1 ms */ - + /* encode msm signal data */ i=encode_msm_psrng (rtcm,i,psrng,ncell); /* fine pseudorange */ i=encode_msm_phrng (rtcm,i,phrng,ncell); /* fine phase-range */ @@ -2475,9 +2475,9 @@ static int encode_msm5(rtcm_t *rtcm, int sys, int sync) float cnr[64]; uint8_t info[64],half[64]; int i,nsat,ncell; - + trace(3,"encode_msm5: sys=%d sync=%d\n",sys,sync); - + /* encode msm header */ if (!(i=encode_msm_head(5,rtcm,sys,sync,&nsat,&ncell,rrng,rrate,info,psrng, phrng,rate,lock,half,cnr))) { @@ -2488,7 +2488,7 @@ static int encode_msm5(rtcm_t *rtcm, int sys, int sync) i=encode_msm_info (rtcm,i,info ,nsat ); /* extended satellite info */ i=encode_msm_mod_rrng(rtcm,i,rrng ,nsat ); /* rough range modulo 1 ms */ i=encode_msm_rrate (rtcm,i,rrate,nsat ); /* rough phase-range-rate */ - + /* encode msm signal data */ i=encode_msm_psrng (rtcm,i,psrng,ncell); /* fine pseudorange */ i=encode_msm_phrng (rtcm,i,phrng,ncell); /* fine phase-range */ @@ -2506,9 +2506,9 @@ static int encode_msm6(rtcm_t *rtcm, int sys, int sync) float cnr[64]; uint8_t half[64]; int i,nsat,ncell; - + trace(3,"encode_msm6: sys=%d sync=%d\n",sys,sync); - + /* encode msm header */ if (!(i=encode_msm_head(6,rtcm,sys,sync,&nsat,&ncell,rrng,rrate,NULL,psrng, phrng,NULL,lock,half,cnr))) { @@ -2517,7 +2517,7 @@ static int encode_msm6(rtcm_t *rtcm, int sys, int sync) /* encode msm satellite data */ i=encode_msm_int_rrng(rtcm,i,rrng ,nsat ); /* rough range integer ms */ i=encode_msm_mod_rrng(rtcm,i,rrng ,nsat ); /* rough range modulo 1 ms */ - + /* encode msm signal data */ i=encode_msm_psrng_ex(rtcm,i,psrng,ncell); /* fine pseudorange ext */ i=encode_msm_phrng_ex(rtcm,i,phrng,ncell); /* fine phase-range ext */ @@ -2534,9 +2534,9 @@ static int encode_msm7(rtcm_t *rtcm, int sys, int sync) float cnr[64]; uint8_t info[64],half[64]; int i,nsat,ncell; - + trace(3,"encode_msm7: sys=%d sync=%d\n",sys,sync); - + /* encode msm header */ if (!(i=encode_msm_head(7,rtcm,sys,sync,&nsat,&ncell,rrng,rrate,info,psrng, phrng,rate,lock,half,cnr))) { @@ -2547,7 +2547,7 @@ static int encode_msm7(rtcm_t *rtcm, int sys, int sync) i=encode_msm_info (rtcm,i,info ,nsat ); /* extended satellite info */ i=encode_msm_mod_rrng(rtcm,i,rrng ,nsat ); /* rough range modulo 1 ms */ i=encode_msm_rrate (rtcm,i,rrate,nsat ); /* rough phase-range-rate */ - + /* encode msm signal data */ i=encode_msm_psrng_ex(rtcm,i,psrng,ncell); /* fine pseudorange ext */ i=encode_msm_phrng_ex(rtcm,i,phrng,ncell); /* fine phase-range ext */ @@ -2562,11 +2562,11 @@ static int encode_msm7(rtcm_t *rtcm, int sys, int sync) static int encode_type1230(rtcm_t *rtcm, int sync) { int i=24,j,align,mask=15,bias[4]; - + trace(3,"encode_type1230: sync=%d\n",sync); - + align=rtcm->sta.glo_cp_align; - + for (j=0;j<4;j++) { bias[j]=ROUND(rtcm->sta.glo_cp_bias[j]/0.02); if (bias[j]<=-32768||bias[j]>32767) { @@ -2645,9 +2645,9 @@ static int encode_type4076(rtcm_t *rtcm, int subtype, int sync) extern int encode_rtcm3(rtcm_t *rtcm, int type, int subtype, int sync) { int ret=0; - + trace(3,"encode_rtcm3: type=%d subtype=%d sync=%d\n",type,subtype,sync); - + switch (type) { case 1001: ret=encode_type1001(rtcm,sync); break; case 1002: ret=encode_type1002(rtcm,sync); break; diff --git a/util/testeph/dumpssr.c b/util/testeph/dumpssr.c index 4630ac71b..fc4891af8 100644 --- a/util/testeph/dumpssr.c +++ b/util/testeph/dumpssr.c @@ -10,9 +10,9 @@ static void printhead(int topt, int mopt) { int i; - + printf("%% %s SAT ",topt?" DAY TIME ":" GPST "); - + if (mopt&1) { printf(" UDI IOD URA REF "); } @@ -33,7 +33,7 @@ static void printssrmsg(int sat, const ssr_t *ssr, int topt, int mopt) double tow; int week; char tstr[40],id[8]; - + if (topt) { time2str(ssr->t0,tstr,0); printf("%s ",tstr); @@ -44,7 +44,7 @@ static void printssrmsg(int sat, const ssr_t *ssr, int topt, int mopt) } satno2id(sat,id); printf("%4s ",id); - + if (mopt&1) { printf("%4.0f %3d %3d %3d ",ssr->udint,ssr->iode,ssr->ura,ssr->refd); } @@ -70,17 +70,17 @@ static void dumpssrmsg(FILE *fp, int sat, int topt, int mopt) static rtcm_t rtcm; static gtime_t t0[MAXSAT]={{0}}; int i,stat; - + init_rtcm(&rtcm); - + while ((stat=input_rtcm3f(&rtcm,fp))>=0) { - + if (stat!=10) continue; /* ssr message */ - + for (i=0;i Date: Mon, 21 Jul 2025 10:44:57 +0200 Subject: [PATCH 2/6] feat: add coding/decoding of IGS IONO SSR (4076/201 RTCM message) A tool iono2rtcm (in utils) to encode BNC VTEC blocks to RTCM 4076/201 message has been added (including tests) The utils tool dumpssr includes also the dumping of IONO SSR corrections in BNC SSR text format (see section 2.8 at https://software.rtcm-ntrip.org/export/HEAD/ntrip/trunk/BNC/src/bnchelp.html (cherry picked from commit 942b034ac226c97be78cb1b5426d8baf4702bffe) --- src/rtcm3.c | 81 ++++++++ src/rtcm3e.c | 70 +++++++ src/rtklib.h | 39 ++++ util/data/IONO00IGS1_BNC.rtcm3 | Bin 0 -> 531 bytes util/data/IONO00IGS1_BNC.txt | 34 ++++ util/data/IONO00IGS1_ICD_example.rtcm3 | Bin 0 -> 47 bytes util/data/IONO00IGS1_ICD_example.txt | 10 + util/iono2rtcm/main.c | 252 +++++++++++++++++++++++++ util/iono2rtcm/makefile | 60 ++++++ util/testeph/dumpssr.c | 148 ++++++++++++++- util/testeph/makefile | 44 ++++- 11 files changed, 726 insertions(+), 12 deletions(-) create mode 100644 util/data/IONO00IGS1_BNC.rtcm3 create mode 100644 util/data/IONO00IGS1_BNC.txt create mode 100644 util/data/IONO00IGS1_ICD_example.rtcm3 create mode 100644 util/data/IONO00IGS1_ICD_example.txt create mode 100644 util/iono2rtcm/main.c create mode 100644 util/iono2rtcm/makefile diff --git a/src/rtcm3.c b/src/rtcm3.c index 62dddc576..ce7e7721c 100644 --- a/src/rtcm3.c +++ b/src/rtcm3.c @@ -1953,6 +1953,86 @@ static int decode_ssr7(rtcm_t *rtcm, int sys, int subtype) } return 20; } + +/* decode IGS SSR message subtype 201: ionospheric spherical harmonics -------*/ +/* Follows ICD description found in https://files.igs.org/pub/data/format/igs_ssr_v1.pdf */ +static int decode_ssr_iono_sph_harmonics(rtcm_t *rtcm, int sys, int subtype) { + + int ret = -1; + int pos = 24; /* Start position after header */ + int DF002, IDF001, IDF002, IDF003, IDF004, IDF005, IDF007, IDF008, IDF009; + int IDF041, IDF035; + + /* header, Table 24*/ + + /* RTCM Message Number ust be 4076 */ + DF002 = getbitu(rtcm->buff, pos, 12); pos += 12; + if (DF002 != 4076) { + trace(2, "rtcm3 %d: invalid message number: %d\n", DF002, 4076); + goto end; + } + + IDF001 = getbitu(rtcm->buff, pos, IDF001_NBITS); pos += IDF001_NBITS; + + IDF002 = getbitu(rtcm->buff, pos, IDF002_NBITS); pos += IDF002_NBITS; + if (IDF002 != subtype) { + trace(2, "rtcm3 %d: invalid IDF002, got: %d, expected: %d\n", DF002, IDF002, 201); + goto end; + } + + IDF003 = getbitu(rtcm->buff, pos, IDF003_NBITS); pos += IDF003_NBITS; /* SSR Epoch Time 1s */ + IDF004 = getbitu(rtcm->buff, pos, IDF004_NBITS); pos += IDF004_NBITS; /* SSR Update Interval */ + IDF005 = getbitu(rtcm->buff, pos, IDF005_NBITS); pos += IDF005_NBITS; /* SSR Multiple Message Indicator */ + IDF007 = getbitu(rtcm->buff, pos, IDF007_NBITS); pos += IDF007_NBITS; /* IOD SSR */ + IDF008 = getbitu(rtcm->buff, pos, IDF008_NBITS); pos += IDF008_NBITS; /* SSR Provider ID */ + IDF009 = getbitu(rtcm->buff, pos, IDF009_NBITS); pos += IDF009_NBITS; /* SSR Solution ID */ + IDF041 = getbitu(rtcm->buff, pos, IDF041_NBITS); pos += IDF041_NBITS; /* VTEC Quality Indicator */ + IDF035 = getbitu(rtcm->buff, pos, IDF035_NBITS); pos += IDF035_NBITS; /* Number of Ionospheric Layers */ + + adjweek(rtcm, IDF003); + rtcm->ssr_ion.update_interval_s = IDF004; + + rtcm->ssr_ion.n_layers = IDF035 + 1; + + for (int i = 0; i < rtcm->ssr_ion.n_layers; i++) { + + ionlayer_sphharm_t* layer = &rtcm->ssr_ion.ionlayers_sph_harm[i]; + + /* Table 25: Header Part for an Ionospheric Layer IM201*/ + + int IDF036 = getbitu(rtcm->buff, pos, IDF036_NBITS); pos += IDF036_NBITS; /* Ionospheric Layer Height */ + int IDF037 = getbitu(rtcm->buff, pos, IDF037_NBITS); pos += IDF037_NBITS; /* Ionospheric Layer degree (N) */ + int IDF038 = getbitu(rtcm->buff, pos, IDF038_NBITS); pos += IDF038_NBITS; /* Ionospheric Layer order (M) */ + int n = IDF037 + 1; + int m = IDF038 + 1; + layer->height_km = IDF036 * IDF036_SCALE; + layer->n_degree = n; + layer->m_order = m; + + trace(3, "Layer %d: Height=%d Degree=%d Order=%d\n", + i+1, layer->height_km, layer->n_degree, layer->m_order); + + /* Table 26: Spherical Harmonics Cosine Coefficients */ + int n_cos_coeffs = (n + 1) * (n + 2) / 2 - (n - m) * (n - m + 1) / 2; + for (int j = 0; j < n_cos_coeffs; j++) { + int IDF039 = getbits(rtcm->buff, pos, IDF039_NBITS); pos += IDF039_NBITS; /* Ionospheric Layer order (M) */ + layer->cos_coeff[j] = IDF039 * IDF039_SCALE; + } + + /* Table 27: Spherical Harmonics Sine Coefficients */ + int n_sin_coeffs = n_cos_coeffs - (n + 1); + for (int j = 0; j < n_sin_coeffs; j++) { + int IDF040 = getbits(rtcm->buff, pos, IDF040_NBITS); pos += IDF040_NBITS; /* Ionospheric Layer order (M) */ + layer->sin_coeff[j] = IDF040 * IDF040_SCALE; + } + } + + ret = 10; /* default return value for SSR message */ +end: + return ret; +} + + /* get signal index ----------------------------------------------------------*/ static void sigindex(int sys, const uint8_t *code, int n, const char *opt, int *idx) @@ -2549,6 +2629,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_ssr_iono_sph_harmonics(rtcm,SYS_NONE,subtype); } trace(2,"rtcm3 4076: unsupported message subtype=%d\n",subtype); return 0; diff --git a/src/rtcm3e.c b/src/rtcm3e.c index 4841351ea..555c73361 100644 --- a/src/rtcm3e.c +++ b/src/rtcm3e.c @@ -1884,6 +1884,75 @@ static int encode_ssr7(rtcm_t *rtcm, int sys, int subtype, int sync) rtcm->nbit=i; return 1; } + +static int encode_ssr_iono_sph_harmonics(rtcm_t *rtcm, int sys, int subtype, int sync) { + + static const int TYPE = 4076; /* IGS SSR Proprietary message type */ + static const int SUBTYPE = 201; /* Spherical harmonics subtype */ + double tow; + int offset=24, week, epoch; + int ret = 0; // if error, return 0 + + trace(4,"encode_head: type=%d subtype=%d sync=%d\n", TYPE, SUBTYPE); + + const ssr_ion_t* ion = &rtcm->ssr_ion; + + /* Table 24 of ICD*/ + + setbitu(rtcm->buff, offset, 12, TYPE); offset+=12; /* DF002 */ + setbitu(rtcm->buff, offset, IDF001_NBITS, 1); offset+=IDF001_NBITS; /* IDF001 IGS SSR Version */ + setbitu(rtcm->buff, offset, IDF002_NBITS, SUBTYPE); offset+=IDF002_NBITS; /* IDF002 Subtype */ + + /* IDF003 SSR Epoch Time */ + tow=time2gpst(rtcm->time,&week); + epoch=ROUND(tow/0.001); + setbitu(rtcm->buff, offset, IDF003_NBITS, epoch); offset+=IDF003_NBITS; + + setbitu(rtcm->buff, offset, IDF004_NBITS, ion->update_interval_s); offset+=IDF004_NBITS; /* IDF004 SSR Update interval */ + setbitu(rtcm->buff, offset, IDF005_NBITS, 0); offset+=IDF005_NBITS; /* IDF005 Multiple message indicator */ + setbitu(rtcm->buff, offset, IDF007_NBITS, 0); offset+=IDF007_NBITS; /* IDF007 IOD SSR */ + setbitu(rtcm->buff, offset, IDF008_NBITS, 0); offset+=IDF008_NBITS; /* IDF008 SSR Provider ID */ + setbitu(rtcm->buff, offset, IDF009_NBITS, 0); offset+=IDF009_NBITS; /* IDF009 SSR Solution ID */ + setbitu(rtcm->buff, offset, IDF041_NBITS, 0); offset+=IDF041_NBITS; /* IDF041 VTEC Quality Indicator */ + + setbitu(rtcm->buff, offset, IDF035_NBITS, ion->n_layers - 1); offset+=IDF035_NBITS; /* IDF035 Number of ionospheric layers */ + + for (int i = 0; i < ion->n_layers; i++) { + + const ionlayer_sphharm_t *layer = &ion->ionlayers_sph_harm[i]; + + /* Table 25 of ICD */ + setbitu(rtcm->buff, offset, IDF036_NBITS, (uint8_t)(layer->height_km / IDF036_SCALE)); offset+=IDF036_NBITS; /* IDF036 Layer height */ + setbitu(rtcm->buff, offset, IDF037_NBITS, layer->n_degree - 1); offset+=IDF037_NBITS; /* IDF037 Degree */ + setbitu(rtcm->buff, offset, IDF038_NBITS, layer->m_order - 1); offset+=4; /* IDF037 Order */ + + /* Table 26 cosine terms */ + int n = layer->n_degree; + int m = layer->m_order; + int n_cos_coeffs = (n + 1) * (n + 2) / 2 - (n - m) * (n - m + 1) / 2; + for (int j = 0; j < n_cos_coeffs; j++) { + float coeff = layer->cos_coeff[j]; + setbits(rtcm->buff, offset, IDF039_NBITS, (int16_t)(coeff / IDF039_SCALE)); offset+=IDF039_NBITS; /* IDF039 cos coeff */ + } + + /* Table 27 sine terms */ + int n_sin_coeffs = n_cos_coeffs - (n + 1); + for (int j = 0; j < n_sin_coeffs; j++) { + float coeff = layer->sin_coeff[j]; + setbits(rtcm->buff, offset, IDF040_NBITS, (int16_t)(coeff / IDF040_SCALE)); offset+=IDF040_NBITS; /* IDF040 sin coeff */ + } + } + + rtcm->nbit=offset; + + trace(3,"encode_ssr_iono_sph_harmonics: sys=%d subtype=%d sync=%d\n",sys,subtype,sync); + + ret = 1; // if successful, return 1 +end: + return ret; + +} + /* satellite no to MSM satellite ID ------------------------------------------*/ static int to_satid(int sys, int sat) { @@ -2637,6 +2706,7 @@ static int encode_type4076(rtcm_t *rtcm, int subtype, int sync) case 125: return encode_ssr3(rtcm,SYS_SBS,subtype,sync); case 126: return encode_ssr7(rtcm,SYS_SBS,subtype,sync); case 127: return encode_ssr5(rtcm,SYS_SBS,subtype,sync); + case 201: return encode_ssr_iono_sph_harmonics(rtcm,SYS_NONE,subtype,sync); } trace(2,"rtcm3 4076: unsupported message subtype=%d\n",subtype); return 0; diff --git a/src/rtklib.h b/src/rtklib.h index c1ff8c16f..8c4859e7b 100644 --- a/src/rtklib.h +++ b/src/rtklib.h @@ -292,6 +292,9 @@ extern "C" { #define MAXRCVCMD 4096 /* max length of receiver commands */ #define MAX_CODE_BIASES 3 /* max # of different code biases per freq */ #define MAX_CODE_BIAS_FREQS 2 /* max # of freqs supported for code biases */ +#define MAX_ION_LAYERS 4 /* max number of ionosphere layers */ +#define MAX_SPH_HARM_COS_COEFF 156 /* max number of cos spherical harmonics coeffs ((n + 1) * (n + 2) / 2 - (n - m) * (n - m + 1) / 2) n=16, m=16 */ +#define MAX_SPH_HARM_SIN_COEFF 136 /* max number of sin spherical harmonics coeffs ((n + 1) * (n + 2) / 2 - (n - m) * (n - m + 1) / 2) - (n+1) n=16, m=16*/ #define RNX2VER 2.10 /* RINEX ver.2 default output version */ #define RNX3VER 3.00 /* RINEX ver.3 default output version */ @@ -552,6 +555,27 @@ extern "C" { #define P2_50 8.881784197001252E-16 /* 2^-50 */ #define P2_55 2.775557561562891E-17 /* 2^-55 */ +#define IDF036_SCALE 10.0f /* SSR IONO Height scale factor */ +#define IDF039_SCALE 0.005f /* Scale factor for the sph. harmonic cosine coefficients */ +#define IDF040_SCALE 0.005f /* Scale factor for the sph. harmonic sine coefficients */ + +#define IDF001_NBITS 3 +#define IDF002_NBITS 8 +#define IDF003_NBITS 20 +#define IDF004_NBITS 4 +#define IDF005_NBITS 1 +#define IDF007_NBITS 4 +#define IDF008_NBITS 16 +#define IDF009_NBITS 4 +#define IDF041_NBITS 9 +#define IDF035_NBITS 2 +#define IDF036_NBITS 8 +#define IDF037_NBITS 4 +#define IDF038_NBITS 4 +#define IDF039_NBITS 16 +#define IDF040_NBITS 16 + + #ifdef WIN32 #define rtklib_thread_t HANDLE #define rtklib_lock_t CRITICAL_SECTION @@ -821,6 +845,20 @@ typedef struct { /* DGPS/GNSS correction type */ double udre; /* UDRE */ } dgps_t; +typedef struct { /* ionosphere spherical harmonics type */ + float height_km; /* height of spherical harmonics */ + int n_degree; /* degree of spherical harmonics */ + int m_order; /* order of spherical harmonics */ + double cos_coeff[MAX_SPH_HARM_COS_COEFF]; /* cosine coefficients, ordered as C00,C10,C20,C30;C11,C21,C31;C22,C32 for N=3 and M=2*/ + double sin_coeff[MAX_SPH_HARM_SIN_COEFF]; /* sine coefficients, ordered as S11,S21,S31;S22,S32 for N=3 and M=2*/ +} ionlayer_sphharm_t; + +typedef struct { + int update_interval_s; /* update interval*/ + int n_layers; /* number of layers */ + ionlayer_sphharm_t ionlayers_sph_harm[MAX_ION_LAYERS]; /* ionosphere spherical harmonics */ +} ssr_ion_t; + typedef struct { /* SSR correction type */ gtime_t t0[6]; /* epoch time (GPST) {eph,clk,hrclk,ura,bias,pbias} */ double udi[6]; /* SSR update interval (s) */ @@ -963,6 +1001,7 @@ typedef struct { /* RTCM control struct type */ sta_t sta; /* station parameters */ dgps_t *dgps; /* output of dgps corrections */ ssr_t ssr[MAXSAT]; /* output of ssr corrections */ + ssr_ion_t ssr_ion; /* ionosphere SSR correction */ char msg[128]; /* special message */ char msgtype[256]; /* last message type */ char msmtype[7][128]; /* msm signal types */ diff --git a/util/data/IONO00IGS1_BNC.rtcm3 b/util/data/IONO00IGS1_BNC.rtcm3 new file mode 100644 index 0000000000000000000000000000000000000000..cf8459cb02be2c5709abe6ac7409ea3a73dcc85b GIT binary patch literal 531 zcmWlWJ&qJH6ojAqhY`c>l59d`O9TX@#MlQQC&LyOai4(;*cU-=!~qBhfCM{0B7tDV zYI=80UU&Y$x+!}S&tXz``|rG&n>n1t!$m9BYHA^}h&TdwLc3ES{bHgohf zJVcrg7Nxg!eiVly@=1{zGZJ)7Foc$ytu#jAS6$;Ujh4@)V`VLZgokB+*E^B6i5wLo zij-Ocbs4^*21F{tZmv^83A4}S+rDbpT_+rAeu-198M;ZtKGn!dGbC{G7w$&>wwH3N zFW3@hfn}z0jdME1-rr4E_P&gIbu`_Cnl4z)=r?SDlNR|830W=5?pEze2NI#uF)UaW z*@wJMd!{9wsWZLb%BT*iPfIw*F816(OO0tqV&;H}5;8NpK~_hBcN$I@5h>K&PXimX siOk%V^`3vkbfwc(wp4R-+z@+ICFa3z`@Q3~Bto^E7_i)|dHv=6e-0lGJpcdz literal 0 HcmV?d00001 diff --git a/util/data/IONO00IGS1_BNC.txt b/util/data/IONO00IGS1_BNC.txt new file mode 100644 index 000000000..be2b1d65d --- /dev/null +++ b/util/data/IONO00IGS1_BNC.txt @@ -0,0 +1,34 @@ +> VTEC 2025 07 15 10 00 45.0 2 1 IONO00IGS1 + 1 15 15 450000.0 + 18.9200 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 3.0400 8.6100 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -6.4700 0.1000 1.7200 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -1.5750 -1.3400 -1.1150 -0.0500 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 2.1000 -0.5850 0.0450 -0.1150 0.1300 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 1.4350 0.8100 -0.4700 0.3600 -0.0900 0.0350 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -0.7550 0.4400 -0.0200 0.0500 0.2650 -0.0500 -0.0850 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -0.1950 -0.3400 0.0500 -0.0450 0.1950 -0.0700 0.1150 0.0650 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.1600 0.0950 0.1450 -0.5200 0.0700 -0.1200 0.0700 -0.1400 0.1100 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -0.2050 -0.2500 0.4850 -0.0700 -0.3950 0.2250 -0.0850 0.0150 -0.0250 0.0850 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0300 -0.2500 -0.0900 0.4150 0.0050 0.0550 0.0400 0.0350 -0.0850 0.0450 0.0850 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.6250 -0.0250 -0.1550 -0.1000 0.2300 -0.1950 0.0550 0.0000 -0.1000 0.0500 0.0450 0.0050 0.0000 0.0000 0.0000 0.0000 + -0.1950 0.1300 0.0650 -0.2950 -0.0150 -0.0950 -0.0300 0.0250 -0.0150 0.0850 -0.0300 -0.0300 -0.0350 0.0000 0.0000 0.0000 + -0.3800 -0.0400 -0.0050 0.0300 -0.0350 0.1800 -0.1550 0.0250 0.0350 0.0450 -0.0550 0.0100 -0.0800 -0.0950 0.0000 0.0000 + 0.2000 -0.0050 -0.0350 0.2900 -0.0950 0.2100 0.0550 -0.0550 0.0500 0.0150 -0.0400 0.1000 -0.0150 0.0500 0.1850 0.0000 + 0.3000 -0.0550 0.0000 -0.0050 -0.0300 -0.0700 0.1350 -0.0400 -0.1450 -0.0700 0.0550 -0.1700 0.1450 0.0900 -0.0050 0.0300 + 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 1.5750 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.5600 -0.5800 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.9900 0.4900 1.2500 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.2050 -0.4450 -0.0100 0.4200 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.8050 0.2050 -0.8000 -0.3250 -0.1250 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.2050 -0.1850 0.0400 -0.1950 0.2050 -0.0700 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.1350 0.0200 -0.1100 0.0100 0.2850 0.2000 -0.1050 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.0100 0.5800 -0.0050 -0.2900 0.1450 0.0150 0.0200 -0.1950 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.0750 -0.0950 0.0900 -0.1400 -0.3400 -0.0600 -0.0550 0.0050 -0.0600 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.3300 -0.2050 0.1850 0.1800 -0.2550 -0.0200 -0.0350 -0.0050 -0.0250 -0.0550 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.0350 -0.1900 -0.1100 0.2100 0.1050 0.0250 -0.0250 -0.1300 -0.0050 -0.0900 -0.1150 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.5550 0.2150 -0.2000 -0.3350 0.1250 -0.0100 0.0500 0.0350 -0.0700 -0.0650 0.0700 -0.0250 0.0000 0.0000 0.0000 + 0.0000 -0.2550 -0.1250 0.2800 -0.2800 -0.0950 -0.0550 0.0650 0.0700 0.1150 0.0150 0.1750 0.0250 -0.1100 0.0000 0.0000 + 0.0000 -0.3750 -0.1250 0.0600 0.0900 -0.0550 -0.0250 0.0200 -0.0050 0.0850 0.0200 -0.1400 0.0000 -0.0950 -0.0300 0.0000 + 0.0000 0.1350 -0.0250 -0.1200 0.1350 -0.0650 0.1600 -0.0850 0.0100 -0.1200 0.0250 -0.0250 -0.0200 0.0750 -0.0850 0.0900 diff --git a/util/data/IONO00IGS1_ICD_example.rtcm3 b/util/data/IONO00IGS1_ICD_example.rtcm3 new file mode 100644 index 0000000000000000000000000000000000000000..6b383443b7e575c8fbc9b603272f0e06bc715821 GIT binary patch literal 47 zcmcc2p!x6cq;N?t1|VQvq9Dm&#Nf@4!cZ%p#ju{?1j9pyUyQ VTEC 2025 07 15 10 00 45.0 2 1 IONO00IGS1 + 1 3 2 450000.0 + 1.0000 0.0000 0.0000 + 2.0000 5.0000 0.0000 + 3.0000 -6.0000 8.0000 + 4.0000 7.0000 9.0000 + 0.0000 0.0000 0.0000 + 0.0000 10.0000 0.0000 + 0.0000 11.0000 -13.0000 + 0.0000 12.0000 14.0000 diff --git a/util/iono2rtcm/main.c b/util/iono2rtcm/main.c new file mode 100644 index 000000000..9468f168b --- /dev/null +++ b/util/iono2rtcm/main.c @@ -0,0 +1,252 @@ +//----------------------------------------------------------------------------- +// rnx2rtcm.c : RINEX to RTCM converter +// +// Copyright (C) 2012 by T.TAKASU, All rights reserved. +// +// Version : $Revision: 1.1 $ $Date: 2008/07/17 21:55:16 $ +// History : 2012/12/12 1.0 new +//----------------------------------------------------------------------------- +#include "rtklib.h" + +#define TRACEFILE "iono2rtcm.trace" // Debug trace file +#define MIN(x,y) ((x)<(y)?(x):(y)) + +// Print usage --------------------------------------------------------------- +static const char *help[] = {"", + "usage: iono2rtcm [options] [infile]", + "", + "options:", + " -out outfile output RTCM file", + " -x level debug trace level", + ""}; +static void print_help(void) { + for (int i = 0; i < sizeof(help) / sizeof(*help); i++) fprintf(stderr, "%s\n", help[i]); + exit(0); +} + +static int compute_index(const int n, const int m, const int n_degree) { + // Compute the index for spherical harmonics coefficients + return n - m + m * (n_degree + 1) - (m * (m - 1)) / 2; +} + +static int compute_update_interval(const int update_interval_code) { + + // Convert update interval code to seconds + switch (update_interval_code) { + case 0: return 1; + case 1: return 2; + case 2: return 5; + case 3: return 10; + case 4: return 15; + case 5: return 30; + case 6: return 60; + case 7: return 120; + case 8: return 240; + case 9: return 300; + case 10: return 600; + case 11: return 900; + case 12: return 1800; + case 13: return 3600; + case 14: return 7200; + case 15: return 10800; + default: return 60; // Default to 1 minute if unknown + } + +} + +/* This function parses the BNC VTEC file and fills the rtcm structure. + +The format of the BNC VTEC file for an example of N(degree)=3 and M(order)=2 + +C_00 0 0 +C_10 C_11 0 +C_20 C_21 C_22 +C_30 C_31 C_32 + 0 0 0 + 0 S_11 0 + 0 S_21 S_22 + 0 S_31 S_32 + +Note however than in the RTCM3 the coefficients are transposed +*/ +static int parse_bnc_vtec(const char *filename, rtcm_t *rtcm) { + + FILE *fp; + char line[1024], mountpoint[128], *p; + int n_degree = 0, m_order = 0; + int year, month, day, hour, min, n_blocks, layer_number; + int update_interval_code; + double sec, height_m; + gtime_t time; + + trace(3, "parse_bnc_vtec: file=%s\n", filename); + + if (!(fp = fopen(filename, "r"))) { + trace(1, "parse_bnc_vtec: file open error %s\n", filename); + return -1; + } + + // Initialize ionospheric SSR structure + memset(&rtcm->ssr_ion, 0, sizeof(ssr_ion_t)); + + while (fgets(line, sizeof(line), fp)) { + if (line[0] != '>') continue; // Skip non-header lines initially + + if (strstr(line, "VTEC") == NULL) { + continue; // Skip non-VTEC lines + } + // Parse header line: > VTEC 2025 07 15 10 00 45.0 2 1 IONO00IGS1 + if (sscanf(line, "> VTEC %d %d %d %d %d %lf %d %d %127s", + &year, &month, &day, &hour, &min, &sec, &update_interval_code, &n_blocks, mountpoint) != 9) { + trace(1, "parse_bnc_vtec: invalid header format\n"); + continue; + } + + time = epoch2time((double[]){year, month, day, hour, min, sec}); + rtcm->time = time; + + rtcm->ssr_ion.update_interval_s = compute_update_interval(update_interval_code);; + + trace(4, "VTEC header: %04d/%02d/%02d %02d:%02d:%05.2f update_interval=%d n_blocks=%d\n", + year, month, day, hour, min, sec, + rtcm->ssr_ion.update_interval_s, n_blocks); + + // Read next line with degree and order info + if (fgets(line, sizeof(line), fp) == NULL) { + trace(1, "parse_bnc_vtec: invalid BNC format\n"); + continue; + } + + if (sscanf(line, " %d %d %d %lf", &layer_number, &n_degree, &m_order, &height_m) != 4) { + trace(4, "Invalid format for layer header [ %s ]\n", line); + continue; + + } else { + trace(4, "Spherical harmonics: degree=%d order=%d height(m)=%.1f\n", + n_degree, m_order, height_m); + } + + // Validate degree and order limits + if (n_degree > 15) n_degree = 15; + if (m_order > 15) m_order = 15; + + // Set up ionospheric layer + ssr_ion_t* ssr_ion = &rtcm->ssr_ion; + ionlayer_sphharm_t *layer = &ssr_ion->ionlayers_sph_harm[0]; + + ssr_ion->n_layers = 1; // Assuming only one layer + layer->height_km = height_m / 1000.0; // Convert height from meters to kilometers + layer->n_degree = n_degree; + layer->m_order = m_order; + + // Initialize coefficients + memset(layer->cos_coeff, 0, sizeof(layer->cos_coeff)); + memset(layer->sin_coeff, 0, sizeof(layer->sin_coeff)); + + // Read cosine coefficients (C_nm) + for (int n = 0; n <= layer->n_degree; n++) { + if (!fgets(line, sizeof(line), fp)) break; + + p = line; + + for (int m = 0; m <= MIN(n, layer->m_order); m++) { + double coeff; + if (sscanf(p, "%lf", &coeff) == 1) { + int index = compute_index(n, m, layer->n_degree); + layer->cos_coeff[index] = coeff; + // Skip to next coefficient (11 characters for each coefficient) + p += 11; + } else { + break; + } + } + + } + + // Read sine coefficients (S_nm) - skip n=0 terms (they don't exist) + for (int n = 0; n <= layer->n_degree; n++) { + if (!fgets(line, sizeof(line), fp)) break; + + if (n == 0) continue; // Skip S_n0 (doesn't exist), start from m=1 + + p = line + 11; // Skip first coefficient (S_0m) + for (int m = 1; m <= MIN(n, layer->m_order); m++) { + + double coeff; + if (sscanf(p, "%lf", &coeff) == 1) { + int index = compute_index(n - 1, m - 1, layer->n_degree - 1); + layer->sin_coeff[index] = coeff; + // Skip to next number + p += 11; + } else { + break; + } + } + } + } + + fclose(fp); + return 0; +} + +// Main ---------------------------------------------------------------------- +int main(int argc, char **argv) { + + int ret = 1; + char *infile = NULL, *outfile = ""; + int trlevel = 0; + + rtcm_t rtcm = {0}; + + for (int i = 1; i < argc; i++) { + if (!strcmp(argv[i], "-out") && i + 1 < argc) + outfile = argv[++i]; + + else if (!strcmp(argv[i], "-x") && i + 1 < argc) + trlevel = atoi(argv[++i]); + + else if (!strncmp(argv[i], "-", 1)) + print_help(); + + else + infile = argv[i]; + } + + if (trlevel > 0) { + traceopen(TRACEFILE); + tracelevel(trlevel); + } + + if (infile == NULL) { + fprintf(stderr, "Input file is required.\n"); + print_help(); + return 2; + } + + FILE *fp = stdout; + if (*outfile && !(fp = fopen(outfile, "wb"))) { + fprintf(stderr, "file open error: %s\n", outfile); + return 0; + } + + if (parse_bnc_vtec(infile, &rtcm) != 0) { + fprintf(stderr, "Error parsing VTEC file: %s\n", infile); + goto end; + } + + if (!gen_rtcm3(&rtcm, 4076, 201, 0)) + goto end; + + if (fwrite(rtcm.buff, rtcm.nbyte, 1, fp) < 1) + goto end; + + ret = 0; +end: + + fclose(fp); + + if (trlevel > 0) { + traceclose(); + } + return ret; +} diff --git a/util/iono2rtcm/makefile b/util/iono2rtcm/makefile new file mode 100644 index 000000000..18a864646 --- /dev/null +++ b/util/iono2rtcm/makefile @@ -0,0 +1,60 @@ +# makefile for iono2rtcm + +BINDIR = /usr/local/bin +SRC = ../../src + + +# Build configuration options +DEBUG ?= 0 +PROFILE ?= 0 + +# Base compiler and linker flags +CFLAGS = -std=c99 -Wall -pedantic -I$(SRC) -DTRACE -DENAQZS -DENAGLO -DEXTLEX +LDFLAGS = -lm + +# Add optimization/debug flags based on configuration +ifeq ($(DEBUG),1) + CFLAGS += -g -O0 -DDEBUG + LDFLAGS += -g +else + CFLAGS += -O3 -DNDEBUG +endif + +# Add profiling flags if requested +ifeq ($(PROFILE),1) + CFLAGS += -pg + LDFLAGS += -pg +endif + +iono2rtcm : trace.o rtkcmn.o rinex.o rtcm.o rtcm2.o rtcm3.o rtcm3e.o main.c + $(CC) $(CFLAGS) $^ -o $@ $(LDFLAGS) + +rtkcmn.o : $(SRC)/rtkcmn.c + $(CC) -c $(CFLAGS) $(SRC)/rtkcmn.c +trace.o : $(SRC)/rtklib.h $(SRC)/trace.c + $(CC) -c $(CFLAGS) $(SRC)/trace.c +rinex.o : $(SRC)/rinex.c + $(CC) -c $(CFLAGS) $(SRC)/rinex.c +rtcm.o : $(SRC)/rtcm.c + $(CC) -c $(CFLAGS) $(SRC)/rtcm.c +rtcm2.o : $(SRC)/rtcm2.c + $(CC) -c $(CFLAGS) $(SRC)/rtcm2.c +rtcm3.o : $(SRC)/rtcm3.c + $(CC) -c $(CFLAGS) $(SRC)/rtcm3.c +rtcm3e.o : $(SRC)/rtcm3e.c + $(CC) -c $(CFLAGS) $(SRC)/rtcm3e.c + +install: + cp iono2rtcm $(BINDIR) + +clean: + rm -f iono2rtcm iono2rtcm.exe *.o + +test: test1 test2 + +test1: + ./iono2rtcm ../data/IONO00IGS1_ICD_example.txt > test.rtcm3 + diff -u test.rtcm3 ../data/IONO00IGS1_ICD_example.rtcm3 +test2: + ./iono2rtcm ../data/IONO00IGS1_BNC.txt > test.rtcm3 + diff -u test.rtcm3 ../data/IONO00IGS1_BNC.rtcm3 diff --git a/util/testeph/dumpssr.c b/util/testeph/dumpssr.c index fc4891af8..2cfc96215 100644 --- a/util/testeph/dumpssr.c +++ b/util/testeph/dumpssr.c @@ -35,18 +35,18 @@ static void printssrmsg(int sat, const ssr_t *ssr, int topt, int mopt) char tstr[40],id[8]; if (topt) { - time2str(ssr->t0,tstr,0); + time2str(ssr->t0[0],tstr,0); printf("%s ",tstr); } else { - tow=time2gpst(ssr->t0,&week); + tow=time2gpst(ssr->t0[0],&week); printf("%4d %6.0f ",week,tow); } satno2id(sat,id); printf("%4s ",id); if (mopt&1) { - printf("%4.0f %3d %3d %3d ",ssr->udint,ssr->iode,ssr->ura,ssr->refd); + printf("%4.0f %3d %3d %3d ", ssr->udi[0], ssr->iode, ssr->ura, ssr->refd); } if (mopt&2) { printf("%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f ",ssr->deph[0],ssr->deph[1], @@ -62,10 +62,136 @@ static void printssrmsg(int sat, const ssr_t *ssr, int topt, int mopt) ssr->cbias[4],ssr->cbias[5],ssr->cbias[6],ssr->cbias[7], ssr->cbias[8],ssr->cbias[9],ssr->cbias[10],ssr->cbias[11]); } + printf("\n"); } + +static int compute_index(const int n, const int m, const int n_degree) { + // Compute the index for spherical harmonics coefficients + return n - m + m * (n_degree + 1) - (m * (m - 1)) / 2; +} + +static int compute_update_interval_code(const int update_interval_s) { + // Convert update interval code to seconds + switch (update_interval_s) { + case 1: return 0; + case 2: return 1; + case 5: return 2; + case 10: return 3; + case 15: return 4; + case 30: return 5; + case 60: return 6; + case 120: return 7; + case 240: return 8; + case 300: return 9; + case 600: return 10; + case 900: return 11; + case 1800: return 12; + case 3600: return 13; + case 7200: return 14; + case 10800: return 15; + default: return 6; + } +} + +/* Print the ION SSR using BNC format + +See format description in Section 2.8 of +https://software.rtcm-ntrip.org/export/HEAD/ntrip/trunk/BNC/src/bnchelp.html#correct + +Sample format: + +> VTEC 2022 10 01 23 59 45.0 6 1 UNKNOWN + 1 3 2 450000.0 + 22.6900 0.0000 0.0000 + 0.2000 11.6050 0.0000 + -8.7500 -0.0200 2.0250 + -0.2900 -1.5900 0.3700 + 0.0000 0.0000 0.0000 + 0.0000 1.5350 0.0000 + 0.0000 0.1150 -1.3850 + 0.0000 -2.0150 0.2500 +*/ +static void print_vtec_ssr_bnc_format( + FILE *fp, const rtcm_t* rtcm, const char *stream) { + + double ep[6]; + int year, month, day, hour, min; + int update_interval_code; + double sec; + int idx; + double coeff; + + const ssr_ion_t *ssr_ion = &rtcm->ssr_ion; + const gtime_t *gtime = &rtcm->time; + + if (!ssr_ion || ssr_ion->n_layers == 0) { + return; /* No ionospheric data available */ + } + + /* Convert GPS time to calendar time */ + time2epoch(*gtime, ep); + year = (int)ep[0]; + month = (int)ep[1]; + day = (int)ep[2]; + hour = (int)ep[3]; + min = (int)ep[4]; + sec = ep[5]; + + update_interval_code = compute_update_interval_code(ssr_ion->update_interval_s); + + /* Print BNC VTEC header */ + fprintf(fp, "> VTEC %04d %02d %02d %02d %02d %04.1f %d %d %s\n", + year, month, day, hour, min, sec, update_interval_code, ssr_ion->n_layers, stream); + + for (int i = 0; i < ssr_ion->n_layers; i++) { + const ionlayer_sphharm_t *layer = &ssr_ion->ionlayers_sph_harm[i]; + + if (layer->n_degree <= 0 || layer->m_order <= 0) { + continue; /* Skip invalid layers */ + } + + /* Print layer header: layer_id, n_degree, m_order, height_km */ + fprintf(fp, "%2d %2d %2d %10.1f\n", + i + 1, layer->n_degree, layer->m_order, layer->height_km * 1000.0); + + /* Print cosine coefficients */ + for (int n = 0; n <= layer->n_degree; n++) { + for (int m = 0; m <= layer->m_order; m++) { + double coeff = 0.0; + if (m <= n) { + idx = compute_index(n, m, layer->n_degree); + coeff = layer->cos_coeff[idx]; + } + if (m > 0) + fprintf(fp, " "); + fprintf(fp, "%10.4f", coeff); + } + fprintf(fp, "\n"); + } + + /* Print sine coefficients */ + for (int n = 0; n <= layer->n_degree; n++) { + for (int m = 0; m <= layer->m_order; m++) { + coeff = 0.0; + + if (n > 0 && m > 0 && m <= n) { // First row and col always zero + idx = compute_index(n - 1, m - 1, layer->n_degree - 1); + coeff = layer->sin_coeff[idx]; + } + if (m > 0) + fprintf(fp, " "); + + fprintf(fp, "%10.4f", coeff); + } + fprintf(fp, "\n"); + } + } +} + + /* dump ssr messages ---------------------------------------------------------*/ -static void dumpssrmsg(FILE *fp, int sat, int topt, int mopt) +static void dumpssrmsg(FILE *fp, int sat, int topt, int mopt, const char *stream) { static rtcm_t rtcm; static gtime_t t0[MAXSAT]={{0}}; @@ -78,22 +204,26 @@ static void dumpssrmsg(FILE *fp, int sat, int topt, int mopt) if (stat!=10) continue; /* ssr message */ for (i=0;i ${TMPFILE1} + tail -8 ../data/$@.txt > ${TMPFILE2} + diff ${TMPFILE1} ${TMPFILE2} + +clean: + rm -f dumpssr diffeph *.o *.stackdump *.trace *.out *.exe $(TMPFILE1) $(TMPFILE2) From 69395186ebc7768ee741c3739518b43219651b0f Mon Sep 17 00:00:00 2001 From: Miquel Garcia Date: Wed, 23 Jul 2025 10:33:42 +0200 Subject: [PATCH 3/6] feat: add compute Legendre Function method Used for IONO SSR spherical harmonics corrections --- src/rtklib.h | 3 ++ src/ssr_ion.c | 83 +++++++++++++++++++++++++++++++++++++++++ test/utest/makefile | 8 +++- test/utest/t_ssr_ion.c | 84 ++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 177 insertions(+), 1 deletion(-) create mode 100644 src/ssr_ion.c create mode 100644 test/utest/t_ssr_ion.c diff --git a/src/rtklib.h b/src/rtklib.h index 8c4859e7b..9e650cf6e 100644 --- a/src/rtklib.h +++ b/src/rtklib.h @@ -1473,6 +1473,9 @@ EXPORT void matfprint(const double *A, int n, int m, int p, int q, FILE *fp); EXPORT void add_fatal(fatalfunc_t *func); +/* math functions ------------------------------------------------------------*/ +EXPORT int compute_legendre_polynomial(const int n, const int m, const double u, double* p); + /* time and string functions -------------------------------------------------*/ EXPORT void setstr(char *dst, const char *src, int n); EXPORT double str2num(const char *s, int i, int n); diff --git a/src/ssr_ion.c b/src/ssr_ion.c new file mode 100644 index 000000000..85f738b9f --- /dev/null +++ b/src/ssr_ion.c @@ -0,0 +1,83 @@ +/*------------------------------------------------------------------------------ +* ionex.c : ionospheric SSR corrections +* +* references: +* [1] IGS, "IGS State Space Representation (SSR)", IGS SSR v1.00 2020-10-05, +* Format Version 1, October 5th 2000 +* Available at: https://files.igs.org/pub/data/format/igs_ssr_v1.pdf" +*-----------------------------------------------------------------------------*/ + +#include "rtklib.h" + +/** \brief Compute the fully normalized Legendre functions (\$\tilde P_n^m(x)\$) + +Compute the Legendre polynomial, based on Section 6.7 of [1] (page 294). This +is equivalent to the function sph_harm_y in the scipy library (\$Y_n^m(\theta, \phi)\$), but +accepting a normalized input to the range [-1, 1]. + +References: +[1] Press, William H., Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. +"Numerical recipes 3rd edition." Cambridge: New York (2007). +[2] Holmes, S. A. and W. E. Featherstone, 2002. "A unified approach to the Clenshaw +summation and the recursive computation of very high degree and order normalised +associated Legendre functions". Journal of Geodesy, 76(5), pp. 279-299. +[3] https://mitgcm.org/~mlosch/geoidcookbook/node11.html +[4] Heiskanen and Moritz, Physical Geodesy, 1967, eq. 1-62 + +* args : int n I degree +* int m I order +* double u I input value, between -1 and 1 +* double p O \tilde P_n^m(u) +* return : >= 0 if successful +*/ +extern int compute_legendre_polynomial(const int n, const int m, const double u, double* p) { + + int ret = -1; + double fact, pll; + double pmm = 1.0; // P_m^m + + if (m < 0 || m > n || fabs(u) > 1.0 || p == NULL) { + goto end; + } + + if (m > 0) { + double omx2 = (1.0 - u) * (1.0 + u); + fact = 1.0; + for (int i = 1; i <= m; i++) { + pmm *= omx2 * fact / (fact + 1.0); + fact += 2.0; + } + } + + pmm = sqrt((2 * m + 1) * pmm / (4.0 * PI)); + + if (m & 1) { // If m is odd, negate pmm, mimics (-1)^m behavior + pmm = -pmm; + } + + if (n == m) { + *p = pmm; // P_m^m(u) = pmm + + } else { + double pmmp1 = u * sqrt(2.0 * m + 3.0) * pmm; + if (n == (m+1)) { + *p = pmmp1; // P_(m+1)^m(u) = pmmp1 + + } else { + double oldfact = sqrt(2.0 * m + 3.0); + + for (int i = m + 2; i <= n; i++) { + fact = sqrt((4.0 * i * i - 1.0)/(i * i - m * m)); + pll = (u * pmmp1 - pmm / oldfact) * fact; + oldfact = fact; + pmm = pmmp1; + pmmp1 = pll; + } + *p = pll; + } + } + + ret = 0; // Success +end: + return ret; +} diff --git a/test/utest/makefile b/test/utest/makefile index 7e505cbcc..641aaa3a0 100644 --- a/test/utest/makefile +++ b/test/utest/makefile @@ -7,7 +7,7 @@ LDLIBS = -lm -llapack -lblas CC = gcc BIN = t_matrix t_time t_coord t_rinex t_lambda t_atmos t_misc t_preceph t_gloeph \ -t_geoid t_ppp t_ionex t_tle +t_geoid t_ppp t_ionex t_tle t_ssr_ion all : $(BIN) t_matrix : t_matrix.o rtkcmn.o trace.o preceph.o @@ -24,6 +24,7 @@ t_ppp : t_ppp.o rtkcmn.o trace.o ephemeris.o preceph.o sbas.o ionex.o pntpo t_ppp : lambda.o tides.o t_ionex : t_ionex.o rtkcmn.o trace.o preceph.o ionex.o t_tle : t_tle.o rtkcmn.o trace.o rinex.o ephemeris.o sbas.o preceph.o tle.o +t_ssr_ion : t_ssr_ion.o rtkcmn.o trace.o preceph.o ephemeris.o sbas.o ionex.o pntpos.o ssr_ion.o rtkcmn.o : $(SRC)/rtklib.h $(SRC)/rtkcmn.c $(CC) -c $(CFLAGS) $(SRC)/rtkcmn.c @@ -55,9 +56,12 @@ tle.o : $(SRC)/rtklib.h $(SRC)/tle.c $(CC) -c $(CFLAGS) $(SRC)/tle.c tides.o : $(SRC)/rtklib.h $(SRC)/tides.c $(CC) -c $(CFLAGS) $(SRC)/tides.c +ssr_ion.o : $(SRC)/rtklib.h $(SRC)/ssr_ion.c + $(CC) -c $(CFLAGS) $(SRC)/ssr_ion.c utest : utest1 utest2 utest3 utest4 utest5 utest6 utest7 utest8 utest : utest9 utest10 utest11 utest12 utest14 +utest : utest_ssr_ion utest1 : ./t_matrix > utest1.out @@ -85,6 +89,8 @@ utest12 : ./t_ionex > utest12.out utest14 : ./t_tle > utest14.out +utest_ssr_ion : + ./t_ssr_ion > utest_ssr_ion.out clean : rm -f *.o *.out *.exe $(BIN) *.stackdump gmon.out diff --git a/test/utest/t_ssr_ion.c b/test/utest/t_ssr_ion.c new file mode 100644 index 000000000..aed92bee1 --- /dev/null +++ b/test/utest/t_ssr_ion.c @@ -0,0 +1,84 @@ +/*------------------------------------------------------------------------------ +* rtklib unit test driver : SSR ionospheric functions +*-----------------------------------------------------------------------------*/ +#include +#include +#include +#include "../../src/rtklib.h" + +/* Test compute_legendre_polynomial function */ +void test_legendre_polynomial(void) +{ + double result; + double tolerance = 1e-10; + + fprintf(stderr, "Testing compute_legendre_polynomial function...\n"); + + /* Test basic Legendre polynomials P_n^m(x) */ + + /* Various test cases for sectorial */ + struct legendre_test_cases { + int n; + int m; + double u; + double result; + } TEST_CASES[] = { + /** test cases built using the following Python code + + \verbatim + from math import acos + from scipy.special import sph_harm_y + + test_cases = [ + ( 1, 0, acos(0.5 )), + ( 2, 1, acos(0.5 )), + ( 3, 2, acos(0.1 )), + ( 4, 1, acos(-0.3 )), + ( 5, 3, acos(0.7 )), + (10, 5, acos(0.9 )), + (15, 15, acos(0.7071)), + (20, 10, acos(0.1 )), + (25, 5, acos(-0.9 )) + ] + + for t in test_cases: + result = sph_harm_y(*t, 0) + print(f"{t} {result}") + \endverbatim + */ + {.n = 1, .m = 0, .u = 0.5 , .result = 0.24430125595145988}, + {.n = 2, .m = 1, .u = 0.5 , .result = -0.33452327177864455}, + {.n = 3, .m = 2, .u = 0.1 , .result = 0.10117656216689495}, + {.n = 4, .m = 1, .u = -0.3 , .result = -0.3208718590203563}, + {.n = 5, .m = 3, .u = 0.7 , .result = -0.429650274134896}, + {.n = 10, .m = 5, .u = 0.9 , .result = -0.2643099927773962}, + {.n = 15, .m = 15, .u = 0.7071, .result = -0.0032983285307971377}, + {.n = 20, .m = 10, .u = 0.1 , .result = 0.07522105324402778}, + {.n = 25, .m = 5, .u = -0.9 , .result = 0.3756544620315938}, + }; + size_t n_test_cases = sizeof(TEST_CASES) / sizeof(TEST_CASES[0]); + + for (int i = 0; i < n_test_cases; i ++) { + const struct legendre_test_cases* tc = &TEST_CASES[i]; + assert (compute_legendre_polynomial(tc->n, tc->m, tc->u, &result) == 0);; + assert (fabs(result - tc->result) < tolerance); + } + + fprintf(stderr, "Testing compute_legendre_polynomial function (degraded cases)...\n"); + assert(compute_legendre_polynomial(-1, 0, 0.0, &result) < 0); + assert(compute_legendre_polynomial(0, 1, 0.0, &result) < 0); + assert(compute_legendre_polynomial(15, 15, 2.0, &result) < 0); + assert(compute_legendre_polynomial(15, 15, 0.0, NULL) < 0); + +} + +/* Main test function */ +int main(void) +{ + printf("=== Test Legendre Polynomial ===\n"); + + test_legendre_polynomial(); + + printf("\n=== All SSR ion tests completed ===\n"); + return 0; +} From 301957b06ae51b62f1a91508ff3038f5085a8706 Mon Sep 17 00:00:00 2001 From: Miquel Garcia Date: Thu, 24 Jul 2025 13:15:11 +0200 Subject: [PATCH 4/6] feat: add compute_stec_from_spherical_harmonics This will be required in order to use the ionospheric model from SSR corrections Also, added a test to exercise this new feature (data files that have been used for the tests have been included in this commit as well) --- src/rtklib.h | 2 + src/ssr_ion.c | 118 +++++++++++++++ .../bnc/IONO00CAS1_S_20252040000_01D_ION.ssr | 34 +++++ .../bnc/IONO01IGS1_S_20252040000_01D_ION.ssr | 34 +++++ test/utest/t_ssr_ion.c | 138 +++++++++++++++++- 5 files changed, 324 insertions(+), 2 deletions(-) create mode 100644 test/data/bnc/IONO00CAS1_S_20252040000_01D_ION.ssr create mode 100644 test/data/bnc/IONO01IGS1_S_20252040000_01D_ION.ssr diff --git a/src/rtklib.h b/src/rtklib.h index 9e650cf6e..db13759f4 100644 --- a/src/rtklib.h +++ b/src/rtklib.h @@ -1615,6 +1615,8 @@ EXPORT int ionocorr(gtime_t time, const nav_t *nav, int sat, const double *pos, EXPORT int tropcorr(gtime_t time, const nav_t *nav, const double *pos, const double *azel, int tropopt, double *trp, double *var); EXPORT int seliflc(int optnf, int sys); +EXPORT int compute_stec_from_spherical_harmonics( const ssr_ion_t* ssr_ion, const int epoch_s, + const double* pos, const double* azel, double re_km, double* stec); /* antenna models ------------------------------------------------------------*/ EXPORT int readpcv(const char *file, pcvs_t *pcvs); diff --git a/src/ssr_ion.c b/src/ssr_ion.c index 85f738b9f..e3e8fef9a 100644 --- a/src/ssr_ion.c +++ b/src/ssr_ion.c @@ -9,6 +9,13 @@ #include "rtklib.h" +#define MIN(x,y) ((x)<(y)?(x):(y)) + +/* Compute the index for spherical harmonics coefficients */ +extern int compute_index_ssr_iono_coeff_arrays(const int n, const int m, const int n_degree) { + return n - m + m * (n_degree + 1) - (m * (m - 1)) / 2; +} + /** \brief Compute the fully normalized Legendre functions (\$\tilde P_n^m(x)\$) Compute the Legendre polynomial, based on Section 6.7 of [1] (page 294). This @@ -81,3 +88,114 @@ extern int compute_legendre_polynomial(const int n, const int m, const double u, end: return ret; } + +/** \brief Compute VTEC using spherical harmonics coefficients */ +static int compute_vtec_for_layer( + const ionlayer_sphharm_t* layer, + const int epoch_s, + const double* pos, + double* vtec) { + + // Offset to approximate TEC maximum at 14:00 local time + // Note: the ICD for the IGS SSR (end of page 16) specifies 50400 s (14h) + // instead of 7200 s (2h). However the computation of the VTEC indicates + // that the IGS SSR products are computed based on these 2h (not 14h). + // This will likely require revision in the future + static const double EPOCH_OFFSET_S = 7200.0; // 50400.0; + double t, C_nm, S_nm; // Coefficients for spherical harmonics + double pnm; // Legendre function value + double lambda_s; // Longitude at sun fixed reference frame + + double ret = -1; + + if (pos == NULL || vtec == NULL) { + goto end; + } + + *vtec = 0.0; // Initialize VTEC for this layer + + /* mean sun fixed and phase shifted longitude of ionospheric pierce point + The mean sun fixed longitude phase shifted by 2 h to the approximate TEC + maximum at 14:00 local time (resp. 50400 s) */ + t = (epoch_s + 86400) % 86400; + lambda_s = pos[1] + (t - EPOCH_OFFSET_S) / 43200.0 * PI; + lambda_s = fmod(lambda_s + 2 * PI, 2 * PI); // Normalize to [0, 2*PI) + + for (int n = 0; n <= layer->n_degree; n++) { + for (int m = 0; m <= MIN(layer->m_order, n); m++) { + + int idx_cos = compute_index_ssr_iono_coeff_arrays(n, m, layer->n_degree); + C_nm = layer->cos_coeff[idx_cos]; + + S_nm = 0.0; + + if (n > 0 && m > 0) { + int idx_sin = compute_index_ssr_iono_coeff_arrays(n - 1, m - 1, layer->n_degree -1); + S_nm = layer->sin_coeff[idx_sin]; + } + + ret = compute_legendre_polynomial(n, m, sin(pos[0]), &pnm); + if (ret < 0) { + goto end; + } + + *vtec += (C_nm * cos(m * lambda_s) + S_nm * sin(m * lambda_s)) * pnm; + } + } + + ret = 0; +end: + return ret; + +} + +/** \brief Compute STEC from spherical harmonics + * + * compute ionospheric pierce point (ipp) position and slant factor +* args : ssr_ion_t* ssr_ion I SSR ionospheric correction data +* int epoch_s I Seconds of the day [0-86399) +* double *pos I receiver position {lat,lon,h} (rad,m) +* double *azel I azimuth/elevation angle {az,el} (rad) +* double re_km I earth radius (km) +* double* stec O Slant Total Electron Content +* return : 0 or positive if successful, negative otherwise +* notes : The pierce point can be computed with methods such as ionppp() +*-----------------------------------------------------------------------------*/ +extern int compute_stec_from_spherical_harmonics( + const ssr_ion_t* ssr_ion, + const int epoch_s, + const double* pos, + const double* azel, + double re_km, + double* stec) { + + double vtec; + double pospp[3]; + double psi_pp; + + int ret = -1; + + *stec = 0.0; + + /* aggregate all STEC values */ + for (int l = 0; l < ssr_ion->n_layers; l ++) { + const ionlayer_sphharm_t* layer = &ssr_ion->ionlayers_sph_harm[l]; + + ionppp(pos, azel, re_km, layer->height_km, pospp); + + ret = compute_vtec_for_layer(layer, epoch_s, pospp, &vtec); + if (ret < 0) { + goto end; + } + psi_pp = PI / 2.0 - azel[1] - asin((re_km + pospp[2]/1000.0) / (re_km + layer->height_km) * cos(azel[1])); + + *stec += vtec / sin(azel[1] + psi_pp); + + // *stec += vtec[l] / cos(asin((re_km + pospp[2]/1000.0) / (re_km + layer->height_km) * cos(azel[1]))); + + } + + ret = 0; +end: + return ret; +} diff --git a/test/data/bnc/IONO00CAS1_S_20252040000_01D_ION.ssr b/test/data/bnc/IONO00CAS1_S_20252040000_01D_ION.ssr new file mode 100644 index 000000000..14d97e7e9 --- /dev/null +++ b/test/data/bnc/IONO00CAS1_S_20252040000_01D_ION.ssr @@ -0,0 +1,34 @@ +> VTEC 2025 07 23 15 29 00.0 6 1 IONO00CAS1 + 1 15 15 450000.0 + 22.1750 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 4.0100 8.8750 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -7.8550 -0.2950 1.2150 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -2.0250 -2.2650 -0.3700 -0.4000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 1.8050 -0.5050 0.4450 0.7300 -0.0700 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 1.7350 1.0400 -0.0450 0.0250 -0.2450 0.1450 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -0.2850 0.2400 0.2450 -0.7050 -0.1300 -0.0100 -0.1850 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -0.3900 -0.2900 -0.0700 0.4750 -0.3150 -0.7950 0.0100 0.2750 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.1850 -0.1250 0.6250 -0.0600 -0.2800 -0.0700 -0.2650 0.2250 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.5600 -0.4600 0.1900 -0.0800 0.0250 -0.1700 -0.1600 -0.3300 0.3950 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -0.3200 0.0150 -0.3850 0.1250 0.0300 -0.0450 -0.0600 -0.0850 -0.1850 -0.0100 -0.0750 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.2150 0.3850 -0.3350 0.1550 0.2950 -0.1100 0.3550 -0.0850 0.1450 -0.6200 -0.1150 -0.0100 0.0000 0.0000 0.0000 0.0000 + 0.3350 -0.3550 -0.2150 0.1200 0.2350 -0.0500 0.3600 -0.1200 -0.0200 0.0850 0.0600 -0.1750 -0.0250 0.0000 0.0000 0.0000 + -0.3000 0.0800 0.2000 0.0800 -0.3550 -0.2550 -0.0200 0.1800 -0.4950 0.3200 0.0500 0.2000 0.0000 -0.2400 0.0000 0.0000 + -0.5000 0.3100 0.1600 0.2000 0.0950 -0.0800 0.0750 0.0300 0.1500 -0.1450 0.0950 -0.3350 0.0600 0.0450 0.3000 0.0000 + 0.2850 0.3050 0.1500 -0.1050 0.1700 0.1800 -0.1500 -0.1250 0.2050 -0.2800 0.0900 -0.0850 -0.0750 0.0550 0.3350 -0.0200 + 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.3200 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 2.9050 -1.7500 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.4450 0.7150 0.7200 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -1.9150 0.6400 0.4650 0.3150 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.0350 -1.0850 0.4850 -0.3650 -0.5900 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 1.0500 -0.4750 -0.2800 0.0500 -0.5500 0.3750 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.0800 0.0650 -0.2200 -0.2150 0.1000 -0.1000 0.1650 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.0000 -0.3150 -0.1500 -0.0750 0.4600 -0.2550 -0.2400 0.0600 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.5200 0.0500 0.3100 0.4650 0.0100 0.1000 0.1400 0.3500 0.3700 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.1100 -0.1000 -0.1350 0.4000 -0.3150 0.0600 0.3300 0.2700 0.0300 0.3900 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.7000 -0.0800 -0.3750 0.0600 -0.2300 -0.0900 0.3000 -0.0250 0.1000 0.0350 -0.1500 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.0550 -0.1300 0.2800 0.0300 -0.0200 -0.2600 0.0250 -0.0950 0.1300 0.0550 0.1850 0.0050 0.0000 0.0000 0.0000 + 0.0000 0.2300 0.1500 -0.2200 -0.0850 -0.0950 -0.0550 -0.1500 -0.1100 0.0900 -0.2600 0.0200 -0.1200 0.1250 0.0000 0.0000 + 0.0000 -0.6200 -0.0300 -0.0950 0.2600 -0.2800 -0.2350 -0.0250 0.1550 -0.2200 0.1850 0.0000 0.0100 -0.0800 0.1950 0.0000 + 0.0000 -0.3200 0.1150 0.2150 0.2050 0.3600 -0.0400 0.0500 0.0550 0.2400 -0.0300 0.2350 0.0400 0.3800 -0.5300 -0.3000 diff --git a/test/data/bnc/IONO01IGS1_S_20252040000_01D_ION.ssr b/test/data/bnc/IONO01IGS1_S_20252040000_01D_ION.ssr new file mode 100644 index 000000000..58b5bd988 --- /dev/null +++ b/test/data/bnc/IONO01IGS1_S_20252040000_01D_ION.ssr @@ -0,0 +1,34 @@ +> VTEC 2025 07 23 15 29 00.0 6 1 IONO01IGS1 + 1 15 15 450000.0 + 21.3050 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 4.1100 9.0850 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -7.3450 -0.3350 1.1900 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -1.7800 -2.1150 -0.6700 -0.6050 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 1.8000 -0.3750 0.2300 0.3850 0.1600 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 1.2500 0.9250 -0.3500 0.2900 0.0000 0.0200 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -0.3400 0.3250 0.0300 -0.3400 0.0550 -0.0850 -0.0800 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -0.2600 -0.3900 0.0750 0.1750 -0.0500 -0.2050 0.0650 0.1650 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.1900 -0.0800 0.0150 0.2050 -0.0600 -0.0850 -0.0800 -0.0500 0.0450 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -0.0850 -0.1600 -0.1350 0.0950 -0.0850 -0.0400 -0.0400 -0.0400 -0.1500 0.0900 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + -0.2400 -0.0850 -0.0050 0.0150 -0.0100 0.0050 -0.0300 -0.0100 -0.0350 -0.0150 -0.0050 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.2000 0.0600 -0.0350 0.0900 0.0950 0.0350 0.1550 -0.0550 -0.0050 -0.1600 0.0200 0.0100 0.0000 0.0000 0.0000 0.0000 + 0.3100 -0.0350 -0.1000 -0.0800 0.0600 -0.0800 0.1050 -0.0350 -0.0150 0.0550 0.0100 -0.0550 -0.0200 0.0000 0.0000 0.0000 + -0.2000 0.0350 0.0350 0.0350 -0.1200 -0.1000 -0.0350 0.0450 -0.1250 0.1100 -0.0100 0.0450 -0.0100 -0.0650 0.0000 0.0000 + -0.2000 0.1600 0.0600 0.0750 0.0500 -0.0200 0.0000 0.0050 0.0650 -0.0550 0.0150 -0.0950 0.0200 0.0250 0.0650 0.0000 + 0.1700 0.0700 0.0350 -0.0300 0.0550 0.0950 -0.0300 -0.0550 0.0550 -0.0800 0.0300 0.0000 -0.0300 0.0200 0.0550 -0.0250 + 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.6000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 2.1700 -1.1400 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.5950 0.7650 0.8400 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -1.4100 0.1500 0.1200 0.3200 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.1300 -0.5700 0.0900 -0.3800 -0.4450 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.7100 -0.2500 -0.1650 0.0150 -0.1850 0.1200 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.1550 0.2950 -0.1550 -0.1200 0.1550 0.1050 0.1500 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.1050 -0.0700 -0.0300 -0.0500 0.2150 -0.0750 -0.0050 -0.0600 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.2850 -0.0150 0.1000 0.1000 -0.0150 0.0400 0.0000 0.0650 0.0650 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.0100 -0.0850 -0.0950 0.1850 -0.1800 0.0000 0.0400 0.0700 -0.0150 0.1100 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 -0.3000 -0.0750 -0.1050 0.0000 -0.0200 -0.0350 0.1050 -0.0350 0.0200 -0.0300 -0.0150 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.2100 0.0000 0.1350 -0.0050 -0.0050 -0.0800 0.0150 0.0000 0.0300 0.0000 0.0950 -0.0200 0.0000 0.0000 0.0000 + 0.0000 0.2250 0.0300 -0.1550 -0.0650 0.0050 -0.0300 -0.0450 0.0300 0.0000 -0.0550 0.0050 -0.0100 0.0350 0.0000 0.0000 + 0.0000 -0.2500 -0.0450 -0.0400 0.0400 -0.0850 -0.0500 0.0200 0.0250 -0.0300 0.0400 0.0150 -0.0050 -0.0150 0.0150 0.0000 + 0.0000 -0.1750 0.0550 0.1550 0.0700 0.0900 0.0000 0.0300 -0.0200 0.0550 -0.0050 0.0550 0.0200 0.0900 -0.1750 -0.0850 diff --git a/test/utest/t_ssr_ion.c b/test/utest/t_ssr_ion.c index aed92bee1..92479b3d4 100644 --- a/test/utest/t_ssr_ion.c +++ b/test/utest/t_ssr_ion.c @@ -7,7 +7,7 @@ #include "../../src/rtklib.h" /* Test compute_legendre_polynomial function */ -void test_legendre_polynomial(void) +static void test_legendre_polynomial(void) { double result; double tolerance = 1e-10; @@ -72,12 +72,146 @@ void test_legendre_polynomial(void) } +static void test_compute_stec_method_to_compute_vtec(void) { + + /* spherical harmonic coefficients extracted from IONO00IGS1_BNC.txt, 2025 07 15 10 00 45.0 */ + const ssr_ion_t SSR_ION_IONO01IGS1 = { + .ionlayers_sph_harm = { + { .height_km = 450.0, .n_degree = 15, .m_order = 15, + .cos_coeff = {2.1305e+01, 4.1100e+00, -7.3450e+00, -1.7800e+00, 1.8000e+00, + 1.2500e+00, -3.4000e-01, -2.6000e-01, 1.9000e-01, -8.5000e-02, + -2.4000e-01, 2.0000e-01, 3.1000e-01, -2.0000e-01, -2.0000e-01, + 1.7000e-01, 9.0850e+00, -3.3500e-01, -2.1150e+00, -3.7500e-01, + 9.2500e-01, 3.2500e-01, -3.9000e-01, -8.0000e-02, -1.6000e-01, + -8.5000e-02, 6.0000e-02, -3.5000e-02, 3.5000e-02, 1.6000e-01, + 7.0000e-02, 1.1900e+00, -6.7000e-01, 2.3000e-01, -3.5000e-01, + 3.0000e-02, 7.5000e-02, 1.5000e-02, -1.3500e-01, -5.0000e-03, + -3.5000e-02, -1.0000e-01, 3.5000e-02, 6.0000e-02, 3.5000e-02, + -6.0500e-01, 3.8500e-01, 2.9000e-01, -3.4000e-01, 1.7500e-01, + 2.0500e-01, 9.5000e-02, 1.5000e-02, 9.0000e-02, -8.0000e-02, + 3.5000e-02, 7.5000e-02, -3.0000e-02, 1.6000e-01, 0.0000e+00, + 5.5000e-02, -5.0000e-02, -6.0000e-02, -8.5000e-02, -1.0000e-02, + 9.5000e-02, 6.0000e-02, -1.2000e-01, 5.0000e-02, 5.5000e-02, + 2.0000e-02, -8.5000e-02, -2.0500e-01, -8.5000e-02, -4.0000e-02, + 5.0000e-03, 3.5000e-02, -8.0000e-02, -1.0000e-01, -2.0000e-02, + 9.5000e-02, -8.0000e-02, 6.5000e-02, -8.0000e-02, -4.0000e-02, + -3.0000e-02, 1.5500e-01, 1.0500e-01, -3.5000e-02, 0.0000e+00, + -3.0000e-02, 1.6500e-01, -5.0000e-02, -4.0000e-02, -1.0000e-02, + -5.5000e-02, -3.5000e-02, 4.5000e-02, 5.0000e-03, -5.5000e-02, + 4.5000e-02, -1.5000e-01, -3.5000e-02, -5.0000e-03, -1.5000e-02, + -1.2500e-01, 6.5000e-02, 5.5000e-02, 9.0000e-02, -1.5000e-02, + -1.6000e-01, 5.5000e-02, 1.1000e-01, -5.5000e-02, -8.0000e-02, + -5.0000e-03, 2.0000e-02, 1.0000e-02, -1.0000e-02, 1.5000e-02, + 3.0000e-02, 1.0000e-02, -5.5000e-02, 4.5000e-02, -9.5000e-02, + 0.0000e+00, -2.0000e-02, -1.0000e-02, 2.0000e-02, -3.0000e-02, + -6.5000e-02, 2.5000e-02, 2.0000e-02, 6.5000e-02, 5.5000e-02, + -2.5000e-02}, + .sin_coeff = { 0.6 , 2.17 , -0.595, -1.41 , 0.13 , 0.71 , -0.155, -0.105, + 0.285, 0.01 , -0.3 , 0.21 , 0.225, -0.25 , -0.175, -1.14 , + 0.765, 0.15 , -0.57 , -0.25 , 0.295, -0.07 , -0.015, -0.085, + -0.075, 0. , 0.03 , -0.045, 0.055, 0.84 , 0.12 , 0.09 , + -0.165, -0.155, -0.03 , 0.1 , -0.095, -0.105, 0.135, -0.155, + -0.04 , 0.155, 0.32 , -0.38 , 0.015, -0.12 , -0.05 , 0.1 , + 0.185, 0. , -0.005, -0.065, 0.04 , 0.07 , -0.445, -0.185, + 0.155, 0.215, -0.015, -0.18 , -0.02 , -0.005, 0.005, -0.085, + 0.09 , 0.12 , 0.105, -0.075, 0.04 , 0. , -0.035, -0.08 , + -0.03 , -0.05 , 0. , 0.15 , -0.005, 0. , 0.04 , 0.105, + 0.015, -0.045, 0.02 , 0.03 , -0.06 , 0.065, 0.07 , -0.035, + 0. , 0.03 , 0.025, -0.02 , 0.065, -0.015, 0.02 , 0.03 , + 0. , -0.03 , 0.055, 0.11 , -0.03 , 0. , -0.055, 0.04 , + -0.005, -0.015, 0.095, 0.005, 0.015, 0.055, -0.02 , -0.01 , + -0.005, 0.02 , 0.035, -0.015, 0.09 , 0.015, -0.175, -0.085} } + }, + .update_interval_s = 30, + .n_layers = 1 + }; + + const ssr_ion_t SSR_ION_IONO00CAS1 = { + .ionlayers_sph_harm = { + { .height_km = 450.0, .n_degree = 15, .m_order = 15, + .cos_coeff = {2.2175e+01, 4.0100e+00, -7.8550e+00, -2.0250e+00, 1.8050e+00, + 1.7350e+00, -2.8500e-01, -3.9000e-01, 0.0000e+00, 0.0000e+00, + -3.2000e-01, 2.1500e-01, 3.3500e-01, -3.0000e-01, -5.0000e-01, + 2.8500e-01, 8.8750e+00, -2.9500e-01, -2.2650e+00, -5.0500e-01, + 1.0400e+00, 2.4000e-01, -2.9000e-01, -1.8500e-01, -5.6000e-01, + 1.5000e-02, 3.8500e-01, -3.5500e-01, 8.0000e-02, 3.1000e-01, + 3.0500e-01, 1.2150e+00, -3.7000e-01, 4.4500e-01, -4.5000e-02, + 2.4500e-01, -7.0000e-02, -1.2500e-01, -4.6000e-01, -3.8500e-01, + -3.3500e-01, -2.1500e-01, 2.0000e-01, 1.6000e-01, 1.5000e-01, + -4.0000e-01, 7.3000e-01, 2.5000e-02, -7.0500e-01, 4.7500e-01, + 6.2500e-01, 1.9000e-01, 1.2500e-01, 1.5500e-01, 1.2000e-01, + 8.0000e-02, 2.0000e-01, -1.0500e-01, -7.0000e-02, -2.4500e-01, + -1.3000e-01, -3.1500e-01, -6.0000e-02, -8.0000e-02, 3.0000e-02, + 2.9500e-01, 2.3500e-01, -3.5500e-01, 9.5000e-02, 1.7000e-01, + 1.4500e-01, -1.0000e-02, -7.9500e-01, -2.8000e-01, 2.5000e-02, + -4.5000e-02, -1.1000e-01, -5.0000e-02, -2.5500e-01, -8.0000e-02, + 1.8000e-01, -1.8500e-01, 1.0000e-02, -7.0000e-02, -1.7000e-01, + -6.0000e-02, 3.5500e-01, 3.6000e-01, -2.0000e-02, 7.5000e-02, + -1.5000e-01, 2.7500e-01, -2.6500e-01, -1.6000e-01, -8.5000e-02, + -8.5000e-02, -1.2000e-01, 1.8000e-01, 3.0000e-02, -1.2500e-01, + 2.2500e-01, -3.3000e-01, -1.8500e-01, 1.4500e-01, -2.0000e-02, + -4.9500e-01, 1.5000e-01, 2.0500e-01, 3.9500e-01, -1.0000e-02, + -6.2000e-01, 8.5000e-02, 3.2000e-01, -1.4500e-01, -2.8000e-01, + -7.5000e-02, -1.1500e-01, 6.0000e-02, 5.0000e-02, 9.5000e-02, + 9.0000e-02, -1.0000e-02, -1.7500e-01, 2.0000e-01, -3.3500e-01, + -8.5000e-02, -2.5000e-02, 0.0000e+00, 6.0000e-02, -7.5000e-02, + -2.4000e-01, 4.5000e-02, 5.5000e-02, 3.0000e-01, 3.3500e-01, + -2.0000e-02}, + .sin_coeff = {0.32 , 2.905, -0.445, -1.915, 0.035, 1.05 , -0.08 , 0. , + 0.52 , 0.11 , -0.7 , 0.055, 0.23 , -0.62 , -0.32 , -1.75 , + 0.715, 0.64 , -1.085, -0.475, 0.065, -0.315, 0.05 , -0.1 , + -0.08 , -0.13 , 0.15 , -0.03 , 0.115, 0.72 , 0.465, 0.485, + -0.28 , -0.22 , -0.15 , 0.31 , -0.135, -0.375, 0.28 , -0.22 , + -0.095, 0.215, 0.315, -0.365, 0.05 , -0.215, -0.075, 0.465, + 0.4 , 0.06 , 0.03 , -0.085, 0.26 , 0.205, -0.59 , -0.55 , + 0.1 , 0.46 , 0.01 , -0.315, -0.23 , -0.02 , -0.095, -0.28 , + 0.36 , 0.375, -0.1 , -0.255, 0.1 , 0.06 , -0.09 , -0.26 , + -0.055, -0.235, -0.04 , 0.165, -0.24 , 0.14 , 0.33 , 0.3 , + 0.025, -0.15 , -0.025, 0.05 , 0.06 , 0.35 , 0.27 , -0.025, + -0.095, -0.11 , 0.155, 0.055, 0.37 , 0.03 , 0.1 , 0.13 , + 0.09 , -0.22 , 0.24 , 0.39 , 0.035, 0.055, -0.26 , 0.185, + -0.03 , -0.15 , 0.185, 0.02 , 0. , 0.235, 0.005, -0.12 , + 0.01 , 0.04 , 0.125, -0.08 , 0.38 , 0.195, -0.53 , -0.3} } + }, + .update_interval_s = 30, + .n_layers = 1 + }; + + int EPOCH_S = 15 * 3600 + 29 * 60 + 00; // 2025 07 23 15 29 00.0 + + double pos[3] = {0.0, 0.0, 0.0}; // Dummy position + double azel[2] = {0.0, 90.0 * D2R}; // Zenith looking direction + double stec_igs = 0.0; + double stec_cas = 0.0; + double re_km = 6371.0; // Earth radius in km + + fprintf(stderr, "Testing compute_stec_from_spherical_harmonics function (VTEC computation)...\n"); + + // Initialize SSR ionospheric data with dummy values + for (int lon = -180; lon < +180; lon += 10) { + pos[1] = lon * D2R; + for (int lat = -90; lat <= +90; lat += 10) { + pos[0] = lat * D2R; + assert(compute_stec_from_spherical_harmonics(&SSR_ION_IONO01IGS1, EPOCH_S, pos, azel, re_km, &stec_igs) == 0); + assert(compute_stec_from_spherical_harmonics(&SSR_ION_IONO00CAS1, EPOCH_S, pos, azel, re_km, &stec_cas) == 0); + assert(fabs(stec_igs - stec_cas) < 2.5); // Check if the STEC values are close enough + // fprintf(stdout, "VTEC %i %lf %lf IONO01IGS1: %lf IONO00CAS1: %lf\n", + // EPOCH_S, pos[0] * R2D, pos[1] * R2D, stec_igs, stec_cas); + } + } + + + fprintf(stderr, "compute_stec_from_spherical_harmonics passed.\n"); + +} + + /* Main test function */ int main(void) { - printf("=== Test Legendre Polynomial ===\n"); test_legendre_polynomial(); + test_compute_stec_method_to_compute_vtec(); printf("\n=== All SSR ion tests completed ===\n"); return 0; From d950c6332cac9ab40e515d2d5a8b8127cbdf4e60 Mon Sep 17 00:00:00 2001 From: Miquel Garcia Date: Thu, 31 Jul 2025 09:53:05 +0200 Subject: [PATCH 5/6] chore(rtksvr): remove trailing spaces --- src/options.c | 66 ++++++++++----------- src/pntpos.c | 110 +++++++++++++++++----------------- src/rtksvr.c | 160 +++++++++++++++++++++++++------------------------- 3 files changed, 168 insertions(+), 168 deletions(-) diff --git a/src/options.c b/src/options.c index 5bfecba87..46d8038e0 100644 --- a/src/options.c +++ b/src/options.c @@ -88,7 +88,7 @@ EXPORT opt_t sysopts[]={ {"pos1-posopt6", 3, (void *)&prcopt_.posopt[5], SWTOPT }, {"pos1-exclsats", 2, (void *)exsats_, "prn ..."}, {"pos1-navsys", 0, (void *)&prcopt_.navsys, NAVOPT }, - + {"pos2-armode", 3, (void *)&prcopt_.modear, ARMOPT }, {"pos2-gloarmode", 3, (void *)&prcopt_.glomodear, GAROPT }, {"pos2-bdsarmode", 3, (void *)&prcopt_.bdsmodear, SWTOPT }, @@ -120,7 +120,7 @@ EXPORT opt_t sysopts[]={ {"pos2-niter", 0, (void *)&prcopt_.niter, "" }, {"pos2-baselen", 1, (void *)&prcopt_.baseline[0],"m" }, {"pos2-basesig", 1, (void *)&prcopt_.baseline[1],"m" }, - + {"out-solformat", 3, (void *)&solopt_.posf, SOLOPT }, {"out-outhead", 3, (void *)&solopt_.outhead, SWTOPT }, {"out-outopt", 3, (void *)&solopt_.outopt, SWTOPT }, @@ -159,7 +159,7 @@ EXPORT opt_t sysopts[]={ {"stats-prntrop", 1, (void *)&prcopt_.prn[2], "m" }, {"stats-prnpos", 1, (void *)&prcopt_.prn[5], "m" }, {"stats-clkstab", 1, (void *)&prcopt_.sclkstab, "s/s" }, - + {"ant1-postype", 3, (void *)&prcopt_.rovpos, POSOPT }, {"ant1-pos1", 1, (void *)&antpos_[0][0], "deg|m"}, {"ant1-pos2", 1, (void *)&antpos_[0][1], "deg|m"}, @@ -168,7 +168,7 @@ EXPORT opt_t sysopts[]={ {"ant1-antdele", 1, (void *)&prcopt_.antdel[0][0],"m" }, {"ant1-antdeln", 1, (void *)&prcopt_.antdel[0][1],"m" }, {"ant1-antdelu", 1, (void *)&prcopt_.antdel[0][2],"m" }, - + {"ant2-postype", 3, (void *)&prcopt_.refpos, POSOPT }, {"ant2-pos1", 1, (void *)&antpos_[1][0], "deg|m"}, {"ant2-pos2", 1, (void *)&antpos_[1][1], "deg|m"}, @@ -179,13 +179,13 @@ EXPORT opt_t sysopts[]={ {"ant2-antdelu", 1, (void *)&prcopt_.antdel[1][2],"m" }, {"ant2-maxaveep", 0, (void *)&prcopt_.maxaveep ,"" }, {"ant2-initrst", 3, (void *)&prcopt_.initrst, SWTOPT }, - + {"misc-timeinterp", 3, (void *)&prcopt_.intpref, SWTOPT }, {"misc-sbasatsel", 0, (void *)&prcopt_.sbassatsel, "0:all"}, {"misc-rnxopt1", 2, (void *)prcopt_.rnxopt[0], "" }, {"misc-rnxopt2", 2, (void *)prcopt_.rnxopt[1], "" }, {"misc-pppopt", 2, (void *)prcopt_.pppopt, "" }, - + {"file-satantfile", 2, (void *)&filopt_.satantp, "" }, {"file-rcvantfile", 2, (void *)&filopt_.rcvantp, "" }, {"file-staposfile", 2, (void *)&filopt_.stapos, "" }, @@ -198,7 +198,7 @@ EXPORT opt_t sysopts[]={ {"file-geexefile", 2, (void *)&filopt_.geexe, "" }, {"file-solstatfile",2, (void *)&filopt_.solstat, "" }, {"file-tracefile", 2, (void *)&filopt_.trace, "" }, - + {"",0,NULL,""} /* terminator */ }; /* discard space characters at tail ------------------------------------------*/ @@ -213,7 +213,7 @@ static int enum2str(char *s, const char *comment, int val) { char str[32],*p,*q; int n; - + n=sprintf(str,"%d:",val); if (!(p=strstr(comment,str))) { return sprintf(s,"%d",val); @@ -263,9 +263,9 @@ static int str2enum(const char *str, const char *comment, int *val) { extern opt_t *searchopt(const char *name, const opt_t *opts) { int i; - + trace(3,"searchopt: name=%s\n",name); - + for (i=0;*opts[i].name;i++) { if (strstr(opts[i].name,name)) return (opt_t *)(opts+i); } @@ -297,9 +297,9 @@ extern int str2opt(opt_t *opt, const char *str) extern int opt2str(const opt_t *opt, char *str) { char *p=str; - + trace(3,"opt2str : name=%s\n",opt->name); - + switch (opt->format) { case 0: p+=sprintf(p,"%d" ,*(int *)opt->var); break; case 1: p+=sprintf(p,"%.15g",*(double*)opt->var); break; @@ -318,9 +318,9 @@ extern int opt2buf(const opt_t *opt, char *buff) { char *p=buff; int n; - + trace(3,"opt2buf : name=%s\n",opt->name); - + p+=sprintf(p,"%-18s =",opt->name); p+=opt2str(opt,p); if (*opt->comment) { @@ -342,9 +342,9 @@ extern int loadopts(const char *file, opt_t *opts) opt_t *opt; char buff[2048],*p; int n=0; - + trace(3,"loadopts: file=%s\n",file); - + if (!(fp=fopen(file,"r"))) { trace(1,"loadopts: options file open error (%s)\n",file); return 0; @@ -352,9 +352,9 @@ extern int loadopts(const char *file, opt_t *opts) while (fgets(buff,sizeof(buff),fp)) { n++; chop(buff); - + if (buff[0]=='\0') continue; - + if (!(p=strstr(buff,"="))) { fprintf(stderr,"invalid option %s (%s:%d)\n",buff,file,n); continue; @@ -362,14 +362,14 @@ extern int loadopts(const char *file, opt_t *opts) *p++='\0'; chop(buff); if (!(opt=searchopt(buff,opts))) continue; - + if (!str2opt(opt,p)) { fprintf(stderr,"invalid option value %s (%s:%d)\n",buff,file,n); continue; } } fclose(fp); - + return 1; } /* save options to file -------------------------------------------------------- @@ -387,15 +387,15 @@ extern int saveopts(const char *file, const char *mode, const char *comment, FILE *fp; char buff[2048]; int i; - + trace(3,"saveopts: file=%s mode=%s\n",file,mode); - + if (!(fp=fopen(file,mode))) { trace(1,"saveopts: options file open error (%s)\n",file); return 0; } if (comment) fprintf(fp,"# %s\n\n",comment); - + for (i=0;*opts[i].name;i++) { opt2buf(opts+i,buff); fprintf(fp,"%s\n",buff); @@ -409,15 +409,15 @@ static void buff2sysopts(void) double pos[3],*rr; char buff[1024],*p,*id; int i,j,sat,ps; - + prcopt_.elmin =elmask_ *D2R; prcopt_.elmaskar =elmaskar_ *D2R; prcopt_.elmaskhold=elmaskhold_*D2R; - + for (i=0;i<2;i++) { ps=i==0?prcopt_.rovpos:prcopt_.refpos; rr=i==0?prcopt_.ru:prcopt_.rb; - + if (ps==POSOPT_POS_LLH) { /* lat/lon/hgt */ pos[0]=antpos_[i][0]*D2R; pos[1]=antpos_[i][1]*D2R; @@ -467,15 +467,15 @@ static void sysopts2buff(void) double pos[3],*rr; char id[8],*p; int i,j,sat,ps; - + elmask_ =prcopt_.elmin *R2D; elmaskar_ =prcopt_.elmaskar *R2D; elmaskhold_=prcopt_.elmaskhold*R2D; - + for (i=0;i<2;i++) { ps=i==0?prcopt_.rovpos:prcopt_.refpos; rr=i==0?prcopt_.ru:prcopt_.rb; - + if (ps==POSOPT_POS_LLH) { ecef2pos(rr,pos); antpos_[i][0]=pos[0]*R2D; @@ -518,9 +518,9 @@ static void sysopts2buff(void) extern void resetsysopts(void) { int i,j; - + trace(3,"resetsysopts:\n"); - + prcopt_=prcopt_default; solopt_=solopt_default; filopt_.satantp[0]='\0'; @@ -550,7 +550,7 @@ extern void resetsysopts(void) extern void getsysopts(prcopt_t *popt, solopt_t *sopt, filopt_t *fopt) { trace(3,"getsysopts:\n"); - + buff2sysopts(); if (popt) *popt=prcopt_; if (sopt) *sopt=solopt_; @@ -568,7 +568,7 @@ extern void setsysopts(const prcopt_t *prcopt, const solopt_t *solopt, const filopt_t *filopt) { trace(3,"setsysopts:\n"); - + resetsysopts(); if (prcopt) prcopt_=*prcopt; if (solopt) solopt_=*solopt; diff --git a/src/pntpos.c b/src/pntpos.c index 27708fc0d..6e56b1b47 100644 --- a/src/pntpos.c +++ b/src/pntpos.c @@ -79,7 +79,7 @@ static double varerr(const prcopt_t *opt, const ssat_t *ssat, const obsd_t *obs, static double gettgd(int sat, const nav_t *nav, int type) { int i,sys=satsys(sat,NULL); - + if (sys==SYS_GLO) { for (i=0;ing;i++) { if (nav->geph[i].sat==sat) break; @@ -120,7 +120,7 @@ static double prange(const obsd_t *obs, const nav_t *nav, const prcopt_t *opt, f2=seliflc(opt->nf,satsys(obs->sat,NULL)); P2=obs->P[f2]; *var=0.0; - + if (P1==0.0||(opt->ionoopt==IONOOPT_IFLC&&P2==0.0)) return 0.0; bias_ix=code2bias_ix(sys,obs->code[0]); /* L1 code bias */ if (bias_ix>0) { /* 0=ref code */ @@ -137,7 +137,7 @@ static double prange(const obsd_t *obs, const nav_t *nav, const prcopt_t *opt, } } if (opt->ionoopt==IONOOPT_IFLC) { /* dual-frequency */ - + if (sys==SYS_GPS||sys==SYS_QZS) { /* L1-L2 or L1-L5 */ gamma=f2==1?SQR(FREQL1/FREQL2):SQR(FREQL1/FREQL5); return (P2-gamma*P1)/(1.0-gamma); @@ -168,7 +168,7 @@ static double prange(const obsd_t *obs, const nav_t *nav, const prcopt_t *opt, } else { /* single-freq (L1/E1/B1) */ *var=SQR(ERR_CBIAS); - + if (sys==SYS_GPS||sys==SYS_QZS) { /* L1 */ b1=gettgd(sat,nav,0); /* TGD (m) */ return P1-b1; @@ -218,7 +218,7 @@ extern int ionocorr(gtime_t time, const nav_t *nav, int sat, const double *pos, trace(4,"ionocorr: time=%s opt=%d sat=%2d pos=%.3f %.3f azel=%.3f %.3f\n", time2str(time,tstr,3),ionoopt,sat,pos[0]*R2D,pos[1]*R2D,azel[0]*R2D, azel[1]*R2D); - + /* SBAS ionosphere model */ if (ionoopt==IONOOPT_SBAS) { if (sbsioncorr(time,nav,pos,azel,ion,var)) return 1; @@ -263,7 +263,7 @@ extern int tropcorr(gtime_t time, const nav_t *nav, const double *pos, trace(4,"tropcorr: time=%s opt=%d pos=%.3f %.3f azel=%.3f %.3f\n", time2str(time,tstr,3),tropopt,pos[0]*R2D,pos[1]*R2D,azel[0]*R2D, azel[1]*R2D); - + /* Saastamoinen model */ if (tropopt==TROPOPT_SAAS||tropopt==TROPOPT_EST||tropopt==TROPOPT_ESTG) { *trp=tropmodel(time,pos,azel,REL_HUMI); @@ -293,16 +293,16 @@ static int rescode(int iter, const obsd_t *obs, int n, const double *rs, for (i=0;i<3;i++) rr[i]=x[i]; dtr=x[3]; - + ecef2pos(rr,pos); trace(3,"rescode: rr=%.3f %.3f %.3f\n",rr[0], rr[1], rr[2]); - + for (i=*ns=0;ielmin) continue; - + if (iter>0) { /* test SNR mask */ if (!snrmask(obs+i,azel+i*2,opt)) continue; - + /* ionospheric correction */ if (!ionocorr(time,nav,sat,pos,azel+i*2,opt->ionoopt,&dion,&vion)) { continue; @@ -329,7 +329,7 @@ static int rescode(int iter, const obsd_t *obs, int n, const double *rs, /* Convert from FREQL1 to freq */ dion*=SQR(FREQL1/freq); vion*=SQR(SQR(FREQL1/freq)); - + /* tropospheric correction */ if (!tropcorr(time,nav,pos,azel+i*2,opt->tropopt,&dtrp,&vtrp)) { continue; @@ -337,12 +337,12 @@ static int rescode(int iter, const obsd_t *obs, int n, const double *rs, } /* pseudorange with code bias correction */ if ((P=prange(obs+i,nav,opt,&vmeas))==0.0) continue; - + /* pseudorange residual */ v[nv]=P-(r+dtr-CLIGHT*dts[i*2]+dion+dtrp); trace(4,"sat=%d: v=%.3f P=%.3f r=%.3f dtr=%.6f dts=%.6f dion=%.3f dtrp=%.3f\n", sat,v[nv],P,r,dtr,dts[i*2],dion,dtrp); - + /* design matrix */ for (j=0;jnx&&vv>chisqr[nv-nx-1]) { @@ -415,11 +415,11 @@ static int estpos(const obsd_t *obs, int n, const double *rs, const double *dts, { double x[NX]={0},dx[NX],Q[NX*NX],*v,*H,*var,sig; int i,j,k,info,stat,nv,ns; - + trace(3,"estpos : n=%d\n",n); - + v=mat(n+NX-3,1); H=mat(NX,n+NX-3); var=mat(n+NX-3,1); - + for (i=0;i<3;i++) x[i]=sol->rr[i]; for (i=0;iqr[5]=(float)Q[2]; /* cov zx */ sol->ns=(uint8_t)ns; sol->age=sol->ratio=0.0; - + /* validate solution */ if ((stat=valsol(azel,vsat,n,opt,v,nv,NX,msg))) { sol->stat=opt->sateph==EPHOPT_SBAS?SOLQ_SBAS:SOLQ_SINGLE; @@ -474,14 +474,14 @@ static int estpos(const obsd_t *obs, int n, const double *rs, const double *dts, } } if (i>=MAXITR) sprintf(msg,"iteration divergent i=%d",i); - + free(v); free(H); free(var); return 0; } /* RAIM FDE (failure detection and exclusion) -------------------------------*/ static int raim_fde(const obsd_t *obs, int n, const double *rs, const double *dts, const double *vare, const int *svh, - const nav_t *nav, const prcopt_t *opt, const ssat_t *ssat, + const nav_t *nav, const prcopt_t *opt, const ssat_t *ssat, sol_t *sol, double *azel, int *vsat, double *resp, char *msg) { obsd_t *obs_e; @@ -489,15 +489,15 @@ static int raim_fde(const obsd_t *obs, int n, const double *rs, char tstr[40],name[8],msg_e[128]; double *rs_e,*dts_e,*vare_e,*azel_e,*resp_e,rms_e,rms=100.0; int i,j,k,nvsat,stat=0,*svh_e,*vsat_e,sat=0; - + trace(3,"raim_fde: %s n=%2d\n",time2str(obs[0].time,tstr,0),n); - + if (!(obs_e=(obsd_t *)malloc(sizeof(obsd_t)*n))) return 0; rs_e = mat(6,n); dts_e = mat(2,n); vare_e=mat(1,n); azel_e=zeros(2,n); - svh_e=imat(1,n); vsat_e=imat(1,n); resp_e=mat(1,n); - + svh_e=imat(1,n); vsat_e=imat(1,n); resp_e=mat(1,n); + for (i=0;irms) continue; - + /* save result */ for (j=k=0;jerr[4]; /* Doppler error (Hz) */ int i,j,nv; - + v=mat(n,1); H=mat(4,n); - + for (i=0;irr,x,azel,vsat,err,v,H))<4) { break; } /* least square estimation */ if (lsq(H,v,4,nv,dx,Q)) break; - + for (j=0;j<4;j++) x[j]+=dx[j]; - + if (norm(dx,4)<1E-6) { trace(3,"estvel : vx=%.3f vy=%.3f vz=%.3f, n=%d\n",x[0],x[1],x[2],n); matcpy(sol->rr+3,x,3,1); @@ -660,12 +660,12 @@ extern int pntpos(const obsd_t *obs, int n, const nav_t *nav, prcopt_t opt_=*opt; double *rs,*dts,*var,*azel_,*resp; int i,stat,vsat[MAXOBS]={0},svh[MAXOBS]; - + char tstr[40]; trace(3,"pntpos : tobs=%s n=%d\n",time2str(obs[0].time,tstr,3),n); - + sol->stat=SOLQ_NONE; - + if (n<=0) { strcpy(msg,"no observation data"); return 0; @@ -673,9 +673,9 @@ extern int pntpos(const obsd_t *obs, int n, const nav_t *nav, sol->time=obs[0].time; msg[0]='\0'; sol->eventime = obs[0].eventime; - + rs=mat(6,n); dts=mat(2,n); var=mat(1,n); azel_=zeros(2,n); resp=mat(1,n); - + if (ssat) { for (i=0;itime,obs,n,nav,opt_.sateph,rs,dts,var,svh); - + /* estimate receiver position and time with pseudorange */ stat=estpos(obs,n,rs,dts,var,svh,nav,&opt_,ssat,sol,azel_,vsat,resp,msg); - + /* RAIM FDE */ if (!stat&&n>=6&&opt->posopt[4]) { stat=raim_fde(obs,n,rs,dts,var,svh,nav,&opt_,ssat,sol,azel_,vsat,resp,msg); diff --git a/src/rtksvr.c b/src/rtksvr.c index 5a5d4648e..f552386bf 100644 --- a/src/rtksvr.c +++ b/src/rtksvr.c @@ -92,11 +92,11 @@ static void writesol(rtksvr_t *svr, int index) solopt_t solopt=solopt_default; uint8_t buff[MAXSOLMSG+1]; int i,n; - + tracet(4,"writesol: index=%d\n",index); - + for (i=0;i<2;i++) { - + if (svr->solopt[i].posf==SOLF_STAT) { /* output solution status */ rtksvrlock(svr); @@ -117,7 +117,7 @@ static void writesol(rtksvr_t *svr, int index) /* output extended solution */ n=outsolexs(buff,&svr->rtk.sol,svr->rtk.ssat,svr->solopt+i); strwrite(svr->stream+i+3,buff,n); - + /* save output buffer */ rtksvrlock(svr); saveoutbuf(svr,buff,n,i); @@ -139,16 +139,16 @@ static void writesol(rtksvr_t *svr, int index) static void update_glofcn(rtksvr_t *svr) { int i,j,sat,frq; - + for (i=0;iraw[j].nav.geph[i].sat!=sat) continue; frq=svr->raw[j].nav.geph[i].frq; } if (frq<-7||frq>6) continue; - + for (j=0;j<3;j++) { if (svr->raw[j].nav.geph[i].sat==sat) continue; svr->raw[j].nav.geph[i].sat=sat; @@ -160,7 +160,7 @@ static void update_glofcn(rtksvr_t *svr) static void update_obs(rtksvr_t *svr, obs_t *obs, int index, int iobs) { int i,n=0,sat,sys; - + if (iobsn;i++) { sat=obs->data[i].sat; @@ -183,7 +183,7 @@ static void update_eph(rtksvr_t *svr, nav_t *nav, int ephsat, int ephset, eph_t *eph1,*eph2,*eph3; geph_t *geph1,*geph2,*geph3; int prn; - + if (satsys(ephsat,&prn)!=SYS_GLO) { if (!svr->navsel||svr->navsel==index+1) { /* svr->nav.eph={current_set1,current_set2,prev_set1,prev_set2} */ @@ -221,7 +221,7 @@ static void update_eph(rtksvr_t *svr, nav_t *nav, int ephsat, int ephset, static void update_sbs(rtksvr_t *svr, sbsmsg_t *sbsmsg, int index) { int i,sbssat=svr->rtk.opt.sbassatsel; - + if (sbsmsg&&(sbssat==sbsmsg->prn||sbssat==0)) { sbsmsg->rcv=index+1; if (svr->nsbsrtcm[index].ssr[i].update) continue; - + /* check consistency between iods of orbit and clock */ if (svr->rtcm[index].ssr[i].iod[0]!=svr->rtcm[index].ssr[i].iod[1]) { continue; } svr->rtcm[index].ssr[i].update=0; - + iode=svr->rtcm[index].ssr[i].iode; sys=satsys(i+1,&prn); - + /* check corresponding ephemeris exists */ if (sys==SYS_GPS||sys==SYS_GAL||sys==SYS_QZS) { if (svr->nav.eph[i ].iode!=iode&& @@ -350,7 +350,7 @@ static void update_svr(rtksvr_t *svr, int ret, obs_t *obs, nav_t *nav, { tracet(4,"updatesvr: ret=%d ephsat=%d ephset=%d index=%d\n",ret,ephsat, ephset,index); - + if (ret==1) { /* observation data */ update_obs(svr,obs,index,iobs); } @@ -383,13 +383,13 @@ static int decoderaw(rtksvr_t *svr, int index) nav_t *nav; sbsmsg_t *sbsmsg=NULL; int i,ret,ephsat,ephset,fobs=0; - + tracet(4,"decoderaw: index=%d\n",index); - + rtksvrlock(svr); - + for (i=0;inb[index];i++) { - + /* input rtcm/receiver raw data from stream */ if (svr->format[index]==STRFMT_RTCM2) { ret=input_rtcm2(svr->rtcm+index,svr->buff[index][i]); @@ -436,9 +436,9 @@ static int decoderaw(rtksvr_t *svr, int index) } } svr->nb[index]=0; - + rtksvrunlock(svr); - + return fobs; } /* decode download file ------------------------------------------------------*/ @@ -453,7 +453,7 @@ static void decodefile(rtksvr_t *svr, int index) } rtksvrlock(svr); - + /* check file path completed */ int nb = svr->nb[index]; if (nb<=2|| @@ -465,11 +465,11 @@ static void decodefile(rtksvr_t *svr, int index) char file[1024]; strncpy(file,(char *)svr->buff[index],nb-2); file[nb-2]='\0'; svr->nb[index]=0; - + rtksvrunlock(svr); - + if (svr->format[index]==STRFMT_SP3) { /* precise ephemeris */ - + /* read sp3 precise ephemeris */ readsp3(file,nav,0); if (nav->ne<=0) { @@ -479,18 +479,18 @@ static void decodefile(rtksvr_t *svr, int index) } /* update precise ephemeris */ rtksvrlock(svr); - + if (svr->nav.peph) free(svr->nav.peph); svr->nav.ne = nav->ne; svr->nav.nemax = nav->nemax; svr->nav.peph=nav->peph; svr->ftime[index]=utc2gpst(timeget()); strcpy(svr->files[index],file); - + rtksvrunlock(svr); } else if (svr->format[index]==STRFMT_RNXCLK) { /* precise clock */ - + /* read rinex clock */ if (readrnxc(file,nav)<=0) { tracet(1,"rinex clock file read error: %s\n",file); @@ -499,14 +499,14 @@ static void decodefile(rtksvr_t *svr, int index) } /* update precise clock */ rtksvrlock(svr); - + if (svr->nav.pclk) free(svr->nav.pclk); svr->nav.nc = nav->nc; svr->nav.ncmax = nav->ncmax; svr->nav.pclk=nav->pclk; svr->ftime[index]=utc2gpst(timeget()); strcpy(svr->files[index],file); - + rtksvrunlock(svr); } free(nav); @@ -517,11 +517,11 @@ static void corr_phase_bias(obsd_t *obs, int n, const nav_t *nav) double freq; uint8_t code; int i,j; - + for (i=0;issr[obs[i].sat-1].pbias[code-1]*freq/CLIGHT; } @@ -532,11 +532,11 @@ static void periodic_cmd(int cycle, const char *cmd, stream_t *stream) const char *p=cmd,*q; char msg[1024],*r; int n,period; - + for (p=cmd;;p=q+1) { for (q=p;;q++) if (*q=='\r'||*q=='\n'||*q=='\0') break; n=(int)(q-p); strncpy(msg,p,n); msg[n]='\0'; - + period=0; if ((r=strrchr(msg,'#'))) { sscanf(r,"# %d",&period); @@ -593,7 +593,7 @@ static void send_nmea(rtksvr_t *svr, uint32_t *tickreset) bl=baseline_len(&svr->rtk); if (bl>=svr->bl_reset&&(int)(tick-*tickreset)>MIN_INT_RESET) { strsendcmd(svr->stream+1,svr->cmd_reset); - + tracet(2,"send reset: bl=%.3f rr=%.3f %.3f %.3f rb=%.3f %.3f %.3f\n", bl,svr->rtk.sol.rr[0],svr->rtk.sol.rr[1],svr->rtk.sol.rr[2], svr->rtk.rb[0],svr->rtk.rb[1],svr->rtk.rb[2]); @@ -631,9 +631,9 @@ static void *rtksvrthread(void *arg) uint8_t *p,*q; char msg[128]; int i,j,n,cycle,cputime; - + tracet(3,"rtksvrthread:\n"); - + obsd_t *data = (obsd_t *)calloc(MAXOBS * 2, sizeof(obsd_t)); if (data == NULL) { trace(1, "rtksvrthread: obsd_t alloc failed\n"); @@ -652,7 +652,7 @@ static void *rtksvrthread(void *arg) tick=tickget(); for (i=0;i<3;i++) { p=svr->buff[i]+svr->nb[i]; q=svr->buff[i]+svr->buffsize; - + /* read receiver raw/rtcm data from input stream */ if ((n=strread(svr->stream+i,p,q-p))<=0) { continue; @@ -660,7 +660,7 @@ static void *rtksvrthread(void *arg) /* write receiver raw/rtcm data to log stream */ strwrite(svr->stream+i+5,p,n); svr->nb[i]+=n; - + /* save peek buffer */ rtksvrlock(svr); n=nbuffsize-svr->npb[i]?n:svr->buffsize-svr->npb[i]; @@ -677,7 +677,7 @@ static void *rtksvrthread(void *arg) else { /* decode receiver raw/rtcm data */ fobs[i]=decoderaw(svr,i); - if (1==i&&svr->rtcm[1].staid>0) sol.refstationid=svr->rtcm[1].staid; + if (1==i&&svr->rtcm[1].staid>0) sol.refstationid=svr->rtcm[1].staid; } } /* averaging single base pos */ @@ -708,13 +708,13 @@ static void *rtksvrthread(void *arg) rtksvrlock(svr); rtkpos(&svr->rtk,obs.data,obs.n,&svr->nav); rtksvrunlock(svr); - + if (svr->rtk.sol.stat!=SOLQ_NONE) { - + /* adjust current time */ tt=(int)(tickget()-tick)/1000.0+DTTOL; timeset(gpst2utc(timeadd(svr->rtk.sol.time,tt))); - + /* write solution */ writesol(svr,i); } @@ -738,7 +738,7 @@ static void *rtksvrthread(void *arg) ticknmea=tick; } if ((cputime=(int)(tickget()-tick))>0) svr->cputime=cputime; - + /* sleep until next cycle */ sleepms(svr->cycle-cputime); } @@ -770,9 +770,9 @@ extern int rtksvrinit(rtksvr_t *svr) geph_t geph0={0,-1}; seph_t seph0={0}; int i,j; - + tracet(3,"rtksvrinit:\n"); - + svr->state=svr->cycle=svr->nmeacycle=svr->nmeareq=0; for (i=0;i<3;i++) svr->nmeapos[i]=0.0; svr->buffsize=0; @@ -795,7 +795,7 @@ extern int rtksvrinit(rtksvr_t *svr) svr->thread=0; svr->cputime=svr->prcout=svr->nave=0; for (i=0;i<3;i++) svr->rb_ave[i]=0.0; - + memset(&svr->nav,0,sizeof(nav_t)); memset(&svr->obs,0,sizeof(svr->obs)); if (!(svr->nav.eph =(eph_t *)malloc(sizeof(eph_t )*MAXSAT*4 ))|| @@ -811,7 +811,7 @@ extern int rtksvrinit(rtksvr_t *svr) svr->nav.n =svr->nav.nmax =MAXSAT *4; svr->nav.ng=svr->nav.ngmax=NSATGLO*2; svr->nav.ns=svr->nav.nsmax=NSATSBS*2; - + for (i=0;i<3;i++) for (j=0;jobs[i][j].data=(obsd_t *)malloc(sizeof(obsd_t)*MAXOBS))) { tracet(1,"rtksvrinit: malloc error\n"); @@ -824,12 +824,12 @@ extern int rtksvrinit(rtksvr_t *svr) memset(svr->rtcm+i,0,sizeof(rtcm_t)); } for (i=0;istream+i); - + for (i=0;i<3;i++) *svr->cmds_periodic[i]='\0'; *svr->cmd_reset='\0'; svr->bl_reset=10.0; rtklib_initlock(&svr->lock); - + return 1; } /* free rtk server ------------------------------------------------------------- @@ -840,7 +840,7 @@ extern int rtksvrinit(rtksvr_t *svr) extern void rtksvrfree(rtksvr_t *svr) { int i,j; - + free(svr->nav.eph ); free(svr->nav.geph); free(svr->nav.seph); @@ -910,10 +910,10 @@ extern int rtksvrstart(rtksvr_t *svr, int cycle, int buffsize, int *strs, { gtime_t time,time0={0}; int i,j,rw; - + tracet(3,"rtksvrstart: cycle=%d buffsize=%d navsel=%d nmeacycle=%d nmeareq=%d\n", cycle,buffsize,navsel,nmeacycle,nmeareq); - + if (svr->state) { sprintf(errmsg,"server already started"); return 0; @@ -931,7 +931,7 @@ extern int rtksvrstart(rtksvr_t *svr, int cycle, int buffsize, int *strs, svr->prcout=0; rtkfree(&svr->rtk); rtkinit(&svr->rtk,prcopt); - + if (prcopt->initrst) { /* init averaging pos by restart */ svr->nave=0; for (i=0;i<3;i++) svr->rb_ave[i]=0.0; @@ -947,15 +947,15 @@ extern int rtksvrstart(rtksvr_t *svr, int cycle, int buffsize, int *strs, for (j=0;j<10;j++) svr->nmsg[i][j]=0; for (j=0;jobs[i][j].n=0; strcpy(svr->cmds_periodic[i],!cmds_periodic[i]?"":cmds_periodic[i]); - + /* initialize receiver raw and rtcm control */ init_raw(svr->raw+i,formats[i]); init_rtcm(svr->rtcm+i); - + /* set receiver and rtcm option */ strcpy(svr->raw [i].opt,rcvopts[i]); strcpy(svr->rtcm[i].opt,rcvopts[i]); - + /* connect dgps corrections */ svr->rtcm[i].dgps=svr->nav.dgps; } @@ -980,10 +980,10 @@ extern int rtksvrstart(rtksvr_t *svr, int cycle, int buffsize, int *strs, for (i=0;inav.eph [i].ttr=time0; for (i=0;inav.geph[i].tof=time0; for (i=0;inav.seph[i].tof=time0; - + /* set monitor stream */ svr->moni=moni; - + /* open input streams */ for (i=0;i<8;i++) { rw=i<3?STR_MODE_R:STR_MODE_W; @@ -1003,7 +1003,7 @@ extern int rtksvrstart(rtksvr_t *svr, int cycle, int buffsize, int *strs, /* sync input streams */ strsync(svr->stream,svr->stream+1); strsync(svr->stream,svr->stream+2); - + /* write start commands to input streams */ for (i=0;i<3;i++) { if (!cmds[i]) continue; @@ -1039,19 +1039,19 @@ extern int rtksvrstart(rtksvr_t *svr, int cycle, int buffsize, int *strs, extern void rtksvrstop(rtksvr_t *svr, const char **cmds) { int i; - + tracet(3,"rtksvrstop:\n"); - + /* write stop commands to input streams */ rtksvrlock(svr); for (i=0;i<3;i++) { if (cmds[i]) strsendcmd(svr->stream+i,cmds[i]); } rtksvrunlock(svr); - + /* stop rtk server */ svr->state=0; - + /* free rtk server thread */ #ifdef WIN32 WaitForSingleObject(svr->thread,10000); @@ -1075,11 +1075,11 @@ extern int rtksvropenstr(rtksvr_t *svr, int index, int str, const char *path, const solopt_t *solopt, const prcopt_t *prcopt) { tracet(3,"rtksvropenstr: index=%d str=%d path=%s\n",index,str,path); - + if (index<3||index>7||!svr->state) return 0; - + rtksvrlock(svr); - + if (svr->stream[index].state>0) { rtksvrunlock(svr); return 0; @@ -1091,7 +1091,7 @@ extern int rtksvropenstr(rtksvr_t *svr, int index, int str, const char *path, } if (index<=4) { svr->solopt[index-3]=*solopt; - + /* write solution header to solution stream */ writesolhead(svr->stream+index,svr->solopt+(index-3),prcopt); } @@ -1109,13 +1109,13 @@ extern int rtksvropenstr(rtksvr_t *svr, int index, int str, const char *path, extern void rtksvrclosestr(rtksvr_t *svr, int index) { tracet(3,"rtksvrclosestr: index=%d\n",index); - + if (index<3||index>7||!svr->state) return; - + rtksvrlock(svr); - + strclose(svr->stream+index); - + rtksvrunlock(svr); } /* get observation data status ------------------------------------------------- @@ -1135,9 +1135,9 @@ extern int rtksvrostat(rtksvr_t *svr, int rcv, gtime_t *time, int *sat, double *az, double *el, int **snr, int *vsat) { int i,j,ns; - + tracet(4,"rtksvrostat: rcv=%d\n",rcv); - + if (!svr->state) return 0; rtksvrlock(svr); ns=svr->obs[rcv][0].n; @@ -1172,9 +1172,9 @@ extern void rtksvrsstat(rtksvr_t *svr, int *sstat, char *msg) { int i; char s[MAXSTRMSG],*p=msg; - + tracet(4,"rtksvrsstat:\n"); - + rtksvrlock(svr); for (i=0;istream+i,s); @@ -1194,17 +1194,17 @@ extern int rtksvrmark(rtksvr_t *svr, const char *name, const char *comment) char buff[MAXSOLMSG+1],tstr[40],*p,*q; double tow,pos[3]; int i,sum,week; - + tracet(4,"rtksvrmark:name=%s comment=%s\n",name,comment); - + if (!svr->state) return 0; - + rtksvrlock(svr); - + time2str(svr->rtk.sol.time,tstr,3); tow=time2gpst(svr->rtk.sol.time,&week); ecef2pos(svr->rtk.sol.rr,pos); - + for (i=0;i<2;i++) { p=buff; if (svr->solopt[i].posf==SOLF_STAT) { From 432d7c777a2a34c90ddeff1c06b178f91a00edcc Mon Sep 17 00:00:00 2001 From: Miquel Garcia Date: Mon, 28 Jul 2025 14:18:09 +0200 Subject: [PATCH 6/6] feat(rtkrcv): add iono model from RTCM SSR These corrections are delivered as Spherical Harmonics. --- app/consapp/rtkrcv/conf/rtk.conf | 4 +- app/consapp/rtkrcv/conf/single.conf | 4 +- app/consapp/rtkrcv/conf/ssr_iono.conf | 113 ++++++++++++++++++++++++++ app/consapp/rtkrcv/gcc/makefile | 8 +- data/config/demo5_m8n_1hz.conf | 2 +- data/config/demo5_m8t_5hz.conf | 2 +- data/config/f9p_ppk.conf | 2 +- data/config/rtknavi_example.conf | 2 +- src/options.c | 2 +- src/pntpos.c | 5 ++ src/ppp.c | 3 + src/rtcm3.c | 3 +- src/rtklib.h | 6 ++ src/rtksvr.c | 9 ++ src/ssr_ion.c | 56 +++++++++++++ 15 files changed, 210 insertions(+), 11 deletions(-) create mode 100644 app/consapp/rtkrcv/conf/ssr_iono.conf diff --git a/app/consapp/rtkrcv/conf/rtk.conf b/app/consapp/rtkrcv/conf/rtk.conf index 90a4e67e8..82c99bb4a 100644 --- a/app/consapp/rtkrcv/conf/rtk.conf +++ b/app/consapp/rtkrcv/conf/rtk.conf @@ -16,7 +16,7 @@ pos1-snrmask_L2 =0,0,0,0,0,0,0,0,0 pos1-snrmask_L5 =0,0,0,0,0,0,0,0,0 pos1-dynamics =off # (0:off,1:on) pos1-tidecorr =off # (0:off,1:on,2:otl) -pos1-ionoopt =brdc # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc,7:qzs-lex,8:stec) +pos1-ionoopt =brdc # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc,7:qzs-lex,8:stec,9:ssr) pos1-tropopt =saas # (0:off,1:saas,2:sbas,3:est-ztd,4:est-ztdgrad,5:ztd) pos1-sateph =brdc # (0:brdc,1:precise,2:brdc+sbas,3:brdc+ssrapc,4:brdc+ssrcom) pos1-posopt1 =on # (0:off,1:on) @@ -116,7 +116,7 @@ file-tempdir = file-geexefile = file-solstatfile = file-tracefile = -# +# inpstr1-type =tcpcli # (0:off,1:serial,2:file,3:tcpsvr,4:tcpcli,7:ntripcli,8:ftp,9:http) inpstr2-type =tcpcli # (0:off,1:serial,2:file,3:tcpsvr,4:tcpcli,7:ntripcli,8:ftp,9:http) diff --git a/app/consapp/rtkrcv/conf/single.conf b/app/consapp/rtkrcv/conf/single.conf index 5fc9c9224..ce9a9f14e 100644 --- a/app/consapp/rtkrcv/conf/single.conf +++ b/app/consapp/rtkrcv/conf/single.conf @@ -16,7 +16,7 @@ pos1-snrmask_L2 =0,0,0,0,0,0,0,0,0 pos1-snrmask_L5 =0,0,0,0,0,0,0,0,0 pos1-dynamics =off # (0:off,1:on) pos1-tidecorr =off # (0:off,1:on,2:otl) -pos1-ionoopt =brdc # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc,7:qzs-lex,8:stec) +pos1-ionoopt =brdc # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc,7:qzs-lex,8:stec,9:ssr) pos1-tropopt =saas # (0:off,1:saas,2:sbas,3:est-ztd,4:est-ztdgrad,5:ztd) pos1-sateph =brdc # (0:brdc,1:precise,2:brdc+sbas,3:brdc+ssrapc,4:brdc+ssrcom) pos1-posopt1 =on # (0:off,1:on) @@ -116,7 +116,7 @@ file-tempdir = file-geexefile = file-solstatfile = file-tracefile = -# +# inpstr1-type =tcpcli # (0:off,1:serial,2:file,3:tcpsvr,4:tcpcli,7:ntripcli,8:ftp,9:http) inpstr2-type =off # (0:off,1:serial,2:file,3:tcpsvr,4:tcpcli,7:ntripcli,8:ftp,9:http) diff --git a/app/consapp/rtkrcv/conf/ssr_iono.conf b/app/consapp/rtkrcv/conf/ssr_iono.conf new file mode 100644 index 000000000..c95622bb1 --- /dev/null +++ b/app/consapp/rtkrcv/conf/ssr_iono.conf @@ -0,0 +1,113 @@ +# This is a sample configuration file to test SSR ionospheric corrections +inpstr1-format =rtcm3 +inpstr1-type =ntripcli +inpstr1-path =igs_username:igs_password@igs-ip.net:2101/ALAC00ESP0 # Contains OBS+EPH + + +outstr1-type =file +outstr1-path = +outstr1-format =4 + +pos1-posmode =single +# pos1-posmode =ppp-kine # Does not seem to converge with the current settings + +pos1-sateph =brdc + +inpstr2-type =ntripcli +inpstr2-format =rtcm3 +# Stream with SSR corrections +# inpstr2-path =igs_username:igs_password@products.igs-ip.net:2101/SSRA00BKG0 +# Stream with ephemeris +inpstr2-path =igs_username:igs_password@products.igs-ip.net:2101/BCEP00BKG0 + +# External Broadcast orbits and clocks to speed-up TTFF +inpstr3-format =rtcm3 +inpstr3-type =ntripcli +# Stream with SSR ionospheric corrections +inpstr3-path =igs_username:igs_password@products.igs-ip.net:2101/IONO00IGS1 + +pos1-ionoopt =ssr +pos1-tropopt =est-ztd + +pos1-frequency =l1+l2+l5 + +pos1-soltype =forward +pos1-elmask =10 +pos1-snrmask_r =off +pos1-snrmask_b =off +pos1-snrmask_L1 =0,0,0,0,0,0,0,0,0 +pos1-snrmask_L2 =0,0,0,0,0,0,0,0,0 +pos1-snrmask_L5 =0,0,0,0,0,0,0,0,0 +pos1-dynamics =on +pos1-tidecorr =on + +pos1-posopt1 =on +pos1-posopt2 =on +pos1-posopt3 =on +pos1-posopt4 =on +pos1-posopt5 =off +pos1-exclsats = +pos1-navsys =63 +pos2-armode =off +pos2-gloarmode =off +pos2-arthres =3 +pos2-arlockcnt =0 +pos2-arelmask =20 +pos2-arminfix =0 +pos2-elmaskhold =0 +pos2-aroutcnt =5 +pos2-maxage =30 +pos2-slipthres =0.05 +pos2-rejionno =30 +pos2-rejgdop =30 +pos2-niter =1 +pos2-baselen =0 +pos2-basesig =0 +out-solformat =xyz +out-outhead =off +out-outopt =off +out-timesys =gpst +out-timeform =hms +out-timendec =3 +out-degform =deg +out-fieldsep = +out-height =geodetic +out-geoid =internal +out-solstatic =all +out-nmeaintv1 =1 +out-nmeaintv2 =1 +out-outstat =1 +stats-eratio1 =300 +stats-eratio2 =300 +stats-errphase =0.003 +stats-errphaseel =0.003 +stats-errphasebl =0 +stats-errdoppler =1 +stats-stdbias =30 +stats-stdiono =0.03 +stats-stdtrop =0.3 +stats-prnaccelh =10 +stats-prnaccelv =10 +stats-prnbias =0.0001 +stats-prniono =0.001 +stats-prntrop =0.0001 +stats-clkstab =5e-12 +ant1-postype =xyz +ant1-antdele =0 +ant1-antdeln =0 +ant1-antdelu =0 +misc-timeinterp =off +misc-sbasatsel =0 +misc-rnxopt1 = +misc-rnxopt2 = + +logstr1-type =off +logstr1-path = +misc-svrcycle =10 +misc-timeout =30000 +misc-reconnect =10000 +misc-nmeacycle =5000 +misc-buffsize =32768 +misc-navmsgsel =all +misc-proxyaddr = +misc-fswapmargin =30 diff --git a/app/consapp/rtkrcv/gcc/makefile b/app/consapp/rtkrcv/gcc/makefile index fdde9a9c7..c970d8dc0 100644 --- a/app/consapp/rtkrcv/gcc/makefile +++ b/app/consapp/rtkrcv/gcc/makefile @@ -16,7 +16,7 @@ rtkrcv : sbas.o stream.o rcvraw.o rtcm.o preceph.o options.o pntpos.o ppp.o rtkrcv : novatel.o ublox.o crescent.o skytraq.o javad.o nvs.o binex.o rtkrcv : rt17.o ephemeris.o rinex.o ionex.o rtcm2.o rtcm3.o rtcm3e.o rtkrcv : tides.o septentrio.o swiftnav.o unicore.o -rtkrcv : sofa.o +rtkrcv : sofa.o ssr_ion.o rtkrcv.o : ../rtkrcv.c @@ -93,6 +93,8 @@ unicore.o: $(SRC)/rcv/unicore.c $(CC) -c $(CFLAGS) $(SRC)/rcv/unicore.c sofa.o : $(SRC)/sofa.c $(CC) -c $(CFLAGS) $(SRC)/sofa.c +ssr_ion.o : $(SRC)/ssr_ion.c + $(CC) -c $(CFLAGS) $(SRC)/ssr_ion.c rtkrcv.o : $(SRC)/rtklib.h ../vt.h vt.o : ../vt.h @@ -143,5 +145,9 @@ test2: test3: ./rtkrcv -o ../rtk_pb.conf +test_ion_ssr: + ./rtkrcv -t 4 -nc -o ../conf/ssr_iono.conf + + clean: rm -f rtkrcv rtkrcv.exe rtkrcv.nav *.o *.out *.trace diff --git a/data/config/demo5_m8n_1hz.conf b/data/config/demo5_m8n_1hz.conf index bed2128d4..fb72d6f83 100644 --- a/data/config/demo5_m8n_1hz.conf +++ b/data/config/demo5_m8n_1hz.conf @@ -9,7 +9,7 @@ pos1-snrmask_b =off # (0:off,1:on) pos1-snrmask_L1 =38,38,38,38,38,38,38,38,38 pos1-dynamics =on # (0:off,1:on) pos1-tidecorr =0 # (1:solid+2:otl+4:spole) -pos1-ionoopt =brdc # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc,7:qzs-lex,8:stec) +pos1-ionoopt =brdc # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc,7:qzs-lex,8:stec,9:ssr) pos1-tropopt =saas # (0:off,1:saas,2:sbas,3:est-ztd,4:est-ztdgrad,5:ztd) pos1-sateph =brdc # (0:brdc,1:precise,2:brdc+sbas,3:brdc+ssrapc,4:brdc+ssrcom) pos1-exclsats = # (prn ...) diff --git a/data/config/demo5_m8t_5hz.conf b/data/config/demo5_m8t_5hz.conf index f325d266b..f6cfb04cc 100644 --- a/data/config/demo5_m8t_5hz.conf +++ b/data/config/demo5_m8t_5hz.conf @@ -9,7 +9,7 @@ pos1-snrmask_b =off # (0:off,1:on) pos1-snrmask_L1 =38,38,38,38,38,38,38,38,38 pos1-dynamics =on # (0:off,1:on) pos1-tidecorr =0 # (1:solid+2:otl+4:spole) -pos1-ionoopt =brdc # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc) +pos1-ionoopt =brdc # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc,9:ssr) pos1-tropopt =saas # (0:off,1:saas,2:sbas,3:est-ztd,4:est-ztdgrad) pos1-sateph =brdc # (0:brdc,1:precise,2:brdc+sbas,3:brdc+ssrapc,4:brdc+ssrcom) pos1-exclsats = # (prn ...) diff --git a/data/config/f9p_ppk.conf b/data/config/f9p_ppk.conf index 464ad142c..dc2c33908 100644 --- a/data/config/f9p_ppk.conf +++ b/data/config/f9p_ppk.conf @@ -11,7 +11,7 @@ pos1-snrmask_L2 =35,35,35,35,35,35,35,35,35 pos1-snrmask_L5 =35,35,35,35,35,35,35,35,35 pos1-dynamics =on # (0:off,1:on) pos1-tidecorr =off # (0:off,1:on,2:otl) -pos1-ionoopt =brdc # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc) +pos1-ionoopt =brdc # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc,9:ssr) pos1-tropopt =saas # (0:off,1:saas,2:sbas,3:est-ztd,4:est-ztdgrad) pos1-sateph =brdc # (0:brdc,1:precise,2:brdc+sbas,3:brdc+ssrapc,4:brdc+ssrcom) pos1-posopt1 =off # (0:off,1:on) diff --git a/data/config/rtknavi_example.conf b/data/config/rtknavi_example.conf index e0b08d63b..47f940677 100644 --- a/data/config/rtknavi_example.conf +++ b/data/config/rtknavi_example.conf @@ -11,7 +11,7 @@ pos1-snrmask_L2 =0,0,0,0,0,0,0,0,0 pos1-snrmask_L5 =0,0,0,0,0,0,0,0,0 pos1-dynamics =on # (0:off,1:on) pos1-tidecorr =0 # (1:solid+2:otl+4:spole) -pos1-ionoopt =brdc # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc) +pos1-ionoopt =brdc # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc,9:ssr) pos1-tropopt =saas # (0:off,1:saas,2:sbas,3:est-ztd,4:est-ztdgrad) pos1-sateph =brdc # (0:brdc,1:precise,2:brdc+sbas,3:brdc+ssrapc,4:brdc+ssrcom) pos1-posopt1 =off # (0:off,1:on) diff --git a/src/options.c b/src/options.c index 46d8038e0..ce3fd9291 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,9:ssr" #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 6e56b1b47..21ebaf370 100644 --- a/src/pntpos.c +++ b/src/pntpos.c @@ -235,6 +235,11 @@ extern int ionocorr(gtime_t time, const nav_t *nav, int sat, const double *pos, *var=SQR(*ion*ERR_BRDCI); return 1; } + /* SSR ionospheric model (spherical harmonics) */ + if (ionoopt==IONOOPT_SSR) { + if (ionssr(time,nav,pos,azel,ion,var)) return 1; + err=1; + } /* GPS broadcast ionosphere model */ if (ionoopt==IONOOPT_BRDC||err==1) { *ion=ionmodel(time,nav->ion_gps,pos,azel); diff --git a/src/ppp.c b/src/ppp.c index ccda4dcb5..3735a045b 100644 --- a/src/ppp.c +++ b/src/ppp.c @@ -909,6 +909,9 @@ static int model_iono(gtime_t time, const double *pos, const double *azel, if (opt->ionoopt==IONOOPT_TEC) { return iontec(time,nav,pos,azel,1,dion,var); } + if (opt->ionoopt==IONOOPT_SSR) { + return ionssr(time,nav,pos,azel,dion,var); + } if (opt->ionoopt==IONOOPT_BRDC) { *dion=ionmodel(time,nav->ion_gps,pos,azel); *var=SQR(*dion*ERR_BRDCI); diff --git a/src/rtcm3.c b/src/rtcm3.c index ce7e7721c..e92a41d93 100644 --- a/src/rtcm3.c +++ b/src/rtcm3.c @@ -1990,6 +1990,7 @@ static int decode_ssr_iono_sph_harmonics(rtcm_t *rtcm, int sys, int subtype) { IDF035 = getbitu(rtcm->buff, pos, IDF035_NBITS); pos += IDF035_NBITS; /* Number of Ionospheric Layers */ adjweek(rtcm, IDF003); + rtcm->ssr_ion.ref_epoch_s = IDF003; rtcm->ssr_ion.update_interval_s = IDF004; rtcm->ssr_ion.n_layers = IDF035 + 1; @@ -2027,7 +2028,7 @@ static int decode_ssr_iono_sph_harmonics(rtcm_t *rtcm, int sys, int subtype) { } } - ret = 10; /* default return value for SSR message */ + ret = RETURN_CODE_DECODE_RTCM_SSR_ION; end: return ret; } diff --git a/src/rtklib.h b/src/rtklib.h index db13759f4..9e3d75b4a 100644 --- a/src/rtklib.h +++ b/src/rtklib.h @@ -432,6 +432,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_SSR 9 /* ionosphere option: SSR model */ #define TROPOPT_OFF 0 /* troposphere option: correction off */ #define TROPOPT_SAAS 1 /* troposphere option: Saastamoinen model */ @@ -575,6 +576,7 @@ extern "C" { #define IDF039_NBITS 16 #define IDF040_NBITS 16 +#define RETURN_CODE_DECODE_RTCM_SSR_ION 11 #ifdef WIN32 #define rtklib_thread_t HANDLE @@ -854,6 +856,7 @@ typedef struct { /* ionosphere spherical harmonics type */ } ionlayer_sphharm_t; typedef struct { + int ref_epoch_s; /* reference epoch time (s) */ int update_interval_s; /* update interval*/ int n_layers; /* number of layers */ ionlayer_sphharm_t ionlayers_sph_harm[MAX_ION_LAYERS]; /* ionosphere spherical harmonics */ @@ -914,6 +917,7 @@ typedef struct { /* navigation data type */ sbsion_t sbsion[MAXBAND+1]; /* SBAS ionosphere corrections */ dgps_t dgps[MAXSAT]; /* DGPS corrections */ ssr_t ssr[MAXSAT]; /* SSR corrections */ + ssr_ion_t ssr_ion; /* SSR ionosphere corrections */ } nav_t; typedef struct { /* station parameter type */ @@ -1617,6 +1621,8 @@ EXPORT int tropcorr(gtime_t time, const nav_t *nav, const double *pos, EXPORT int seliflc(int optnf, int sys); EXPORT int compute_stec_from_spherical_harmonics( const ssr_ion_t* ssr_ion, const int epoch_s, const double* pos, const double* azel, double re_km, double* stec); +EXPORT int ionssr(gtime_t time, const nav_t *nav, const double *pos, + const double *azel, double *delay, double *var); /* antenna models ------------------------------------------------------------*/ EXPORT int readpcv(const char *file, pcvs_t *pcvs); diff --git a/src/rtksvr.c b/src/rtksvr.c index f552386bf..5eb32e54d 100644 --- a/src/rtksvr.c +++ b/src/rtksvr.c @@ -343,6 +343,12 @@ static void update_ssr(rtksvr_t *svr, int index) } svr->nmsg[index][7]++; } + +/* update ssr corrections ----------------------------------------------------*/ +static void update_ssr_ion(rtksvr_t *svr, int index) { + svr->nav.ssr_ion = svr->rtcm[index].ssr_ion; +} + /* update rtk server struct --------------------------------------------------*/ static void update_svr(rtksvr_t *svr, int ret, obs_t *obs, nav_t *nav, int ephsat, int ephset, sbsmsg_t *sbsmsg, int index, @@ -372,6 +378,9 @@ static void update_svr(rtksvr_t *svr, int ret, obs_t *obs, nav_t *nav, else if (ret==10) { /* ssr message */ update_ssr(svr,index); } + else if (ret == RETURN_CODE_DECODE_RTCM_SSR_ION) { + update_ssr_ion(svr, index); + } else if (ret==-1) { /* error */ svr->nmsg[index][9]++; } diff --git a/src/ssr_ion.c b/src/ssr_ion.c index e3e8fef9a..c31324aa9 100644 --- a/src/ssr_ion.c +++ b/src/ssr_ion.c @@ -10,6 +10,11 @@ #include "rtklib.h" #define MIN(x,y) ((x)<(y)?(x):(y)) +#define SQR(x) ((x)*(x)) +#define VAR_NOTEC SQR(30.0) /* variance of no tec */ +#define MIN_EL 0.0 /* min elevation angle (rad) */ +#define MIN_HGT -1000.0 /* min user height (m) */ +#define ERR_SSRION 0.5 /* error scaling of SSR ionosphere corrections (m) */ /* Compute the index for spherical harmonics coefficients */ extern int compute_index_ssr_iono_coeff_arrays(const int n, const int m, const int n_degree) { @@ -199,3 +204,54 @@ extern int compute_stec_from_spherical_harmonics( end: return ret; } + +/* ionosphere model by SSR spherical harmonics data ---------------------------- +* +* args : gtime_t time I time (gpst) +* nav_t *nav I navigation data +* double *pos I receiver position {lat,lon,h} (rad,m) +* double *azel I azimuth/elevation angle {az,el} (rad) +* double *delay O ionospheric delay (L1) (m) +* double *var O ionospheric delay (L1) variance (m^2) +* +* return : status (1:ok,0:error) +* notes : return ok with delay = 0 and var = VAR_NOTEC if el < MIN_EL or h < MIN_HGT +*-----------------------------------------------------------------------------*/ +extern int ionssr(gtime_t time, const nav_t *nav, const double *pos, + const double *azel, double *delay, double *var) +{ + static const double ALPHA = 40.30E16 / FREQL1 / FREQL1; /* tecu->L1 iono (m) */ + + int ret = 0; // Default return value, assume failure + int epoch_s, ref_epoch_s; + double stec_tecu; + double stec_var_tecu; + + char tstr[40]; + trace(3,"ionssr : time=%s pos=%.1f %.1f azel=%.1f %.1f\n",time2str(time, tstr, 0), + pos[0] * R2D, pos[1] * R2D, azel[0] * R2D, azel[1] * R2D); + + *delay = 0.0; + *var = VAR_NOTEC; + + epoch_s = (int)time2gpst(time, NULL); + ref_epoch_s = nav->ssr_ion.ref_epoch_s; + + if (azel[1] < MIN_EL || pos[2] < MIN_HGT || abs(epoch_s - ref_epoch_s) > 7200) { + goto end; + } + + if (compute_stec_from_spherical_harmonics( + &nav->ssr_ion, epoch_s, pos, azel, RE_WGS84 / 1000.0, &stec_tecu) < 0) { + goto end; + } + + // Variance of the STEC, unknown, but set to half the square of the STEC, as done with the BRDC corrections + *delay = ALPHA * stec_tecu; + *var = SQR((*delay) * ERR_SSRION); + + ret = 1; // Successful completion +end: + trace(3,"ionssr : delay=%5.2f std=%5.2f\n", *delay, sqrt(*var)); + return ret; +}