Skip to content

Ssr vtec misc #851

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ourairquality wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Ssr vtec misc #851

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions app/qtapp/rtknavi_qt/mondlg.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -421,7 +421,7 @@ void MonitorDialog::showRtk()
unsigned int nmsg[3][10] = {{0}};
char tstr[40], id[8], s1[40] = "-", s2[40] = "-", s3[40] = "-";
char file[1024] = "";
const QString ionoopt[] = {tr("OFF"), tr("Broadcast"), tr("SBAS"), tr("Dual-Frequency"), tr("Estimate STEC"), tr("IONEX TEC"), tr("QZSS LEX")};
const QString ionoopt[] = {tr("OFF"), tr("Broadcast"), tr("SBAS"), tr("Dual-Frequency"), tr("Estimate STEC"), tr("IONEX TEC"), tr("QZSS LEX"), tr("SSR VTEC")};
const QString tropopt[] = {tr("OFF"), tr("Saastamoinen"), tr("SBAS"), tr("Estimate ZTD"), tr("Estimate ZTD+Grad")};
const QString ephopt [] = {tr("Broadcast"), tr("Precise"), tr("Broadcast+SBAS"), tr("Broadcast+SSR APC"), tr("Broadcast+SSR CoM"), tr("QZSS LEX")};

Expand Down Expand Up @@ -530,7 +530,7 @@ void MonitorDialog::showRtk()
ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1, %2").arg(rtk->opt.dynamics ? tr("ON") : tr("OFF"), rtk->opt.tidecorr ? tr("ON") : tr("OFF")));

ui->tWConsole->item(row, 0)->setText(tr("Ionosphere/Troposphere Model"));
ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1, %2").arg(ionoopt[rtk->opt.ionoopt < 0 || rtk->opt.ionoopt > 6 ? 0 : rtk->opt.ionoopt],
ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1, %2").arg(ionoopt[rtk->opt.ionoopt < 0 || rtk->opt.ionoopt > 7 ? 0 : rtk->opt.ionoopt],
tropopt[rtk->opt.tropopt < 0 || rtk->opt.tropopt > 4 ? 0 : rtk->opt.tropopt]));

ui->tWConsole->item(row, 0)->setText(tr("Satellite Ephemeris"));
Expand Down
2 changes: 1 addition & 1 deletion app/winapp/rtknavi/mondlg.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -372,7 +372,7 @@ void __fastcall TMonitorDialog::ShowRtk(void)
int cputime,nb[3]={0},nmsg[3][10]={{0}},ne;
char tstr[40],*ant,id[8],s1[40]="-",s2[40]="-",s3[40]="-";
char file[1024]="";
const char *ionoopt[]={"OFF","Broadcast","SBAS","Dual-Frequency","Estimate STEC","IONEX TEC","QZSS LEX",""};
const char *ionoopt[]={"OFF","Broadcast","SBAS","Dual-Frequency","Estimate STEC","IONEX TEC","QZSS LEX","SSR VTEC",""};
const char *tropopt[]={"OFF","Saastamoinen","SBAS","Estimate ZTD","Estimate ZTD+Grad",""};
const char *ephopt []={"Broadcast","Precise","Broadcast+SBAS","Broadcat+SSR APC","Broadcast+SSR CoM","QZSS LEX",""};

Expand Down
3 changes: 2 additions & 1 deletion app/winapp/rtknavi/naviopt.dfm
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -238,7 +238,8 @@ object OptDialog: TOptDialog
'Iono-Free LC'
'Estimate STEC'
'IONEX TEC'
'QZSS Broadcast')
'QZSS Broadcast'
'SSR VTEC')
end
object TropOpt: TComboBox
Left = 248
Expand Down
3 changes: 2 additions & 1 deletion app/winapp/rtkpost/postopt.dfm
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -257,7 +257,8 @@ object OptDialog: TOptDialog
'Iono-Free LC'
'Estimate STEC'
'IONEX TEC'
'QZSS Broadcast')
'QZSS Broadcast'
'SSR VTEC')
end
object TropOpt: TComboBox
Left = 248
Expand Down
93 changes: 50 additions & 43 deletions src/ionex.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
#define SQR(x) ((x)*(x))
#define MIN(x,y) ((x)<=(y)?(x):(y))
#define VAR_NOTEC SQR(30.0) /* variance of no tec */
#define VAR_SSR_VTEC SQR(10.0); /* variance of SSR VTEC */
#define VAR_SSR_VTEC SQR(2.0); /* variance of SSR VTEC */
#define MIN_EL 0.0 /* min elevation angle (rad) */
#define MIN_HGT -1000.0 /* min user height (m) */

Expand Down Expand Up @@ -500,9 +500,6 @@ extern int ionvtec(gtime_t time, const nav_t *nav, const double *pos,
const double *azel, double freq, double *delay, double *var)
{
const double re=6371.0; /* fallback single-layer height (km) */
double pppos[2],fs,vtec,P[17][17],cosl[17],sinl[17];
double sinp,cosp,x;
int i,j,k,nmax,mmax;
char tstr[40];

trace(4,"ionvtec: time=%s pos=%.3f %.3f azel=%.3f %.3f\n",
Expand All @@ -516,52 +513,62 @@ extern int ionvtec(gtime_t time, const nav_t *nav, const double *pos,
}
if (pos[2]<MIN_HGT||azel[1]<=0.0) return 0;

vtec=0.0;
for (i=0;i<nav->vtec.nlay;i++) {
/* ionospheric pierce point position */
ionppp(pos,azel,re,nav->vtec.hgt[i],pppos);
double tow = time2gpst(time, NULL), tod = fmod(tow, 86400);
double vtec = 0.0;
for (int i = 0; i < nav->vtec.nlay; i++) {
// Ionospheric pierce point position.
double pppos[3];
double fs = ionppp(pos, azel, re, nav->vtec.hgt[i], pppos);

nmax=nav->vtec.nmax[i];
mmax=nav->vtec.mmax[i];
sinp=sin(pppos[0]);
cosp=cos(pppos[0]);
// Mean sun fixed longitude pierce point phase shifted to the local
// time of maximum TEC, 14:00, and modulo 2 * PI.
double pplon = pppos[1] + (tod - 50400) * PI / 43200;
pplon -= floor(pplon / (2 * PI)) * 2 * PI;

/* fully normalized associated Legendre polynomials */
P[0][0]=1.0;
for (j=1;j<=nmax;j++) {
/* diagonal term */
P[j][j]=sqrt((2*j+1)/(2.0*j))*cosp*P[j-1][j-1];
/* one above diagonal */
P[j][j-1]=sqrt(2*j+1)*sinp*P[j-1][j-1];
}
for (j=2;j<=nmax;j++) {
for (k=0;k<=j-2;k++) {
P[j][k]=sqrt((2*j+1)*(2*j-1)/((double)(j-k)*(j+k)))*sinp*P[j-1][k]
-sqrt((2*j+1)*(j+k-1)*(j-k-1)/((double)(j-k)*(j+k)*(2*j-3)))*P[j-2][k];
int nmax = nav->vtec.nmax[i];
int mmax = nav->vtec.mmax[i];

// Fully normalized associated Legendre polynomials.
double sinp = sin(pppos[0]), cosp = cos(pppos[0]);
double P[17][17];
P[0][0] = 1.0;
P[1][0] = sqrt(3) * sinp;
P[1][1] = sqrt(3) * cosp;
// Diagonal term.
for (int n = 2; n <= nmax; n++) {
for (int m = 0; m < n; m++) {
double anm = sqrt((2 * n + 1) * (2 * n - 1) / ((double)(n - m) * (n + m)));
double bnm = sqrt((2 * n + 1) * (n + m - 1) * (n - m - 1) /
((double)(n - m) * (n + m) * (2 * n - 3)));
P[n][m] = anm * sinp * P[n - 1][m] - bnm * P[n - 2][m];
}
P[n][n] = cosp * sqrt((2 * n + 1) / (2.0 * n)) * P[n - 1][n - 1];
}
/* cosine and sine of longitude multiples */
cosl[0]=1.0; sinl[0]=0.0;
for (j=1;j<=mmax;j++) {
cosl[j]=cos(j*pppos[1]);
sinl[j]=sin(j*pppos[1]);

// Cosine and sine of longitude multiples.
double cosl[17], sinl[17];
cosl[0] = 1.0; sinl[0] = 0.0;
for (int m = 1; m <= mmax; m++) {
cosl[m] = cos(m * pplon);
sinl[m] = sin(m * pplon);
}
/* evaluate spherical harmonic expansion (TECU) */
x=0.0;
for (j=0;j<=nmax;j++) {
x+=nav->vtec.cosC[i][j][0]*P[j][0];
for (k=1;k<=MIN(j,mmax);k++) {
x+=nav->vtec.cosC[i][j][k]*P[j][k]*cosl[k]+
nav->vtec.sinC[i][j][k]*P[j][k]*sinl[k];
// Evaluate spherical harmonic expansion (TECU).
double x = 0.0;
for (int n = 0; n <= nmax; n++) {
x += nav->vtec.cosC[i][n][0] * P[n][0];
for (int m = 1; m <= MIN(n, mmax); m++) {
x += nav->vtec.cosC[i][n][m] * P[n][m] * cosl[m] +
nav->vtec.sinC[i][n][m] * P[n][m] * sinl[m];
}
}
vtec+=x;
// Mapping function.
vtec += fs * x;
*var += SQR(fs) * VAR_SSR_VTEC; // Use fixed variance for now.
}
/* convert TECU to metres on L1, then scale to signal frequency */
fs=ionmapf(pos,azel); /* mapping function */
*delay=40.3E16/SQR(FREQL1)*vtec*SQR(FREQL1/freq)*fs;
*var=VAR_SSR_VTEC; /* use fixed variance for now */
trace(4,"ionvtec: pppos=%.3f %.3f vtec=%.1f delay=%.3f var=%.2f, freq=%.1f\n",
pppos[0]*R2D,pppos[1]*R2D,vtec,*delay,*var,freq);
// Guard against negative vtec.
if (vtec < 0) return 0;
// Convert TECU to metres on L1, then scale to signal frequency.
*delay = vtec * 40.3E16 / SQR(freq);
trace(4,"ionvtec: vtec=%.1f delay=%.3f var=%.2f, freq=%.1f\n", vtec, *delay, *var, freq);
return 1;
}
2 changes: 1 addition & 1 deletion src/options.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ static char snrmask_[NFREQ][1024];
#define MODOPT "0:single,1:dgps,2:kinematic,3:static,4:static-start,5:movingbase,6:fixed,7:ppp-kine,8:ppp-static,9:ppp-fixed"
#define FRQOPT "1:l1,2:l1+l2,3:l1+l2+l5,4:l1+l2+l5+l6"
#define TYPOPT "0:forward,1:backward,2:combined,3:combined-nophasereset"
#define IONOPT "0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc"
#define IONOPT "0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc,7:ssr-vtec"
#define TRPOPT "0:off,1:saas,2:sbas,3:est-ztd,4:est-ztdgrad"
#define EPHOPT "0:brdc,1:precise,2:brdc+sbas,3:brdc+ssrapc,4:brdc+ssrcom"
#define NAVOPT "1:gps+2:sbas+4:glo+8:gal+16:qzs+32:bds+64:navic"
Expand Down
5 changes: 5 additions & 0 deletions src/pntpos.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -221,6 +221,11 @@ extern int ionocorr(gtime_t time, const nav_t *nav, int sat, const double *pos,
if (iontec(time,nav,pos,azel,1,ion,var)) return 1;
err=1;
}
// IONEX VTEC model.
if (ionoopt == IONOOPT_VTEC) {
if (ionvtec(time, nav, pos, azel, FREQL1, ion, var)) return 1;
err=1;
}
/* QZSS broadcast ionosphere model */
if (ionoopt==IONOOPT_QZS&&norm(nav->ion_qzs,8)>0.0) {
*ion=ionmodel(time,nav->ion_qzs,pos,azel);
Expand Down
2 changes: 1 addition & 1 deletion src/ppp.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -713,7 +713,7 @@ static void udiono_ppp(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
int sys=satsys(sat,NULL);
double P0_corr=obs[i].P[0];
double Pf_corr=obs[i].P[f2];
if (&rtk->opt.sateph==EPHOPT_SSRAPC||&rtk->opt.sateph==EPHOPT_SSRCOM) {
if (rtk->opt.sateph == EPHOPT_SSRAPC || rtk->opt.sateph == EPHOPT_SSRCOM) {
/* apply SSR correction */
P0_corr-=nav->ssr[obs->sat-1].cbias[obs[i].code[0]-1];
Pf_corr-=nav->ssr[obs->sat-1].cbias[obs[i].code[f2]-1];
Expand Down
19 changes: 10 additions & 9 deletions src/rtcm3.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2059,8 +2059,8 @@ static int decode_ssr7(rtcm_t *rtcm, int sys, int subtype)
/* decode SSR 8: VTEC ionosphere --------------------------------------------- */
static int decode_ssr8(rtcm_t *rtcm, int subtype)
{
double udint,hgt,cosC,sinC,tow;
int i=0,j,k,l,type,sync,iod,udi,nlay,nmax,mmax,np,offp,qi;
double udint,hgt,cosC,sinC;
int i=0,type,sync,iod,udi,nlay,nmax,mmax,qi;
char *msg,tstr[40];

type=getbitu(rtcm->buff,24,12);
Expand All @@ -2086,7 +2086,7 @@ static int decode_ssr8(rtcm_t *rtcm, int subtype)
return -1;
}

for (j=0;j<nlay&&i+16<=rtcm->len*8;j++) {
for (int j=0;j<nlay&&i+16<=rtcm->len*8;j++) {
hgt =getbitu(rtcm->buff,i,8)*10; i+=8; /* height (km) DF472 */
nmax =getbitu(rtcm->buff,i, 4)+1; i+= 4; /* degree DF473 */
mmax =getbitu(rtcm->buff,i, 4)+1; i+= 4; /* order DF474 */
Expand All @@ -2099,17 +2099,17 @@ static int decode_ssr8(rtcm_t *rtcm, int subtype)
memset(rtcm->nav.vtec.cosC, 0, sizeof(rtcm->nav.vtec.cosC));
memset(rtcm->nav.vtec.sinC, 0, sizeof(rtcm->nav.vtec.sinC));
/* cosine coefficients: for order o=0..mmax, degree n=o..nmax */
for (k=0;k<=mmax;k++) {
for (l=k;l<=nmax;l++) {
for (int m=0;m<=mmax;m++) {
for (int n=m;n<=nmax&&i+16<=rtcm->len*8;n++) {
cosC=getbits(rtcm->buff,i,16); i+=16;
rtcm->nav.vtec.cosC[j][l][k]=cosC/200.0;
rtcm->nav.vtec.cosC[j][n][m]=cosC/200.0;
}
}
/* sine coefficients: for order o=1..mmax, degree n=o..nmax */
for (k=1;k<=mmax;k++) {
for (l=k;l<=nmax;l++) {
for (int m=1;m<=mmax;m++) {
for (int n=m;n<=nmax&&i+16<=rtcm->len*8;n++) {
sinC=getbits(rtcm->buff,i,16); i+=16;
rtcm->nav.vtec.sinC[j][l][k]=sinC/200.0;
rtcm->nav.vtec.sinC[j][n][m]=sinC/200.0;
}
}
}
Expand Down Expand Up @@ -2717,6 +2717,7 @@ static int decode_type4076(rtcm_t *rtcm)
case 125: return decode_ssr3(rtcm,SYS_SBS,subtype);
case 126: return decode_ssr7(rtcm,SYS_SBS,subtype);
case 127: return decode_ssr5(rtcm,SYS_SBS,subtype);
case 201: return decode_ssr8(rtcm, subtype);
}
trace(2,"rtcm3 4076: unsupported message subtype=%d\n",subtype);
return 0;
Expand Down
1 change: 1 addition & 0 deletions src/rtklib.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -431,6 +431,7 @@ extern "C" {
#define IONOOPT_EST 4 /* ionosphere option: estimation */
#define IONOOPT_TEC 5 /* ionosphere option: IONEX TEC model */
#define IONOOPT_QZS 6 /* ionosphere option: QZSS broadcast model */
#define IONOOPT_VTEC 7 /* ionosphere option: IONEX VTEC model */

#define TROPOPT_OFF 0 /* troposphere option: correction off */
#define TROPOPT_SAAS 1 /* troposphere option: Saastamoinen model */
Expand Down