From 6bf6d0b4015e77bee254d90b14553662a93e1de9 Mon Sep 17 00:00:00 2001 From: oaq Date: Wed, 20 Nov 2024 20:26:14 +1100 Subject: [PATCH] ANTEX: PCV with azimuth, full freq range, ANTEX 2, SINEX metadata Improve the ANTEX support. * Support PCV with azimuth. Was just NOAZI. Use linear interpolation between azimuth points, as well as the zenith. Satellite PCV still only supports NOAZI as that is most common, but all the support is there to extend this to vary with azimuth * Support ZEN1, ZEN2, DZEN, and DAZI. The PHV size is also now variable so the smallest size array is allocated. There was a fixed azimuth length of 19. It now correctly handles variations of these parameters. Require some more memory management, to free these. * Support the full range of frequencies, across systems. Was just GPS and NFREQ. The usage of the antmodel() etc functions needed to change. There were returning a NFREQ array with the offsets, but that did not work with multiple systems with different frequencies mapped to the frequency same indices - these functions are now single frequency functions. * Support PCO and PCV interpolation between frequency entries, with heuristics to avoid extrapolation of noisy data, and avoid interpolation across frequency bands (which may well be implemented with separate antenna elements e.g. stacked elements). * Add support for ANTEX 2 which avoids much redundancy. This required the PRN to SVN mapping in the SINEX satellite meta data. * Add minimal SINEX satellite meta data support, just for the PRN to SVN mappings. Adding the file-satmetafile option. * Support for the console apps. * Support for the qt apps, adding the satellite meta data file to the GUI. * Some windows app support, untested, and not yet the satellite meta data file for ANTEX 2 which needs adding to the GUI. --- app/consapp/rtkrcv/rtkrcv.c | 84 +- app/qtapp/appcmn_qt/navi_post_opt.cpp | 60 +- app/qtapp/appcmn_qt/navi_post_opt.h | 2 + app/qtapp/appcmn_qt/navi_post_opt.ui | 153 +-- app/qtapp/rtknavi_qt/mondlg.cpp | 22 +- app/qtapp/rtknavi_qt/navimain.cpp | 19 +- app/winapp/rtknavi/navimain.cpp | 33 +- app/winapp/rtknavi/naviopt.cpp | 6 +- app/winapp/rtkpost/postopt.cpp | 6 +- src/options.c | 1 + src/postpos.c | 102 +- src/ppp.c | 166 +-- src/preceph.c | 179 ++-- src/rtkcmn.c | 1350 +++++++++++++++++++++---- src/rtklib.h | 75 +- src/rtkpos.c | 231 +++-- 16 files changed, 1794 insertions(+), 695 deletions(-) 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 97c1768dd..72d30f14f 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 48b65c966..e5d29d560 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 ebced0c37..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, @@ -935,8 +951,7 @@ 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],dtrp=0.0,dion=0.0,vart=0.0,vari=0.0,dcb,freq; - double dantr[NFREQ]={0},dants[NFREQ]={0}; + double var[MAXOBS*2*NFREQ],dtrp=0.0,dion=0.0,vart=0.0,vari=0.0,dcb,freq; 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 ee38090b5..4aeef7569 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..8832cc000 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,359 @@ 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 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 offsets for the respective frequency (m) + * + * Note: this applies the PCO and PCV when available, and for the PCO only see + * antpco(). + * + *-----------------------------------------------------------------------------*/ +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 e[3], cosel = cos(azel[1]); + e[0] = sin(azel[0]) * cosel; + e[1] = cos(azel[0]) * cosel; + e[2] = sin(azel[1]); + + double freq1, freq2; + int idx2, idx = antpcoidx(pcv, freq, &freq1, &freq2, &idx2); + if (idx < 0) { + double dant = -dot3(del, e); + trace(3, "antmodel: no pco freq=%.0lf, dant=%6.3f\n", freq, dant); + return dant; + } + + if (!(pcv->init[idx] & PCV_PCO)) + trace(2, "antmodel: error pcv %d not init\n", idx); + + // Apply the PCO. + double off[3]; + for (int i = 0; i < 3; i++) off[i] = del[i] + pcv->off[idx][i]; + + if (idx2 >= 0 && freq2 > 0) { + // PCO with interpolation 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++) + off[j] += -a * pcv->off[idx][j] + a * pcv->off[idx2][j]; + } + + double dant = -dot3(off, e); + if (opt && pcv->init[idx] & (PCV_PHV | PCV_NOAZI) && pcv->dzen[idx] > 0.01) { + if (pcv->zen_len[idx] != (pcv->zen2[idx] - pcv->zen1[idx]) / pcv->dzen[idx] + 1) + trace(2, "antmodel: unexpected zen_len\n"); + 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, "antmodel azi %d >= len %d\n", i + 1, pcv->azi_len[idx] - 1); + if (i < 0) { + i = 0; + r = 0; + trace(2, "antmodel azi < 0 clipped\n"); + } else if (i >= pcv->azi_len[idx] - 1) { + i = pcv->azi_len[idx] - 2; + r = 0; + trace(2, "antmodel azi >= len clipped\n"); + } + // Interpolate for the zenith at each azimuth step. + double dant1 = interpvar(90.0 - azel[1] * R2D, pcv->zen1[idx], pcv->zen2[idx], pcv->dzen[idx], + pcv->var[idx] + (1 + i) * pcv->zen_len[idx]); + pcv1 += dant1; + if (i + 1 < pcv->azi_len[idx] - 1) { + // Interpolate for the zenith at each azimuth step. + double dant2 = interpvar(90.0 - azel[1] * R2D, pcv->zen1[idx], pcv->zen2[idx], pcv->dzen[idx], + pcv->var[idx] + (1 + i + 1) * pcv->zen_len[idx]); + // Interpolate for the azimuth. + pcv1 += -r * dant1 + r * dant2; + } + } + dant += pcv1; + + if (idx2 >= 0 && freq2 > 0 && pcv->init[idx2] & (PCV_PHV | PCV_NOAZI) && pcv->dzen[idx2] > 0.01 && pcv->dazi[idx2] > 0.01) { + // PCV with interpolation between frequencies. + if (pcv->zen_len[idx2] != (pcv->zen2[idx2] - pcv->zen1[idx2]) / pcv->dzen[idx2] + 1) + trace(2, "antmodel: 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, "antmodel azi %d >= len %d\n", i + 1, pcv->azi_len[idx2] - 1); + if (i < 0) { + i = 0; + r = 0; + trace(2, "antmodel azi < 0 clipped\n"); + } else if (i >= pcv->azi_len[idx2] - 1) { + i = pcv->azi_len[idx2] - 2; + r = 0; + trace(2, "antmodel 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, "antmodel: dant=%6.3f\n", dant); + return dant; + } - dant[i]=-dot3(off,e)+(opt?interpvar(90.0-azel[1]*R2D,pcv->var[i]):0.0); + // NOAZI + if (pcv->init[idx] & PCV_PHV) { + double pcv1 = interpvar(90.0 - azel[1] * R2D, pcv->zen1[idx], pcv->zen2[idx], pcv->dzen[idx], pcv->var[idx]); + dant += pcv1; + if (idx2 >= 0 && freq2 > 0 && pcv->init[idx2] & (PCV_PHV | PCV_NOAZI) && pcv->dzen[idx2] > 0.01) { + // PCV with interpolation between frequencies. + double pcv2 = interpvar(90.0 - azel[1] * R2D, pcv->zen1[idx2], pcv->zen2[idx2], pcv->dzen[idx2], pcv->var[idx2]); + double lam = CLIGHT / freq, lam1 = CLIGHT / freq1, lam2 = CLIGHT / freq2, dlam = lam2 - lam1; + double a = (lam - lam1) / dlam; + dant += -a * pcv1 + a * pcv2; + } } - trace(4,"antmodel: dant=%6.3f %6.3f\n",dant[0],dant[1]); + } + + trace(4, "antmodel: dant=%6.3f\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; + } - for (i=0;ivar[i]); + double freq1, freq2; + int idx2, idx = antpcoidx(pcv, freq, &freq1, &freq2, &idx2); + if (idx < 0 || !(pcv->init[idx] & PCV_NOAZI)) { + trace(2, "antmodel_s: NOAZI PHV not found for sys=%d svn=%d freq=%.0lf idx=%d\n", pcv->satsys, pcv->svn, freq, idx); + return 0.0; + } + + if (pcv->dzen[idx] < 0.01) { + trace(2, "antmodel_s: NOAZI PHV has no zenith entries for sys=%d svn=%d freq=%.0lf idx=%d\n", pcv->satsys, pcv->svn, freq, idx); + return 0.0; + } + + double dant = interpvar(nadir * R2D, pcv->zen1[idx], pcv->zen2[idx], pcv->dzen[idx], pcv->var[idx]); + trace(4, "antmodel_s: dant=%6.4f\n", dant); + return dant; +} +// Returns the antenna PCO +// For a receiver this is enu, for a satellite xyz. +extern void antpco(const pcv_t *pcv, double freq, double pco[3]) { + double freq1, freq2; + int idx2, idx = antpcoidx(pcv, freq, &freq1, &freq2, &idx2); + if (idx < 0) { + pco[0] = pco[1] = pco[2] = 0.0; + trace(4, "antpco PCO not found for sys=%d svn=%d freq=%.0lf idx=%d\n", pcv->satsys, pcv->svn, freq, idx); + return; + } + + 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: PCO sys=%d svn=%d type=%s %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 46716827d..50df4de53 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 */ @@ -1556,12 +1594,19 @@ 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 antcode2freq(int code); +EXPORT const char *antcode2id(int code); /* earth tide models ---------------------------------------------------------*/ EXPORT void sunmoonpos(gtime_t tutc, const double *erpv, double *rsun, @@ -1629,7 +1674,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 7ae6dc72d..574bdd8ba 100644 --- a/src/rtkpos.c +++ b/src/rtkpos.c @@ -929,55 +929,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; - - C1= SQR(freq1)/(SQR(freq1)-SQR(freq2)); - C2=-SQR(freq2)/(SQR(freq1)-SQR(freq2)); - dant_if=C1*dant[0]+C2*dant[f2]; - - if (obs->L[0]!=0.0&&obs->L[f2]!=0.0) { - y[0]=C1*obs->L[0]*CLIGHT/freq1+C2*obs->L[f2]*CLIGHT/freq2-r-dant_if; - } - if (obs->P[0]!=0.0&&obs->P[f2]!=0.0) { - y[nf]=C1*obs->P[0]+C2*obs->P[f2]-r-dant_if; - } - freq[0]=1.0; - } - else { - for (i=0;isat,obs->code[i],nav))==0.0) continue; - - /* check SNR mask */ - if (testsnr(base,i,azel[1],obs->SNR[i]*SNR_UNIT,&opt->snrmask)) { - continue; - } - /* residuals = observable - estimated range */ - if (obs->L[i]!=0.0) y[i ]=obs->L[i]*CLIGHT/freq[i]-r-dant[i]; - if (obs->P[i]!=0.0) y[i+nf]=obs->P[i] -r-dant[i]; - trace(4,"zdres_sat: %d: L=%.6f P=%.6f r=%.6f f=%.0f\n",obs->sat,obs->L[i], - obs->P[i],r,freq[i]); - } - } +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; + } + + 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; + + /* Check SNR mask */ + if (testsnr(base, i, azel[1], obs->SNR[i] * SNR_UNIT, &opt->snrmask)) { + continue; + } + /* Calc receiver antenna phase center correction */ + double dant = antmodel(opt->pcvr + base, opt->antdel[base], azel, opt->posopt[1], 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) @@ -990,70 +1001,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); - - trace(3,"zdres : n=%d rr=%.2f %.2f %.2f\n",n,rr[0], rr[1], rr[2]); - - /* 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; } /* test valid observation data -----------------------------------------------*/ static int validobs(int i, int j, int f, int nf, double *y)