diff --git a/app/consapp/rtkrcv/rtkrcv.c b/app/consapp/rtkrcv/rtkrcv.c index e909f5107..14286985a 100644 --- a/app/consapp/rtkrcv/rtkrcv.c +++ b/app/consapp/rtkrcv/rtkrcv.c @@ -362,41 +362,54 @@ 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) { + trace(3, "readant:\n"); + + const pcv_t pcv0 = {0}; + opt->pcvr[0] = opt->pcvr[1] = pcv0; + + if (*filopt.rcvantp) { + gtime_t time = timeget(); + pcvs_t pcvr = {0}; + if (readpcv(filopt.rcvantp, 2, &pcvr)) { + for (int i = 0; i < 2; i++) { + if (!*opt->anttype[i]) continue; + pcv_t *pcv = searchpcv(0, opt->anttype[i], time, NULL, &pcvr); + 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); + } + free_pcvs(&pcvr); + } 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,7 +453,7 @@ 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); @@ -523,6 +536,9 @@ 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]); + for (int i = 0; i < MAXSAT; i++) free_pcv(&svr.nav.pcvs[i]); + vt_printf(vt,"stop rtk server\n"); } /* print time ----------------------------------------------------------------*/ 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 80679b77a..22d4cadcd 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,10 +660,11 @@ 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', 3).arg(pco[1], 0, 'f', 3).arg(pco[2], 0, 'f', 3)); } del = rtk.opt.antdel[0]; @@ -673,10 +674,11 @@ void MonitorDialog::showRtk() 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', 3).arg(pco[1], 0, 'f', 3).arg(pco[2], 0, 'f', 3)); } del = rtk.opt.antdel[1]; diff --git a/app/qtapp/rtknavi_qt/navimain.cpp b/app/qtapp/rtknavi_qt/navimain.cpp index 42aa2b65a..439706eba 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,24 @@ 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 (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)); 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 diff --git a/app/winapp/rtknavi/navimain.cpp b/app/winapp/rtknavi/navimain.cpp index c005beb83..89579ced1 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"); @@ -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); @@ -1203,19 +1205,28 @@ void __fastcall TMainForm::SvrStart(void) for (i=0;i<3;i++) PrcOpt.antdel[1][i]=RefAntDel[i]; } if (RovAntPcvF||RefAntPcvF) { - free(pcvr.pcv); + free_pcvs(&pcvr); } - 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;iAdd(""); 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/rtkpost/postopt.cpp b/app/winapp/rtkpost/postopt.cpp index 350ae1c25..b97481c11 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/src/options.c b/src/options.c index 9f6f40075..55ec42e3d 100644 --- a/src/options.c +++ b/src/options.c @@ -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/postpos.c b/src/postpos.c index c1dffe76e..724488165 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 */ @@ -897,22 +898,30 @@ static int antpos(prcopt_t *opt, int rcvno, const obs_t *obs, const nav_t *nav, } /* 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)) { @@ -922,14 +931,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(); @@ -941,26 +955,28 @@ 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[64]; - /* 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 */ @@ -976,13 +992,14 @@ static void setpcv(gtime_t time, prcopt_t *popt, nav_t *nav, const pcvs_t *pcvs, 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 ----------------------------------------*/ @@ -1084,12 +1101,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 */ @@ -1111,7 +1129,7 @@ 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, + setpcv(obss.n>0?obss.data[0].time:timeget(),&popt_,&navs,&satsvns,&pcvss,&pcvsr, stas); } /* read ocean tide loading parameters */ @@ -1121,21 +1139,18 @@ static int execses(gtime_t ts, gtime_t te, double ti, const prcopt_t *popt, /* 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; + err = 1; + goto done; } if (!antpos(&popt_,2,&obss,&navs,stas,fopt->stapos)) { - freeobsnav(&obss,&navs); - free(rtk_ptr); - return 0; + 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; + err = 1; + goto done; } } /* open solution statistics */ @@ -1147,9 +1162,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); @@ -1224,10 +1238,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, @@ -1396,18 +1412,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; } } @@ -1490,7 +1506,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 9965e1e54..1e5e32aeb 100644 --- a/src/ppp.c +++ b/src/ppp.c @@ -401,55 +401,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) @@ -716,7 +733,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); @@ -751,7 +768,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; @@ -828,24 +845,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, @@ -936,7 +952,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[32]; int ne=0,obsi[MAXOBS*2*NFREQ]={0},frqi[MAXOBS*2*NFREQ],maxobs,maxfrq,rej; @@ -967,17 +982,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 1b8f4e292..c1e3e0f34 100644 --- a/src/preceph.c +++ b/src/preceph.c @@ -353,22 +353,23 @@ 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; +extern int readsap(const char *file, const gtime_t time, const satsvns_t *satsvns, nav_t *nav) { + trace(3, "readsap : file=%s time=%s\n", file, time_str(time,0)); - trace(3,"readsap : file=%s time=%s\n",file,time_str(time,0)); + pcvs_t pcvs = {0}; + if (!readpcv(file, 1, &pcvs)) return 0; - if (!readpcv(file,&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 @@ -683,80 +684,84 @@ 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; - - trace(4,"satantoff: time=%s sat=%2d\n",time_str(time,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) { + trace(4, "satantoff: time=%s sat=%2d\n", time_str(time,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/rtkcmn.c b/src/rtkcmn.c index fe6aa5beb..a66581ec8 100644 --- a/src/rtkcmn.c +++ b/src/rtkcmn.c @@ -2383,6 +2383,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 @@ -2485,206 +2613,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 @@ -3816,61 +4434,373 @@ extern double tropmapf(gtime_t time, const double pos[], const double azel[], return nmf(time,pos,azel,mapfw); /* NMF */ #endif } -/* interpolate antenna phase center variation --------------------------------*/ -static double interpvar(double ang, const double *var) -{ - double a=ang/5.0; /* ang=0-90 */ - int i=(int)a; - if (i<0) return var[0]; else if (i>=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; + } + } + } - trace(4,"antmodel: azel=%6.1f %4.1f opt=%d\n",azel[0]*R2D,azel[1]*R2D,opt); + 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; + } + + 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; + } - e[0]=sin(azel[0])*cosel; - e[1]=cos(azel[0])*cosel; - e[2]=sin(azel[1]); + // 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; + } - for (i=0;ioff[i][j]+del[j]; + // 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; + } - dant[i]=-dot3(off,e)+(opt?interpvar(90.0-azel[1]*R2D,pcv->var[i]):0.0); + 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; + } } - trace(4,"antmodel: dant=%6.3f %6.3f\n",dant[0],dant[1]); + 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; + } + + // 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, "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; + } - for (i=0;ivar[i]); + 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]; + } + + 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 8681cf56e..61d4ec772 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 */ @@ -1557,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, @@ -1630,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);*/ diff --git a/src/rtkpos.c b/src/rtkpos.c index 6558c8071..e1ac730c5 100644 --- a/src/rtkpos.c +++ b/src/rtkpos.c @@ -953,55 +953,66 @@ static void udstate(rtk_t *rtk, const obsd_t *obs, const int *sat, } } /* 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; +static void zdres_sat(int base, double r, const obsd_t *obs, const nav_t *nav, const double *azel, + const prcopt_t *opt, double *y, double *freq) { + int nf = NF(opt); + if (opt->ionoopt == IONOOPT_IFLC) { /* Iono-free linear combination */ + int sat = obs->sat; + int sys = satsys(sat, NULL); + int code1 = obs->code[0]; + double freq1 = sat2freq(sat, code1, nav); + int f2 = seliflc(opt->nf, sys); + int code2 = obs->code[f2]; + double freq2 = sat2freq(sat, code2, 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; + + double C1 = SQR(freq1) / (SQR(freq1) - SQR(freq2)); + double C2 = -SQR(freq2) / (SQR(freq1) - SQR(freq2)); + + /* Calc receiver antenna phase center correction */ + double dant1 = antmodel(opt->pcvr + base, opt->antdel[base], azel, opt->posopt[1], freq1); + double dant2 = antmodel(opt->pcvr + base, opt->antdel[base], azel, opt->posopt[1], freq2); + double dant_if = C1 * dant1 + C2 * dant2; + + 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; + return; + } - C1= SQR(freq1)/(SQR(freq1)-SQR(freq2)); - C2=-SQR(freq2)/(SQR(freq1)-SQR(freq2)); - dant_if=C1*dant[0]+C2*dant[f2]; + int sat = obs->sat; + for (int i = 0; i < nf; i++) { + int code = obs->code[i]; + freq[i] = sat2freq(sat, code, nav); + if (freq[i] == 0.0) continue; - 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; + /* Check SNR mask */ + if (testsnr(base, i, azel[1], obs->SNR[i] * SNR_UNIT, &opt->snrmask)) { + continue; } - else { - for (i=0;isat,obs->code[i],nav))==0.0) continue; + /* Calc receiver antenna phase center correction */ + double dant = antmodel(opt->pcvr + base, opt->antdel[base], azel, opt->posopt[1], freq[i]); - /* 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]); - } - } + /* Residuals = observable - estimated range */ + if (obs->L[i] != 0.0) y[i] = obs->L[i] * CLIGHT / freq[i] - r - dant; + if (obs->P[i] != 0.0) y[i + nf] = obs->P[i] - r - dant; + 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,70 +1025,68 @@ 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; - - /* calc receiver antenna phase center correction */ - antmodel(opt->pcvr+base,opt->antdel[base],azel+i*2,opt->posopt[1], - dant); - - /* 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); - } - 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;itidecorr) { + 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++) { + /* Compute geometric-range and azimuth/elevation angle */ + double r = geodist(rs + i * 6, rr_, e + i * 3); + if (r <= 0.0) continue; + if (satazel(pos, e + i * 3, azel + i * 2) < opt->elmin) 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) */ + double zazel[] = {0.0, 90.0 * D2R}; + double zhd = tropmodel(obs[0].time, pos, zazel, 0.0); + double mapfh = tropmapf(obs[i].time,pos, azel + i * 2, NULL); + r += mapfh * zhd; + trace(4, "sat=%d r=%.6f c*dts=%.6f zhd=%.6f map=%.6f\n", obs[i].sat, r, + CLIGHT * dts[i * 2], zhd, mapfh); + + /* Calc undifferenced phase/code residual for satellite */ + zdres_sat(base, r, obs + i, nav, azel + i * 2, opt, y + i * nf * 2, freq + i * nf); + } + 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)