diff --git a/app/consapp/rtkrcv/rtkrcv.c b/app/consapp/rtkrcv/rtkrcv.c index 8dd375ec2..e18c2b0c8 100644 --- a/app/consapp/rtkrcv/rtkrcv.c +++ b/app/consapp/rtkrcv/rtkrcv.c @@ -362,41 +362,52 @@ static int readcmd(const char *file, char *cmd, int type) fclose(fp); return 1; } -/* read antenna file ---------------------------------------------------------*/ -static void readant(vt_t *vt, prcopt_t *opt, nav_t *nav) -{ - const pcv_t pcv0={0}; - pcvs_t pcvr={0},pcvs={0}; - pcv_t *pcv; - gtime_t time=timeget(); - int i; - - trace(3,"readant:\n"); - - opt->pcvr[0]=opt->pcvr[1]=pcv0; - if (!*filopt.rcvantp) return; - - if (readpcv(filopt.rcvantp,&pcvr)) { - for (i=0;i<2;i++) { - if (!*opt->anttype[i]) continue; - if (!(pcv=searchpcv(0,opt->anttype[i],time,&pcvr))) { - vt_printf(vt,"no antenna %s in %s",opt->anttype[i],filopt.rcvantp); - continue; - } - opt->pcvr[i]=*pcv; +/* Read antenna file ---------------------------------------------------------*/ +static void readant(vt_t *vt, prcopt_t *opt, nav_t *nav, pcvs_t *pcvsr) { + trace(3, "readant:\n"); + + const pcv_t pcv0 = {0}; + opt->pcvr[0] = opt->pcvr[1] = pcv0; + + if (*filopt.rcvantp) { + gtime_t time = timeget(); + if (readpcv(filopt.rcvantp, 2, pcvsr)) { + for (int i = 0; i < 2; i++) { + if (!*opt->anttype[i]) continue; + pcv_t *pcv = searchpcv(0, opt->anttype[i], time, NULL, pcvsr); + if (!pcv) { + vt_printf(vt, "no antenna %s in %s", opt->anttype[i], filopt.rcvantp); + continue; } - } - else vt_printf(vt,"antenna file open error %s",filopt.rcvantp); - - if (readpcv(filopt.satantp,&pcvs)) { - for (i=0;ipcvs[i]=*pcv; + copy_pcv(&opt->pcvr[i], pcv); + } + } else + vt_printf(vt, "antenna file open error %s", filopt.rcvantp); + } + + if (*filopt.satantp) { + satsvns_t satsvns = {0}; + if (*filopt.satmeta && !readsinex(filopt.satmeta, &satsvns)) + vt_printf(vt, "error : reading sat meta sinex %s", filopt.satmeta); + + pcvs_t pcvs = {0}; + if (readpcv(filopt.satantp, 1, &pcvs)) { + gtime_t time = timeget(); + for (int i = 0; i < MAXSAT; i++) { + pcv_t *pcv = searchpcv(i + 1, "", time, &satsvns, &pcvs); + if (!pcv) { + char id[8]; + satno2id(i + 1, id); + vt_printf(vt, "no satellite %s pcv in %s", id, filopt.satantp); + continue; } - } - else vt_printf(vt,"antenna file open error %s",filopt.satantp); - - free(pcvr.pcv); free(pcvs.pcv); + copy_pcv(&nav->pcvs[i], pcv); + } + free_pcvs(&pcvs); + } else + vt_printf(vt, "antenna file open error %s", filopt.satantp); + free(satsvns.satsvn); + } } /* start rtk server ----------------------------------------------------------*/ static int startsvr(vt_t *vt) @@ -440,9 +451,9 @@ static int startsvr(vt_t *vt) pos[1]=nmeapos[1]*D2R; pos[2]=nmeapos[2]; pos2ecef(pos,npos); - + /* read antenna file */ - readant(vt,&prcopt,&svr.nav); + readant(vt,&prcopt,&svr.nav,&svr.pcvsr); /* read dcb file */ if (*filopt.dcb) { @@ -523,6 +534,10 @@ static void stopsvr(vt_t *vt) #endif if (solopt[0].geoid>0) closegeoid(); + for (int i = 0; i < 2; i++) free_pcv(&prcopt.pcvr[i]); + free_pcvs(&svr.pcvsr); + for (int i = 0; i < MAXSAT; i++) free_pcv(&svr.nav.pcvs[i]); + vt_printf(vt,"stop rtk server\n"); } /* print time ----------------------------------------------------------------*/ @@ -771,10 +786,10 @@ static void prstatus(vt_t *vt) vt_printf(vt,"%-28s: %d\n","# of average single pos base",nave); vt_printf(vt,"%-28s: %s\n","ant type rover",rtk.opt.pcvr[0].type); del=rtk.opt.antdel[0]; - vt_printf(vt,"%-28s: %.3f %.3f %.3f\n","ant delta rover",del[0],del[1],del[2]); + vt_printf(vt,"%-28s: %.4f %.4f %.4f\n","ant delta rover",del[0],del[1],del[2]); vt_printf(vt,"%-28s: %s\n","ant type base" ,rtk.opt.pcvr[1].type); del=rtk.opt.antdel[1]; - vt_printf(vt,"%-28s: %.3f %.3f %.3f\n","ant delta base",del[0],del[1],del[2]); + vt_printf(vt,"%-28s: %.4f %.4f %.4f\n","ant delta base",del[0],del[1],del[2]); ecef2enu(pos,rtk.rb+3,vel); vt_printf(vt,"%-28s: %.3f,%.3f,%.3f\n","vel enu (m/s) base", vel[0],vel[1],vel[2]); diff --git a/app/qtapp/appcmn_qt/navi_post_opt.cpp b/app/qtapp/appcmn_qt/navi_post_opt.cpp index 356ed0898..9283333ce 100644 --- a/app/qtapp/appcmn_qt/navi_post_opt.cpp +++ b/app/qtapp/appcmn_qt/navi_post_opt.cpp @@ -145,6 +145,7 @@ OptDialog::OptDialog(QWidget *parent, int opts) fileOptions.geoid[0] = '\0'; fileOptions.iono[0] = '\0'; fileOptions.rcvantp[0] = '\0'; + fileOptions.satmeta[0] = '\0'; fileOptions.satantp[0] = '\0'; fileOptions.solstat[0] = '\0'; fileOptions.stapos[0] = '\0'; @@ -200,6 +201,7 @@ OptDialog::OptDialog(QWidget *parent, int opts) fileCompleter->setModel(fileModel); ui->lEStationPositionFile->setCompleter(fileCompleter); ui->lEAntennaPcvFile->setCompleter(fileCompleter); + ui->lESatelliteMetaFile->setCompleter(fileCompleter); ui->lESatellitePcvFile->setCompleter(fileCompleter); ui->lEDCBFile->setCompleter(fileCompleter); ui->lEGeoidDataFile->setCompleter(fileCompleter); @@ -220,7 +222,14 @@ OptDialog::OptDialog(QWidget *parent, int opts) acStationPositionFileView->setToolTip(tr("View File")); acStationPositionFileView->setEnabled(false); - // satllite PCV line edit actions + // Satellite meta data line edit actions + QAction *acSatelliteMetaFileSelect = ui->lESatelliteMetaFile->addAction(QIcon(":/buttons/folder"), QLineEdit::TrailingPosition); + acSatelliteMetaFileSelect->setToolTip(tr("Select File")); + QAction *acSatelliteMetaFileView = ui->lESatelliteMetaFile->addAction(QIcon(":/buttons/doc"), QLineEdit::TrailingPosition); + acSatelliteMetaFileView->setToolTip(tr("View File")); + acSatelliteMetaFileView->setEnabled(false); + + // satellite PCV line edit actions QAction *acSatellitePcvFileSelect = ui->lESatellitePcvFile->addAction(QIcon(":/buttons/folder"), QLineEdit::TrailingPosition); acSatellitePcvFileSelect->setToolTip(tr("Select File")); QAction *acSatellitePcvFileView = ui->lESatellitePcvFile->addAction(QIcon(":/buttons/doc"), QLineEdit::TrailingPosition); @@ -271,6 +280,8 @@ OptDialog::OptDialog(QWidget *parent, int opts) connect(ui->buttonBox, &QDialogButtonBox::rejected, this, &OptDialog::reject); connect(ui->btnLoad, &QPushButton::clicked, this, &OptDialog::loadSettings); connect(ui->btnSave, &QPushButton::clicked, this, &OptDialog::saveSettings); + connect(acSatelliteMetaFileSelect, &QAction::triggered, this, &OptDialog::selectSatelliteMetaFile); + connect(acSatelliteMetaFileView, &QAction::triggered, this, &OptDialog::viewSatelliteMetaFile); connect(acAntennaPcvFileSelect, &QAction::triggered, this, &OptDialog::selectAntennaPcvFile); connect(acAntennaPcvFileView, &QAction::triggered, this, &OptDialog::viewAntennaPcvFile); connect(ui->lEAntennaPcvFile, &QLineEdit::textChanged, this, [acAntennaPcvFileView, this]() @@ -504,6 +515,23 @@ void OptDialog::selectReferencePosition() setPosition(ui->cBReferencePositionType->currentIndex(), edit, p); } //--------------------------------------------------------------------------- +void OptDialog::viewSatelliteMetaFile() +{ + if (ui->lESatelliteMetaFile->text().isEmpty()) return; + + textViewer->read(ui->lESatelliteMetaFile->text()); + + textViewer->show(); +} +//--------------------------------------------------------------------------- +void OptDialog::selectSatelliteMetaFile() +{ + QString filename = QFileDialog::getOpenFileName(this, tr("Satellite Meta Data Sinex File"), ui->lESatelliteMetaFile->text(), tr("Sinex File (*.snx);All (*.*)")); + + if (!filename.isEmpty()) + ui->lESatelliteMetaFile->setText(QDir::toNativeSeparators(filename)); +} +//--------------------------------------------------------------------------- void OptDialog::viewSatellitePcvFile() { if (ui->lESatellitePcvFile->text().isEmpty()) return; @@ -650,13 +678,14 @@ void OptDialog::updateOptions() QLineEdit *editu[] = {ui->lERoverPosition1, ui->lERoverPosition2, ui->lERoverPosition3 }; QLineEdit *editr[] = {ui->lEReferencePosition1, ui->lEReferencePosition2, ui->lEReferencePosition3 }; pcvs_t pcvr; - pcv_t *pcv, pcv0; + pcv_t pcv0; gtime_t time = timeget(); memset(&pcvr, 0, sizeof(pcvs_t)); memset(&pcv0, 0, sizeof(pcv_t)); // file options + strncpy(fileOptions.satmeta, qPrintable(ui->lESatelliteMetaFile->text()), MAXSTRPATH-1); strncpy(fileOptions.satantp, qPrintable(ui->lESatellitePcvFile->text()), MAXSTRPATH-1); strncpy(fileOptions.rcvantp, qPrintable(ui->lEAntennaPcvFile->text()), MAXSTRPATH-1); strncpy(fileOptions.stapos, qPrintable(ui->lEStationPositionFile->text()), MAXSTRPATH-1); @@ -797,24 +826,27 @@ void OptDialog::updateOptions() } processingOptions.pcvr[0] = processingOptions.pcvr[1] = pcv0; // initialize antenna PCV - if ((ui->cBRoverAntennaPcv->isChecked() || ui->cBReferenceAntennaPcv->isChecked()) && !readpcv(fileOptions.rcvantp, &pcvr)) { + if ((ui->cBRoverAntennaPcv->isChecked() || ui->cBReferenceAntennaPcv->isChecked()) && !readpcv(fileOptions.rcvantp, 2, &pcvr)) { QMessageBox::warning(this, tr("Error"), tr("Antenna file read error: \"%1\"").arg(fileOptions.rcvantp)); return; } if (ui->cBRoverAntennaPcv->isChecked() && (processingOptions.anttype[0] != QStringLiteral("*"))) { - if ((pcv = searchpcv(0, processingOptions.anttype[0], time, &pcvr))) - processingOptions.pcvr[0] = *pcv; + pcv_t *pcv = searchpcv(0, processingOptions.anttype[0], time, NULL, &pcvr); + if (pcv) + copy_pcv(&processingOptions.pcvr[0], pcv); else QMessageBox::warning(this, tr("Error"), tr("No rover antenna PCV: \"%1\"").arg(processingOptions.anttype[0])); } if (ui->cBReferenceAntennaPcv->isChecked()&& (processingOptions.anttype[1] != QStringLiteral("*"))) { - if ((pcv = searchpcv(0, processingOptions.anttype[1], time, &pcvr))) - processingOptions.pcvr[1] = *pcv; + pcv_t *pcv = searchpcv(0, processingOptions.anttype[1], time, NULL, &pcvr); + if (pcv) + copy_pcv(&processingOptions.pcvr[1], pcv); else QMessageBox::warning(this, tr("Error"), tr("No reference station antenna PCV: \"%1\"").arg(processingOptions.anttype[1])); } - if (ui->cBRoverAntennaPcv->isChecked() || ui->cBReferenceAntennaPcv->isChecked()) - free(pcvr.pcv); + if (ui->cBRoverAntennaPcv->isChecked() || ui->cBReferenceAntennaPcv->isChecked()) { + free_pcvs(&pcvr); + } fillExcludedSatellites(&processingOptions, ui->lEExcludedSatellites->text()); processingOptions.maxaveep = ui->sBMaxAveEp->value(); processingOptions.initrst = ui->cBInitRestart->isChecked(); @@ -1041,6 +1073,7 @@ void OptDialog::load(const QString &file) ui->sBMaxSolutionStd->setValue(solopt.maxsolstd); // file options + ui->lESatelliteMetaFile->setText(filopt.satmeta); ui->lESatellitePcvFile->setText(filopt.satantp); ui->lEAntennaPcvFile->setText(filopt.rcvantp); ui->lEStationPositionFile->setText(filopt.stapos); @@ -1215,6 +1248,7 @@ void OptDialog::save(const QString &file) strncpy(solOpts.sep, qPrintable(ui->lEFieldSeperator->text()), 63); solOpts.maxsolstd = ui->sBMaxSolutionStd->value(); + strncpy(filopt.satmeta, qPrintable(ui->lESatelliteMetaFile->text()), MAXSTRPATH-1); strncpy(filopt.satantp, qPrintable(ui->lESatellitePcvFile->text()), MAXSTRPATH-1); strncpy(filopt.rcvantp, qPrintable(ui->lEAntennaPcvFile->text()), MAXSTRPATH-1); strncpy(filopt.stapos, qPrintable(ui->lEStationPositionFile->text()), MAXSTRPATH-1); @@ -1388,6 +1422,7 @@ void OptDialog::saveOptions(QSettings &settings) settings.setValue("prcopt/outsingle", processingOptions.outsingle); settings.setValue("solopt/maxsolstd", ui->sBMaxSolutionStd->value()); + settings.setValue("setting/satmetafile", ui->lESatelliteMetaFile->text()); settings.setValue("setting/satpcvfile", ui->lESatellitePcvFile->text()); settings.setValue("setting/antpcvfile", ui->lEAntennaPcvFile->text()); settings.setValue("setting/staposfile", ui->lEStationPositionFile->text()); @@ -1533,6 +1568,7 @@ void OptDialog::loadOptions(QSettings &settings) ui->sBBaselineLen->setValue(settings.value("prcopt/baseline1", 0.0).toDouble()); ui->sBBaselineSig->setValue(settings.value("prcopt/baseline2", 0.0).toDouble()); + ui->lESatelliteMetaFile->setText(settings.value("setting/satmetafile", "").toString()); ui->lESatellitePcvFile->setText(settings.value("setting/satpcvfile", "").toString()); ui->lEAntennaPcvFile->setText(settings.value("setting/antpcvfile", "").toString()); readAntennaList(); @@ -1906,7 +1942,7 @@ void OptDialog::readAntennaList() QString currentRoverAntenna, currentReferenceAntenna; int i; - if (!readpcv(qPrintable(ui->lEAntennaPcvFile->text()), &pcvs)) return; + if (!readpcv(qPrintable(ui->lEAntennaPcvFile->text()), 2, &pcvs)) return; /* Save currently defined antennas */ currentRoverAntenna = ui->cBRoverAntenna->currentText(); @@ -1920,7 +1956,7 @@ void OptDialog::readAntennaList() ui->cBRoverAntenna->addItem("*"); ui->cBReferenceAntenna->addItem("*"); for (int i = 0; i < pcvs.n; i++) { - if (pcvs.pcv[i].sat) continue; + if (pcvs.pcv[i].sat || pcvs.pcv[i].svn) continue; if ((p = strchr(pcvs.pcv[i].type, ' '))) *p = '\0'; if (i > 0 && !strcmp(pcvs.pcv[i].type, pcvs.pcv[i - 1].type)) continue; ui->cBRoverAntenna->addItem(pcvs.pcv[i].type); @@ -1933,7 +1969,7 @@ void OptDialog::readAntennaList() i = ui->cBReferenceAntenna->findText(currentReferenceAntenna); ui->cBReferenceAntenna->setCurrentIndex(i == -1 ? 0 : i); - free(pcvs.pcv); + free_pcvs(&pcvs); } //--------------------------------------------------------------------------- void OptDialog::showKeyDialog() diff --git a/app/qtapp/appcmn_qt/navi_post_opt.h b/app/qtapp/appcmn_qt/navi_post_opt.h index 0903271e8..6979928e2 100644 --- a/app/qtapp/appcmn_qt/navi_post_opt.h +++ b/app/qtapp/appcmn_qt/navi_post_opt.h @@ -79,6 +79,8 @@ protected slots: void selectPanelFont(); void selectSolutionFont(); void selectGeoidDataFile(); + void viewSatelliteMetaFile(); + void selectSatelliteMetaFile(); void viewSatellitePcvFile(); void selectSatellitePcvFile(); void selectLocalDirectory(); diff --git a/app/qtapp/appcmn_qt/navi_post_opt.ui b/app/qtapp/appcmn_qt/navi_post_opt.ui index 12139a7c2..d885f749c 100644 --- a/app/qtapp/appcmn_qt/navi_post_opt.ui +++ b/app/qtapp/appcmn_qt/navi_post_opt.ui @@ -2534,27 +2534,55 @@ &Files - - + + - DCB Data File + Satellite Meta Data File SINEX - lEDCBFile + lESatelliteMetaFile - - + + + + <html><head/><body><p>Satellite meta data, needed for satellite Antex version 2. Antenna parameters file path (option: <span style=" font-style:italic;">file-satmetafile</span>).</p><p>An example of the satellite meta data file is found at &quot;data\igs_satellite_metadata.snx&quot;.</p></body></html> + + + + + - Ocean Loading BL&Q Format + Satellite Antenna PCV File ANTEX/NGS PCV - lEBLQFile + lESatellitePcvFile + + + + + + + <html><head/><body><p>If you use the precise ephemeris or SSR correction, input the ANTEX antenna parameters file path for the satellite antenna PCV (phase center variation) correction (option: <span style=" font-style:italic;">file-satantfile</span>). </p><p>Usually use latest igs08.atx file provided by IGS.</p><p>An example of the ANTEX file is found in the rtklib directory &quot;data\igs08.atx&quot;</p></body></html> + + + Receiver Antenna PCV File ANTEX/NGS PCV + + + + + + + <html><head/><body><p>If you apply the receiver antenna phase center offset and PCV correction, input ANTEX or NGS type antenna parameters file path (option: <span style=" font-style:italic;">file-rcvantfile</span>).</p><p>An example of the antenna parameter file is found at &quot;data\igs08.atx&quot; or &quot;data\ngs_abs.pcv&quot;.</p></body></html> + + + + Qt::Orientation::Vertical @@ -2567,54 +2595,41 @@ - - - - <html><head/><body><p>Input the file path of an ionospheric corrections file (*.yyi)</p></body></html> - - - - - + + - Satellite Antenna PCV File ANTEX/NGS PCV + Geoid Data File - lESatellitePcvFile + lEGeoidDataFile - - - - Ionosphere Data File - - - lEIonosphereFile + + + + <html><head/><body><p>Input the file path of the geoid data file if selecting the external model as Geoid Model (option: <span style=" font-style:italic;">file-geoidfile</span>)</p></body></html> - - + + - FTP/HTTP Local Directory + DCB Data File - lELocalDirectory + lEDCBFile - - - - Geoid Data File - - - lEGeoidDataFile + + + + <html><head/><body><p>Input the file path of DCB correction for PPP in CODE format (option: <span style=" font-style:italic;">file-dcbfile</span>)</p></body></html> - + EOP Data File @@ -2624,59 +2639,61 @@ - - + + + + <html><head/><body><p>Input the file path of an EOP data file (option: <span style=" font-style:italic;">file-eopfile</span>). </p><p>The format of the EOP data file shall be IGS ERP format version 2. </p></body></html> + + + + + - Receiver Antenna PCV File ANTEX/NGS PCV + Ocean Loading BL&Q Format + + + lEBLQFile - + <html><head/><body><p>Input the file path of an OTL coefficients file (option: <span style=" font-style:italic;">file-blqfile</span>). </p><p>The format of the OTL coefficients file is BLQ format.</p></body></html> - - - - Local directory for downloaded data + + + + Ionosphere Data File - - - - - - <html><head/><body><p>Input the file path of an EOP data file (option: <span style=" font-style:italic;">file-eopfile</span>). </p><p>The format of the EOP data file shall be IGS ERP format version 2. </p></body></html> + + lEIonosphereFile - - + + - <html><head/><body><p>Input the file path of DCB correction for PPP in CODE format (option: <span style=" font-style:italic;">file-dcbfile</span>)</p></body></html> + <html><head/><body><p>Input the file path of an ionospheric corrections file (*.yyi)</p></body></html> - - - - <html><head/><body><p>Input the file path of the geoid data file if selecting the external model as Geoid Model (option: <span style=" font-style:italic;">file-geoidfile</span>)</p></body></html> + + + + FTP/HTTP Local Directory - - - - - - <html><head/><body><p>If you apply the receiver antenna phase center offset and PCV correction, input ANTEX or NGS type antenna parameters file path (option: <span style=" font-style:italic;">file-rcvantfile</span>).</p><p>An example of the antenna parameter file is found at &quot;data\igs08.atx&quot; or &quot;data\ngs_abs.pcv&quot;.</p></body></html> + + lELocalDirectory - - + + - <html><head/><body><p>If you use the precise ephemeris or SSR correction, input the ANTEX antenna parameters file path for the satellite antenna PCV (phase center variation) correction (option: <span style=" font-style:italic;">file-satantfile</span>). </p><p>Usually use latest igs08.atx file provided by IGS.</p><p>An example of the ANTEX file is found in the rtklib directory &quot;data\igs08.atx&quot;</p></body></html> + Local directory for downloaded data diff --git a/app/qtapp/rtknavi_qt/mondlg.cpp b/app/qtapp/rtknavi_qt/mondlg.cpp index 589866b31..abc8b9613 100644 --- a/app/qtapp/rtknavi_qt/mondlg.cpp +++ b/app/qtapp/rtknavi_qt/mondlg.cpp @@ -386,7 +386,7 @@ void MonitorDialog::setRtk() int width[] = {500, 500}; ui->tWConsole->setColumnCount(2); - ui->tWConsole->setRowCount(52 + NFREQ*2); + ui->tWConsole->setRowCount(49 + NFREQ + ANTNFREQ * 2); ui->tWConsole->setHorizontalHeaderLabels(header); for (int i = 0; (i < ui->tWConsole->columnCount()) && (i < 2); i++) @@ -471,7 +471,7 @@ void MonitorDialog::showRtk() if (rtk.opt.navsys & SYS_IRN) navsys = navsys + tr("NavIC "); if (rtk.opt.navsys & SYS_SBS) navsys = navsys + tr("SBAS "); - if (ui->tWConsole->rowCount() < 55) return; + if (ui->tWConsole->rowCount() < 49 + NFREQ + ANTNFREQ * 2) return; row = 0; @@ -660,28 +660,30 @@ void MonitorDialog::showRtk() ui->tWConsole->item(row, 0)->setText(tr("Antenna Type Rover")); ui->tWConsole->item(row++, 1)->setText(rtk.opt.pcvr[0].type); - for (j = 0; j < NFREQ; j++) { - off = rtk.opt.pcvr[0].off[j]; - ui->tWConsole->item(row, 0)->setText(tr("Antenna Phase Center L%1 E/N/U Rover").arg(j+1)); - ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1 m, %2 m, %3 m").arg(off[0], 0, 'f', 3).arg(off[1], 0, 'f', 3).arg(off[2], 0, 'f', 3)); + for (j = 0; j < ANTNFREQ; j++) { + double pco[3]; + antpco(&rtk.opt.pcvr[0], antcode2freq(j), pco); + ui->tWConsole->item(row, 0)->setText(tr("Antenna Phase Center Ant %1 E/N/U Rover").arg(antcode2id(j))); + ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1 m, %2 m, %3 m").arg(pco[0], 0, 'f', 4).arg(pco[1], 0, 'f', 4).arg(pco[2], 0, 'f', 4)); } del = rtk.opt.antdel[0]; ui->tWConsole->item(row, 0)->setText(tr("Antenna Delta E/N/U Rover")); - ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1 m, %2 m, %3 m").arg(del[0], 0, 'f', 3).arg(del[1], 0, 'f', 3).arg(del[2], 0, 'f', 3)); + ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1 m, %2 m, %3 m").arg(del[0], 0, 'f', 4).arg(del[1], 0, 'f', 4).arg(del[2], 0, 'f', 4)); ui->tWConsole->item(row, 0)->setText(tr("Antenna Type Base Station")); ui->tWConsole->item(row++, 1)->setText(rtk.opt.pcvr[1].type); - for (j = 0; j < NFREQ; j++) { - off = rtk.opt.pcvr[1].off[0]; - ui->tWConsole->item(row, 0)->setText(tr("Antenna Phase Center L%1 E/N/U Base Station").arg(j+1)); - ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1 m, %2 m, %3 m").arg(off[0], 0, 'f', 3).arg(off[1], 0, 'f', 3).arg(off[2], 0, 'f', 3)); + for (j = 0; j < ANTNFREQ; j++) { + double pco[3]; + antpco(&rtk.opt.pcvr[1], antcode2freq(j), pco); + ui->tWConsole->item(row, 0)->setText(tr("Antenna Phase Center Ant %1 E/N/U Base Station").arg(antcode2id(j))); + ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1 m, %2 m, %3 m").arg(pco[0], 0, 'f', 4).arg(pco[1], 0, 'f', 4).arg(pco[2], 0, 'f', 4)); } del = rtk.opt.antdel[1]; ui->tWConsole->item(row, 0)->setText(tr("Antenna Delta E/N/U Base Station")); - ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1 m, %2 m, %3 m").arg(del[0], 0, 'f', 3).arg(del[1], 0, 'f', 3).arg(del[2], 0, 'f', 3)); + ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1 m, %2 m, %3 m").arg(del[0], 0, 'f', 4).arg(del[1], 0, 'f', 4).arg(del[2], 0, 'f', 4)); ui->tWConsole->item(row, 0)->setText(tr("Precise Ephemeris Time/# of Epoch")); ui->tWConsole->item(row++, 1)->setText(QStringLiteral("%1-%2 (%3)").arg(s1, s2).arg(ne)); @@ -817,7 +819,7 @@ void MonitorDialog::showSat() ui->tWConsole->item(n, j++)->setText(fix == 1 ? tr("FLOAT") : (fix == 2 ? tr("FIX") : (fix == 3 ? tr("HOLD") : tr("-")))); } for (k = 0; k < nfreq; k++) - ui->tWConsole->item(n, j++)->setText(QString::number(ssat->resp[k], 'f', 2)); + ui->tWConsole->item(n, j++)->setText(QString::number(ssat->resp[k], 'f', 3)); for (k = 0; k < nfreq; k++) ui->tWConsole->item(n, j++)->setText(QString::number(ssat->resc[k], 'f', 4)); for (k = 0; k < nfreq; k++) diff --git a/app/qtapp/rtknavi_qt/navimain.cpp b/app/qtapp/rtknavi_qt/navimain.cpp index eb2571c16..90581852f 100644 --- a/app/qtapp/rtknavi_qt/navimain.cpp +++ b/app/qtapp/rtknavi_qt/navimain.cpp @@ -1118,7 +1118,6 @@ void MainWindow::serverStart() char errmsg[20148]; gtime_t time = timeget(); pcvs_t pcvs; - pcv_t *pcv; trace(3, "serverStart\n"); @@ -1131,16 +1130,30 @@ void MainWindow::serverStart() tracelevel(optDialog->solutionOptions.trace); } - if (optDialog->processingOptions.sateph == EPHOPT_PREC || optDialog->processingOptions.sateph == EPHOPT_SSRCOM) { - if (!readpcv(optDialog->fileOptions.satantp, &pcvs)) { + if (!readpcv(optDialog->fileOptions.rcvantp, 2, &rtksvr->pcvsr)) { + ui->lblMessage->setText(tr("Receiver antenna file read error: %1").arg(optDialog->fileOptions.rcvantp)); + return; + } + + if (optDialog->processingOptions.sateph == EPHOPT_PREC || + optDialog->processingOptions.sateph == EPHOPT_SSRCOM || + optDialog->processingOptions.mode >= PMODE_PPP_KINEMA) { + satsvns_t satsvns = {0}; + if (!readsinex(optDialog->fileOptions.satmeta, &satsvns)) + ui->lblMessage->setText(tr("Satellite meta data sinex file read error: %1").arg(optDialog->fileOptions.satmeta)); + + if (!readpcv(optDialog->fileOptions.satantp, 1, &pcvs)) { ui->lblMessage->setText(tr("Satellite antenna file read error: %1").arg(optDialog->fileOptions.satantp)); + free_pcvs(&rtksvr->pcvsr); return; } for (i = 0; i < MAXSAT; i++) { - if (!(pcv = searchpcv(i + 1, "", time, &pcvs))) continue; - rtksvr->nav.pcvs[i] = *pcv; + pcv_t *pcv = searchpcv(i + 1, "", time, &satsvns, &pcvs); + if (!pcv) continue; + copy_pcv(&rtksvr->nav.pcvs[i], pcv); } - free(pcvs.pcv); + free_pcvs(&pcvs); + free(satsvns.satsvn); } for (i = 0; i < 3; i++) streamTypes[i] = streamEnabled[i] ? itype[streamType[i]] : STR_NONE; // input stream @@ -1290,6 +1303,9 @@ void MainWindow::serverStop() } rtksvrstop(rtksvr, (const char **)cmds); + free_pcvs(&rtksvr->pcvsr); + for (int i = 0; i < MAXSAT; i++) free_pcv(&rtksvr->nav.pcvs[i]); + for (i = 0; i < 3; i++) delete[] cmds[i]; ui->btnStart->setVisible(true); diff --git a/app/qtapp/rtkplot_qt/plotdata.cpp b/app/qtapp/rtkplot_qt/plotdata.cpp index a4b872880..2140d39ba 100644 --- a/app/qtapp/rtkplot_qt/plotdata.cpp +++ b/app/qtapp/rtkplot_qt/plotdata.cpp @@ -1213,7 +1213,7 @@ void Plot::updateObservation(int nobs) if (plotOptDialog->getReceiverPosition() == 0) { // single point position opt.err[0] = 900.0; - pntpos(observation.data + i, j - i, &navigation, &opt, &sol, azel, NULL, msg); + pntpos(observation.data + i, j - i, &navigation, &opt, 0, &sol, azel, NULL, msg); matcpy(rr, sol.rr, 3, 1); ecef2pos(rr, pos); } else if (plotOptDialog->getReceiverPosition() == 1) { // lat/lon/height diff --git a/app/winapp/rtknavi/navimain.cpp b/app/winapp/rtknavi/navimain.cpp index e29305fb2..f78d9de12 100644 --- a/app/winapp/rtknavi/navimain.cpp +++ b/app/winapp/rtknavi/navimain.cpp @@ -1128,7 +1128,7 @@ void __fastcall TMainForm::SvrStart(void) FILE *fp; gtime_t time=timeget(); pcvs_t pcvr={0},pcvs={0}; - pcv_t *pcv,pcv0={0}; + pcv_t pcv0={0}; trace(3,"SvrStart\n"); @@ -1173,7 +1173,7 @@ void __fastcall TMainForm::SvrStart(void) PrcOpt.exsats[sat-1]=ex; } } - if ((RovAntPcvF||RefAntPcvF)&&!readpcv(AntPcvFileF.c_str(),&pcvr)) { + if ((RovAntPcvF||RefAntPcvF)&&!readpcv(AntPcvFileF.c_str(),2,&rtksvr.pcvr)) { Message->Caption=s.sprintf("rcv ant file read error %s",AntPcvFileF.c_str()); Message->Hint=Message->Caption; return; @@ -1182,8 +1182,9 @@ void __fastcall TMainForm::SvrStart(void) if (RovAntPcvF) { type=RovAntF.c_str(); - if ((pcv=searchpcv(0,type,time,&pcvr))) { - PrcOpt.pcvr[0]=*pcv; + pcv_t *pcv=searchpcv(0,type,time,NULL,&pcvr); + if (pcv) { + copy_pcv(&PrcOpt.pcvr[0],pcv); } else { Message->Caption=s.sprintf("no antenna pcv %s",type); @@ -1193,8 +1194,9 @@ void __fastcall TMainForm::SvrStart(void) } if (RefAntPcvF) { type=RefAntF.c_str(); - if ((pcv=searchpcv(0,type,time,&pcvr))) { - PrcOpt.pcvr[1]=*pcv; + pcv_t *pcv=searchpcv(0,type,time,NULL,&pcvr); + if (pcv) { + copy_pcv(&PrcOpt.pcvr[1],pcv); } else { Message->Caption=s.sprintf("no antenna pcv %s",type); @@ -1202,20 +1204,26 @@ void __fastcall TMainForm::SvrStart(void) } for (i=0;i<3;i++) PrcOpt.antdel[1][i]=RefAntDel[i]; } - if (RovAntPcvF||RefAntPcvF) { - free(pcvr.pcv); - } - if (PrcOpt.sateph==EPHOPT_PREC||PrcOpt.sateph==EPHOPT_SSRCOM) { - if (!readpcv(SatPcvFileF.c_str(),&pcvs)) { + if (PrcOpt.sateph==EPHOPT_PREC||PrcOpt.sateph==EPHOPT_SSRCOM||PrcOpt.mode>=PMODE_PPP_KINEMA) { + satsvns_t satsvns = {0}; +#if 0 // TODO add the sat meta data file option, for Antex 2 sat support + if (!readsinex(SatMetaFileF.c_str(),&satsvns)) { + Message->Caption=s.sprintf("sat meta data file read error %s",SatMetaFileF.c_str()); + Message->Hint=Message->Caption; + } +#endif + if (!readpcv(SatPcvFileF.c_str(),1,&pcvs)) { Message->Caption=s.sprintf("sat ant file read error %s",SatPcvFileF.c_str()); Message->Hint=Message->Caption; return; } for (i=0;iVisible=true; BtnOpt ->Enabled=true; BtnExit ->Enabled=true; diff --git a/app/winapp/rtknavi/naviopt.cpp b/app/winapp/rtknavi/naviopt.cpp index 68a5d2447..11a22032d 100644 --- a/app/winapp/rtknavi/naviopt.cpp +++ b/app/winapp/rtknavi/naviopt.cpp @@ -1247,14 +1247,14 @@ void __fastcall TOptDialog::ReadAntList(void) pcvs_t pcvs={0}; char *p; - if (!readpcv(AntPcvFile_Text.c_str(),&pcvs)) return; + if (!readpcv(AntPcvFile_Text.c_str(),2,&pcvs)) return; list=new TStringList; list->Add(""); list->Add("*"); for (int i=0;i0&&!strcmp(pcvs.pcv[i].type,pcvs.pcv[i-1].type)) continue; list->Add(pcvs.pcv[i].type); @@ -1262,7 +1262,7 @@ void __fastcall TOptDialog::ReadAntList(void) RovAnt->Items=list; RefAnt->Items=list; - free(pcvs.pcv); + free_pcvs(&pcvs); } //--------------------------------------------------------------------------- void __fastcall TOptDialog::NavSys6Click(TObject *Sender) diff --git a/app/winapp/rtkplot/plotdata.cpp b/app/winapp/rtkplot/plotdata.cpp index c9eaebed0..895689285 100644 --- a/app/winapp/rtkplot/plotdata.cpp +++ b/app/winapp/rtkplot/plotdata.cpp @@ -1223,7 +1223,7 @@ void __fastcall TPlot::UpdateObs(int nobs) char msg[128]; opt.err[0]=900.0; - pntpos(Obs.data+i,j-i,&Nav,&opt,&sol,azel,NULL,msg); + pntpos(Obs.data+i,j-i,&Nav,&opt,0,&sol,azel,NULL,msg); matcpy(rr,sol.rr,3,1); ecef2pos(rr,pos); } diff --git a/app/winapp/rtkpost/postopt.cpp b/app/winapp/rtkpost/postopt.cpp index 34b9a5db6..e3f87f1a2 100644 --- a/app/winapp/rtkpost/postopt.cpp +++ b/app/winapp/rtkpost/postopt.cpp @@ -1079,14 +1079,14 @@ void __fastcall TOptDialog::ReadAntList(void) pcvs_t pcvs={0}; char *p; - if (!readpcv(AntPcvFile_Text.c_str(),&pcvs)) return; + if (!readpcv(AntPcvFile_Text.c_str(),2,&pcvs)) return; list=new TStringList; list->Add(""); list->Add("*"); for (int i=0;i0&&!strcmp(pcvs.pcv[i].type,pcvs.pcv[i-1].type)) continue; list->Add(pcvs.pcv[i].type); @@ -1094,7 +1094,7 @@ void __fastcall TOptDialog::ReadAntList(void) RovAnt->Items=list; RefAnt->Items=list; - free(pcvs.pcv); + free_pcvs(&pcvs); } //--------------------------------------------------------------------------- void __fastcall TOptDialog::BtnHelpClick(TObject *Sender) diff --git a/data/config/demo5_m8n_1hz.conf b/data/config/demo5_m8n_1hz.conf index 70d39364b..1d7a2ae56 100644 --- a/data/config/demo5_m8n_1hz.conf +++ b/data/config/demo5_m8n_1hz.conf @@ -84,7 +84,7 @@ stats-prniono =0.001 # (m) stats-prntrop =0.0001 # (m) stats-prnpos =0 # (m) stats-clkstab =5e-12 # (s/s) -ant1-postype =rinexhead # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm,6:raw) +ant1-postype =rinexhead # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm) ant1-pos1 =0 # (deg|m) ant1-pos2 =0 # (deg|m) ant1-pos3 =0 # (m|m) @@ -92,7 +92,7 @@ ant1-anttype = ant1-antdele =0 # (m) ant1-antdeln =0 # (m) ant1-antdelu =0 # (m) -ant2-postype =rinexhead # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm,6:raw) +ant2-postype =rinexhead # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm) ant2-pos1 =0 # (deg|m) ant2-pos2 =0 # (deg|m) ant2-pos3 =0 # (m|m) diff --git a/data/config/demo5_m8t_5hz.conf b/data/config/demo5_m8t_5hz.conf index dd95e8793..de19477e9 100644 --- a/data/config/demo5_m8t_5hz.conf +++ b/data/config/demo5_m8t_5hz.conf @@ -83,7 +83,7 @@ stats-prniono =0.001 # (m) stats-prntrop =0.0001 # (m) stats-prnpos =0 # (m) stats-clkstab =5e-12 # (s/s) -ant1-postype =rinexhead # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm,6:raw) +ant1-postype =rinexhead # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm ant1-pos1 =0 # (deg|m) ant1-pos2 =0 # (deg|m) ant1-pos3 =0 # (m|m) @@ -91,7 +91,7 @@ ant1-anttype = ant1-antdele =0 # (m) ant1-antdeln =0 # (m) ant1-antdelu =0 # (m) -ant2-postype =rinexhead # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm,6:raw) +ant2-postype =rinexhead # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm) ant2-pos1 =0 # (deg|m) ant2-pos2 =0 # (deg|m) ant2-pos3 =0 # (m|m) diff --git a/data/config/f9p_ppk.conf b/data/config/f9p_ppk.conf index 0b0d14924..464ad142c 100644 --- a/data/config/f9p_ppk.conf +++ b/data/config/f9p_ppk.conf @@ -90,7 +90,7 @@ stats-prniono =0.001 # (m) stats-prntrop =0.0001 # (m) stats-prnpos =0 # (m) stats-clkstab =5e-12 # (s/s) -ant1-postype =llh # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm,6:raw) +ant1-postype =llh # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm) ant1-pos1 =90 # (deg|m) ant1-pos2 =0 # (deg|m) ant1-pos3 =-6335367.6285 # (m|m) @@ -98,7 +98,7 @@ ant1-anttype = ant1-antdele =0 # (m) ant1-antdeln =0 # (m) ant1-antdelu =0 # (m) -ant2-postype =rinexhead # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm,6:raw) +ant2-postype =rinexhead # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm) ant2-pos1 =0 # (deg|m) ant2-pos2 =0 # (deg|m) ant2-pos3 =0 # (m|m) diff --git a/data/config/rtknavi_example.conf b/data/config/rtknavi_example.conf index 7457190fb..9c717f334 100644 --- a/data/config/rtknavi_example.conf +++ b/data/config/rtknavi_example.conf @@ -90,7 +90,7 @@ stats-prniono =0.001 # (m) stats-prntrop =0.0001 # (m) stats-prnpos =0 # (m) stats-clkstab =5e-12 # (s/s) -ant1-postype =llh # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm,6:raw) +ant1-postype =llh # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm) ant1-pos1 =0 # (deg|m) ant1-pos2 =0 # (deg|m) ant1-pos3 =0 # (m|m) @@ -98,7 +98,7 @@ ant1-anttype = ant1-antdele =0 # (m) ant1-antdeln =0 # (m) ant1-antdelu =0 # (m) -ant2-postype =single # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm,6:raw) +ant2-postype =single # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm) ant2-pos1 =0 # (deg|m) ant2-pos2 =0 # (deg|m) ant2-pos3 =0 # (m|m) diff --git a/src/convrnx.c b/src/convrnx.c index 503b45159..3f2af2bc1 100644 --- a/src/convrnx.c +++ b/src/convrnx.c @@ -1218,7 +1218,7 @@ static void setopt_apppos(strfile_t *str, rnxopt_t *opt) prcopt.navsys=opt->navsys; /* point positioning with last obs data */ - if (!pntpos(str->obs->data,str->obs->n,str->nav,&prcopt,&sol,NULL,NULL, + if (!pntpos(str->obs->data,str->obs->n,str->nav,&prcopt,0,&sol,NULL,NULL, msg)) { trace(2,"point position error (%s)\n",msg); return; diff --git a/src/options.c b/src/options.c index 9f6f40075..02ac61d24 100644 --- a/src/options.c +++ b/src/options.c @@ -60,7 +60,7 @@ static char snrmask_[NFREQ][1024]; #define STAOPT "0:all,1:single" #define STSOPT "0:off,1:state,2:residual" #define ARMOPT "0:off,1:continuous,2:instantaneous,3:fix-and-hold" -#define POSOPT "0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm,6:raw" +#define POSOPT "0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm" #define TIDEOPT "0:off,1:on,2:otl" #define PHWOPT "0:off,1:on,2:precise" @@ -184,6 +184,7 @@ EXPORT opt_t sysopts[]={ {"misc-rnxopt2", 2, (void *)prcopt_.rnxopt[1], "" }, {"misc-pppopt", 2, (void *)prcopt_.pppopt, "" }, + {"file-satmetafile", 2, (void *)&filopt_.satmeta, "" }, {"file-satantfile", 2, (void *)&filopt_.satantp, "" }, {"file-rcvantfile", 2, (void *)&filopt_.rcvantp, "" }, {"file-staposfile", 2, (void *)&filopt_.stapos, "" }, diff --git a/src/pntpos.c b/src/pntpos.c index 68f8dab77..db4bc280e 100644 --- a/src/pntpos.c +++ b/src/pntpos.c @@ -94,16 +94,16 @@ static double gettgd(int sat, const nav_t *nav, int type) } } /* test SNR mask -------------------------------------------------------------*/ -static int snrmask(const obsd_t *obs, const double *azel, const prcopt_t *opt) +static int snrmask(const obsd_t *obs, const double *azel, const prcopt_t *opt, int base) { int f2; - if (testsnr(0,0,azel[1],obs->SNR[0]*SNR_UNIT,&opt->snrmask)) { + if (testsnr(base,0,azel[1],obs->SNR[0]*SNR_UNIT,&opt->snrmask)) { return 0; } if (opt->ionoopt==IONOOPT_IFLC) { f2=seliflc(opt->nf,satsys(obs->sat,NULL)); - if (testsnr(0,f2,azel[1],obs->SNR[f2]*SNR_UNIT,&opt->snrmask)) return 0; + if (testsnr(base,f2,azel[1],obs->SNR[f2]*SNR_UNIT,&opt->snrmask)) return 0; } return 1; } @@ -283,68 +283,114 @@ extern int tropcorr(gtime_t time, const nav_t *nav, const double *pos, /* pseudorange residuals -----------------------------------------------------*/ static int rescode(int iter, const obsd_t *obs, int n, const double *rs, const double *dts, const double *vare, const int *svh, - const nav_t *nav, const double *x, const prcopt_t *opt, + const nav_t *nav, const double *x, const prcopt_t *opt, int base, const ssat_t *ssat, double *v, double *H, double *var, double *azel, int *vsat, double *resp, int *ns) { - gtime_t time; - double r,freq,dion=0.0,dtrp=0.0,vmeas,vion=0.0,vtrp=0.0,rr[3],pos[3],dtr,e[3],P; - int i,j,nv=0,sat,sys,mask[NX-3]={0}; + double rr[3]; + for (int i=0;i<3;i++) rr[i]=x[i]; + double dtr=x[3]; - for (i=0;i<3;i++) rr[i]=x[i]; - dtr=x[3]; - + // Geodetic position of the marker. + double pos[3]; ecef2pos(rr,pos); - trace(3,"rescode: rr=%.3f %.3f %.3f\n",rr[0], rr[1], rr[2]); - - for (i=*ns=0;ierr[5]; + *ns=0; + for (int i=0;ielmin) continue; - + double freq=sat2freq(sat,obs[i].code[0],nav); + if (freq==0.0) continue; + + // Receiver antenna phase center offset. + double pco[3]; // ENU + double freq2 = 0.0, C1 = 0.0, C2 = 0.0; // Also needed for PCV + if (opt->ionoopt==IONOOPT_IFLC) { + int f2 = seliflc(opt->nf, sys); + freq2 = sat2freq(sat, obs[i].code[f2], nav); + if (freq2 == 0.0) continue; + // Iono-free LC + C1 = SQR(freq) / (SQR(freq) - SQR(freq2)); + C2 = -SQR(freq2) / (SQR(freq) - SQR(freq2)); + double pco1[3], pco2[3]; + antpco(opt->pcvr + base, freq, pco1); + antpco(opt->pcvr + base, freq2, pco2); + for (int j = 0; j < 3; j++) pco[j] = C1 * pco1[j] + C2 * pco2[j]; + } else { + antpco(opt->pcvr + base, freq, pco); + } + + // Add in the antenna delta, which includes the offset from the marker + // to the antenna reference point plus eccentricities. + for (int j = 0; j < 3; j++) pco[j] += opt->antdel[base][j]; + double dr[3]; + enu2ecef(pos, pco, dr); + double rpc[3]; + for (int j = 0; j < 3; j++) rpc[j] = rr[j] + dr[j]; + // Recalculate the geodetic position, now for the phase center. + double rpc_pos[3]; + ecef2pos(rpc, rpc_pos); + + // Compute geometric distance and azimuth/elevation angle + double e[3], r = geodist(rs + i * 6, rpc, e); + if (r <= 0.0) continue; + if (satazel(rpc_pos,e,azel+i*2)elmin) continue; + + double dion = 0.0, vion = 0.0, dtrp = 0.0, vtrp = 0.0, dant = 0.0; 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)) { + /* Test SNR mask */ + if (!snrmask(obs+i,azel+i*2,opt,base)) continue; + + if (opt->ionoopt!=IONOOPT_IFLC) { + // Ionospheric correction. + if (!ionocorr(time,nav,sat,rpc_pos,azel+i*2,opt->ionoopt,&dion,&vion)) continue; + // Convert from FREQL1 to freq. + dion*=SQR(FREQL1/freq); + vion*=SQR(SQR(FREQL1/freq)); } - if ((freq=sat2freq(sat,obs[i].code[0],nav))==0.0) continue; - /* 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)) { + + /* Tropospheric correction */ + if (!tropcorr(time,nav,rpc_pos,azel+i*2,opt->tropopt,&dtrp,&vtrp)) continue; + + // Receiver antenna phase center variation. + if (opt->posopt[1]) { + if (opt->ionoopt==IONOOPT_IFLC) { + dant = C1 * antpcv(opt->pcvr + base, azel + i * 2, freq) + + C2 * antpcv(opt->pcvr + base, azel + i * 2, freq2); + } else { + dant = antpcv(opt->pcvr + base, azel + i * 2, freq); + } } } - /* 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); - + /* Pseudorange with code bias correction */ + double vmeas, P = prange(obs+i,nav,opt,&vmeas); + if (P==0.0) continue; + + /* Pseudorange residual */ + v[nv]=P-(r+dant+dtr-CLIGHT*dts[i*2]+dion+dtrp); + trace(4,"sat=%d: v=%.3f P=%.3f r=%.3f dant=%.4f dtr=%.6f dts=%.6f dion=%.3f dtrp=%.3f\n", + sat,v[nv],P,r,dant,dtr,dts[i*2],dion,dtrp); + /* design matrix */ - for (j=0;jstat=SOLQ_NONE; @@ -693,11 +740,11 @@ extern int pntpos(const obsd_t *obs, int n, const nav_t *nav, satposs(sol->time,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); + stat=estpos(obs,n,rs,dts,var,svh,nav,&opt_,base,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); + stat=raim_fde(obs,n,rs,dts,var,svh,nav,&opt_,base,ssat,sol,azel_,vsat,resp,msg); } /* estimate receiver velocity with Doppler */ if (stat) { diff --git a/src/postpos.c b/src/postpos.c index 77ddaf883..74e46714d 100644 --- a/src/postpos.c +++ b/src/postpos.c @@ -55,6 +55,7 @@ /* constants/global variables ------------------------------------------------*/ +static satsvns_t satsvns = {0}; // Satellite SINEX meta data SVN to PRN mapping static pcvs_t pcvss={0}; /* satellite antenna parameters */ static pcvs_t pcvsr={0}; /* receiver antenna parameters */ static obs_t obss={0}; /* observation data */ @@ -796,7 +797,7 @@ static void freeobsnav(obs_t *obs, nav_t *nav) free(nav->seph); nav->seph=NULL; nav->ns=nav->nsmax=0; } /* average of single position ------------------------------------------------*/ -static int avepos(double *ra, int rcv, const obs_t *obs, const nav_t *nav, +static int avepos(double *ra, int base, const obs_t *obs, const nav_t *nav, const prcopt_t *opt) { obsd_t data[MAXOBS]; @@ -805,11 +806,11 @@ static int avepos(double *ra, int rcv, const obs_t *obs, const nav_t *nav, int i,j,n=0,m,iobs; char msg[128]; - trace(3,"avepos: rcv=%d obs.n=%d\n",rcv,obs->n); + trace(3,"avepos: base=%d obs.n=%d\n",base,obs->n); for (i=0;i<3;i++) ra[i]=0.0; - for (iobs=0;(m=nextobsf(obs,&iobs,rcv))>0;iobs+=m) { + for (iobs=0;(m=nextobsf(obs,&iobs,base+1))>0;iobs+=m) { for (i=j=0;idata[iobs+i]; @@ -818,7 +819,7 @@ static int avepos(double *ra, int rcv, const obs_t *obs, const nav_t *nav, } if (j<=0||!screent(data[0].time,ts,ts,1.0)) continue; /* only 1 hz */ - if (!pntpos(data,j,nav,opt,&sol,NULL,NULL,msg)) continue; + if (!pntpos(data,j,nav,opt,base,&sol,NULL,NULL,msg)) continue; for (i=0;i<3;i++) ra[i]+=sol.rr[i]; n++; @@ -864,68 +865,65 @@ static int getstapos(const char *file, const char *name, double *r) trace(1,"no station position: %s %s\n",name,file); return 0; } -/* antenna phase center position ---------------------------------------------*/ -static int antpos(prcopt_t *opt, int rcvno, const obs_t *obs, const nav_t *nav, +// Antenna marker position -------------------------------------------------- +// The delta from the marker to the antenna ARP is in opt->antdel[] and this +// comes from an antdel config option or from RINEX or RTCM station info. +static int antpos(prcopt_t *opt, int base, const obs_t *obs, const nav_t *nav, const sta_t *sta, const char *posfile) { - double *rr=rcvno==1?opt->ru:opt->rb,del[3],pos[3],dr[3]={0}; - int i,postype=rcvno==1?opt->rovpos:opt->refpos; - char *name; + trace(3,"antpos : base=%d\n",base); - trace(3,"antpos : rcvno=%d\n",rcvno); - - if (postype==POSOPT_SINGLE) { /* average of single position */ - if (!avepos(rr,rcvno,obs,nav,opt)) { + double *rr=base==0?opt->ru:opt->rb; + int postype=base==0?opt->rovpos:opt->refpos; + if (postype==POSOPT_SINGLE) { // Average of single position. + if (!avepos(rr,base,obs,nav,opt)) { showmsg("error : station pos computation"); return 0; } } - else if (postype==POSOPT_FILE) { /* read from position file */ - name=stas[rcvno==1?0:1].name; + else if (postype==POSOPT_FILE) { // Read from position file. + char *name=stas[base].name; if (!getstapos(posfile,name,rr)) { showmsg("error : no position of %s in %s",name,posfile); return 0; } } - else if (postype==POSOPT_RINEX) { /* get from rinex header */ - if (norm(stas[rcvno==1?0:1].pos,3)<=0.0) { + else if (postype==POSOPT_RINEX) { // Get from rinex header. + if (norm(stas[base].pos,3)<=0.0) { showmsg("error : no position in rinex header"); trace(1,"no position in rinex header\n"); return 0; } - /* add antenna delta unless already done in antpcv() */ - if (!strcmp(opt->anttype[rcvno],"*")) { - if (stas[rcvno==1?0:1].deltype==0) { /* enu */ - for (i=0;i<3;i++) del[i]=stas[rcvno==1?0:1].del[i]; - del[2]+=stas[rcvno==1?0:1].hgt; - ecef2pos(stas[rcvno==1?0:1].pos,pos); - enu2ecef(pos,del,dr); - } else { /* xyz */ - for (i=0;i<3;i++) dr[i]=stas[rcvno==1?0:1].del[i]; - } - } - for (i=0;i<3;i++) rr[i]=stas[rcvno==1?0:1].pos[i]+dr[i]; + for (int i=0;i<3;i++) rr[i]=stas[base].pos[i]; } return 1; } /* open processing session ----------------------------------------------------*/ static int openses(const prcopt_t *popt, const solopt_t *sopt, - const filopt_t *fopt, nav_t *nav, pcvs_t *pcvs, pcvs_t *pcvr) + const filopt_t *fopt, nav_t *nav, satsvns_t *satsvns, pcvs_t *pcvs, pcvs_t *pcvr) { trace(3,"openses :\n"); + // Read satellite meta data, svn to prn mapping. + if (*fopt->satmeta && !readsinex(fopt->satmeta, satsvns)) { + showmsg("error : reading sat meta sinex %s", fopt->satmeta); + trace(1,"sat meta sinex read error: %s\n", fopt->satmeta); + } + /* read satellite antenna parameters */ - if (*fopt->satantp&&!(readpcv(fopt->satantp,pcvs))) { - showmsg("error : no sat ant pcv in %s",fopt->satantp); - trace(1,"sat antenna pcv read error: %s\n",fopt->satantp); - return 0; + if (*fopt->satantp && !readpcv(fopt->satantp,1,pcvs)) { + showmsg("error : no sat ant pcv in %s",fopt->satantp); + trace(1,"sat antenna pcv read error: %s\n",fopt->satantp); + return 0; } + /* read receiver antenna parameters */ - if (*fopt->rcvantp&&!(readpcv(fopt->rcvantp,pcvr))) { + if (*fopt->rcvantp&&!(readpcv(fopt->rcvantp,2,pcvr))) { showmsg("error : no rec ant pcv in %s",fopt->rcvantp); trace(1,"rec antenna pcv read error: %s\n",fopt->rcvantp); return 0; } + /* open geoid data */ if (sopt->geoid>0&&*fopt->geoid) { if (!opengeoid(sopt->geoid,fopt->geoid)) { @@ -935,14 +933,19 @@ static int openses(const prcopt_t *popt, const solopt_t *sopt, } return 1; } -/* close processing session ---------------------------------------------------*/ -static void closeses(nav_t *nav, pcvs_t *pcvs, pcvs_t *pcvr) +/* Close processing session ---------------------------------------------------*/ +static void closeses(nav_t *nav, satsvns_t *satsvns, pcvs_t *pcvs, pcvs_t *pcvr) { trace(3,"closeses:\n"); - /* free antenna parameters */ - free(pcvs->pcv); pcvs->pcv=NULL; pcvs->n=pcvs->nmax=0; - free(pcvr->pcv); pcvr->pcv=NULL; pcvr->n=pcvr->nmax=0; + // Free SINEX SVN to PRN mappings. + free(satsvns->satsvn); + + // Free antenna parameters. + // Clear copies in nav. + for (int i = 0; i < MAXSAT; i++) free_pcv(&nav->pcvs[i]); + free_pcvs(pcvs); + free_pcvs(pcvr); /* close geoid data */ closegeoid(); @@ -954,48 +957,56 @@ static void closeses(nav_t *nav, pcvs_t *pcvs, pcvs_t *pcvr) rtkclosestat(); traceclose(); } -/* set antenna parameters ----------------------------------------------------*/ -static void setpcv(gtime_t time, prcopt_t *popt, nav_t *nav, const pcvs_t *pcvs, - const pcvs_t *pcvr, const sta_t *sta) +/* Set antenna parameters ----------------------------------------------------*/ +static void setpcv(gtime_t time, prcopt_t *popt, nav_t *nav, const satsvns_t *satsvns, + const pcvs_t *pcvs, const pcvs_t *pcvr, const sta_t *sta) { - pcv_t *pcv,pcv0={0}; + pcv_t pcv0={0}; double pos[3],del[3]; int i,j,mode=PMODE_DGPS<=popt->mode&&popt->mode<=PMODE_FIXED; char id[8]; - /* set satellite antenna parameters */ + /* Set satellite antenna parameters */ for (i=0;ipcvs[i]=pcv0; if (!(satsys(i+1,NULL)&popt->navsys)) continue; - if (!(pcv=searchpcv(i+1,"",time,pcvs))) { + const pcv_t *pcv = searchpcv(i + 1, "", time, satsvns, pcvs); + if (!pcv) { satno2id(i+1,id); trace(4,"no satellite antenna pcv: %s\n",id); continue; } - nav->pcvs[i]=*pcv; + copy_pcv(&nav->pcvs[i], pcv); } + for (i=0;i<(mode?2:1);i++) { popt->pcvr[i]=pcv0; - if (!strcmp(popt->anttype[i],"*")) { /* set by station parameters */ + if (!strcmp(popt->anttype[i],"*")) { + // Set by station parameters. This overrides the config antenna + // type and antenna delta, which are ignored in this case. The + // delta is applied even if the antenna type does not map to a + // known PCV, to at least get to then ARP. strcpy(popt->anttype[i],sta[i].antdes); - if (sta[i].deltype==1) { /* xyz */ + if (sta[i].deltype==1) { // XYZ + // Need at least an approx position to map the delta. if (norm(sta[i].pos,3)>0.0) { ecef2pos(sta[i].pos,pos); ecef2enu(pos,sta[i].del,del); for (j=0;j<3;j++) popt->antdel[i][j]=del[j]; } } - else { /* enu */ + else { // ENU for (j=0;j<3;j++) popt->antdel[i][j]=stas[i].del[j]; } } - if (!(pcv=searchpcv(0,popt->anttype[i],time,pcvr))) { + const pcv_t *pcv = searchpcv(0, popt->anttype[i], time, NULL, pcvr); + if (!pcv) { trace(2,"no receiver antenna pcv: %s\n",popt->anttype[i]); *popt->anttype[i]='\0'; continue; } strcpy(popt->anttype[i],pcv->type); - popt->pcvr[i]=*pcv; + copy_pcv(&popt->pcvr[i], pcv); } } /* read ocean tide loading parameters ----------------------------------------*/ @@ -1097,12 +1108,13 @@ static int execses(gtime_t ts, gtime_t te, double ti, const prcopt_t *popt, trace(2,"no erp data %s\n",path); } } - /* read obs and nav data */ + + int err = 0; + + /* Read obs and nav data */ if (!readobsnav(ts,te,ti,infile,index,n,&popt_,&obss,&navs,stas)) { - /* free obs and nav data */ - freeobsnav(&obss, &navs); - free(rtk_ptr); - return 0; + err = 1; + goto done; } /* read dcb parameters from DCB, BIA, BSX files */ @@ -1123,32 +1135,26 @@ static int execses(gtime_t ts, gtime_t te, double ti, const prcopt_t *popt, } /* set antenna parameters */ - if (popt_.mode!=PMODE_SINGLE) { - setpcv(obss.n>0?obss.data[0].time:timeget(),&popt_,&navs,&pcvss,&pcvsr, - stas); - } + setpcv(obss.n>0?obss.data[0].time:timeget(),&popt_,&navs,&satsvns,&pcvss,&pcvsr,stas); /* read ocean tide loading parameters */ if (popt_.mode>PMODE_SINGLE&&*fopt->blq) { readotl(&popt_,fopt->blq,stas); } /* rover/reference fixed position */ if (popt_.mode==PMODE_FIXED) { - if (!antpos(&popt_,1,&obss,&navs,stas,fopt->stapos)) { - freeobsnav(&obss,&navs); - free(rtk_ptr); - return 0; + if (!antpos(&popt_,0,&obss,&navs,stas,fopt->stapos)) { + err = 1; + goto done; } - if (!antpos(&popt_,2,&obss,&navs,stas,fopt->stapos)) { - freeobsnav(&obss,&navs); - free(rtk_ptr); - return 0; + if (!antpos(&popt_,1,&obss,&navs,stas,fopt->stapos)) { + err = 1; + goto done; } } else if (PMODE_DGPS<=popt_.mode&&popt_.mode<=PMODE_STATIC_START) { - if (!antpos(&popt_,2,&obss,&navs,stas,fopt->stapos)) { - freeobsnav(&obss,&navs); - free(rtk_ptr); - return 0; + if (!antpos(&popt_,1,&obss,&navs,stas,fopt->stapos)) { + err = 1; + goto done; } } /* open solution statistics */ @@ -1160,9 +1166,8 @@ static int execses(gtime_t ts, gtime_t te, double ti, const prcopt_t *popt, } /* write header to output file */ if (flag&&!outhead(outfile,infile,n,&popt_,sopt)) { - freeobsnav(&obss,&navs); - free(rtk_ptr); - return 0; + err = 1; + goto done; } /* name time events file */ namefiletm(outfiletm,outfile); @@ -1237,10 +1242,12 @@ static int execses(gtime_t ts, gtime_t te, double ti, const prcopt_t *popt, free(rbb); } /* free rtk, obs and nav data */ +done: + for (int i = 0; i < 2; i++) free_pcv(&popt_.pcvr[i]); free(rtk_ptr); freeobsnav(&obss,&navs); - return aborts?1:0; + return (!err && aborts)?1:0; } /* execute processing session for each rover ---------------------------------*/ static int execses_r(gtime_t ts, gtime_t te, double ti, const prcopt_t *popt, @@ -1409,18 +1416,18 @@ extern int postpos(gtime_t ts, gtime_t te, double ti, double tu, trace(3,"postpos : ti=%.0f tu=%.0f n=%d outfile=%s\n",ti,tu,n,outfile); /* open processing session */ - if (!openses(popt,sopt,fopt,&navs,&pcvss,&pcvsr)) return -1; + if (!openses(popt,sopt,fopt,&navs,&satsvns,&pcvss,&pcvsr)) return -1; if (ts.time!=0&&te.time!=0&&tu>=0.0) { if (timediff(te,ts)<0.0) { showmsg("error : no period"); - closeses(&navs,&pcvss,&pcvsr); + closeses(&navs,&satsvns,&pcvss,&pcvsr); return 0; } for (i=0;i=0;i--) free(ifile[i]); - closeses(&navs,&pcvss,&pcvsr); + closeses(&navs,&satsvns,&pcvss,&pcvsr); return -1; } } @@ -1504,7 +1511,7 @@ extern int postpos(gtime_t ts, gtime_t te, double ti, double tu, base); } /* close processing session */ - closeses(&navs,&pcvss,&pcvsr); + closeses(&navs,&satsvns,&pcvss,&pcvsr); return stat; } diff --git a/src/ppp.c b/src/ppp.c index b65ea556c..5459859e8 100644 --- a/src/ppp.c +++ b/src/ppp.c @@ -402,55 +402,72 @@ static double mwmeas(const obsd_t *obs, const nav_t *nav) (freq1*obs->P[0]+freq2*obs->P[1])/(freq1+freq2); } /* antenna corrected measurements --------------------------------------------*/ -static void corr_meas(const obsd_t *obs, const nav_t *nav, const double *azel, - const prcopt_t *opt, const double *dantr, - const double *dants, double phw, double *L, double *P, - double *Lc, double *Pc) -{ - double freq[NFREQ]={0},C1,C2; - int i,ix=0,frq,frq2,bias_ix,sys=satsys(obs->sat,NULL); - - for (i=0;inf;i++) { - L[i]=P[i]=0.0; - /* skip if low SNR or missing observations */ - freq[i]=sat2freq(obs->sat,obs->code[i],nav); - if (freq[i]==0.0||obs->L[i]==0.0||obs->P[i]==0.0) continue; - if (testsnr(0,0,azel[1],obs->SNR[i]*SNR_UNIT,&opt->snrmask)) continue; - - /* antenna phase center and phase windup correction */ - L[i]=obs->L[i]*CLIGHT/freq[i]-dants[i]-dantr[i]-phw*CLIGHT/freq[i]; - P[i]=obs->P[i] -dants[i]-dantr[i]; - - if (opt->sateph==EPHOPT_SSRAPC||opt->sateph==EPHOPT_SSRCOM) { - /* select SSR code correction based on code */ - if (sys==SYS_GPS) - ix=(i==0?CODE_L1W-1:CODE_L2W-1); - else if (sys==SYS_GLO) - ix=(i==0?CODE_L1P-1:CODE_L2P-1); - else if (sys==SYS_GAL) - ix=(i==0?CODE_L1X-1:CODE_L7X-1); - /* apply SSR correction */ - P[i]+=(nav->ssr[obs->sat-1].cbias[obs->code[i]-1]-nav->ssr[obs->sat-1].cbias[ix]); - } - else { /* apply code bias corrections from file */ - if (sys==SYS_GAL&&(i==1||i==2)) frq=3-i; /* GAL biases are L1/L5 */ - else frq=i; /* other biases are L1/L2 */ - if (frq>=MAX_CODE_BIAS_FREQS) continue; /* only 2 freqs per system supported in code bias table */ - bias_ix=code2bias_ix(sys,obs->code[i]); /* look up bias index in table */ - if (bias_ix>0) { /* 0=ref code */ - P[i]+=nav->cbias[obs->sat-1][frq][bias_ix-1]; /* code bias */ - } - } - } - /* choose freqs for iono-free LC */ - *Lc=*Pc=0.0; - frq2=L[1]==0?2:1; /* if L[1]==0, try L[2] */ - if (freq[0]==0.0||freq[frq2]==0.0) return; - C1= SQR(freq[0])/(SQR(freq[0])-SQR(freq[frq2])); - C2=-SQR(freq[frq2])/(SQR(freq[0])-SQR(freq[frq2])); - - if (L[0]!=0.0&&L[frq2]!=0.0) *Lc=C1*L[0]+C2*L[frq2]; - if (P[0]!=0.0&&P[frq2]!=0.0) *Pc=C1*P[0]+C2*P[frq2]; +static void corr_meas(const obsd_t *obs, const nav_t *nav, const double *azel, const prcopt_t *opt, + int dantrp, int dantsp, double nadir, double phw, double *L, double *P, + double *Lc, double *Pc) { + double freq[NFREQ] = {0}; + int sys = satsys(obs->sat, NULL); + + for (int i = 0; i < opt->nf; i++) { + L[i] = P[i] = 0.0; + /* Skip if low SNR or missing observations */ + int code = obs->code[i]; + freq[i] = sat2freq(obs->sat, code, nav); + if (freq[i] == 0.0 || obs->L[i] == 0.0 || obs->P[i] == 0.0) continue; + if (testsnr(0, 0, azel[1], obs->SNR[i] * SNR_UNIT, &opt->snrmask)) continue; + + /* Satellite and receiver antenna model */ + double dantr = 0.0, dants = 0.0; + if (dantsp) { + const pcv_t *pcv = nav->pcvs + obs->sat - 1; + { + char id[8]; + satno2id(obs->sat, id); + if (pcv->sat != obs->sat) + fprintf(stderr, "XX pcv sat mismatch sat=%d %s expected pcv sat=%d\n", obs->sat, id, pcv->sat); + } + dants = antmodel_s(pcv, nadir, freq[i]); + } + if (dantrp) dantr = antmodel(opt->pcvr, opt->antdel[0], azel, opt->posopt[1], freq[i]); + + /* Antenna phase center and phase windup correction */ + L[i] = obs->L[i] * CLIGHT / freq[i] - dants - dantr - phw * CLIGHT / freq[i]; + P[i] = obs->P[i] - dants - dantr; + + if (opt->sateph == EPHOPT_SSRAPC || opt->sateph == EPHOPT_SSRCOM) { + /* Select SSR code correction based on code */ + int ix = 0; + if (sys == SYS_GPS) + ix = (i == 0 ? CODE_L1W - 1 : CODE_L2W - 1); + else if (sys == SYS_GLO) + ix = (i == 0 ? CODE_L1P - 1 : CODE_L2P - 1); + else if (sys == SYS_GAL) + ix = (i == 0 ? CODE_L1X - 1 : CODE_L7X - 1); + /* Apply SSR correction */ + P[i] += (nav->ssr[obs->sat - 1].cbias[code - 1] - nav->ssr[obs->sat - 1].cbias[ix]); + } else { /* Apply code bias corrections from file */ + int frq; + if (sys == SYS_GAL && (i == 1 || i == 2)) + frq = 3 - i; /* GAL biases are L1/L5 */ + else + frq = i; /* Other biases are L1/L2 */ + if (frq >= MAX_CODE_BIAS_FREQS) + continue; /* Only 2 freqs per system supported in code bias table */ + int bias_ix = code2bias_ix(sys, code); /* Look up bias index in table */ + if (bias_ix > 0) { /* 0=ref code */ + P[i] += nav->cbias[obs->sat - 1][frq][bias_ix - 1]; /* Code bias */ + } + } + } + /* Choose freqs for iono-free LC */ + *Lc = *Pc = 0.0; + int frq2 = L[1] == 0 ? 2 : 1; /* If L[1]==0, try L[2] */ + if (freq[0] == 0.0 || freq[frq2] == 0.0) return; + double C1 = SQR(freq[0]) / (SQR(freq[0]) - SQR(freq[frq2])); + double C2 = -SQR(freq[frq2]) / (SQR(freq[0]) - SQR(freq[frq2])); + + if (L[0] != 0.0 && L[frq2] != 0.0) *Lc = C1 * L[0] + C2 * L[frq2]; + if (P[0] != 0.0 && P[frq2] != 0.0) *Pc = C1 * P[0] + C2 * P[frq2]; } /* detect cycle slip by LLI --------------------------------------------------*/ static void detslp_ll(rtk_t *rtk, const obsd_t *obs, int n) @@ -717,7 +734,7 @@ static void uddcb_ppp(rtk_t *rtk) static void udbias_ppp(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav) { double L[NFREQ],P[NFREQ],Lc,Pc,bias[MAXOBS],offset=0.0,pos[3]={0}; - double freq1,freq2,ion,dantr[NFREQ]={0},dants[NFREQ]={0}; + double freq1,freq2,ion; int i,j,k,f,sat,slip[MAXOBS]={0},clk_jump=0; trace(3,"udbias : n=%d\n",n); @@ -752,7 +769,7 @@ static void udbias_ppp(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav) for (i=k=0;iopt); - corr_meas(obs+i,nav,rtk->ssat[sat-1].azel,&rtk->opt,dantr,dants, + corr_meas(obs+i,nav,rtk->ssat[sat-1].azel,&rtk->opt,0,0,0.0, 0.0,L,P,&Lc,&Pc); bias[i]=0.0; @@ -830,24 +847,23 @@ static void udstate_ppp(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav) /* temporal update of phase-bias */ udbias_ppp(rtk,obs,n,nav); } -/* satellite antenna phase center variation ----------------------------------*/ -static void satantpcv(const double *rs, const double *rr, const pcv_t *pcv, - double *dant) -{ - double ru[3],rz[3],eu[3],ez[3],nadir,cosa; - int i; - - for (i=0;i<3;i++) { - ru[i]=rr[i]-rs[i]; - rz[i]=-rs[i]; - } - if (!normv3(ru,eu)||!normv3(rz,ez)) return; - - cosa=dot3(eu,ez); - cosa=cosa<-1.0?-1.0:(cosa>1.0?1.0:cosa); - nadir=acos(cosa); +/* Satellite antenna phase center variation setup ----------------------------*/ +static int satantpcv(const double *rs, const double *rr, double *nadir) { + double ru[3], rz[3]; + for (int i = 0; i < 3; i++) { + ru[i] = rr[i] - rs[i]; + rz[i] = -rs[i]; + } + double eu[3], ez[3]; + if (!normv3(ru, eu) || !normv3(rz, ez)) { + *nadir = 0.0; + return 0; + } - antmodel_s(pcv,nadir,dant); + double cosa = dot3(eu, ez); + cosa = cosa < -1.0 ? -1.0 : (cosa > 1.0 ? 1.0 : cosa); + *nadir = acos(cosa); + return 1; } /* precise tropospheric model ------------------------------------------------*/ static double trop_model_prec(gtime_t time, const double *pos, @@ -938,7 +954,6 @@ static int ppp_res(int post, const obsd_t *obs, int n, const double *rs, prcopt_t *opt=&rtk->opt; double y,r,cdtr,bias,rr[3],pos[3],e[3],dtdx[3],L[NFREQ],P[NFREQ],Lc,Pc; double var[MAXOBS*2*NFREQ],dtrp=0.0,dion=0.0,vart=0.0,vari=0.0,dcb,freq; - double dantr[NFREQ]={0},dants[NFREQ]={0}; double ve[MAXOBS*2*NFREQ]={0},vmax=0; char str[40]; int ne=0,obsi[MAXOBS*2*NFREQ]={0},frqi[MAXOBS*2*NFREQ],maxobs,maxfrq,rej; @@ -969,17 +984,18 @@ static int ppp_res(int post, const obsd_t *obs, int n, const double *rs, !model_iono(obs[i].time,pos,azel+i*2,opt,sat,x,nav,&dion,&vari)) { continue; } - /* satellite and receiver antenna model */ - if (opt->posopt[0]) satantpcv(rs+i*6,rr,nav->pcvs+sat-1,dants); - antmodel(opt->pcvr,opt->antdel[0],azel+i*2,opt->posopt[1],dantr); + /* Satellite and receiver antenna model */ + double nadir = 0.0; + int dantsp = 0; + if (opt->posopt[0]) dantsp = satantpcv(rs + i * 6, rr, &nadir); /* phase windup model */ if (!model_phw(rtk->sol.time,sat,nav->pcvs[sat-1].type, opt->posopt[2]?2:0,rs+i*6,rr,&rtk->ssat[sat-1].phw)) { continue; } - /* corrected phase and code measurements */ - corr_meas(obs+i,nav,azel+i*2,&rtk->opt,dantr,dants, + /* Corrected phase and code measurements */ + corr_meas(obs+i,nav,azel+i*2,&rtk->opt,1,dantsp,nadir, rtk->ssat[sat-1].phw,L,P,&Lc,&Pc); /* stack phase and code residuals {L1,P1,L2,P2,...} */ diff --git a/src/preceph.c b/src/preceph.c index 0ff919545..9e5fb19be 100644 --- a/src/preceph.c +++ b/src/preceph.c @@ -353,23 +353,24 @@ extern void readsp3(const char *file, nav_t *nav, int opt) * return : status (1:ok,0:error) * notes : only support antex format for the antenna parameter file *-----------------------------------------------------------------------------*/ -extern int readsap(const char *file, gtime_t time, nav_t *nav) -{ - pcvs_t pcvs={0}; - pcv_t pcv0={0},*pcv; - int i; - - char tstr[40]; - trace(3,"readsap : file=%s time=%s\n",file,time2str(time,tstr,0)); +extern int readsap(const char *file, const gtime_t time, const satsvns_t *satsvns, nav_t *nav) { + char tstr[40]; + trace(3, "readsap : file=%s time=%s\n", file, time2str(time, tstr, 0)); - if (!readpcv(file,&pcvs)) return 0; + pcvs_t pcvs = {0}; + if (!readpcv(file, 1, &pcvs)) return 0; - for (i=0;ipcvs[i]=pcv?*pcv:pcv0; + for (int i = 0; i < MAXSAT; i++) { + pcv_t *pcv = searchpcv(i + 1, "", time, satsvns, &pcvs); + if (pcv) + copy_pcv(&nav->pcvs[i], pcv); + else { + pcv_t pcv0 = {0}; + nav->pcvs[i] = pcv0; } - free(pcvs.pcv); - return 1; + } + free_pcvs(&pcvs); + return 1; } /* read DCB parameters from DCB file ------------------------------------------- * - supports satellite and receiver biases @@ -686,81 +687,85 @@ static int pephclk(gtime_t time, int sat, const nav_t *nav, double *dts, if (varc) *varc=SQR(std); return 1; } -/* satellite antenna phase center offset --------------------------------------- -* compute satellite antenna phase center offset in ecef -* args : gtime_t time I time (gpst) -* double *rs I satellite position and velocity (ecef) -* {x,y,z,vx,vy,vz} (m|m/s) -* int sat I satellite number -* nav_t *nav I navigation data -* double *dant O satellite antenna phase center offset (ecef) -* {dx,dy,dz} (m) (iono-free LC value) -* return : none -* notes : iono-free LC frequencies defined as follows: -* GPS/QZSS : L1-L2 -* GLONASS : G1-G2 -* Galileo : E1-E5b -* BDS : B1I-B2I -* NavIC : L5-S -*-----------------------------------------------------------------------------*/ -extern void satantoff(gtime_t time, const double *rs, int sat, const nav_t *nav, - double *dant) -{ - const pcv_t *pcv=nav->pcvs+sat-1; - double ex[3],ey[3],ez[3],es[3],r[3],rsun[3],gmst,erpv[5]={0},freq[2]; - double C1,C2,dant1,dant2; - int i,sys; - - char tstr[40]; - trace(4,"satantoff: time=%s sat=%2d\n",time2str(time,tstr,3),sat); - - dant[0]=dant[1]=dant[2]=0.0; - - /* sun position in ecef */ - sunmoonpos(gpst2utc(time),erpv,rsun,NULL,&gmst); - - /* unit vectors of satellite fixed coordinates */ - for (i=0;i<3;i++) r[i]=-rs[i]; - if (!normv3(r,ez)) return; - for (i=0;i<3;i++) r[i]=rsun[i]-rs[i]; - if (!normv3(r,es)) return; - cross3(ez,es,r); - if (!normv3(r,ey)) return; - cross3(ey,ez,ex); - - /* iono-free LC coefficients */ - sys=satsys(sat,NULL); - if (sys==SYS_GPS||sys==SYS_QZS) { /* L1-L2 */ - freq[0]=FREQL1; - freq[1]=FREQL2; - } - else if (sys==SYS_GLO) { /* G1-G2 */ - freq[0]=sat2freq(sat,CODE_L1C,nav); - freq[1]=sat2freq(sat,CODE_L2C,nav); - } - else if (sys==SYS_GAL) { /* E1-E5b */ - freq[0]=FREQL1; - freq[1]=FREQE5b; - } - else if (sys==SYS_CMP) { /* B1I-B2I */ - freq[0]=FREQ1_CMP; - freq[1]=FREQ2_CMP; - } - else if (sys==SYS_IRN) { /* L5-S */ - freq[0]=FREQL5; - freq[1]=FREQs; - } - else return; - - C1= SQR(freq[0])/(SQR(freq[0])-SQR(freq[1])); - C2=-SQR(freq[1])/(SQR(freq[0])-SQR(freq[1])); - - /* iono-free LC */ - for (i=0;i<3;i++) { - dant1=pcv->off[0][0]*ex[i]+pcv->off[0][1]*ey[i]+pcv->off[0][2]*ez[i]; - dant2=pcv->off[1][0]*ex[i]+pcv->off[1][1]*ey[i]+pcv->off[1][2]*ez[i]; - dant[i]=C1*dant1+C2*dant2; - } +/* Satellite antenna phase center offset --------------------------------------- + * Compute satellite antenna phase center offset in ECEF + * Args : gtime_t time I time (GPST) + * double *rs I satellite position and velocity (ECEF) + * {x,y,z,vx,vy,vz} (m|m/s) + * int sat I satellite number + * nav_t *nav I navigation data + * double *dant O satellite antenna phase center offset (ECEF) + * {dx,dy,dz} (m) (iono-free LC value) + * Return : none + * Notes : iono-free LC frequencies defined as follows: + * GPS/QZSS : L1-L2 + * GLONASS : G1-G2 + * Galileo : E1-E5b + * BDS : B1I-B2I + * NavIC : L5-S + *----------------------------------------------------------------------------*/ +extern void satantoff(gtime_t time, const double *rs, int sat, const nav_t *nav, double *dant) { + char tstr[40]; + trace(4, "satantoff: time=%s sat=%2d\n", time2str(time, tstr, 3), sat); + + dant[0] = dant[1] = dant[2] = 0.0; + + /* Sun position in ECEF */ + double rsun[3], gmst; + const double erpv[5] = {0}; + sunmoonpos(gpst2utc(time), erpv, rsun, NULL, &gmst); + + /* Unit vectors of satellite fixed coordinates */ + double r[3]; + for (int i = 0; i < 3; i++) r[i] = -rs[i]; + double ez[3]; + if (!normv3(r, ez)) return; + for (int i = 0; i < 3; i++) r[i] = rsun[i] - rs[i]; + double es[3]; + if (!normv3(r, es)) return; + cross3(ez, es, r); + double ey[3]; + if (!normv3(r, ey)) return; + double ex[3]; + cross3(ey, ez, ex); + + /* Iono-free LC coefficients */ + int sys = satsys(sat, NULL); + int code1, code2; + if (sys == SYS_GPS || sys == SYS_QZS) { /* L1-L2 */ + code1 = CODE_L1C; + code2 = CODE_L2W; + } else if (sys == SYS_GLO) { /* G1-G2 */ + code1 = CODE_L1C; + code2 = CODE_L2C; + } else if (sys == SYS_GAL) { /* E1-E5b */ + code1 = CODE_L1C; + code2 = CODE_L7I; + } else if (sys == SYS_CMP) { /* B1I-B2I */ + code1 = CODE_L2I; + code2 = CODE_L7I; + } else if (sys == SYS_IRN) { /* L5-S */ + code1 = CODE_L5A; + code2 = CODE_L9A; + } else + return; + + double freq1 = sat2freq(sat, code1, nav); + double freq2 = sat2freq(sat, code2, nav); + double C1 = SQR(freq1) / (SQR(freq1) - SQR(freq2)); + double C2 = -SQR(freq2) / (SQR(freq1) - SQR(freq2)); + + /* Iono-free LC */ + double pco1[3], pco2[3]; + const pcv_t *pcv = nav->pcvs + sat - 1; + antpco(pcv, freq1, pco1); + antpco(pcv, freq2, pco2); + + for (int i = 0; i < 3; i++) { + double dant1 = pco1[0] * ex[i] + pco1[1] * ey[i] + pco1[2] * ez[i]; + double dant2 = pco2[0] * ex[i] + pco2[1] * ey[i] + pco2[2] * ez[i]; + dant[i] = C1 * dant1 + C2 * dant2; + } } /* satellite position/clock by precise ephemeris/clock ------------------------- * compute satellite position/clock with precise ephemeris/clock diff --git a/src/rinex.c b/src/rinex.c index f43004813..09f65dcb8 100644 --- a/src/rinex.c +++ b/src/rinex.c @@ -157,15 +157,6 @@ typedef struct { /* signal index type */ double shift[MAXOBSTYPE]; /* phase shift (cycle) */ } sigind_t; -/* set string without tail space ---------------------------------------------*/ -static void setstr(char *dst, const char *src, int n) -{ - char *p=dst; - const char *q=src; - while (*q&&q=dst&&*p==' ') *p--='\0'; -} /* adjust time considering week handover -------------------------------------*/ static gtime_t adjweek(gtime_t t, gtime_t t0) { diff --git a/src/rtkcmn.c b/src/rtkcmn.c index f3de7c8c0..6b16d8230 100644 --- a/src/rtkcmn.c +++ b/src/rtkcmn.c @@ -1553,6 +1553,15 @@ extern void matprint(const double A[], int n, int m, int p, int q) { matfprint(A,n,m,p,q,stdout); } +/* set string without tail space ---------------------------------------------*/ +extern void setstr(char *dst, const char *src, int n) +{ + char *p=dst; + const char *q=src; + while (*q&&q=dst&&*p==' ') *p--='\0'; +} /* string to number ------------------------------------------------------------ * convert substring in string to number * args : char *s I string ("... nnn.nnn ...") @@ -2365,6 +2374,134 @@ static void nut_iau1980(double t, const double *f, double *dpsi, double *deps) *dpsi*=1E-4*AS2R; /* 0.1 mas -> rad */ *deps*=1E-4*AS2R; } + +static int code2sys(char code) { + switch (code) { + case 'G': return SYS_GPS; + case 'R': return SYS_GLO; + case 'E': return SYS_GAL; + case 'J': return SYS_QZS; + case 'C': return SYS_CMP; + case 'I': return SYS_IRN; + case 'L': return SYS_LEO; + case 'S': return SYS_SBS; + } + return 0; +} + +// Add satellite svn to prn mapping ------------------------------------------- +static void addsatsvn(const satsvn_t *satsvn, satsvns_t *satsvns) +{ + satsvn_t *data; + if (satsvns->nmax <= satsvns->n) { + satsvns->nmax += 256; + data = (satsvn_t *)realloc(satsvns->satsvn, sizeof(satsvn_t) * satsvns->nmax); + if (!data) { + trace(1,"addsatsvn: memory allocation error\n"); + free(satsvns->satsvn); + satsvns->satsvn = NULL; + satsvns->n = satsvns->nmax = 0; + return; + } + satsvns->satsvn = data; + } + satsvns->satsvn[satsvns->n++] = *satsvn; +} +/* Read SINEX meta data file --------------------------------------------------*/ +extern int readsinex(const char *file, satsvns_t *satsvns) { + trace(3, "readsinex: file=%s\n", file); + + FILE *fp = fopen(file, "r"); + if (!fp) { + trace(2, "sinex file open error: %s\n", file); + return 0; + } + int state = 0; + char buff[512]; + while (fgets(buff, sizeof(buff), fp)) { + if (buff[0] == '*') continue; // Comment + if (strncmp(buff, "+SATELLITE/PRN", 14) == 0) { + state = 1; + continue; + } + if (state == 1) { + if (strncmp(buff, "-SATELLITE/PRN", 14) == 0) { + state = 0; + continue; + } + if (strlen(buff) < 40) { + trace(0, "Sinex: unexpect line '%s'\n", buff); + continue; + } + satsvn_t satsvn = {0}; + int svnsys = code2sys(buff[1]); + satsvn.svn = str2num(buff, 2, 3); + char satid[8]=""; + snprintf(satid, 8, "%.3s", buff + 36); + satsvn.sat = satid2no(satid); + if (satsvn.sat == 0) { + trace(1, "Sinex: no mapping to sat for '%s'\n", satid); + continue; // Give up on this mapping + } + if (svnsys != satsys(satsvn.sat, NULL)) { + trace(0, "Error: sinex svn sys '%c' and prn sys '%c' not consistent.\n", buff[1], buff[36]); + continue; + } + + double start_year = str2num(buff, 6, 4); + double start_doy = str2num(buff, 11, 3); + double start_sod = str2num(buff, 15, 5); + if (start_year != 0 || start_doy != 0 || start_sod > 0.5) { + double ep[6] = {2000, 1, 1, 0, 0, 0}; + ep[0] = start_year; + satsvn.ts = timeadd(epoch2time(ep), (start_doy - 1.0) * 86400.0 + start_sod); + } + double end_year = str2num(buff, 21, 4); + double end_doy = str2num(buff, 26, 3); + double end_sod = str2num(buff, 30, 5); + if (end_year != 0 || end_doy != 0 || end_sod > 0.5) { + double ep[6] = {2000, 1, 1, 0, 0, 0}; + ep[0] = end_year; + satsvn.te = timeadd(epoch2time(ep), (end_doy - 1.0) * 86400.0 + end_sod); + } + addsatsvn(&satsvn, satsvns); + continue; + } + // Ignored + } + + return 1; +} + +// Search for a mapping from the SVN (system/num) to satellite number (system/PRN). +// Return: the satellite number (0: error) +extern int searchsvnsat(int sys, int svn, const gtime_t time, const satsvns_t *satsvns) { + trace(4, "searchsvnsat: sys=%2d svn=%s\n", sys, svn); + for (int i = 0; i < satsvns->n; i++) { + satsvn_t *satsvn = satsvns->satsvn + i; + if (satsvn->svn != svn) continue; + if (satsys(satsvn->sat, NULL) != sys) continue; + if (satsvn->ts.time != 0 && timediff(satsvn->ts, time) > 0.0) continue; + if (satsvn->te.time != 0 && timediff(satsvn->te, time) < 0.0) continue; + return satsvn->sat; + } + return 0; +} +// Search for a mapping from the satellite number (system/PRN) to and SVN number. +// Return: the satellite prn number (0: error) +extern int searchsatsvn(int sat, const gtime_t time, const satsvns_t *satsvns) { + trace(4, "searchsatsvn: sat=%d\n", sat); + for (int i = 0; i < satsvns->n; i++) { + satsvn_t *satsvn = satsvns->satsvn + i; + if (satsvn->sat != sat) continue; + if (satsvn->ts.time != 0 && timediff(satsvn->ts, time) > 0.0) continue; + if (satsvn->te.time != 0 && timediff(satsvn->te, time) < 0.0) continue; + return satsvn->svn; + } + return 0; +} + + /* eci to ecef transformation matrix ------------------------------------------- * compute eci to ecef transformation matrix * args : gtime_t tutc I time in utc @@ -2468,206 +2605,696 @@ static void addpcv(const pcv_t *pcv, pcvs_t *pcvs) } pcvs->pcv[pcvs->n++]=*pcv; } +// The pcv init bit encoding. +#define PCV_PCO 1 +#define PCV_NOAZI 2 +#define PCV_PHV 4 /* read ngs antenna parameter file -------------------------------------------*/ -static int readngspcv(const char *file, pcvs_t *pcvs) -{ - FILE *fp; - static const pcv_t pcv0={0}; - pcv_t pcv; - double neu[3]; - int n=0; - char buff[256]; - - if (!(fp=fopen(file,"r"))) { - trace(2,"ngs pcv file open error: %s\n",file); - return 0; +static int readngspcv(const char *file, pcvs_t *pcvs) { + FILE *fp; + static const pcv_t pcv0 = {0}; + pcv_t pcv; + double neu[3]; + int n = 0; + char buff[256]; + + if (!(fp = fopen(file, "r"))) { + trace(2, "ngs pcv file open error: %s\n", file); + return 0; + } + while (fgets(buff, sizeof(buff), fp)) { + if (strlen(buff) >= 62 && buff[61] == '|') continue; + + if (buff[0] != ' ') n = 0; /* start line */ + if (++n == 1) { + pcv = pcv0; + setstr(pcv.type, buff, 61); + pcv.zen1[0] = pcv.zen1[1] = 0; + pcv.zen2[0] = pcv.zen2[1] = 90; + pcv.dzen[0] = pcv.dzen[1] = 5; + pcv.zen_len[0] = pcv.zen_len[1] = 19; + pcv.azi_len[0] = pcv.azi_len[1] = 1; + pcv.var[0] = calloc(pcv.azi_len[0] * pcv.zen_len[0], sizeof(double)); + pcv.var[1] = calloc(pcv.azi_len[1] * pcv.zen_len[1], sizeof(double)); + } else if (n == 2) { + if (decodef(buff, 3, neu) < 3) continue; + pcv.off[0][0] = neu[1]; + pcv.off[0][1] = neu[0]; + pcv.off[0][2] = neu[2]; + } else if (n == 3) + decodef(buff, 10, pcv.var[0] + 0); + else if (n == 4) + decodef(buff, 9, pcv.var[0] + 10); + else if (n == 5) { + if (decodef(buff, 3, neu) < 3) continue; + pcv.off[1][0] = neu[1]; + pcv.off[1][1] = neu[0]; + pcv.off[1][2] = neu[2]; + } else if (n == 6) + decodef(buff, 10, pcv.var[1] + 0); + else if (n == 7) { + decodef(buff, 9, pcv.var[1] + 10); + pcv.init[0] = pcv.init[1] = PCV_PCO | PCV_NOAZI; + addpcv(&pcv, pcvs); } - while (fgets(buff,sizeof(buff),fp)) { + } + fclose(fp); + + return 1; +} + +// Antex frequencyindex ---------------------------------------------- +// +// A frequency index is arbitrarily assigned to each disinct frequency +// used in Antex definitions, as a common key between the observation +// codes and the Antex codes. +static int antfreqidx(int sys, int frq) { + switch (sys) { + case SYS_GPS: + switch (frq) { + case 1: + return ANTFREQL1; + case 2: + return ANTFREQL2; + case 5: + return ANTFREQL5; + } + break; + case SYS_GLO: + switch (frq) { + // For receivers, R01 and R02 do not differentiate by satellite. + // For a satellite this could well be the particular frequency. + case 1: + return ANTFREQG1; + case 2: + return ANTFREQG2; + case 3: + return ANTFREQG3; + // Note R04 entries often share with R01 (1600.995 vs 1602.0 Mhz) + case 4: + return ANTFREQG4; + // Note R06 entries often share with R02 (1248.06 vs 1246.0 MHz) + case 6: + return ANTFREQG6; + } + break; + case SYS_GAL: + switch (frq) { + case 1: + return ANTFREQL1; + case 5: + return ANTFREQL5; + case 6: + return ANTFREQL6; + case 7: + return ANTFREQL7; + case 8: + return ANTFREQL8; + } + break; + case SYS_QZS: + switch (frq) { + case 1: + return ANTFREQL1; + case 2: + return ANTFREQL2; + case 5: + return ANTFREQL5; + case 6: + return ANTFREQL6; + } + break; + case SYS_CMP: + switch (frq) { + case 1: + return ANTFREQL1; + case 2: + return ANTFREQB1; + case 5: + return ANTFREQL5; + case 6: + return ANTFREQB3; + case 7: + return ANTFREQL7; + case 8: + return ANTFREQL8; + } + break; + case SYS_IRN: + switch (frq) { + case 1: + return ANTFREQL1; + case 5: + return ANTFREQL5; + case 9: + return ANTFREQI9; + } + break; + case SYS_SBS: + switch (frq) { + case 1: + return ANTFREQL1; + case 5: + return ANTFREQL5; + } + break; + } + trace(1, "antfreqidx: unexpected sys=%d frq=%d\n", sys, frq); + return -1; +} +static double antfreq[ANTNFREQ] = {FREQL1, FREQL2, FREQL5, FREQL6, FREQE5b, + FREQE5ab, FREQ1_GLO, FREQ2_GLO, FREQ3_GLO, FREQ1a_GLO, + FREQ2a_GLO, FREQ1_CMP, FREQ3_CMP, FREQs}; - if (strlen(buff)>=62&&buff[61]=='|') continue; +extern double antcode2freq(int code) { + if (code < 0 || code >= ANTNFREQ) return 0.0; + return antfreq[code]; +} - if (buff[0]!=' ') n=0; /* start line */ - if (++n==1) { - pcv=pcv0; - strncpy(pcv.type,buff,61); pcv.type[61]='\0'; - } - else if (n==2) { - if (decodef(buff,3,neu)<3) continue; - pcv.off[0][0]=neu[1]; - pcv.off[0][1]=neu[0]; - pcv.off[0][2]=neu[2]; +static const char *antcodeid[ANTNFREQ] = {"L1", "L2", "L5", "L6", "L7", "L8", + "G1", "G2", "G3", "G4", "G6", + "B1", "B3", "I9"}; +extern const char *antcode2id(int code) { + if (code < 0 || code >= ANTNFREQ) return ""; + return antcodeid[code]; +} + + +/* Read antex file ----------------------------------------------------------*/ +static int readantex(const char *file, int filter, pcvs_t *pcvs) { + trace(3, "readantex: file=%s filter=%x\n", file, filter); + + FILE *fp = fopen(file, "r"); + if (!fp) { + trace(2, "antex pcv file open error: %s\n", file); + return 0; + } + static const pcv_t pcv0 = {0}; + pcv_t pcv; + int state = 0, azi_num = 0, group[ANTNFREQ] = {0}, fdup[ANTNFREQ] = {0}; + int procp = 0, phvp = 0; + double dazi = 0, zen1 = 0, zen2 = 0, dzen = 0; + int azi_len = 0, zen_len = 0; + char buff[512]; + while (fgets(buff, sizeof(buff), fp)) { + if (strlen(buff) < 60) continue; + if (strstr(buff + 60, "COMMENT")) continue; + if (strstr(buff + 60, "ANTEX VERSION / SYST")) continue; + if (strstr(buff + 60, "PCV TYPE / REFANT")) continue; + if (strstr(buff + 60, "ANTENNA TYPES")) { + // ATX2: SATELLITE, MIXED, RECEIVER + continue; + } + if (strstr(buff + 60, "REFERENCE FRAME")) { + // ATX2: e.g. IGS20 + continue; + } + if (strstr(buff + 60, "RELEASE")) { + // ATX2 + continue; + } + if (strstr(buff + 60, "END OF HEADER")) continue; + + if (strstr(buff + 60, "START OF ANTENNA")) { + pcv = pcv0; + state = 1; + procp = 0; + phvp = 0; + azi_len = 0; + continue; + } + if (strstr(buff + 60, "END OF ANTENNA")) { + if (state != 1) { + trace(1, "readantex: unexpected end of antenna\n"); + continue; + } + if (procp) addpcv(&pcv, pcvs); + state = 0; + continue; + } + if (state == 1 && strstr(buff + 60, "TYPE / SERIAL NO")) { + setstr(pcv.type, buff, 20); + setstr(pcv.code, buff + 20, 20); + if (strlen(pcv.code) == 3) { + pcv.sat = satid2no(pcv.code); + } + if (buff[40] != ' ') { + pcv.satsys = code2sys(buff[40]); + if (pcv.satsys == 0) { + trace(1, "readantex: unexpect SVN sys at line '%s'\n", buff); + continue; } - else if (n==3) decodef(buff,10,pcv.var[0]); - else if (n==4) decodef(buff,9,pcv.var[0]+10); - else if (n==5) { - if (decodef(buff,3,neu)<3) continue; - pcv.off[1][0]=neu[1]; - pcv.off[1][1]=neu[0]; - pcv.off[1][2]=neu[2]; + pcv.svn = str2num(buff, 41, 3); + } + procp = (pcv.sat || pcv.svn) ? !(filter & 2) : !(filter & 1); + continue; + } + if (state == 1 && strstr(buff + 60, "TYPE / SN")) { + // ATX2, receiver + setstr(pcv.type, buff, 20); + procp = !(filter & 1); + continue; + } + if (state == 1 && strstr(buff + 60, "TYPE / SVN / SAT ID")) { + // ATX2, satellite + setstr(pcv.type, buff, 20); + pcv.satsys = code2sys(buff[40]); + if (pcv.satsys == 0) { + trace(1, "readantex: unexpect SVN sys at line '%s'\n", buff); + continue; + } + pcv.svn = str2num(buff, 41, 3); + // The mapping to a RTKLIB satellite number needs a time, and with ATX2 + // uses the satellite meta data. + procp = !(filter & 2); + continue; + } + // + // Ignore content if it has been filtered out, if it is not being + // processed, until the end of this antenna entry. + if (!procp) continue; + // + if (state == 1 && strstr(buff + 60, "METH / BY / # / DATE")) { + continue; + } + if (state == 1 && strstr(buff, "DAZI")) { + if (sscanf(buff, "%lf", &dazi) < 1) continue; + // azi_len includes the NOAZI entry. + azi_len = (dazi < 0.01 ? 0 : 360 / dazi + 1) + 1; + continue; + } + if (state == 1 && strstr(buff, "ZEN1 / ZEN2 / DZEN")) { + if (sscanf(buff, "%lf%lf%lf", &zen1, &zen2, &dzen) < 3) continue; + zen_len = dzen < 0.01 ? 0 : (zen2 - zen1) / dzen + 1; + continue; + } + if (state == 1 && strstr(buff, "# OF FREQUENCIES")) { + continue; + } + if (state == 1 && strstr(buff, "# OF PHV")) { + // ATX2 + continue; + } + if (state == 1 && strstr(buff + 60, "VALID FROM")) { + if (phvp) { + // Save the current entry and start a new entry. + addpcv(&pcv, pcvs); + // Clear time period and the frequency data. + memset(&pcv.ts, 0, sizeof(pcv.ts)); + memset(&pcv.te, 0, sizeof(pcv.te)); + for (int idx = 0; idx < ANTNFREQ; idx++) { + pcv.zen1[idx] = pcv.zen2[idx] = pcv.dzen[idx] = pcv.dazi[idx] = 0.0; + pcv.zen_len[idx] = pcv.azi_len[idx] = pcv.init[idx] = 0; + for (int j = 0; j < 3; j ++) pcv.off[idx][j] = 0.0; + pcv.var[idx] = NULL; } - else if (n==6) decodef(buff,10,pcv.var[1]); - else if (n==7) { - decodef(buff,9,pcv.var[1]+10); - addpcv(&pcv,pcvs); + phvp = 0; + } + if (str2time(buff, 0, 43, &pcv.ts) < 0) + trace(1, "readantex: error parsing pcv 'VALID FROM' time at line '%s'\n", buff); + continue; + } + if (state == 1 && strstr(buff + 60, "VALID UNTIL")) { + if (phvp) { + // Save the current entry and start a new entry. + addpcv(&pcv, pcvs); + // Clear time period and the frequency data. + memset(&pcv.ts, 0, sizeof(pcv.ts)); + memset(&pcv.te, 0, sizeof(pcv.te)); + for (int idx = 0; idx < ANTNFREQ; idx++) { + pcv.zen1[idx] = pcv.zen2[idx] = pcv.dzen[idx] = pcv.dazi[idx] = 0.0; + pcv.zen_len[idx] = pcv.azi_len[idx] = pcv.init[idx] = 0; + for (int j = 0; j < 3; j ++) pcv.off[idx][j] = 0.0; + pcv.var[idx] = NULL; } + phvp = 0; + } + if (str2time(buff, 0, 43, &pcv.te) < 0) + trace(1, "readantex: error parsing pcv 'VALID UNTIL' time at line '%s'\n", buff); + continue; } - fclose(fp); - - return 1; -} -/* read antex file ----------------------------------------------------------*/ -static int readantex(const char *file, pcvs_t *pcvs) -{ - FILE *fp; - static const pcv_t pcv0={0}; - pcv_t pcv; - double neu[3]; - int i,f,freq=0,state=0,freqs[]={1,2,5,0}; - char buff[256]; - - trace(3,"readantex: file=%s\n",file); - - if (!(fp=fopen(file,"r"))) { - trace(2,"antex pcv file open error: %s\n",file); - return 0; + if (state == 1 && strstr(buff, "SINEX CODE")) { + continue; } - while (fgets(buff,sizeof(buff),fp)) { - - if (strlen(buff)<60||strstr(buff+60,"COMMENT")) continue; - - if (strstr(buff+60,"START OF ANTENNA")) { - pcv=pcv0; - state=1; + if (state == 1 && (strstr(buff + 60, "START OF FREQUENCY") || + strstr(buff + 60, "START OF PHV"))) { + phvp = 1; + if (strstr(buff + 60, "START OF PHV")) { + // ATX2 - "START OF PHV" can include a group of signals. The systems + // are not distinguished here but an entry is allocated for each + // distinct frequency, e.g. R02 and R06 are often grouped but have + // slightly differency frequencies. + state = 3; + } else { + state = 2; + } + azi_num = 0; + // Note all the distrinct frequencies in the group. + for (int idx = 0; idx < ANTNFREQ; idx++) { + group[idx] = 0; + fdup[idx] = 0; + } + for (int i = 0; i < 10; i++) { + int start = 6 * i + 3; + // TODO should it break here, is this the end of the list? + if (buff[start] == ' ') continue; + int sys = code2sys(buff[start]); + if (sys == 0) { + trace(1, "readantex: unexpected code system='%c'\n", buff[start]); + continue; } - if (strstr(buff+60,"END OF ANTENNA")) { - addpcv(&pcv,pcvs); - state=0; + int frq = -1; + if (sscanf(buff + start + 1, "%d", &frq) < 1) { + trace(1, "readantex: unexpected code at '%s'\n", buff + start); + continue; } - if (!state) continue; - - if (strstr(buff+60,"TYPE / SERIAL NO")) { - strncpy(pcv.type,buff ,20); pcv.type[20]='\0'; - strncpy(pcv.code,buff+20,20); pcv.code[20]='\0'; - if (!strncmp(pcv.code+3," ",8)) { - pcv.sat=satid2no(pcv.code); - } + if (frq < 0 || frq > 9) { + trace(1, "readantex: unexpected code frq=%d\n", frq); + continue; } - else if (strstr(buff+60,"VALID FROM")) { - if (!str2time(buff,0,43,&pcv.ts)) continue; + if (pcv.sat && satsys(pcv.sat, NULL) != sys) { + trace(1, "readantex: unexpected sat sys %d vs %d\n", satsys(pcv.sat, NULL), sys); + continue; } - else if (strstr(buff+60,"VALID UNTIL")) { - if (!str2time(buff,0,43,&pcv.te)) continue; + int idx = antfreqidx(sys, frq); + if (idx < 0) { + trace(1, "readantex: unknown freq idx for sys=%d frq=%d\n", sys, frq); + continue; } - else if (strstr(buff+60,"START OF FREQUENCY")) { - if (!pcv.sat&&buff[3]!='G') continue; /* only read rec ant for GPS */ - if (sscanf(buff+4,"%d",&f)<1) continue; - for (i=0;freqs[i];i++) if (freqs[i]==f) break; - if (freqs[i]) freq=i+1; - /* for Galileo E5b: save to E2, not E7 */ - if (satsys(pcv.sat,NULL)==SYS_GAL&&f==7) freq=2; + group[idx] = 1; + } + // Allocate a var array for each distinct frequency index in the group. + for (int idx = 0; idx < ANTNFREQ; idx++) { + if (group[idx]) { + if (pcv.var[idx] != NULL) { + if (zen_len > 0) { + fdup[idx] = 1; + char code[4]; + memcpy(code, buff + 3, 3); + code[3] = '\0'; + trace(4, "readantex: duplicate entry type='%s' code=%s idx=%d\n", pcv.type, code, idx); + } + } else { + fdup[idx] = 0; + pcv.dazi[idx] = dazi; + pcv.azi_len[idx] = azi_len; + pcv.zen1[idx] = zen1; + pcv.zen2[idx] = zen2; + pcv.dzen[idx] = dzen; + pcv.zen_len[idx] = zen_len; + if (zen_len > 0) + pcv.var[idx] = calloc(azi_len * zen_len, sizeof(double)); + } + trace(4, "readantex: alloc pcv idx=%d azi_len=%d zen_len=%d %p\n", idx, azi_len, zen_len, (void *)pcv.var[idx]); } - else if (strstr(buff+60,"END OF FREQUENCY")) { - freq=0; + } + continue; + } + if ((state == 2 && strstr(buff + 60, "END OF FREQUENCY")) || (state == 3 && strstr(buff + 60, "END OF PHV"))) { + for (int idx = 0; idx < ANTNFREQ; idx++) { + if (group[idx]) { + if (pcv.azi_len[idx] > 1) { + if (azi_num == pcv.azi_len[idx] - 1) { + pcv.init[idx] |= PCV_PHV; + } else { + trace(1, "readantex: %d of %d azimuth elements supplied\n", azi_num, pcv.azi_len[idx]); + } + } + group[idx] = 0; } - else if (strstr(buff+60,"NORTH / EAST / UP")) { - if (freq<1||NFREQ= azi_len - 1) { + trace(1, "readantex: PHV AZI overflow %d expected %d\n", azi_num, azi_len); + continue; + } + double *offs = malloc((zen_len + 1) * sizeof(double)); + if (!offs) { + trace(1, "readantex: PHV AZI offs alloc failed\n"); + continue; + } + int i = decodef(buff + 1, zen_len + 1, offs); + if (i < zen_len) { + trace(1, "readantex: PHV AZI only %d of %d values read\n", i, zen_len); + free(offs); + continue; + } + // Check the azi, the first number. + if (fabs(offs[0] * 1000.0 - azi_num * dazi) > 0.01) { + trace(1, "readantex: PHV AZI error %lf expected %lf azi_num=%d dazi=%lf\n", offs[0], azi_num * dazi, azi_num, dazi); + free(offs); + continue; + } + for (int idx = 0; idx < ANTNFREQ; idx++) { + if (group[idx]) { + if (fdup[idx]) { + // Guard the dup size being the same. + if (azi_num >= pcv.azi_len[idx] - 1) { + trace(1, "readantex PHV AZI dup overflows azi length %d\n", pcv.azi_len[idx]); + continue; + } + if (zen_len != pcv.zen_len[idx]) { + trace(1, "readantex PHV AZI dup zen length mismatch %d expected %d\n", zen_len, pcv.zen_len[idx]); + continue; + } + for (int j = 0; j < zen_len; j++) { + if (offs[1 + j] != pcv.var[idx][(1 + azi_num) * zen_len + j]) + trace(1, "readantex PHV AZI dup data differs j=%d %lf %lf\n", j, offs[j], pcv.var[idx][(1 + azi_num) * zen_len + j]); + } + continue; + } + memcpy(&pcv.var[idx][(1 + azi_num) * zen_len], offs + 1, zen_len * sizeof(double)); + } + } + azi_num++; + free(offs); + continue; } - for (i=0;in;i++) { - pcv=pcvs->pcv+i; - trace(4,"sat=%2d type=%20s code=%s off=%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n", - pcv->sat,pcv->type,pcv->code,pcv->off[0][0],pcv->off[0][1], - pcv->off[0][2],pcv->off[1][0],pcv->off[1][1],pcv->off[1][2]); + if (state == 1 && strstr(buff + 60, "START OF GDV")) { + state = 4; + continue; } - return stat; + if (state == 4 && strstr(buff + 60, "END OF GDV")) { + state = 1; + continue; + } + if (state == 4) { + // GDV Ignored. + continue; + } + if (state == 1 && strstr(buff + 60, "START OF GAIN")) { + state = 5; + continue; + } + if (state == 5 && strstr(buff + 60, "END OF GAIN")) { + state = 1; + continue; + } + if (state == 5) { + // GAIN Ignored. + continue; + } + + trace(1, "readantex: unexpected state=%d at line %s\n", state, buff); + } + fclose(fp); + return 1; +} +/* read antenna parameters ------------------------------------------------------ + * read antenna parameters + * args : char *file I antenna parameter file (antex) + * int filter I bit 0 - exclude rcv, bit 1 - exclude sat + * pcvs_t *pcvs IO antenna parameters + * return : status (1:ok,0:file open error) + * notes : file with the extension .atx or .ATX is recognized as antex + * file except for antex is recognized ngs antenna parameters + * see reference [3] + * only support non-azimuth-depedent parameters + *-----------------------------------------------------------------------------*/ +extern int readpcv(const char *file, int filter, pcvs_t *pcvs) { + trace(3, "readpcv: file=%s\n", file); + + char *ext = strrchr(file, '.'); + if (!ext) ext = ""; + + int stat = 0; + if (!strcmp(ext, ".atx") || !strcmp(ext, ".ATX") || + !strcmp(ext, ".atx2") || !strcmp(ext, ".ATX2")) { + stat = readantex(file, filter, pcvs); + } else if (!(filter & 1)) { + stat = readngspcv(file, pcvs); + } + for (int i = 0; i < pcvs->n; i++) { + pcv_t *pcv = pcvs->pcv + i; + trace(4, "sat=%2d type=%20s code=%s off=%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n", pcv->sat, + pcv->type, pcv->code, pcv->off[0][0], pcv->off[0][1], pcv->off[0][2], pcv->off[1][0], + pcv->off[1][1], pcv->off[1][2]); + } + return stat; } /* search antenna parameter ---------------------------------------------------- -* read satellite antenna phase center position -* args : int sat I satellite number (0: receiver antenna) -* char *type I antenna type for receiver antenna -* gtime_t time I time to search parameters -* pcvs_t *pcvs IO antenna parameters -* return : antenna parameter (NULL: no antenna) -*-----------------------------------------------------------------------------*/ -extern pcv_t *searchpcv(int sat, const char *type, gtime_t time, - const pcvs_t *pcvs) -{ - pcv_t *pcv; - char buff[MAXANT],*types[2],*p; - int i,j,n=0; - - trace(4,"searchpcv: sat=%2d type=%s\n",sat,type); - - if (sat) { /* search satellite antenna */ - for (i=0;in;i++) { - pcv=pcvs->pcv+i; - if (pcv->sat!=sat) continue; - if (pcv->ts.time!=0&&timediff(pcv->ts,time)>0.0) continue; - if (pcv->te.time!=0&&timediff(pcv->te,time)<0.0) continue; - return pcv; + * read satellite antenna phase center position + * args : int sat I satellite number (0: receiver antenna) + * char *type I antenna type for receiver antenna + * gtime_t time I time to search parameters + * satsvns_t satsvns I satellite SVN to PRN mapping + * pcvs_t *pcvs IO antenna parameters + * return : antenna parameter (NULL: no antenna) + *-----------------------------------------------------------------------------*/ +extern pcv_t *searchpcv(int sat, const char *type, gtime_t time, const satsvns_t *satsvns, const pcvs_t *pcvs) { + trace(4, "searchpcv: sat=%2d type=%s\n", sat, type); + + if (sat) { + // Search satellite antenna + int prn, sys = satsys(sat, &prn); + // Map to an SVN if there is satellite meta data. + int svn = 0; + if (satsvns) svn = searchsatsvn(sat, time, satsvns); + for (int i = 0; i < pcvs->n; i++) { + pcv_t *pcv = pcvs->pcv + i; + if (pcv->satsys && pcv->satsys != sys) continue; + if (svn) { + if (pcv->svn != svn) { + if (pcv->sat && pcv->sat == sat) { + trace(1, "searchpcv: matching sat=%d but mismatching svn=%d %d\n", sat, svn, pcv->svn); + // TODO ?? + } + continue; } - } - else { - strcpy(buff,type); - char *q; - for (p=strtok_r(buff," ",&q);p&&n<2;p=strtok_r(NULL," ",&q)) types[n++]=p; - if (n<=0) return NULL; - - /* search receiver antenna with radome at first */ - for (i=0;in;i++) { - pcv=pcvs->pcv+i; - for (j=0;jtype,types[j])) break; - if (j>=n) return pcv; + if (pcv->sat && pcv->sat != sat) { + trace(1, "searchpcv: matching svn=%d but mismatching sat=%d %d\n", svn, sat, pcv->sat); + // TODO ?? + continue; } - /* search receiver antenna without radome */ - for (i=0;in;i++) { - pcv=pcvs->pcv+i; - if (strstr(pcv->type,types[0])!=pcv->type) continue; + } else { + // No SVN info, just check the sat no + if (pcv->sat != sat) continue; + } + if (pcv->ts.time != 0 && timediff(pcv->ts, time) > 0.0) continue; + if (pcv->te.time != 0 && timediff(pcv->te, time) < 0.0) continue; + trace(3, "searchpcv: sat=%2d sys=%2d prn=%3d svn=%3d type=%s\n", sat, sys, prn, pcv->svn, pcv->type); + return pcv; + } + trace(2, "searchpcv: failed for sat=%d sys=%2d prn=%3d\n", sat, sys, prn); + } else { + char buff[MAXANT]; + strcpy(buff, type); + char *types[2]; + int n = 0; + for (char *q, *p = strtok_r(buff, " ", &q); p && n < 2; p = strtok_r(NULL, " ", &q)) + types[n++] = p; + if (n <= 0) return NULL; + + // Search receiver antenna with radome at first + for (int i = 0; i < pcvs->n; i++) { + pcv_t *pcv = pcvs->pcv + i; + int j; + for (j = 0; j < n; j++) + if (!strstr(pcv->type, types[j])) break; + if (j >= n) return pcv; + } + // Search receiver antenna without radome + for (int i = 0; i < pcvs->n; i++) { + pcv_t *pcv = pcvs->pcv + i; + if (strstr(pcv->type, types[0]) != pcv->type) continue; - trace(2,"pcv without radome is used type=%s\n",type); - return pcv; - } + trace(2, "pcv without radome is used type=%s\n", type); + return pcv; } - return NULL; + trace(2, "searchpcv: failed for rcv type='%s'\n", type); + } + return NULL; } /* read station positions ------------------------------------------------------ * read positions from station position file @@ -3669,35 +4296,40 @@ extern int seliflc(int optnf,int sys) /* use L1/L5 for Galileo if L5 is enabled */ return((optnf==2||sys!=SYS_GAL)?1:2); } -/* troposphere model ----------------------------------------------------------- -* compute tropospheric delay by standard atmosphere and saastamoinen model -* args : gtime_t time I time + +/* Troposphere model ----------------------------------------------------------- +* Compute tropospheric delay by standard atmosphere and saastamoinen model +* Args : gtime_t time I time * double *pos I receiver position {lat,lon,h} (rad,m) * double *azel I azimuth/elevation angle {az,el} (rad) * double humi I relative humidity -* return : tropospheric delay (m) -*-----------------------------------------------------------------------------*/ -extern double tropmodel(gtime_t time, const double *pos, const double *azel, - double humi) -{ - const double temp0=15.0; /* temperature at sea level */ - double hgt,pres,temp,e,z,trph,trpw; +* Return : tropospheric delay (m) +*----------------------------------------------------------------------------*/ +extern double tropmodel(gtime_t time, const double *pos, const double *azel, double humi) { + const double temp0 = 15.0; // Temperature at sea level + + if (pos[2] < -100.0 || 1E4 < pos[2] || azel[1] <= 0) return 0.0; - if (pos[2]<-100.0||1E4=18) return var[18]; - return var[i]*(1.0-a+i)+var[i+1]*(a-i); -} -/* receiver antenna model ------------------------------------------------------ -* compute antenna offset by antenna phase center parameters -* args : pcv_t *pcv I antenna phase center parameters -* double *del I antenna delta {e,n,u} (m) -* double *azel I azimuth/elevation for receiver {az,el} (rad) -* int opt I option (0:only offset,1:offset+pcv) -* double *dant O range offsets for each frequency (m) -* return : none -* notes : current version does not support azimuth dependent terms -*-----------------------------------------------------------------------------*/ -extern void antmodel(const pcv_t *pcv, const double *del, const double *azel, - int opt, double *dant) -{ - double e[3],off[3],cosel=cos(azel[1]); - int i,j; +/* Interpolate antenna phase center variation --------------------------------*/ +static double interpvar(double ang, double start, double end, double delta, const double *var) { + double a = (ang - start) / delta; + int i = trunc(a); + int n = (end - start) / delta + 1; + if (i < 0) { + trace(3, "interpvar: %d < 0\n", i); + return var[0]; + } + else if (i >= n) { + trace(3, "interpvar: %d >= %d\n", i, n); + return var[n - 1]; + } + return var[i] * (1.0 - (a - i)) + var[i + 1] * (a - i); +} + +// Find the closest available PCO entry for a given frequency. +// +// Returns the closest available Antex frequency index for the supplied +// frequency. Also the second closest in the same band and on the other side +// of the supplied frequency, for interpolation, unless the the first is +// within tolerance and consider a match. +// +// Different bands can have distinct offsets that appear unrelated to the +// change of the PCO with frequency, such as from stacked antenna elements, so +// interpolation is not attempted between bands. +// +// The relationship between the PCO and frequency is typically not be linear +// even within a band and measurements can have some noise so only interpolate +// if entries are found either side of the respective frequency. While for +// some antennas, with a few calibrated frequencies, the relationship between +// the offset and frequency might suggest fitting a curve for interpolation, +// there are other antenna calibrations that do not have smooth data. +// +static int antpcoidx(const pcv_t *pcv, double freq, double *freq1, double *freq2, int *idx2) { + if (freq2) *freq2 = 0.0; + if (idx2) *idx2 = -1; + + // Closest higher and lower available entries, in the same band. + int idxh = -1, idxl = -1; + double frqh = 0.0, frql = 0.0; + // Closest available entries, in the same band. + int idxa1 = -1, idxa2 = -1; + double frqa1 = 0.0, frqa2 = 0.0; + + for (int i = 0; i < ANTNFREQ; i++) { + // Only consider if there is a PCO entry. + if (pcv->init[i] & PCV_PCO) { + double frq = antcode2freq(i); + // Alternatives in the same band, within 140MHz. + if (fabs(frq - freq) > 145e6) continue; + + // The closest two entries in the same band. + if (idxa1 < 0) { + idxa1 = i; + frqa1 = frq; + } else if (fabs(frq - freq) < fabs(frqa1 - freq)) { + idxa2 = idxa1; + frqa2 = frqa1; + idxa1 = i; + frqa1 = frq; + } else if (idxa2 < 0 || fabs(frq - freq) < fabs(frqa2 - freq)) { + idxa2 = i; + frqa2 = frq; + } + + // Closest higher freq entry + if (frq >= freq && (idxh < 0 || frq < frqh)) { + idxh = i; + frqh = frq; + } + // Closest lower freq entry + if (frq <= freq && (idxl < 0 || frq > frql)) { + idxl = i; + frql = frq; + } + } + } + + if (idxh < 0 && idxl < 0) { + // No entries in the same band, return what was found and issue a warning. + trace(3, "antpcoidx: '%s' close PCO not found for freq=%.0lf\n", pcv->type, freq); + if (freq1) *freq1 = 0.0; + return -1; + } + + if (idxl < 0) { + // Only a higher frequency entry. + if (freq1) *freq1 = frqh; + return idxh; + } - trace(4,"antmodel: azel=%6.1f %4.1f opt=%d\n",azel[0]*R2D,azel[1]*R2D,opt); + if (idxh < 0) { + // Only a lower frequency entry. + if (freq1) *freq1 = frql; + return idxl; + } + + if (idxh == idxl) { + // Only one if there were the same frequency. + if (freq1) *freq1 = frqh; + return idxh; + } + + // Two close alternatives on either side, order them. + if (fabs(frqh - freq) <= fabs(frql - freq)) { + // If the closer higher frequency is within tolerance then consider it a + // match, avoiding interpolation. + if (fabs(frqh - freq) > 0.250E6) { + if (freq2) *freq2 = frql; + if (idx2) *idx2 = idxl; + } + if (freq1) *freq1 = frqh; + return idxh; + } - e[0]=sin(azel[0])*cosel; - e[1]=cos(azel[0])*cosel; - e[2]=sin(azel[1]); + // If the closer lower frequency is within tolerance then consider it a + // match, avoiding interpolation. + if (fabs(frql - freq) > 0.250E6) { + if (freq2) *freq2 = frqh; + if (idx2) *idx2 = idxh; + } + if (freq1) *freq1 = frql; + return idxl; +} + +/* Receiver antenna phase center variation -------------------------------------- + * Compute antenna offset by antenna phase center parameters + * Args : pcv_t *pcv I antenna phase center parameters + * double *azel I azimuth/elevation for receiver {az,el} (rad) + * double freq I frequency + * Return : range offsets for the respective frequency (m) + * + * Note: this applies the PCV when available, and for the PCO only see antpco(). + *-----------------------------------------------------------------------------*/ +extern double antpcv(const pcv_t *pcv, const double *azel, double freq) { + trace(4, "antpcv: azel=%6.1f %4.1f\n", azel[0] * R2D, azel[1] * R2D); + + double freq1, freq2; + int idx2, idx = antpcoidx(pcv, freq, &freq1, &freq2, &idx2); + if (idx < 0 || !(pcv->init[idx] & (PCV_PHV | PCV_NOAZI)) || pcv->dzen[idx] < 0.01) { + trace(3, "antpcv: no pcv freq=%.0lf\n", freq); + return 0.0; + } - for (i=0;ioff[i][j]+del[j]; + if (pcv->zen_len[idx] != (pcv->zen2[idx] - pcv->zen1[idx]) / pcv->dzen[idx] + 1) + trace(2, "antpcv: unexpected zen_len\n"); + + double dant = 0.0; + + if ((pcv->init[idx] & PCV_PHV) && pcv->dazi[idx] > 0.01) { + double pcv1 = 0.0; + { + double a = azel[0] * R2D / pcv->dazi[idx]; + int i = trunc(a); + double r = a - i; + if (i + 1 >= pcv->azi_len[idx] - 1) + trace(2, "antpcv azi %d >= len %d\n", i + 1, pcv->azi_len[idx] - 1); + if (i < 0) { + i = 0; + r = 0; + trace(2, "antpcv azi < 0 clipped\n"); + } else if (i >= pcv->azi_len[idx] - 1) { + i = pcv->azi_len[idx] - 2; + r = 0; + trace(2, "antpcv azi >= len clipped\n"); + } + // Interpolate for the zenith at each azimuth step. + double dant1 = interpvar(90.0 - azel[1] * R2D, pcv->zen1[idx], pcv->zen2[idx], pcv->dzen[idx], + pcv->var[idx] + (1 + i) * pcv->zen_len[idx]); + pcv1 += dant1; + if (i + 1 < pcv->azi_len[idx] - 1) { + // Interpolate for the zenith at each azimuth step. + double dant2 = interpvar(90.0 - azel[1] * R2D, pcv->zen1[idx], pcv->zen2[idx], pcv->dzen[idx], + pcv->var[idx] + (1 + i + 1) * pcv->zen_len[idx]); + // Interpolate for the azimuth. + pcv1 += -r * dant1 + r * dant2; + } + } + dant += pcv1; + + if (idx2 >= 0 && freq2 > 0 && pcv->init[idx2] & (PCV_PHV | PCV_NOAZI) && pcv->dzen[idx2] > 0.01 && pcv->dazi[idx2] > 0.01) { + // PCV with interpolation between frequencies. + if (pcv->zen_len[idx2] != (pcv->zen2[idx2] - pcv->zen1[idx2]) / pcv->dzen[idx2] + 1) + trace(2, "antpcv: unexpected zen_len\n"); + double pcv2 = 0.0; + double a = azel[0] * R2D / pcv->dazi[idx2]; + int i = trunc(a); + double r = a - i; + if (i + 1 >= pcv->azi_len[idx2] - 1) + trace(2, "antpcv azi %d >= len %d\n", i + 1, pcv->azi_len[idx2] - 1); + if (i < 0) { + i = 0; + r = 0; + trace(2, "antpcv azi < 0 clipped\n"); + } else if (i >= pcv->azi_len[idx2] - 1) { + i = pcv->azi_len[idx2] - 2; + r = 0; + trace(2, "antpcv azi >= len clipped\n"); + } + // Interpolate for the zenith at each azimuth step. + double dant1 = interpvar(90.0 - azel[1] * R2D, pcv->zen1[idx2], pcv->zen2[idx2], pcv->dzen[idx2], + pcv->var[idx2] + (1 + i) * pcv->zen_len[idx2]); + pcv2 += dant1; + if (i + 1 < pcv->azi_len[idx2] - 1) { + // Interpolate for the zenith at each azimuth step. + double dant2 = interpvar(90.0 - azel[1] * R2D, pcv->zen1[idx2], pcv->zen2[idx2], pcv->dzen[idx2], + pcv->var[idx2] + (1 + i + 1) * pcv->zen_len[idx2]); + // Interpolate for the azimuth. + pcv2 += -r * dant1 + r * dant2; + } + + // Interpolate between frequencies. + double lam = CLIGHT / freq, lam1 = CLIGHT / freq1, lam2 = CLIGHT / freq2, dlam = lam2 - lam1; + double b = (lam - lam1) / dlam; + dant += -b * pcv1 + b * pcv2; + } + trace(4, "antpcv: dant=%6.3f\n", dant); + return dant; + } - dant[i]=-dot3(off,e)+(opt?interpvar(90.0-azel[1]*R2D,pcv->var[i]):0.0); + // NOAZI + if (pcv->init[idx] & PCV_PHV) { + double pcv1 = interpvar(90.0 - azel[1] * R2D, pcv->zen1[idx], pcv->zen2[idx], pcv->dzen[idx], pcv->var[idx]); + dant += pcv1; + if (idx2 >= 0 && freq2 > 0 && pcv->init[idx2] & (PCV_PHV | PCV_NOAZI) && pcv->dzen[idx2] > 0.01) { + // PCV with interpolation between frequencies. + double pcv2 = interpvar(90.0 - azel[1] * R2D, pcv->zen1[idx2], pcv->zen2[idx2], pcv->dzen[idx2], pcv->var[idx2]); + double lam = CLIGHT / freq, lam1 = CLIGHT / freq1, lam2 = CLIGHT / freq2, dlam = lam2 - lam1; + double a = (lam - lam1) / dlam; + dant += -a * pcv1 + a * pcv2; } - trace(4,"antmodel: dant=%6.3f %6.3f\n",dant[0],dant[1]); + } + + trace(4, "antpcv: dant=%6.4f\n", dant); + return dant; +} + +/* Receiver antenna model ------------------------------------------------------ + * Compute antenna offset by antenna phase center parameters + * Args : pcv_t *pcv I antenna phase center parameters + * double *del I antenna delta {e,n,u} (m) + * double *azel I azimuth/elevation for receiver {az,el} (rad) + * int opt I option (0:only offset,1:offset+pcv) + * double freq I frequency + * Return : range offset for the respective frequency (m) + * + * Note: this applies the PCO and PCV when available. For the PCO only see + * antpco(). The receiver PCO does not depend on the azimuth and elevation to + * the satellite (less the PCV) and might affect these angles so the caller + * may wish to apply the PCO first and then call antpcv(), particularly if the + * antenna delta can be large. + *-----------------------------------------------------------------------------*/ +extern double antmodel(const pcv_t *pcv, const double *del, const double *azel, int opt, double freq) { + trace(4, "antmodel: azel=%6.1f %4.1f opt=%d\n", azel[0] * R2D, azel[1] * R2D, opt); + + double pco[3]; + antpco(pcv, freq, pco); + + double off[3]; + for (int i = 0; i < 3; i++) off[i] = pco[i] + del[i]; + + double e[3], cosel = cos(azel[1]); + e[0] = sin(azel[0]) * cosel; + e[1] = cos(azel[0]) * cosel; + e[2] = sin(azel[1]); + + double dant = -dot3(off, e); + if (opt) dant += antpcv(pcv, azel, freq); + + trace(4, "antmodel: dant=%6.4f\n", dant); + return dant; } -/* satellite antenna model ------------------------------------------------------ -* compute satellite antenna phase center parameters -* args : pcv_t *pcv I antenna phase center parameters -* double nadir I nadir angle for satellite (rad) -* double *dant O range offsets for each frequency (m) -* return : none -*-----------------------------------------------------------------------------*/ -extern void antmodel_s(const pcv_t *pcv, double nadir, double *dant) -{ - int i; - trace(4,"antmodel_s: nadir=%6.1f\n",nadir*R2D); +/* Satellite antenna model ------------------------------------------------------ + * Compute satellite antenna phase center parameters + * Args : pcv_t *pcv I antenna phase center parameters + * double nadir I nadir angle for satellite (rad) + * double freq I frequency + * Return : range offsets for respective frequency (m) + * Note this applies the PCV only, for the PCO see antpco(). + * Most satellites have only NOAZI entries, and only these are supported here. + *-----------------------------------------------------------------------------*/ +extern double antmodel_s(const pcv_t *pcv, double nadir, double freq) { + trace(4, "antmodel_s: sys=%d svn=%d nadir=%6.1f\n", pcv->satsys, pcv->svn, nadir * R2D); + + if (pcv->satsys == 0) { + trace(2, "antmodel_s: no pcv\n"); + return 0.0; + } + + double freq1, freq2; + int idx2, idx = antpcoidx(pcv, freq, &freq1, &freq2, &idx2); + if (idx < 0 || !(pcv->init[idx] & PCV_NOAZI)) { + trace(2, "antmodel_s: NOAZI PHV not found for sys=%d svn=%d freq=%.0lf idx=%d\n", pcv->satsys, pcv->svn, freq, idx); + return 0.0; + } + + if (pcv->dzen[idx] < 0.01) { + trace(2, "antmodel_s: NOAZI PHV has no zenith entries for sys=%d svn=%d freq=%.0lf idx=%d\n", pcv->satsys, pcv->svn, freq, idx); + return 0.0; + } + + double dant = interpvar(nadir * R2D, pcv->zen1[idx], pcv->zen2[idx], pcv->dzen[idx], pcv->var[idx]); + trace(4, "antmodel_s: dant=%6.4f\n", dant); + return dant; +} +// Returns the antenna PCO +// For a receiver this is enu, for a satellite xyz. +extern void antpco(const pcv_t *pcv, double freq, double pco[3]) { + double freq1, freq2; + int idx2, idx = antpcoidx(pcv, freq, &freq1, &freq2, &idx2); + if (idx < 0) { + pco[0] = pco[1] = pco[2] = 0.0; + trace(4, "antpco PCO not found for sys=%d svn=%d freq=%.0lf idx=%d\n", pcv->satsys, pcv->svn, freq, idx); + return; + } + + if (!(pcv->init[idx] & PCV_PCO)) + trace(2, "antpco: error no pco for sys=%d svn=%d freq=%.0lf idx=%d\n", pcv->satsys, pcv->svn, freq, idx); + + for (int i = 0; i < 3; i++) pco[i] = pcv->off[idx][i]; + + if (idx2 >= 0 && freq1 > 0 && freq2 > 0) { + // Interpolate between frequencies. + double lam = CLIGHT / freq, lam1 = CLIGHT / freq1, lam2 = CLIGHT / freq2, dlam = lam2 - lam1; + double a = (lam - lam1) / dlam; + for (int j = 0; j < 3; j++) + pco[j] += -a * pcv->off[idx][j] + a * pcv->off[idx2][j]; + } - for (i=0;ivar[i]); + trace(4, "antpco: sys=%d svn=%d type=%s pco=%6.4f %6.4f %6.4f\n", pcv->satsys, pcv->svn, pcv->type, pco[0], pco[1], pco[2]); + return; +} + +// Copy the pcv structure and the var[] data. +extern void copy_pcv(pcv_t *dst, const pcv_t *src) { + memcpy(dst, src, sizeof(pcv_t)); + for (int i = 0; i < ANTNFREQ; i++) { + if (dst->zen_len[i] > 0 && dst->azi_len[i] > 0) { + size_t size = dst->azi_len[i] * dst->zen_len[i] * sizeof(double); + dst->var[i] = malloc(size); + if (!dst->var[i]) { + trace(1, "copy_pcv alloc failured\n"); + } else { + memcpy(dst->var[i], src->var[i], size); + } + } else if (dst->var[i]) { + trace(2, "copy_pcv: empty PHV, zen_len=%d azi_len=%d var=%p\n", dst->zen_len[i], dst->azi_len[i], (void *)dst->var[i]); } - trace(4,"antmodel_s: dant=%6.3f %6.3f\n",dant[0],dant[1]); + } +} +// Free all pcv var[] data. +extern void free_pcv(pcv_t *pcv) { + for (int i = 0; i < ANTNFREQ; i++) { + free(pcv->var[i]); + pcv->var[i] = NULL; + } +} +// Free all pcvs var[] data, and the pcv array. +void free_pcvs(pcvs_t *pcvs) { + for (int i = 0; i < pcvs->n; i++) free_pcv(&pcvs->pcv[i]); + free(pcvs->pcv); + pcvs->pcv = NULL; + pcvs->n = pcvs->nmax = 0; } + /* sun and moon position in eci (ref [4] 5.1.1, 5.2.1) -----------------------*/ static void sunmoonpos_eci(gtime_t tut, double *rsun, double *rmoon) { diff --git a/src/rtklib.h b/src/rtklib.h index 1a073dab2..297f34da1 100644 --- a/src/rtklib.h +++ b/src/rtklib.h @@ -602,14 +602,51 @@ typedef struct { /* earth rotation parameter type */ erpd_t *data; /* earth rotation parameter data */ } erp_t; -typedef struct { /* antenna parameter type */ - int sat; /* satellite number (0:receiver) */ - char type[MAXANT]; /* antenna type */ - char code[MAXANT]; /* serial number or satellite code */ - gtime_t ts,te; /* valid time start and end */ - double off[NFREQ][ 3]; /* phase center offset e/n/u or x/y/z (m) */ - double var[NFREQ][19]; /* phase center variation (m) */ - /* el=90,85,...,0 or nadir=0,1,2,3,... (deg) */ +#define ANTFREQL1 0 // FREQL1 1.57542E9 +#define ANTFREQL2 1 // FREQL2 1.22760E9 +#define ANTFREQL5 2 // FREQL5 1.17645E9 +#define ANTFREQL6 3 // FREQL6 1.27875E9 +#define ANTFREQL7 4 // FREQE5b 1.20714E9 +#define ANTFREQL8 5 // FREQE5ab 1.191795E9 +#define ANTFREQG1 6 // FREQ1_GLO 1.60200E9 +#define ANTFREQG2 7 // FREQ2_GLO 1.24600E9 +#define ANTFREQG3 8 // FREQ3_GLO 1.202025E9 +#define ANTFREQG4 9 // FREQ1a_GLO 1.600995E9 +#define ANTFREQG6 10 // FREQ2a_GLO 1.248060E9 +#define ANTFREQB1 11 // FREQ1_CMP 1.561098E9 +#define ANTFREQB3 12 // FREQ3_CMP 1.26852E9 +#define ANTFREQI9 13 // FREQs 2.492028E9 +#define ANTNFREQ 14 + +typedef struct { // Satellite meta data, svn to sat mapping + gtime_t ts, te; // Valid time start and end + int sat; // Satellite number; system and PRN. + int svn; // Satellite SVN number, for the sat system. +} satsvn_t; + +typedef struct { // Satellite meta data, svn to sat mapping + int n,nmax; // Number of data/allocated + satsvn_t *satsvn; // Data +} satsvns_t; + +typedef struct { // Antenna parameter type + int sat; // Satellite number (0:receiver) + int satsys; // Satellite system + int svn; // Satellite SVN number (0:receiver) + char type[MAXANT]; // Antenna type + char code[MAXANT]; // Serial number or satellite code + gtime_t ts, te; // Valid time start and end + double zen1[ANTNFREQ]; // Zenith range, start. + double zen2[ANTNFREQ]; // Zenith range, end. + double dzen[ANTNFREQ]; // Zenith range, delta. + int zen_len[ANTNFREQ]; // Number of zenith elements in var[] below. + double dazi[ANTNFREQ]; // Azimuth range delta + int azi_len[ANTNFREQ]; // Number of azimuth elements in var[] below, including the NOAZI entry. + int init[ANTNFREQ]; // Frequency initialization state: bit 0 for the offset; bit 1 for the noazi entry; bit 2 for the azi entries. + double off[ANTNFREQ][3]; // Phase center offset e/n/u or x/y/z (m) + // The azimuth / zenith values. The first azimuth entry is for the NOAZI / zenith. + // el=90,85,...,0 or nadir=0,1,2,3,... (deg) + double *var[ANTNFREQ]; // Phase center variation (m) } pcv_t; typedef struct { /* antenna parameters type */ @@ -1091,6 +1128,7 @@ typedef struct { /* solution options type */ } solopt_t; typedef struct { /* file options type */ + char satmeta[MAXSTRPATH]; /* Satellite meta data sinex file */ char satantp[MAXSTRPATH]; /* satellite antenna parameters file */ char rcvantp[MAXSTRPATH]; /* receiver antenna parameters file */ char stapos [MAXSTRPATH]; /* station positions file */ @@ -1324,6 +1362,7 @@ typedef struct { /* RTK server type */ char cmds_periodic[3][MAXRCVCMD]; /* periodic commands */ char cmd_reset[MAXRCVCMD]; /* reset command */ double bl_reset; /* baseline length to reset (km) */ + pcvs_t pcvsr; // Receiver antenna parameters. rtklib_lock_t lock; /* lock flag */ } rtksvr_t; @@ -1416,6 +1455,7 @@ EXPORT void matfprint(const double *A, int n, int m, int p, int q, FILE *fp); EXPORT void add_fatal(fatalfunc_t *func); /* time and string functions -------------------------------------------------*/ +EXPORT void setstr(char *dst, const char *src, int n); EXPORT double str2num(const char *s, int i, int n); EXPORT int str2time(const char *s, int i, int n, gtime_t *t); EXPORT char *time2str(gtime_t t, char str[40], int n); @@ -1555,12 +1595,20 @@ EXPORT int tropcorr(gtime_t time, const nav_t *nav, const double *pos, EXPORT int seliflc(int optnf, int sys); /* antenna models ------------------------------------------------------------*/ -EXPORT int readpcv(const char *file, pcvs_t *pcvs); -EXPORT pcv_t *searchpcv(int sat, const char *type, gtime_t time, - const pcvs_t *pcvs); -EXPORT void antmodel(const pcv_t *pcv, const double *del, const double *azel, - int opt, double *dant); -EXPORT void antmodel_s(const pcv_t *pcv, double nadir, double *dant); +EXPORT int readsinex(const char *file, satsvns_t *satsvns); +EXPORT int searchsvnsat(int sys, int svn, const gtime_t time, const satsvns_t *satsvns); +EXPORT int searchsatsvn(int sat, const gtime_t time, const satsvns_t *satsvns); +EXPORT int readpcv(const char *file, int filter, pcvs_t *pcvs); +EXPORT pcv_t *searchpcv(int sat, const char *type, gtime_t time, const satsvns_t *satsvns, const pcvs_t *pcvs); +EXPORT void copy_pcv(pcv_t *dst, const pcv_t *src); +EXPORT void free_pcv(pcv_t *pcv); +EXPORT void free_pcvs(pcvs_t *pcvs); +EXPORT double antmodel(const pcv_t *pcv, const double *del, const double *azel, int opt, double freq); +EXPORT double antmodel_s(const pcv_t *pcv, double nadir, double freq); +EXPORT void antpco(const pcv_t *pcv, double freq, double pco[3]); +EXPORT double antpcv(const pcv_t *pcv, const double *azel, double freq); +EXPORT double antcode2freq(int code); +EXPORT const char *antcode2id(int code); /* earth tide models ---------------------------------------------------------*/ EXPORT void sunmoonpos(gtime_t tutc, const double *erpv, double *rsun, @@ -1628,7 +1676,7 @@ EXPORT void satposs(gtime_t time, const obsd_t *obs, int n, const nav_t *nav, EXPORT void setseleph(int sys, int sel); EXPORT int getseleph(int sys); EXPORT void readsp3(const char *file, nav_t *nav, int opt); -EXPORT int readsap(const char *file, gtime_t time, nav_t *nav); +EXPORT int readsap(const char *file, const gtime_t time, const satsvns_t *satsvns, nav_t *nav); EXPORT int readdcb(const char *file, nav_t *nav, const sta_t *sta); EXPORT int code2bias_ix(const int sys,const int code); /*EXPORT int readfcb(const char *file, nav_t *nav);*/ @@ -1816,7 +1864,7 @@ EXPORT int lambda_search(int n, int m, const double *a, const double *Q, /* standard positioning ------------------------------------------------------*/ EXPORT int pntpos(const obsd_t *obs, int n, const nav_t *nav, - const prcopt_t *opt, sol_t *sol, double *azel, + const prcopt_t *opt, int base, sol_t *sol, double *azel, ssat_t *ssat, char *msg); /* precise positioning -------------------------------------------------------*/ diff --git a/src/rtkpos.c b/src/rtkpos.c index 800be2cf9..bd69cfd5a 100644 --- a/src/rtkpos.c +++ b/src/rtkpos.c @@ -54,6 +54,7 @@ #define STD_PREC_VAR_THRESH 0 /* pos variance threshold to skip standard precision */ /* solution: 0 = run every epoch, */ /* 0.5 = skip except for first*/ +#define REL_HUMI 0.7 /* Relative humidity for Saastamoinen model */ /* constants/macros ----------------------------------------------------------*/ @@ -648,7 +649,7 @@ static void udrcvbias(rtk_t *rtk, double tt) } } // Detect a change in the observation code for a given frequency index. -// Only one bias per frequency index per satallite is supported, so if the +// Only one bias per frequency index per satellite is supported, so if the // observation code changes then consider it a slip. static void detslp_code(rtk_t *rtk, const obsd_t *obs, int i, int rcv) { int sat = obs[i].sat; @@ -952,56 +953,10 @@ static void udstate(rtk_t *rtk, const obsd_t *obs, const int *sat, udbias(rtk,tt,obs,sat,iu,ir,ns,nav); } } -/* UD (undifferenced) phase/code residual for satellite ----------------------*/ -static void zdres_sat(int base, double r, const obsd_t *obs, const nav_t *nav, - const double *azel, const double *dant, - const prcopt_t *opt, double *y, double *freq) -{ - double freq1,freq2,C1,C2,dant_if; - int i,nf=NF(opt),f2; - - if (opt->ionoopt==IONOOPT_IFLC) { /* iono-free linear combination */ - freq1=sat2freq(obs->sat,obs->code[0],nav); - f2=seliflc(opt->nf,satsys(obs->sat,NULL)); - freq2=sat2freq(obs->sat,obs->code[f2],nav); - - if (freq1==0.0||freq2==0.0) return; - - if (testsnr(base,0,azel[1],obs->SNR[0]*SNR_UNIT,&opt->snrmask)|| - testsnr(base,f2,azel[1],obs->SNR[f2]*SNR_UNIT,&opt->snrmask)) return; - - C1= SQR(freq1)/(SQR(freq1)-SQR(freq2)); - C2=-SQR(freq2)/(SQR(freq1)-SQR(freq2)); - dant_if=C1*dant[0]+C2*dant[f2]; - - if (obs->L[0]!=0.0&&obs->L[f2]!=0.0) { - y[0]=C1*obs->L[0]*CLIGHT/freq1+C2*obs->L[f2]*CLIGHT/freq2-r-dant_if; - } - if (obs->P[0]!=0.0&&obs->P[f2]!=0.0) { - y[nf]=C1*obs->P[0]+C2*obs->P[f2]-r-dant_if; - } - freq[0]=1.0; - } - else { - for (i=0;isat,obs->code[i],nav))==0.0) continue; - - /* check SNR mask */ - if (testsnr(base,i,azel[1],obs->SNR[i]*SNR_UNIT,&opt->snrmask)) { - continue; - } - /* residuals = observable - estimated range */ - if (obs->L[i]!=0.0) y[i ]=obs->L[i]*CLIGHT/freq[i]-r-dant[i]; - if (obs->P[i]!=0.0) y[i+nf]=obs->P[i] -r-dant[i]; - trace(4,"zdres_sat: %d: L=%.6f P=%.6f r=%.6f f=%.0f\n",obs->sat,obs->L[i], - obs->P[i],r,freq[i]); - } - } -} -/* undifferenced phase/code residuals ---------------------------------------- - calculate zero diff residuals [observed pseudorange - range] +/* Undifferenced phase/code residuals ---------------------------------------- + Calculate zero diff residuals [observed pseudorange - range] output is in y[0:nu-1], only shared input with base is nav - args: I base: 1=base,0=rover + Args: I base: 1=base,0=rover I obs = sat observations I n = # of sats I rs [(0:2)+i*6]= sat position {x,y,z} (m) @@ -1014,73 +969,194 @@ static void zdres_sat(int base, double r, const obsd_t *obs, const nav_t *nav, O y[(0:1)+i*2] = zero diff residuals {phase,code} (m) O e = line of sight unit vectors to sats O azel = [az, el] to sats */ -static int zdres(int base, const obsd_t *obs, int n, const double *rs, - const double *dts, const double *var, const int *svh, - const nav_t *nav, const double *rr, const prcopt_t *opt, - double *y, double *e, double *azel, double *freq) -{ - double r,rr_[3],pos[3],dant[NFREQ]={0},disp[3]; - double mapfh,zhd,zazel[]={0.0,90.0*D2R}; - int i,nf=NF(opt); +static int zdres(int base, const obsd_t *obs, int n, const double *rs, const double *dts, + const double *var, const int *svh, const nav_t *nav, const double *rr, + const prcopt_t *opt, double *y, double *e, double *azel, double *freq) { + trace(3, "zdres : n=%d rr=%.2f %.2f %.2f\n", n, rr[0], rr[1], rr[2]); - trace(3,"zdres : n=%d rr=%.2f %.2f %.2f\n",n,rr[0], rr[1], rr[2]); + int nf = NF(opt); - /* init residuals to zero */ - for (i=0;itidecorr) { - tidedisp(gpst2utc(obs[0].time),rr_,opt->tidecorr,&nav->erp, - opt->odisp[base],disp); - for (i=0;i<3;i++) rr_[i]+=disp[i]; - } - /* translate rcvr pos from ecef to geodetic */ - ecef2pos(rr_,pos); - - /* loop through satellites */ - for (i=0;ielmin) continue; - - /* excluded satellite? */ - if (satexclude(obs[i].sat,var[i],svh[i],opt)) continue; - - /* adjust range for satellite clock-bias */ - r+=-CLIGHT*dts[i*2]; - - /* adjust range for troposphere delay model (hydrostatic) */ - zhd=tropmodel(obs[0].time,pos,zazel,0.0); - mapfh=tropmapf(obs[i].time,pos,azel+i*2,NULL); - r+=mapfh*zhd; + /* Adjust rcvr pos for earth tide correction */ + if (opt->tidecorr) { + double disp[3]; + tidedisp(gpst2utc(obs[0].time), rr_, opt->tidecorr, &nav->erp, opt->odisp[base], disp); + for (int i = 0; i < 3; i++) rr_[i] += disp[i]; + } + /* Translate rcvr pos from ECEF to geodetic */ + double pos[3]; + ecef2pos(rr_, pos); + + /* Loop through satellites */ + for (int i = 0; i < n; i++) { + int sat = obs[i].sat; + + if (opt->ionoopt == IONOOPT_IFLC) { /* Iono-free linear combination */ + int sys = satsys(sat, NULL); + int code1 = obs[i].code[0]; + double freq1 = sat2freq(sat, code1, nav); + int f2 = seliflc(opt->nf, sys); + int code2 = obs[i].code[f2]; + double freq2 = sat2freq(sat, code2, nav); + + if (freq1 == 0.0 || freq2 == 0.0) continue; + double C1 = SQR(freq1) / (SQR(freq1) - SQR(freq2)); + double C2 = -SQR(freq2) / (SQR(freq1) - SQR(freq2)); + + // Receiver antenna phase center offset. + double pco[3]; // ENU + double pco1[3], pco2[3]; + antpco(opt->pcvr + base, freq1, pco1); + antpco(opt->pcvr + base, freq2, pco2); + for (int j = 0; j < 3; j++) pco[j] = C1 * pco1[j] + C2 * pco2[j]; + + // Add in the antenna delta, which includes the offset from the marker + // to the antenna reference point plus eccentricities. + for (int j = 0; j < 3; j++) pco[j] += opt->antdel[base][j]; + double dr[3]; + enu2ecef(pos, pco, dr); + double rpc[3]; + for (int j = 0; j < 3; j++) rpc[j] = rr_[j] + dr[j]; + // Recalculate the position, now for the phase center. + double rpc_pos[3]; + ecef2pos(rpc, rpc_pos); + + // Compute geometric-range and azimuth/elevation angle + double r = geodist(rs + i * 6, rpc, e + i * 3); + if (r <= 0.0) continue; + if (satazel(rpc_pos, e + i * 3, azel + i * 2) < opt->elmin) continue; + + // Check SNR mask + if (testsnr(base, 0, azel[1 + i * 2], obs[i].SNR[0], &opt->snrmask) || + testsnr(base, f2, azel[1 + i * 2], obs[i].SNR[f2], &opt->snrmask)) + continue; + + // Excluded satellite? + if (satexclude(sat, var[i], svh[i], opt)) continue; + + // Adjust range for satellite clock-bias. + r += -CLIGHT * dts[i * 2]; + + // Adjust range for troposphere delay model. + double dtrp = 0.0; + if (opt->tropopt <= TROPOPT_SAAS) { + dtrp = tropmodel(obs[0].time, rpc_pos, azel + i * 2, REL_HUMI); + } else if (opt->tropopt == TROPOPT_SBAS) { + double vart; + dtrp = sbstropcorr(obs[0].time, rpc_pos, azel + i * 2, &vart); + } else if (opt->tropopt >= TROPOPT_EST) { + // Hydrostatic only. + dtrp = tropmodel(obs[0].time, rpc_pos, azel + i * 2, 0.0); + } - /* calc receiver antenna phase center correction */ - antmodel(opt->pcvr+base,opt->antdel[base],azel+i*2,opt->posopt[1], - dant); + // Calc receiver antenna phase center variation, from the phase center + double dant_if = 0.0; + if (opt->posopt[1]) { + double dant1 = antpcv(opt->pcvr + base, azel + i * 2, freq1); + double dant2 = antpcv(opt->pcvr + base, azel + i * 2, freq2); + dant_if = C1 * dant1 + C2 * dant2; + } - /* calc undifferenced phase/code residual for satellite */ - trace(4,"sat=%d r=%.6f c*dts=%.6f zhd=%.6f map=%.6f\n",obs[i].sat,r,CLIGHT*dts[i*2],zhd,mapfh); - zdres_sat(base,r,obs+i,nav,azel+i*2,dant,opt,y+i*nf*2,freq+i*nf); + if (obs[i].L[0] != 0.0 && obs[i].L[f2] != 0.0) + y[0 + i * nf * 2] = C1 * obs[i].L[0] * CLIGHT / freq1 + C2 * obs[i].L[f2] * CLIGHT / freq2 - (r + dant_if + dtrp); + if (obs->P[0] != 0.0 && obs->P[f2] != 0.0) + y[0 + nf + i * nf * 2] = C1 * obs[i].P[0] + C2 * obs[i].P[f2] - (r + dant_if + dtrp); + freq[i * nf] = 1.0; + trace(4, "zdres: sat=%d L=%.6f P=%.6f r=%.6f c*dts=%.6f dtrp=%.6f dant=%.6lf\n", sat, obs[i].L[0], obs[i].P[0], r, CLIGHT * dts[i * 2], dtrp, dant_if); } - trace(4,"rr_=%.3f %.3f %.3f\n",rr_[0],rr_[1],rr_[2]); - trace(4,"pos=%.9f %.9f %.3f\n",pos[0]*R2D,pos[1]*R2D,pos[2]); - for (i=0;ipcvr + base, frq, pco); + + // Add in the antenna delta, which includes the offset from the marker + // to the antenna reference point plus eccentricities. + for (int j = 0; j < 3; j++) pco[j] += opt->antdel[base][j]; + double dr[3]; + enu2ecef(pos, pco, dr); + double rpc[3]; + for (int j = 0; j < 3; j++) rpc[j] = rr_[j] + dr[j]; + // Recalculate the position, now for the phase center. + double rpc_pos[3]; + ecef2pos(rpc, rpc_pos); + + // Compute geometric-range and azimuth/elevation angle + double r = geodist(rs + i * 6, rpc, e + f * 3 + i * nf * 3); + if (r <= 0.0) continue; + if (satazel(rpc_pos, e + f * 3 + i * nf * 3, azel + i * 2) < opt->elmin) continue; + + // Check SNR mask + if (testsnr(base, f, azel[1 + i * 2], obs[i].SNR[f], &opt->snrmask)) continue; + + // Excluded satellite? + if (satexclude(sat, var[i], svh[i], opt)) continue; + + // Adjust range for satellite clock-bias + r += -CLIGHT * dts[i * 2]; + + // Adjust range for troposphere delay model. + double dtrp = 0.0; + if (opt->tropopt <= TROPOPT_SAAS) { + dtrp = tropmodel(obs[0].time, rpc_pos, azel + i * 2, REL_HUMI); + } else if (opt->tropopt == TROPOPT_SBAS) { + double vart; + dtrp = sbstropcorr(obs[0].time, rpc_pos, azel + i * 2, &vart); + } else if (opt->tropopt >= TROPOPT_EST) { + // Hydrostatic only. + dtrp = tropmodel(obs[0].time, rpc_pos, azel + i * 2, 0.0); + } + + // Ionospheric correction. + double dion = 0.0, vion = 0.0; +#ifdef RTK_DISABLED + // TODO can this help for longer baselines? + ionocorr(obs[0].time, nav, sat, rpc_pos, azel + i * 2, opt->ionoopt, &dion, &vion); +#endif + // Scale iono delay for frequency. The slant factor has already been applied. + double C = SQR(FREQL1 / frq); + + // Calc receiver antenna phase center variation, from the phase center. + double dant = 0.0; + if (opt->posopt[1]) dant = antpcv(opt->pcvr + base, azel + i * 2, frq); + + // Calc undifferenced phase/code residual for satellite + // Residuals = observable - estimated range + if (obs[i].L[f] != 0.0) + y[f + i * nf * 2] = obs[i].L[f] * CLIGHT / frq - (r + dant + dtrp - C * dion); + if (obs[i].P[f] != 0.0) + y[f + nf + i * nf * 2] = obs[i].P[f] - (r + dant + dtrp + C * dion); + trace(4, "zdres: sat=%d f=%d frq=%.0f L=%.6f P=%.6f r=%.6f c*dts=%.6f dtrp=%.6f dion=%.6f dant=%.6lf\n", sat, f, frq, obs[i].L[f], obs[i].P[f], r, CLIGHT * dts[i * 2], dtrp, C * dion, dant); + } } - trace(3,"y=\n"); tracemat(3,y,nf*2,n,13,3); + } + trace(4, "rr_=%.3f %.3f %.3f\n", rr_[0], rr_[1], rr_[2]); + trace(4, "pos=%.9f %.9f %.3f\n", pos[0] * R2D, pos[1] * R2D, pos[2]); + for (int i = 0; i < n; i++) { + if ((obs[i].L[0] == 0 && obs[i].L[1] == 0 && obs[i].L[2] == 0) || base == 0) continue; + trace(3, "sat=%2d rs=%13.3f %13.3f %13.3f dts=%13.10f az=%6.1f el=%5.1f\n", obs[i].sat, + rs[i * 6], rs[1 + i * 6], rs[2 + i * 6], dts[i * 2], azel[i * 2] * R2D, + azel[1 + i * 2] * R2D); + } + trace(3, "y=\n"); + tracemat(3, y, nf * 2, n, 13, 3); - return 1; + return 1; } /* test valid observation data -----------------------------------------------*/ -static int validobs(int i, int j, int f, int nf, double *y) +static int validobs(int i, int j, int f, int nf, const double *y) { /* check for valid residuals */ return y[f+i*nf*2]!=0.0&&y[f+j*nf*2]!=0.0; @@ -1197,8 +1273,8 @@ static int test_sys(int sys, int m) I P = error covariance matrix of float states I sat = list of common sats I y = zero diff residuals (code and phase, base and rover) - I e = line of sight unit vectors to sats - I azel = [az, el] to sats + I e = line of sight unit vectors to sats, from phase center + I azel = [az, el] to sats, from phase center I iu,ir = user and ref indices to sats I ns = # of sats O v = double diff innovations (measurement-model) (phase and code) @@ -1206,8 +1282,8 @@ static int test_sys(int sys, int m) O R = measurement error covariances O vflg = bit encoded list of sats used for each double diff */ static int ddres(rtk_t *rtk, const obsd_t *obs, double dt, const double *x, - const double *P, const int *sat, double *y, double *e, - double *azel, double *freq, const int *iu, const int *ir, + const double *P, const int *sat, const double *y, const double *e, + const double *azel, const double *freq, const int *iu, const int *ir, int ns, double *v, double *H, double *R, int *vflg) { prcopt_t *opt=&rtk->opt; @@ -1280,10 +1356,11 @@ static int ddres(rtk_t *rtk, const obsd_t *obs, double dt, const double *x, v[nv]=(y[f+iu[i]*nf*2]-y[f+ir[i]*nf*2])- (y[f+iu[j]*nf*2]-y[f+ir[j]*nf*2]); - /* partial derivatives by rover position, combine unit vectors from two sats */ + /* Partial derivatives by rover position, combine unit vectors from two sats */ if (H) { for (k=0;k<3;k++) { - Hi[k]=-e[k+iu[i]*3]+e[k+iu[j]*3]; /* translation of innovation to position states */ + /* Translation of innovation to position states */ + Hi[k]=-e[k+frq*3+iu[i]*nf*3]+e[k+frq*3+iu[j]*nf*3]; } } if (opt->ionoopt==IONOOPT_EST) { @@ -1433,7 +1510,7 @@ static double intpres(gtime_t time, const obsd_t *obs, int n, const nav_t *nav, { static obsd_t obsb[MAXOBS]; static double yb[MAXOBS*NFREQ*2],rs[MAXOBS*6],dts[MAXOBS*2],var[MAXOBS]; - static double e[MAXOBS*3],azel[MAXOBS*2],freq[MAXOBS*NFREQ]; + static double e[MAXOBS*NFREQ*3],azel[MAXOBS*2],freq[MAXOBS*NFREQ]; static int nb=0,svh[MAXOBS*2]; prcopt_t *opt=&rtk->opt; double tt,ttb,*p,*q; @@ -1948,7 +2025,7 @@ static int relpos(rtk_t *rtk, const obsd_t *obs, int nu, int nr, dts=mat(2,n); /* satellite clock biases */ var=mat(1,n); y=mat(nf*2,n); - e=mat(3,n); + e=mat(3*nf,n); azel=zeros(2,n); /* [az, el] */ freq=zeros(nf,n); @@ -1968,7 +2045,7 @@ static int relpos(rtk_t *rtk, const obsd_t *obs, int nu, int nr, output is in y[nu:nu+nr], see call for rover below for more details */ trace(3,"base station:\n"); if (!zdres(1,obs+nu,nr,rs+nu*6,dts+nu*2,var+nu,svh+nu,nav,rtk->rb,opt, - y+nu*nf*2,e+nu*3,azel+nu*2,freq+nu*nf)) { + y+nu*nf*2,e+nu*nf*3,azel+nu*2,freq+nu*nf)) { errmsg(rtk,"initial base station position error\n"); free(rs); free(dts); free(var); free(y); free(e); free(azel); free(freq); @@ -2322,7 +2399,7 @@ extern int rtkpos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav) /* rover position and time by single point positioning, skip if position variance smaller than threshold */ if (rtk->P[0]==0||rtk->P[0]>STD_PREC_VAR_THRESH) { - if (!pntpos(obs,nu,nav,&rtk->opt,&rtk->sol,NULL,rtk->ssat,msg)) { + if (!pntpos(obs,nu,nav,&rtk->opt,0,&rtk->sol,NULL,rtk->ssat,msg)) { errmsg(rtk,"point pos error (%s)\n",msg); if (!rtk->opt.dynamics) { @@ -2369,7 +2446,7 @@ extern int rtkpos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav) /* estimate position/velocity of base station, skip if position variance below threshold*/ if (rtk->P[0]==0||rtk->P[0]>STD_PREC_VAR_THRESH) { - if (!pntpos(obs+nu,nr,nav,&rtk->opt,&solb,NULL,NULL,msg)) { + if (!pntpos(obs+nu,nr,nav,&rtk->opt,1,&solb,NULL,NULL,msg)) { errmsg(rtk,"base station position error (%s)\n",msg); return 0; } diff --git a/src/rtksvr.c b/src/rtksvr.c index d765bb772..761ab3d7b 100644 --- a/src/rtksvr.c +++ b/src/rtksvr.c @@ -231,42 +231,57 @@ static void update_ionutc(rtksvr_t *svr, nav_t *nav, int index) } svr->nmsg[index][2]++; } -/* update antenna position ---------------------------------------------------*/ -static void update_antpos(rtksvr_t *svr, int index) -{ - sta_t *sta; - double pos[3],del[3]={0},dr[3]; - int i; +// Update antenna position --------------------------------------------------- +static void update_antpos(rtksvr_t *svr, int index) { + sta_t *sta; + if (svr->format[index] == STRFMT_RTCM2 || svr->format[index] == STRFMT_RTCM3) { + sta = &svr->rtcm[index].sta; + } else { + sta = &svr->raw[index].sta; + } + if (index == 1 && svr->rtk.opt.refpos == POSOPT_RTCM) { + // Update base station position. + for (int i = 0; i < 3; i++) svr->rtk.rb[i] = sta->pos[i]; + } + if (index == 0 && svr->rtk.opt.rovpos == POSOPT_RTCM && + (svr->rtk.opt.mode == PMODE_FIXED || svr->rtk.opt.mode == PMODE_PPP_FIXED)) { + // Update rover fixed position. + for (int i = 0; i < 3; i++) svr->rtk.opt.ru[i] = sta->pos[i]; + } - if (svr->rtk.opt.refpos==POSOPT_RTCM&&index==1) { - if (svr->format[1]==STRFMT_RTCM2||svr->format[1]==STRFMT_RTCM3) { - sta=&svr->rtcm[1].sta; - } - else { - sta=&svr->raw[1].sta; - } - /* update base station position */ - for (i=0;i<3;i++) { - svr->rtk.rb[i]=sta->pos[i]; - } - /* antenna delta */ - ecef2pos(svr->rtk.rb,pos); - if (sta->deltype) { /* xyz */ - del[2]=sta->hgt; - enu2ecef(pos,del,dr); - for (i=0;i<3;i++) { - svr->rtk.rb[i]+=sta->del[i]+dr[i]; - } - } - else { /* enu */ - enu2ecef(pos,sta->del,dr); - for (i=0;i<3;i++) { - svr->rtk.rb[i]+=dr[i]; - } - } - } - svr->nmsg[index][4]++; + // Antenna type and delta. These are updated independently of the antenna + // marker position above when the anttype is "*". + if (strcmp(svr->rtk.opt.anttype[index], "*") == 0) { + if (sta->antdes[0] != '\0' && strcmp(svr->rtk.opt.pcvr[index].type, sta->antdes) != 0) { + // Antenna type is to be set from the RTCM stream, and does not match + // the current pcv_t type, so search for this pcv. + pcv_t *pcv = searchpcv(0, sta->antdes, utc2gpst(timeget()), NULL, &svr->pcvsr); + if (!pcv) { + tracet(2, "no index=%d antenna '%s'", index, sta->antdes); + } else { + tracet(2, "updated index=%d antenna to '%s'", index, sta->antdes); + copy_pcv(&svr->rtk.opt.pcvr[index], pcv); + } } + // Update the delta from the marker position to the antenna ARP, taking + // into account the RCTM antenna height. This overrides any config delta + // values which are ignored in this path. + if (sta->deltype == 1) { // XYZ + // Convert to the antdel ENU, adding the height. + // Need at least an approx position to map the delta. + if (norm(sta->pos, 3) > 0.0) { + double pos[3]; + ecef2pos(sta->pos, pos); + ecef2enu(pos, sta->del, svr->rtk.opt.antdel[index]); + svr->rtk.opt.antdel[index][2] += sta->hgt; + } + } else { // ENU + for (int i = 0; i < 3; i++) svr->rtk.opt.antdel[index][i] = sta->del[i]; + svr->rtk.opt.antdel[index][2] += sta->hgt; + } + } + svr->nmsg[index][4]++; +} /* update ssr corrections ----------------------------------------------------*/ static void update_ssr(rtksvr_t *svr, int index) { @@ -618,7 +633,7 @@ static void *rtksvrthread(void *arg) if (fobs[1]>0&&svr->rtk.opt.refpos==POSOPT_SINGLE) { if ((svr->rtk.opt.maxaveep<=0||svr->navertk.opt.maxaveep)&& pntpos(svr->obs[1][0].data,svr->obs[1][0].n,&svr->nav, - &svr->rtk.opt,&sol,NULL,NULL,msg)) { + &svr->rtk.opt,1,&sol,NULL,NULL,msg)) { svr->nave++; for (i=0;i<3;i++) { svr->rb_ave[i]+=(sol.rr[i]-svr->rb_ave[i])/svr->nave;