Changeset 7203 in ntrip for trunk/BNC/src


Ignore:
Timestamp:
Aug 17, 2015, 12:30:54 PM (9 years ago)
Author:
stuerze
Message:

some renaming regarding PPP

Location:
trunk/BNC/src
Files:
15 added
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/BNC/src/PPP/pppClient.cpp

    r6653 r7203  
     1// Part of BNC, a utility for retrieving decoding and
     2// converting GNSS data streams from NTRIP broadcasters.
     3//
     4// Copyright (C) 2007
     5// German Federal Agency for Cartography and Geodesy (BKG)
     6// http://www.bkg.bund.de
     7// Czech Technical University Prague, Department of Geodesy
     8// http://www.fsv.cvut.cz
     9//
     10// Email: euref-ip@bkg.bund.de
     11//
     12// This program is free software; you can redistribute it and/or
     13// modify it under the terms of the GNU General Public License
     14// as published by the Free Software Foundation, version 2.
     15//
     16// This program is distributed in the hope that it will be useful,
     17// but WITHOUT ANY WARRANTY; without even the implied warranty of
     18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     19// GNU General Public License for more details.
     20//
     21// You should have received a copy of the GNU General Public License
     22// along with this program; if not, write to the Free Software
     23// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
     24
    125/* -------------------------------------------------------------------------
    226 * BKG NTRIP Client
     
    529 * Class:      t_pppClient
    630 *
    7  * Purpose:    PPP Client processing starts here
     31 * Purpose:    Precise Point Positioning
    832 *
    933 * Author:     L. Mervart
    1034 *
    11  * Created:    29-Jul-2014
     35 * Created:    21-Nov-2009
    1236 *
    1337 * Changes:   
     
    1539 * -----------------------------------------------------------------------*/
    1640
    17 #include <QThreadStorage>
    18 
    19 #include <iostream>
     41#include <newmatio.h>
    2042#include <iomanip>
    21 #include <stdlib.h>
    22 #include <string.h>
    23 #include <stdexcept>
     43#include <sstream>
    2444
    2545#include "pppClient.h"
    26 #include "pppEphPool.h"
    27 #include "pppObsPool.h"
    28 #include "bncconst.h"
     46#include "bncephuser.h"
    2947#include "bncutils.h"
    30 #include "pppStation.h"
    31 #include "bncantex.h"
    32 #include "pppFilter.h"
    3348
    3449using namespace BNC_PPP;
    3550using namespace std;
    3651
    37 // Global variable holding thread-specific pointers
     52// Constructor
     53////////////////////////////////////////////////////////////////////////////
     54t_pppClient::t_pppClient(const t_pppOptions* opt) {
     55
     56  _opt     = new t_pppOptions(*opt);
     57  _filter  = new t_pppFilter(this);
     58  _epoData = new t_epoData();
     59  _log     = new ostringstream();
     60  _ephUser = new bncEphUser(false);
     61
     62  for (unsigned ii = 0; ii <= t_prn::MAXPRN; ii++) {
     63    _satCodeBiases[ii] = 0;
     64  }
     65
     66}
     67
     68// Destructor
     69////////////////////////////////////////////////////////////////////////////
     70t_pppClient::~t_pppClient() {
     71
     72  for (unsigned ii = 0; ii <= t_prn::MAXPRN; ii++) {
     73    delete _satCodeBiases[ii];
     74  }
     75
     76  delete _filter;
     77  delete _epoData;
     78  delete _opt;
     79  delete _ephUser;
     80  delete _log;
     81}
     82
     83//
     84////////////////////////////////////////////////////////////////////////////
     85void t_pppClient::processEpoch(const vector<t_satObs*>& satObs, t_output* output) {
     86 
     87  // Convert and store observations
     88  // ------------------------------
     89  _epoData->clear();
     90  for (unsigned ii = 0; ii < satObs.size(); ii++) {
     91    const t_satObs* obs     = satObs[ii];
     92    t_prn prn = obs->_prn;
     93    if (prn.system() == 'E') {prn.setFlags(1);} // force I/NAV usage
     94    t_satData*   satData = new t_satData();
     95
     96    if (_epoData->tt.undef()) {
     97      _epoData->tt = obs->_time;
     98    }
     99
     100    satData->tt       = obs->_time;
     101    satData->prn      = QString(prn.toInternalString().c_str());
     102    satData->slipFlag = false;
     103    satData->P1       = 0.0;
     104    satData->P2       = 0.0;
     105    satData->P5       = 0.0;
     106    satData->P7       = 0.0;
     107    satData->L1       = 0.0;
     108    satData->L2       = 0.0;
     109    satData->L5       = 0.0;
     110    satData->L7       = 0.0;
     111    for (unsigned ifrq = 0; ifrq < obs->_obs.size(); ifrq++) {
     112      t_frqObs* frqObs = obs->_obs[ifrq];
     113      double cb = 0.0;
     114      const t_satCodeBias* satCB = satCodeBias(prn);
     115      if (satCB && satCB->_bias.size()) {
     116        for (unsigned ii = 0; ii < satCB->_bias.size(); ii++) {
     117
     118          const t_frqCodeBias& bias = satCB->_bias[ii];
     119          if (frqObs && frqObs->_rnxType2ch == bias._rnxType2ch) {
     120            cb  = bias._value;
     121          }
     122        }
     123      }
     124      if      (frqObs->_rnxType2ch[0] == '1') {
     125        if (frqObs->_codeValid)  satData->P1       = frqObs->_code + cb;
     126        if (frqObs->_phaseValid) satData->L1       = frqObs->_phase;
     127        if (frqObs->_slip)       satData->slipFlag = true;
     128      }
     129      else if (frqObs->_rnxType2ch[0] == '2') {
     130        if (frqObs->_codeValid)  satData->P2       = frqObs->_code + cb;
     131        if (frqObs->_phaseValid) satData->L2       = frqObs->_phase;
     132        if (frqObs->_slip)       satData->slipFlag = true;
     133      }
     134      else if (frqObs->_rnxType2ch[0] == '5') {
     135        if (frqObs->_codeValid)  satData->P5       = frqObs->_code + cb;
     136        if (frqObs->_phaseValid) satData->L5       = frqObs->_phase;
     137        if (frqObs->_slip)       satData->slipFlag = true;
     138      }
     139      else if (frqObs->_rnxType2ch[0] == '7') {
     140        if (frqObs->_codeValid)  satData->P7       = frqObs->_code + cb;
     141        if (frqObs->_phaseValid) satData->L7       = frqObs->_phase;
     142        if (frqObs->_slip)       satData->slipFlag = true;
     143      }
     144    }
     145    putNewObs(satData);
     146  }
     147
     148  // Data Pre-Processing
     149  // -------------------
     150  QMutableMapIterator<QString, t_satData*> it(_epoData->satData);
     151  while (it.hasNext()) {
     152    it.next();
     153    QString    prn     = it.key();
     154    t_satData* satData = it.value();
     155   
     156    if (cmpToT(satData) != success) {
     157      delete satData;
     158      it.remove();
     159      continue;
     160    }
     161
     162  }
     163
     164  // Filter Solution
     165  // ---------------
     166  if (_filter->update(_epoData) == success) {
     167    output->_error = false;
     168    output->_epoTime     = _filter->time();
     169    output->_xyzRover[0] = _filter->x();
     170    output->_xyzRover[1] = _filter->y();
     171    output->_xyzRover[2] = _filter->z();
     172    copy(&_filter->Q().data()[0], &_filter->Q().data()[6], output->_covMatrix);
     173    output->_neu[0]      = _filter->neu()[0];
     174    output->_neu[1]      = _filter->neu()[1];
     175    output->_neu[2]      = _filter->neu()[2];
     176    output->_numSat      = _filter->numSat();
     177    output->_pDop        = _filter->PDOP();
     178    output->_trp0        = _filter->trp0();
     179    output->_trp         = _filter->trp();
     180  }
     181  else {
     182    output->_error = true;
     183  }
     184
     185  output->_log = _log->str();
     186  delete _log; _log = new ostringstream();
     187}
     188
     189//
     190////////////////////////////////////////////////////////////////////////////
     191void t_pppClient::putNewObs(t_satData* satData) {
     192
     193  // Set Observations GPS
     194  // --------------------
     195  if      (satData->system() == 'G') {
     196    if (satData->P1 != 0.0 && satData->P2 != 0.0 &&
     197        satData->L1 != 0.0 && satData->L2 != 0.0 ) {
     198      t_frequency::type fType1 = t_lc::toFreq(satData->system(), t_lc::l1);
     199      t_frequency::type fType2 = t_lc::toFreq(satData->system(), t_lc::l2);
     200      double f1 = t_CST::freq(fType1, 0);
     201      double f2 = t_CST::freq(fType2, 0);
     202      double a1 =   f1 * f1 / (f1 * f1 - f2 * f2);
     203      double a2 = - f2 * f2 / (f1 * f1 - f2 * f2);
     204      satData->L1      = satData->L1 * t_CST::c / f1;
     205      satData->L2      = satData->L2 * t_CST::c / f2;
     206      satData->P3      = a1 * satData->P1 + a2 * satData->P2;
     207      satData->L3      = a1 * satData->L1 + a2 * satData->L2;
     208      satData->lambda3 = a1 * t_CST::c / f1 + a2 * t_CST::c / f2;
     209      satData->lkA     = a1;
     210      satData->lkB     = a2;
     211      _epoData->satData[satData->prn] = satData;
     212    }
     213    else {
     214      delete satData;
     215    }
     216  }
     217
     218  // Set Observations Glonass
     219  // ------------------------
     220  else if (satData->system() == 'R' && _opt->useSystem('R')) {
     221    if (satData->P1 != 0.0 && satData->P2 != 0.0 &&
     222        satData->L1 != 0.0 && satData->L2 != 0.0 ) {
     223
     224      int channel = 0;
     225      if (satData->system() == 'R') {
     226        const t_eph* eph = _ephUser->ephLast(satData->prn);
     227        if (eph) {
     228          channel = eph->slotNum();
     229        }
     230        else {
     231          delete satData;
     232          return;
     233        }
     234      }
     235
     236      t_frequency::type fType1 = t_lc::toFreq(satData->system(), t_lc::l1);
     237      t_frequency::type fType2 = t_lc::toFreq(satData->system(), t_lc::l2);
     238      double f1 = t_CST::freq(fType1, channel);
     239      double f2 = t_CST::freq(fType2, channel);
     240      double a1 =   f1 * f1 / (f1 * f1 - f2 * f2);
     241      double a2 = - f2 * f2 / (f1 * f1 - f2 * f2);
     242      satData->L1      = satData->L1 * t_CST::c / f1;
     243      satData->L2      = satData->L2 * t_CST::c / f2;
     244      satData->P3      = a1 * satData->P1 + a2 * satData->P2;
     245      satData->L3      = a1 * satData->L1 + a2 * satData->L2;
     246      satData->lambda3 = a1 * t_CST::c / f1 + a2 * t_CST::c / f2;
     247      satData->lkA     = a1;
     248      satData->lkB     = a2;
     249      _epoData->satData[satData->prn] = satData;
     250    }
     251    else {
     252      delete satData;
     253    }
     254  }
     255
     256  // Set Observations Galileo
     257  // ------------------------
     258  else if (satData->system() == 'E' && _opt->useSystem('E')) {
     259    if (satData->P1 != 0.0 && satData->P5 != 0.0 &&
     260        satData->L1 != 0.0 && satData->L5 != 0.0 ) {
     261      double f1 = t_CST::freq(t_frequency::E1, 0);
     262      double f5 = t_CST::freq(t_frequency::E5, 0);
     263      double a1 =   f1 * f1 / (f1 * f1 - f5 * f5);
     264      double a5 = - f5 * f5 / (f1 * f1 - f5 * f5);
     265      satData->L1      = satData->L1 * t_CST::c / f1;
     266      satData->L5      = satData->L5 * t_CST::c / f5;
     267      satData->P3      = a1 * satData->P1 + a5 * satData->P5;
     268      satData->L3      = a1 * satData->L1 + a5 * satData->L5;
     269      satData->lambda3 = a1 * t_CST::c / f1 + a5 * t_CST::c / f5;
     270      satData->lkA     = a1;
     271      satData->lkB     = a5;
     272      _epoData->satData[satData->prn] = satData;
     273    }
     274    else {
     275      delete satData;
     276    }
     277  }
     278
     279  // Set Observations BDS
     280  // ---------------------
     281  else if (satData->system() == 'C' && _opt->useSystem('C')) {
     282    if (satData->P2 != 0.0 && satData->P7 != 0.0 &&
     283        satData->L2 != 0.0 && satData->L7 != 0.0 ) {
     284      double f2 = t_CST::freq(t_frequency::C2, 0);
     285      double f7 = t_CST::freq(t_frequency::C7, 0);
     286      double a2 =   f2 * f2 / (f2 * f2 - f7 * f7);
     287      double a7 = - f7 * f7 / (f2 * f2 - f7 * f7);
     288      satData->L2      = satData->L2 * t_CST::c / f2;
     289      satData->L7      = satData->L7 * t_CST::c / f7;
     290      satData->P3      = a2 * satData->P2 + a7 * satData->P7;
     291      satData->L3      = a2 * satData->L2 + a7 * satData->L7;
     292      satData->lambda3 = a2 * t_CST::c / f2 + a7 * t_CST::c / f7;
     293      satData->lkA     = a2;
     294      satData->lkB     = a7;
     295      _epoData->satData[satData->prn] = satData;
     296    }
     297    else {
     298      delete satData;
     299    }
     300  }
     301}
     302
     303//
     304////////////////////////////////////////////////////////////////////////////
     305void t_pppClient::putOrbCorrections(const std::vector<t_orbCorr*>& corr) {
     306  for (unsigned ii = 0; ii < corr.size(); ii++) {
     307    QString prn = QString(corr[ii]->_prn.toInternalString().c_str());
     308    t_eph* eLast = _ephUser->ephLast(prn);
     309    t_eph* ePrev = _ephUser->ephPrev(prn);
     310    if      (eLast && eLast->IOD() == corr[ii]->_iod) {
     311      eLast->setOrbCorr(corr[ii]);
     312    }
     313    else if (ePrev && ePrev->IOD() == corr[ii]->_iod) {
     314      ePrev->setOrbCorr(corr[ii]);
     315    }
     316  }
     317}
     318
     319//
     320////////////////////////////////////////////////////////////////////////////
     321void t_pppClient::putClkCorrections(const std::vector<t_clkCorr*>& corr) {
     322  for (unsigned ii = 0; ii < corr.size(); ii++) {
     323    QString prn = QString(corr[ii]->_prn.toInternalString().c_str());
     324    t_eph* eLast = _ephUser->ephLast(prn);
     325    t_eph* ePrev = _ephUser->ephPrev(prn);
     326    if      (eLast && eLast->IOD() == corr[ii]->_iod) {
     327      eLast->setClkCorr(corr[ii]);
     328    }
     329    else if (ePrev && ePrev->IOD() == corr[ii]->_iod) {
     330      ePrev->setClkCorr(corr[ii]);
     331    }
     332  }
     333}
     334
     335//
    38336//////////////////////////////////////////////////////////////////////////////
    39 QThreadStorage<t_pppClient*> CLIENTS;
    40 
    41 // Static function returning thread-specific pointer
    42 //////////////////////////////////////////////////////////////////////////////
    43 t_pppClient* t_pppClient::instance() {
    44   return CLIENTS.localData();
    45 }
    46 
    47 // Constructor
    48 //////////////////////////////////////////////////////////////////////////////
    49 t_pppClient::t_pppClient(const t_pppOptions* opt) {
    50   _output   = 0;
    51   _opt      = new t_pppOptions(*opt);
    52   _log      = new ostringstream();
    53   _ephPool  = new t_pppEphPool();
    54   _obsPool  = new t_pppObsPool();
    55   _staRover = new t_pppStation();
    56   _filter   = new t_pppFilter();
    57   _tides    = new t_tides();
    58 
    59   if (!_opt->_antexFileName.empty()) {
    60     _antex = new bncAntex(_opt->_antexFileName.c_str());
    61   }
    62   else {
    63     _antex = 0;
    64   }
    65 
    66   CLIENTS.setLocalData(this);  // CLIENTS takes ownership over "this"
    67 }
    68 
    69 // Destructor
    70 //////////////////////////////////////////////////////////////////////////////
    71 t_pppClient::~t_pppClient() {
    72   delete _log;
    73   delete _opt;
    74   delete _ephPool;
    75   delete _obsPool;
    76   delete _staRover;
    77   delete _antex;
    78   delete _filter;
    79   delete _tides;
    80   clearObs();
     337void t_pppClient::putCodeBiases(const std::vector<t_satCodeBias*>& satCodeBias) {
     338  for (unsigned ii = 0; ii < satCodeBias.size(); ii++) {
     339    putCodeBias(new t_satCodeBias(*satCodeBias[ii]));
     340  }
    81341}
    82342
     
    84344//////////////////////////////////////////////////////////////////////////////
    85345void t_pppClient::putEphemeris(const t_eph* eph) {
     346  bool check = _opt->_realTime;
    86347  const t_ephGPS* ephGPS = dynamic_cast<const t_ephGPS*>(eph);
    87348  const t_ephGlo* ephGlo = dynamic_cast<const t_ephGlo*>(eph);
    88349  const t_ephGal* ephGal = dynamic_cast<const t_ephGal*>(eph);
     350  const t_ephBDS* ephBDS = dynamic_cast<const t_ephBDS*>(eph);
    89351  if      (ephGPS) {
    90     _ephPool->putEphemeris(new t_ephGPS(*ephGPS));
     352    _ephUser->putNewEph(new t_ephGPS(*ephGPS), check);
    91353  }
    92354  else if (ephGlo) {
    93     _ephPool->putEphemeris(new t_ephGlo(*ephGlo));
     355    _ephUser->putNewEph(new t_ephGlo(*ephGlo), check);
    94356  }
    95357  else if (ephGal) {
    96     _ephPool->putEphemeris(new t_ephGal(*ephGal));
    97   }
    98 }
    99 
    100 //
    101 //////////////////////////////////////////////////////////////////////////////
    102 void t_pppClient::putOrbCorrections(const vector<t_orbCorr*>& corr) {
    103   for (unsigned ii = 0; ii < corr.size(); ii++) {
    104     _ephPool->putOrbCorrection(new t_orbCorr(*corr[ii]));
    105   }
    106 }
    107 
    108 //
    109 //////////////////////////////////////////////////////////////////////////////
    110 void t_pppClient::putClkCorrections(const vector<t_clkCorr*>& corr) {
    111   for (unsigned ii = 0; ii < corr.size(); ii++) {
    112     _ephPool->putClkCorrection(new t_clkCorr(*corr[ii]));
    113   }
    114 }
    115 
    116 //
    117 //////////////////////////////////////////////////////////////////////////////
    118 void t_pppClient::putCodeBiases(const vector<t_satCodeBias*>& biases) {
    119   for (unsigned ii = 0; ii < biases.size(); ii++) {
    120     _obsPool->putCodeBias(new t_satCodeBias(*biases[ii]));
    121   }
    122 }
    123 
    124 //
    125 //////////////////////////////////////////////////////////////////////////////
    126 t_irc t_pppClient::prepareObs(const vector<t_satObs*>& satObs,
    127                               vector<t_pppSatObs*>& obsVector, bncTime& epoTime) {
    128   // Default
    129   // -------
    130   epoTime.reset();
    131 
    132   // Create vector of valid observations
    133   // -----------------------------------
    134   for (unsigned ii = 0; ii < satObs.size(); ii++) {
    135     char system = satObs[ii]->_prn.system();
    136     if (OPT->useSystem(system)) {
    137       t_pppSatObs* pppSatObs = new t_pppSatObs(*satObs[ii]);
    138       if (pppSatObs->isValid()) {
    139         obsVector.push_back(pppSatObs);
    140       }
    141       else {
    142         delete pppSatObs;
    143       }
    144     }
    145   }
    146 
    147   // Check whether data are synchronized, compute epoTime
    148   // ----------------------------------------------------
    149   const double MAXSYNC = 0.05; // synchronization limit
    150   double meanDt = 0.0;
    151   for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    152     const t_pppSatObs* satObs = obsVector.at(ii);
    153     if (epoTime.undef()) {
    154       epoTime = satObs->time();
    155     }
    156     else {
    157       double dt = satObs->time() - epoTime;
    158       if (fabs(dt) > MAXSYNC) {
    159         LOG << "t_pppClient::prepareObs asynchronous observations" << endl;
    160         return failure;
    161       }
    162       meanDt += dt;
    163     }
    164   }
    165 
    166   if (obsVector.size() > 0) {
    167     epoTime += meanDt / obsVector.size();
    168   }
    169 
    170   return success;
    171 }
    172 
    173 // Compute the Bancroft position, check for blunders
    174 //////////////////////////////////////////////////////////////////////////////
    175 t_irc t_pppClient::cmpBancroft(const bncTime& epoTime,
    176                                   vector<t_pppSatObs*>& obsVector,
    177                                   ColumnVector& xyzc, bool print) {
    178 
    179   t_lc::type tLC = t_lc::dummy;
    180 
    181   while (true) {
    182     Matrix BB(obsVector.size(), 4);
    183     int iObs = -1;
    184     for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    185       const t_pppSatObs* satObs = obsVector.at(ii);
    186       if (satObs->prn().system() == 'G') {
    187         if (tLC == t_lc::dummy) {
    188           tLC = satObs->isValid(t_lc::cIF) ? t_lc::cIF : t_lc::c1;
    189         }
    190         if ( satObs->isValid(tLC) && (!satObs->modelSet() || satObs->eleSat() >= OPT->_minEle) ) {
    191           ++iObs;   
    192           BB[iObs][0] = satObs->xc()[0];
    193           BB[iObs][1] = satObs->xc()[1];
    194           BB[iObs][2] = satObs->xc()[2];
    195           BB[iObs][3] = satObs->obsValue(tLC) - satObs->cmpValueForBanc(tLC);
    196         }
    197       }
    198     }
    199     if (iObs + 1 < OPT->_minObs) {
    200       LOG << "t_pppClient::cmpBancroft not enough observations" << endl;
     358    _ephUser->putNewEph(new t_ephGal(*ephGal), check);
     359  }
     360  else if (ephBDS) {
     361      _ephUser->putNewEph(new t_ephBDS(*ephBDS), check);
     362    }
     363}
     364
     365// Satellite Position
     366////////////////////////////////////////////////////////////////////////////
     367t_irc t_pppClient::getSatPos(const bncTime& tt, const QString& prn,
     368                              ColumnVector& xc, ColumnVector& vv) {
     369
     370  t_eph* eLast = _ephUser->ephLast(prn);
     371  t_eph* ePrev = _ephUser->ephPrev(prn);
     372  if      (eLast && eLast->getCrd(tt, xc, vv, _opt->useOrbClkCorr()) == success) {
     373    return success;
     374  }
     375  else if (ePrev && ePrev->getCrd(tt, xc, vv, _opt->useOrbClkCorr()) == success) {
     376    return success;
     377  }
     378  return failure;
     379}
     380
     381// Correct Time of Transmission
     382////////////////////////////////////////////////////////////////////////////
     383t_irc t_pppClient::cmpToT(t_satData* satData) {
     384
     385  double prange = satData->P3;
     386  if (prange == 0.0) {
     387    return failure;
     388  }
     389
     390  double clkSat = 0.0;
     391  for (int ii = 1; ii <= 10; ii++) {
     392
     393    bncTime ToT = satData->tt - prange / t_CST::c - clkSat;
     394
     395    ColumnVector xc(4);
     396    ColumnVector vv(3);
     397    if (getSatPos(ToT, satData->prn, xc, vv) != success) {
    201398      return failure;
    202399    }
    203     BB = BB.Rows(1,iObs+1);
    204     bancroft(BB, xyzc);
    205 
    206     xyzc[3] /= t_CST::c;
    207 
    208     // Check Blunders
    209     // --------------
    210     const double BLUNDER = 100.0;
    211     double   maxRes      = 0.0;
    212     unsigned maxResIndex = 0;
    213     for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    214       const t_pppSatObs* satObs = obsVector.at(ii);
    215       if ( satObs->isValid() && satObs->prn().system() == 'G' &&
    216            (!satObs->modelSet() || satObs->eleSat() >= OPT->_minEle) ) {
    217         ColumnVector rr = satObs->xc().Rows(1,3) - xyzc.Rows(1,3);
    218         double res = rr.norm_Frobenius() - satObs->obsValue(tLC)
    219           - (satObs->xc()[3] - xyzc[3]) * t_CST::c;
    220         if (fabs(res) > maxRes) {
    221           maxRes      = fabs(res);
    222           maxResIndex = ii;
    223         }
    224       }
    225     }
    226     if (maxRes < BLUNDER) {
    227       if (print) {
    228         LOG.setf(ios::fixed);
    229         LOG << string(epoTime) << " BANCROFT:"        << ' '
    230             << setw(14) << setprecision(3) << xyzc[0] << ' '
    231             << setw(14) << setprecision(3) << xyzc[1] << ' '
    232             << setw(14) << setprecision(3) << xyzc[2] << ' '
    233             << setw(14) << setprecision(3) << xyzc[3] * t_CST::c << endl << endl;
    234       }
    235       break;
    236     }
    237     else {
    238       t_pppSatObs* satObs = obsVector.at(maxResIndex);
    239       LOG << "t_pppClient::cmpBancroft outlier " << satObs->prn().toString()
    240           << " " << maxRes << endl;
    241       delete satObs;
    242       obsVector.erase(obsVector.begin() + maxResIndex);
    243     }
    244   }
    245 
    246   return success;
    247 }
    248 
    249 // Compute A Priori GPS-Glonass Offset
    250 //////////////////////////////////////////////////////////////////////////////
    251 double t_pppClient::cmpOffGG(vector<t_pppSatObs*>& obsVector) {
    252 
    253   t_lc::type tLC   = t_lc::dummy;
    254   double     offGG = 0.0;
    255 
    256   if (OPT->useSystem('R')) {
    257 
    258     while (obsVector.size() > 0) {
    259       offGG = 0.0;
    260       double   maxRes      = 0.0;
    261       int      maxResIndex = -1;
    262       t_prn    maxResPrn;
    263       unsigned nObs        = 0;
    264       for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    265         t_pppSatObs* satObs = obsVector.at(ii);
    266         if (satObs->prn().system() == 'R') {
    267           if (tLC == t_lc::dummy) {
    268             tLC = satObs->isValid(t_lc::cIF) ? t_lc::cIF : t_lc::c1;
    269           }
    270           if (satObs->isValid(tLC) && (!satObs->modelSet() || satObs->eleSat() >= OPT->_minEle)) {
    271             double ll = satObs->obsValue(tLC) - satObs->cmpValue(tLC);
    272             ++nObs;
    273             offGG += ll;
    274             if (fabs(ll) > fabs(maxRes)) {
    275               maxRes      = ll;
    276               maxResIndex = ii;
    277               maxResPrn   = satObs->prn();
    278             }
    279           }
    280         }
    281       }
    282 
    283       if (nObs > 0) {
    284         offGG = offGG / nObs;
    285       }
    286       else {
    287         offGG = 0.0;
    288       }
    289 
    290       if (fabs(maxRes) > 1000.0) {
    291         LOG << "t_pppClient::cmpOffGG outlier " << maxResPrn.toString() << " " << maxRes << endl;
    292         obsVector.erase(obsVector.begin() + maxResIndex);
    293       }
    294       else {
    295         break;
    296       }
    297     }
    298   }
    299 
    300   return offGG;
    301 }
    302 
    303 //
    304 //////////////////////////////////////////////////////////////////////////////
    305 void t_pppClient::initOutput(t_output* output) {
    306   _output = output;
    307   _output->_numSat = 0;
    308   _output->_pDop   = 0.0;
    309   _output->_error  = false;
    310 }
    311 
    312 //
    313 //////////////////////////////////////////////////////////////////////////////
    314 void t_pppClient::clearObs() {
    315   for (unsigned ii = 0; ii < _obsRover.size(); ii++) {
    316     delete _obsRover.at(ii);
    317   }
    318   _obsRover.clear();
    319 }
    320 
    321 //
    322 //////////////////////////////////////////////////////////////////////////////
    323 void t_pppClient::finish(t_irc irc) {
    324 
    325   clearObs();
    326 
    327   _output->_epoTime = _epoTimeRover;
    328 
    329   if (irc == success) {
    330     _output->_xyzRover[0] = _staRover->xyzApr()[0] + _filter->x()[0];
    331     _output->_xyzRover[1] = _staRover->xyzApr()[1] + _filter->x()[1];
    332     _output->_xyzRover[2] = _staRover->xyzApr()[2] + _filter->x()[2];
    333 
    334     xyz2neu(_staRover->ellApr().data(), _filter->x().data(), _output->_neu);
    335 
    336     copy(&_filter->Q().data()[0], &_filter->Q().data()[6], _output->_covMatrix);
    337 
    338     _output->_trp0     = t_tropo::delay_saast(_staRover->xyzApr(), M_PI/2.0);
    339     _output->_trp      = _filter->trp();
    340     _output->_trpStdev = _filter->trpStdev();
    341 
    342     _output->_numSat     = _filter->numSat();
    343     _output->_pDop       = _filter->PDOP();
    344     _output->_error = false;
    345   }
    346   else {
    347     _output->_error = true;
    348   }
    349   _output->_log = _log->str();
    350   delete _log; _log = new ostringstream();
    351 }
    352 
    353 //
    354 //////////////////////////////////////////////////////////////////////////////
    355 t_irc t_pppClient::cmpModel(t_pppStation* station, const ColumnVector& xyzc,
    356                                vector<t_pppSatObs*>& obsVector) {
    357 
    358   bncTime time;
    359   time = _epoTimeRover;
    360   station->setName(OPT->_roverName);
    361   station->setAntName(OPT->_antNameRover);
    362   if (OPT->xyzAprRoverSet()) {
    363     station->setXyzApr(OPT->_xyzAprRover);
    364   }
    365   else {
    366     station->setXyzApr(xyzc.Rows(1,3));
    367   }
    368   station->setNeuEcc(OPT->_neuEccRover);
    369 
    370   // Receiver Clock
    371   // --------------
    372   station->setDClk(xyzc[3]);
    373 
    374   // Tides
    375   // -----
    376   station->setTideDspl( _tides->displacement(time, station->xyzApr()) );
    377  
    378   // Observation model
    379   // -----------------
    380   vector<t_pppSatObs*>::iterator it = obsVector.begin();
    381   while (it != obsVector.end()) {
    382     t_pppSatObs* satObs = *it;
    383     if (satObs->isValid()) {
    384       satObs->cmpModel(station);
    385     }
    386     if (satObs->isValid() && satObs->eleSat() >= OPT->_minEle) {
    387       ++it;
    388     }
    389     else {
    390       delete satObs;
    391       it = obsVector.erase(it);
    392     }
    393   }
    394 
    395   return success;
    396 }
    397 
    398 //
    399 //////////////////////////////////////////////////////////////////////////////
    400 void t_pppClient::processEpoch(const vector<t_satObs*>& satObs, t_output* output) {
    401 
    402   try {
    403     initOutput(output);
    404 
    405     // Prepare Observations of the Rover
    406     // ---------------------------------   
    407     if (prepareObs(satObs, _obsRover, _epoTimeRover) != success) {
    408       return finish(failure);
    409     }
    410 
    411     LOG << "\nResults of Epoch ";
    412     if (!_epoTimeRover.undef()) LOG << string(_epoTimeRover);
    413     LOG << "\n--------------------------------------\n";
    414  
    415     for (int iter = 1; iter <= 2; iter++) {
    416       ColumnVector xyzc(4); xyzc = 0.0;
    417       bool print = (iter == 2);
    418       if (cmpBancroft(_epoTimeRover, _obsRover, xyzc, print) != success) {
    419         return finish(failure);
    420       }
    421       if (cmpModel(_staRover, xyzc, _obsRover) != success) {
    422         return finish(failure);
    423       }
    424     }
    425 
    426     _offGG = cmpOffGG(_obsRover);
    427 
    428     // Store last epoch of data
    429     // ------------------------   
    430     _obsPool->putEpoch(_epoTimeRover, _obsRover);
    431 
    432     // Process Epoch in Filter
    433     // -----------------------
    434     if (_filter->processEpoch(_obsPool) != success) {
    435       return finish(failure);
    436     }
    437   }
    438   catch (Exception& exc) {
    439     LOG << exc.what() << endl;
    440     return finish(failure);
    441   }
    442   catch (t_except msg) {
    443     LOG << msg.what() << endl;
    444     return finish(failure);
    445   }
    446   catch (const char* msg) {
    447     LOG << msg << endl;
    448     return finish(failure);
    449   }
    450   catch (...) {
    451     LOG << "unknown exception" << endl;
    452     return finish(failure);
    453   }
    454 
    455   return finish(success);
    456 }
    457 
    458 //
    459 ////////////////////////////////////////////////////////////////////////////
    460 double lorentz(const ColumnVector& aa, const ColumnVector& bb) {
    461   return aa[0]*bb[0] +  aa[1]*bb[1] +  aa[2]*bb[2] -  aa[3]*bb[3];
    462 }
    463 
    464 //
    465 ////////////////////////////////////////////////////////////////////////////
    466 void t_pppClient::bancroft(const Matrix& BBpass, ColumnVector& pos) {
    467 
    468   if (pos.Nrows() != 4) {
    469     pos.ReSize(4);
    470   }
    471   pos = 0.0;
    472 
    473   for (int iter = 1; iter <= 2; iter++) {
    474     Matrix BB = BBpass;
    475     int mm = BB.Nrows();
    476     for (int ii = 1; ii <= mm; ii++) {
    477       double xx = BB(ii,1);
    478       double yy = BB(ii,2);
    479       double traveltime = 0.072;
    480       if (iter > 1) {
    481         double zz  = BB(ii,3);
    482         double rho = sqrt( (xx-pos(1)) * (xx-pos(1)) +
    483                            (yy-pos(2)) * (yy-pos(2)) +
    484                            (zz-pos(3)) * (zz-pos(3)) );
    485         traveltime = rho / t_CST::c;
    486       }
    487       double angle = traveltime * t_CST::omega;
    488       double cosa  = cos(angle);
    489       double sina  = sin(angle);
    490       BB(ii,1) =  cosa * xx + sina * yy;
    491       BB(ii,2) = -sina * xx + cosa * yy;
    492     }
    493    
    494     Matrix BBB;
    495     if (mm > 4) {
    496       SymmetricMatrix hlp; hlp << BB.t() * BB;
    497       BBB = hlp.i() * BB.t();
    498     }
    499     else {
    500       BBB = BB.i();
    501     }
    502     ColumnVector ee(mm); ee = 1.0;
    503     ColumnVector alpha(mm); alpha = 0.0;
    504     for (int ii = 1; ii <= mm; ii++) {
    505       alpha(ii) = lorentz(BB.Row(ii).t(),BB.Row(ii).t())/2.0;
    506     }
    507     ColumnVector BBBe     = BBB * ee;
    508     ColumnVector BBBalpha = BBB * alpha;
    509     double aa = lorentz(BBBe, BBBe);
    510     double bb = lorentz(BBBe, BBBalpha)-1;
    511     double cc = lorentz(BBBalpha, BBBalpha);
    512     double root = sqrt(bb*bb-aa*cc);
    513 
    514     Matrix hlpPos(4,2);
    515     hlpPos.Column(1) = (-bb-root)/aa * BBBe + BBBalpha;
    516     hlpPos.Column(2) = (-bb+root)/aa * BBBe + BBBalpha;
    517 
    518     ColumnVector omc(2);
    519     for (int pp = 1; pp <= 2; pp++) {
    520       hlpPos(4,pp)      = -hlpPos(4,pp);
    521       omc(pp) = BB(1,4) -
    522                 sqrt( (BB(1,1)-hlpPos(1,pp)) * (BB(1,1)-hlpPos(1,pp)) +
    523                       (BB(1,2)-hlpPos(2,pp)) * (BB(1,2)-hlpPos(2,pp)) +
    524                       (BB(1,3)-hlpPos(3,pp)) * (BB(1,3)-hlpPos(3,pp)) ) -
    525                 hlpPos(4,pp);
    526     }
    527     if ( fabs(omc(1)) > fabs(omc(2)) ) {
    528       pos = hlpPos.Column(2);
    529     }
    530     else {
    531       pos = hlpPos.Column(1);
    532     }
    533   }
    534 }
    535 
     400
     401    double clkSatOld = clkSat;
     402    clkSat = xc(4);
     403
     404    if ( fabs(clkSat-clkSatOld) * t_CST::c < 1.e-4 ) {
     405      satData->xx      = xc.Rows(1,3);
     406      satData->vv      = vv;
     407      satData->clk     = clkSat * t_CST::c;
     408      return success;
     409    }
     410  }
     411
     412  return failure;
     413}
     414
     415void t_pppClient::putCodeBias(t_satCodeBias* satCodeBias) {
     416  int iPrn = satCodeBias->_prn.toInt();
     417  delete _satCodeBiases[iPrn];
     418  _satCodeBiases[iPrn] = satCodeBias;
     419}
  • trunk/BNC/src/PPP/pppClient.h

    r6463 r7203  
     1// Part of BNC, a utility for retrieving decoding and
     2// converting GNSS data streams from NTRIP broadcasters.
     3//
     4// Copyright (C) 2007
     5// German Federal Agency for Cartography and Geodesy (BKG)
     6// http://www.bkg.bund.de
     7// Czech Technical University Prague, Department of Geodesy
     8// http://www.fsv.cvut.cz
     9//
     10// Email: euref-ip@bkg.bund.de
     11//
     12// This program is free software; you can redistribute it and/or
     13// modify it under the terms of the GNU General Public License
     14// as published by the Free Software Foundation, version 2.
     15//
     16// This program is distributed in the hope that it will be useful,
     17// but WITHOUT ANY WARRANTY; without even the implied warranty of
     18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     19// GNU General Public License for more details.
     20//
     21// You should have received a copy of the GNU General Public License
     22// along with this program; if not, write to the Free Software
     23// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
     24
    125#ifndef PPPCLIENT_H
    226#define PPPCLIENT_H
    327
    4 #include <sstream>
    528#include <vector>
     29#include <QtCore>
     30
    631#include "pppInclude.h"
    7 #include "ephemeris.h"
    832#include "pppOptions.h"
    9 #include "pppModel.h"
    10 #include "satObs.h"
     33#include "pppFilter.h"
    1134
    12 class bncAntex;
     35class bncEphUser;
     36class t_eph;
    1337
    1438namespace BNC_PPP {
    15 
    16 class t_pppEphPool;
    17 class t_pppObsPool;
    18 class t_pppSatObs;
    19 class t_pppStation;
    20 class t_pppFilter;
    21 
     39 
    2240class t_pppClient : public interface_pppClient {
    2341 public:
    24   t_pppClient(const t_pppOptions* opt);                                                     
    25   ~t_pppClient();                                                     
     42  t_pppClient(const t_pppOptions* opt);
     43  ~t_pppClient();
     44  void                processEpoch(const std::vector<t_satObs*>& satObs, t_output* output);
     45  void                putEphemeris(const t_eph* eph);                 
     46  void                putOrbCorrections(const std::vector<t_orbCorr*>& corr);
     47  void                putClkCorrections(const std::vector<t_clkCorr*>& corr);
     48  void                putCodeBiases(const std::vector<t_satCodeBias*>& satCodeBias);   
     49  std::ostringstream& log() {return *_log;}
     50  const t_pppOptions* opt() const {return _opt;}
     51  void putCodeBias(t_satCodeBias* satCodeBias);
     52  const t_satCodeBias* satCodeBias(const t_prn& prn) const {
     53    return _satCodeBiases[prn.toInt()];
     54  }
     55 private:
     56  t_irc getSatPos(const bncTime& tt, const QString& prn, ColumnVector& xc, ColumnVector& vv);
     57  void  putNewObs(t_satData* satData);
     58  t_irc cmpToT(t_satData* satData);
    2659
    27   void putEphemeris(const t_eph* eph);                 
    28   void putOrbCorrections(const std::vector<t_orbCorr*>& corr);
    29   void putClkCorrections(const std::vector<t_clkCorr*>& corr);
    30   void putCodeBiases(const std::vector<t_satCodeBias*>& satBias);   
    31   void processEpoch(const std::vector<t_satObs*>& satObs, t_output* output);
    32 
    33   const t_pppEphPool* ephPool() const {return _ephPool;}
    34   const t_pppObsPool* obsPool() const {return _obsPool;}
    35   const bncAntex*  antex() const {return _antex;}
    36   const t_pppStation* staRover() const {return _staRover;}
    37   double           offGG() const {return _offGG;}
    38 
    39   std::ostringstream& log() {return *_log;}
    40   const t_pppOptions*    opt() const {return _opt;}
    41 
    42   static void bancroft(const Matrix& BBpass, ColumnVector& pos);
    43 
    44   static t_pppClient* instance();
    45 
    46  private:
    47   void initOutput(t_output* output);
    48   void finish(t_irc irc);
    49   void clearObs();
    50   t_irc prepareObs(const std::vector<t_satObs*>& satObs,
    51                    std::vector<t_pppSatObs*>& obsVector, bncTime& epoTime);
    52   t_irc cmpModel(t_pppStation* station, const ColumnVector& xyzc,
    53                  std::vector<t_pppSatObs*>& obsVector);
    54   t_irc cmpBancroft(const bncTime& epoTime, std::vector<t_pppSatObs*>& obsVector,
    55                     ColumnVector& xyzc, bool print);
    56   double cmpOffGG(std::vector<t_pppSatObs*>& obsVector);
    57 
    58   t_output*                 _output;
    59   t_pppEphPool*             _ephPool;
    60   t_pppObsPool*             _obsPool;
    61   bncTime                   _epoTimeRover;
    62   t_pppStation*             _staRover;
    63   bncAntex*                 _antex;
    64   t_pppFilter*              _filter;
    65   double                    _offGG;
    66   std::vector<t_pppSatObs*> _obsRover;
    67   std::ostringstream*       _log;
    68   t_pppOptions*             _opt;
    69   t_tides*                  _tides;
     60  t_satCodeBias*       _satCodeBiases[t_prn::MAXPRN+1];
     61  bncEphUser*         _ephUser;
     62  t_pppOptions*       _opt;
     63  t_epoData*          _epoData;
     64  t_pppFilter*        _filter;
     65  std::ostringstream* _log;
    7066};
    7167
    72 }; // namespace BNC_PPP
    73 
    74 #define PPP_CLIENT (BNC_PPP::t_pppClient::instance())
    75 #define LOG        (BNC_PPP::t_pppClient::instance()->log())
    76 #define OPT        (BNC_PPP::t_pppClient::instance()->opt())
     68} // namespace
    7769
    7870#endif
  • trunk/BNC/src/PPP/pppFilter.cpp

    r6997 r7203  
     1// Part of BNC, a utility for retrieving decoding and
     2// converting GNSS data streams from NTRIP broadcasters.
     3//
     4// Copyright (C) 2007
     5// German Federal Agency for Cartography and Geodesy (BKG)
     6// http://www.bkg.bund.de
     7// Czech Technical University Prague, Department of Geodesy
     8// http://www.fsv.cvut.cz
     9//
     10// Email: euref-ip@bkg.bund.de
     11//
     12// This program is free software; you can redistribute it and/or
     13// modify it under the terms of the GNU General Public License
     14// as published by the Free Software Foundation, version 2.
     15//
     16// This program is distributed in the hope that it will be useful,
     17// but WITHOUT ANY WARRANTY; without even the implied warranty of
     18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     19// GNU General Public License for more details.
     20//
     21// You should have received a copy of the GNU General Public License
     22// along with this program; if not, write to the Free Software
     23// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
     24
    125/* -------------------------------------------------------------------------
    226 * BKG NTRIP Client
    327 * -------------------------------------------------------------------------
    428 *
    5  * Class:      t_pppFilter
     29 * Class:      t_pppParam, t_pppFilter
    630 *
    7  * Purpose:    Filter Adjustment
     31 * Purpose:    Model for PPP
    832 *
    933 * Author:     L. Mervart
    1034 *
    11  * Created:    29-Jul-2014
     35 * Created:    01-Dec-2009
    1236 *
    1337 * Changes:   
     
    1539 * -----------------------------------------------------------------------*/
    1640
    17 #include <iostream>
    1841#include <iomanip>
    1942#include <cmath>
    20 #include <newmat.h>
     43#include <sstream>
    2144#include <newmatio.h>
    2245#include <newmatap.h>
    2346
    2447#include "pppFilter.h"
     48#include "pppClient.h"
    2549#include "bncutils.h"
    26 #include "pppParlist.h"
    27 #include "pppObsPool.h"
    28 #include "pppStation.h"
    29 #include "pppClient.h"
     50#include "bncantex.h"
     51#include "pppOptions.h"
     52#include "pppModel.h"
    3053
    3154using namespace BNC_PPP;
    3255using namespace std;
    3356
     57const double   MAXRES_CODE           = 2.98 * 3.0;
     58const double   MAXRES_PHASE_GPS      = 2.98 * 0.03;
     59const double   MAXRES_PHASE_GLONASS  = 2.98 * 0.03;
     60const double   GLONASS_WEIGHT_FACTOR = 1.0;
     61const double   BDS_WEIGHT_FACTOR     = 2.0;
     62
     63#define LOG (_pppClient->log())
     64#define OPT (_pppClient->opt())
     65
    3466// Constructor
    3567////////////////////////////////////////////////////////////////////////////
    36 t_pppFilter::t_pppFilter() {
    37   _parlist = 0;
     68t_pppParam::t_pppParam(t_pppParam::parType typeIn, int indexIn,
     69                   const QString& prnIn) {
     70  type      = typeIn;
     71  index     = indexIn;
     72  prn       = prnIn;
     73  index_old = 0;
     74  xx        = 0.0;
     75  numEpo    = 0;
    3876}
    3977
    4078// Destructor
    4179////////////////////////////////////////////////////////////////////////////
     80t_pppParam::~t_pppParam() {
     81}
     82
     83// Partial
     84////////////////////////////////////////////////////////////////////////////
     85double t_pppParam::partial(t_satData* satData, bool phase) {
     86
     87  Tracer tracer("t_pppParam::partial");
     88
     89  // Coordinates
     90  // -----------
     91  if      (type == CRD_X) {
     92    return (xx - satData->xx(1)) / satData->rho;
     93  }
     94  else if (type == CRD_Y) {
     95    return (xx - satData->xx(2)) / satData->rho;
     96  }
     97  else if (type == CRD_Z) {
     98    return (xx - satData->xx(3)) / satData->rho;
     99  }
     100
     101  // Receiver Clocks
     102  // ---------------
     103  else if (type == RECCLK) {
     104    return 1.0;
     105  }
     106
     107  // Troposphere
     108  // -----------
     109  else if (type == TROPO) {
     110    return 1.0 / sin(satData->eleSat);
     111  }
     112
     113  // Glonass Offset
     114  // --------------
     115  else if (type == GLONASS_OFFSET) {
     116    if (satData->prn[0] == 'R') {
     117      return 1.0;
     118    }
     119    else {
     120      return 0.0;
     121    }
     122  }
     123
     124  // Galileo Offset
     125  // --------------
     126  else if (type == GALILEO_OFFSET) {
     127    if (satData->prn[0] == 'E') {
     128      return 1.0;
     129    }
     130    else {
     131      return 0.0;
     132    }
     133  }
     134
     135  // BDS Offset
     136  // ----------
     137  else if (type == BDS_OFFSET) {
     138    if (satData->prn[0] == 'C') {
     139      return 1.0;
     140    }
     141    else {
     142      return 0.0;
     143    }
     144  }
     145
     146  // Ambiguities
     147  // -----------
     148  else if (type == AMB_L3) {
     149    if (phase && satData->prn == prn) {
     150      return 1.0;
     151    }
     152    else {
     153      return 0.0;
     154    }
     155  }
     156
     157  // Default return
     158  // --------------
     159  return 0.0;
     160}
     161
     162// Constructor
     163////////////////////////////////////////////////////////////////////////////
     164t_pppFilter::t_pppFilter(t_pppClient* pppClient) {
     165
     166  _pppClient = pppClient;
     167  _tides     = new t_tides();
     168
     169  // Antenna Name, ANTEX File
     170  // ------------------------
     171  _antex = 0;
     172  if (!OPT->_antexFileName.empty()) {
     173    _antex = new bncAntex(OPT->_antexFileName.c_str());
     174  }
     175
     176  // Bancroft Coordinates
     177  // --------------------
     178  _xcBanc.ReSize(4);  _xcBanc  = 0.0;
     179  _ellBanc.ReSize(3); _ellBanc = 0.0;
     180
     181  // Save copy of data (used in outlier detection)
     182  // ---------------------------------------------
     183  _epoData_sav = new t_epoData();
     184
     185  // Some statistics
     186  // ---------------
     187  _neu.ReSize(3); _neu = 0.0;
     188  _numSat = 0;
     189  _pDop   = 0.0;
     190}
     191
     192// Destructor
     193////////////////////////////////////////////////////////////////////////////
    42194t_pppFilter::~t_pppFilter() {
    43   delete _parlist;
    44 }
    45 
    46 // Process Single Epoch
    47 ////////////////////////////////////////////////////////////////////////////
    48 t_irc t_pppFilter::processEpoch(t_pppObsPool* obsPool) {
    49 
    50   _numSat     = 0;
    51 
    52   if (!_parlist) {
    53     _parlist = new t_pppParlist();
    54   }
    55 
    56   // Vector of all Observations
    57   // --------------------------
    58   t_pppObsPool::t_epoch* epoch = obsPool->lastEpoch();
    59   if (!epoch) {
     195  delete _tides;
     196  delete _antex;
     197  for (int iPar = 1; iPar <= _params.size(); iPar++) {
     198    delete _params[iPar-1];
     199  }
     200  for (int iPar = 1; iPar <= _params_sav.size(); iPar++) {
     201    delete _params_sav[iPar-1];
     202  }
     203  delete _epoData_sav;
     204}
     205
     206// Reset Parameters and Variance-Covariance Matrix
     207////////////////////////////////////////////////////////////////////////////
     208void t_pppFilter::reset() {
     209
     210  Tracer tracer("t_pppFilter::reset");
     211
     212  double lastTrp = 0.0;
     213  for (int ii = 0; ii < _params.size(); ii++) {
     214    t_pppParam* pp = _params[ii];
     215    if (pp->type == t_pppParam::TROPO) {
     216      lastTrp = pp->xx;
     217    }
     218    delete pp;
     219  }
     220  _params.clear();
     221
     222  int nextPar = 0;
     223  _params.push_back(new t_pppParam(t_pppParam::CRD_X,  ++nextPar, ""));
     224  _params.push_back(new t_pppParam(t_pppParam::CRD_Y,  ++nextPar, ""));
     225  _params.push_back(new t_pppParam(t_pppParam::CRD_Z,  ++nextPar, ""));
     226  _params.push_back(new t_pppParam(t_pppParam::RECCLK, ++nextPar, ""));
     227  if (OPT->estTrp()) {
     228    _params.push_back(new t_pppParam(t_pppParam::TROPO, ++nextPar, ""));
     229  }
     230  if (OPT->useSystem('R')) {
     231    _params.push_back(new t_pppParam(t_pppParam::GLONASS_OFFSET, ++nextPar, ""));
     232  }
     233  if (OPT->useSystem('E')) {
     234    _params.push_back(new t_pppParam(t_pppParam::GALILEO_OFFSET, ++nextPar, ""));
     235  }
     236  if (OPT->useSystem('C')) {
     237    _params.push_back(new t_pppParam(t_pppParam::BDS_OFFSET, ++nextPar, ""));
     238  }
     239
     240  _QQ.ReSize(_params.size());
     241  _QQ = 0.0;
     242  for (int iPar = 1; iPar <= _params.size(); iPar++) {
     243    t_pppParam* pp = _params[iPar-1];
     244    pp->xx = 0.0;
     245    if      (pp->isCrd()) {
     246      _QQ(iPar,iPar) = OPT->_aprSigCrd(1) * OPT->_aprSigCrd(1);
     247    }
     248    else if (pp->type == t_pppParam::RECCLK) {
     249      _QQ(iPar,iPar) = OPT->_noiseClk * OPT->_noiseClk;
     250    }
     251    else if (pp->type == t_pppParam::TROPO) {
     252      _QQ(iPar,iPar) = OPT->_aprSigTrp * OPT->_aprSigTrp;
     253      pp->xx = lastTrp;
     254    }
     255    else if (pp->type == t_pppParam::GLONASS_OFFSET) {
     256      _QQ(iPar,iPar) = 1000.0 * 1000.0;
     257    }
     258    else if (pp->type == t_pppParam::GALILEO_OFFSET) {
     259      _QQ(iPar,iPar) = 1000.0 * 1000.0;
     260    }
     261    else if (pp->type == t_pppParam::BDS_OFFSET) {
     262      _QQ(iPar,iPar) = 1000.0 * 1000.0;
     263    }
     264  }
     265}
     266
     267// Bancroft Solution
     268////////////////////////////////////////////////////////////////////////////
     269t_irc t_pppFilter::cmpBancroft(t_epoData* epoData) {
     270
     271  Tracer tracer("t_pppFilter::cmpBancroft");
     272
     273  if (int(epoData->sizeSys('G')) < OPT->_minObs) {
     274    LOG << "t_pppFilter::cmpBancroft: not enough data\n";
    60275    return failure;
    61276  }
    62   vector<t_pppSatObs*>& allObs = epoch->obsVector();
    63 
    64   // Time of the Epoch
    65   // -----------------
    66   _epoTime = epoch->epoTime();
    67 
    68   if (!_firstEpoTime.valid()) {
    69     _firstEpoTime = _epoTime;
    70   }
    71 
    72   // Set Parameters
    73   // --------------
    74   _parlist->set(_epoTime, allObs);
    75   const vector<t_pppParam*>& params = _parlist->params();
    76 
    77   // Status Vector, Variance-Covariance Matrix
    78   // -----------------------------------------
    79   ColumnVector    xFltOld = _xFlt;
    80   SymmetricMatrix QFltOld = _QFlt;
    81 
    82   _QFlt.ReSize(_parlist->nPar()); _QFlt = 0.0;
    83   _xFlt.ReSize(_parlist->nPar()); _xFlt = 0.0;
    84   _x0.ReSize(_parlist->nPar());   _x0   = 0.0;
     277
     278  Matrix BB(epoData->sizeSys('G'), 4);
     279
     280  QMapIterator<QString, t_satData*> it(epoData->satData);
     281  int iObsBanc = 0;
     282  while (it.hasNext()) {
     283    it.next();
     284    t_satData* satData = it.value();
     285    if (satData->system() == 'G') {
     286      ++iObsBanc;
     287      QString    prn     = it.key();
     288      BB(iObsBanc, 1) = satData->xx(1);
     289      BB(iObsBanc, 2) = satData->xx(2);
     290      BB(iObsBanc, 3) = satData->xx(3);
     291      BB(iObsBanc, 4) = satData->P3 + satData->clk;
     292    }
     293  }
     294
     295  bancroft(BB, _xcBanc);
     296
     297  // Ellipsoidal Coordinates
     298  // ------------------------
     299  xyz2ell(_xcBanc.data(), _ellBanc.data());
     300
     301  // Compute Satellite Elevations
     302  // ----------------------------
     303  QMutableMapIterator<QString, t_satData*> im(epoData->satData);
     304  while (im.hasNext()) {
     305    im.next();
     306    t_satData* satData = im.value();
     307    cmpEle(satData);
     308    if (satData->eleSat < OPT->_minEle) {
     309      delete satData;
     310      im.remove();
     311    }
     312  }
     313
     314  return success;
     315}
     316
     317// Computed Value
     318////////////////////////////////////////////////////////////////////////////
     319double t_pppFilter::cmpValue(t_satData* satData, bool phase) {
     320
     321  Tracer tracer("t_pppFilter::cmpValue");
     322
     323  ColumnVector xRec(3);
     324  xRec(1) = x();
     325  xRec(2) = y();
     326  xRec(3) = z();
     327
     328  double rho0 = (satData->xx - xRec).norm_Frobenius();
     329  double dPhi = t_CST::omega * rho0 / t_CST::c;
     330
     331  xRec(1) = x() * cos(dPhi) - y() * sin(dPhi);
     332  xRec(2) = y() * cos(dPhi) + x() * sin(dPhi);
     333  xRec(3) = z();
     334
     335  xRec += _tides->displacement(_time, xRec);
     336
     337  satData->rho = (satData->xx - xRec).norm_Frobenius();
     338
     339  double tropDelay = delay_saast(satData->eleSat) +
     340                     trp() / sin(satData->eleSat);
     341
     342  double wind = 0.0;
     343  if (phase) {
     344    wind = windUp(satData->prn, satData->xx, xRec) * satData->lambda3;
     345  }
     346
     347  double offset = 0.0;
     348  t_frequency::type frqA = t_frequency::G1;
     349  t_frequency::type frqB = t_frequency::G2;
     350  if      (satData->prn[0] == 'R') {
     351    offset = Glonass_offset();
     352    frqA = t_frequency::R1;
     353    frqB = t_frequency::R2;
     354  }
     355  else if (satData->prn[0] == 'E') {
     356    offset = Galileo_offset();
     357    //frqA = t_frequency::E1; as soon as available
     358    //frqB = t_frequency::E5; -"-
     359  }
     360  else if (satData->prn[0] == 'C') {
     361    offset = Bds_offset();
     362    //frqA = t_frequency::C2; as soon as available
     363    //frqB = t_frequency::C7; -"-
     364  }
     365  double phaseCenter = 0.0;
     366  if (_antex) {
     367    bool found;
     368    phaseCenter = satData->lkA * _antex->rcvCorr(OPT->_antNameRover, frqA,
     369                                                 satData->eleSat, satData->azSat,
     370                                                 found)
     371                + satData->lkB * _antex->rcvCorr(OPT->_antNameRover, frqB,
     372                                                 satData->eleSat, satData->azSat,
     373                                                 found);
     374    if (!found) {
     375      LOG << "ANTEX: antenna >" << OPT->_antNameRover << "< not found\n";
     376    }
     377  }
     378
     379  double antennaOffset = 0.0;
     380  double cosa = cos(satData->azSat);
     381  double sina = sin(satData->azSat);
     382  double cose = cos(satData->eleSat);
     383  double sine = sin(satData->eleSat);
     384  antennaOffset = -OPT->_neuEccRover(1) * cosa*cose
     385                  -OPT->_neuEccRover(2) * sina*cose
     386                  -OPT->_neuEccRover(3) * sine;
     387
     388  return satData->rho + phaseCenter + antennaOffset + clk()
     389                      + offset - satData->clk + tropDelay + wind;
     390}
     391
     392// Tropospheric Model (Saastamoinen)
     393////////////////////////////////////////////////////////////////////////////
     394double t_pppFilter::delay_saast(double Ele) {
     395
     396  Tracer tracer("t_pppFilter::delay_saast");
     397
     398  double xyz[3];
     399  xyz[0] = x();
     400  xyz[1] = y();
     401  xyz[2] = z();
     402  double ell[3];
     403  xyz2ell(xyz, ell);
     404  double height = ell[2];
     405
     406  double pp =  1013.25 * pow(1.0 - 2.26e-5 * height, 5.225);
     407  double TT =  18.0 - height * 0.0065 + 273.15;
     408  double hh =  50.0 * exp(-6.396e-4 * height);
     409  double ee =  hh / 100.0 * exp(-37.2465 + 0.213166*TT - 0.000256908*TT*TT);
     410
     411  double h_km = height / 1000.0;
    85412 
    86   for (unsigned ii = 0; ii < params.size(); ii++) {
    87     const t_pppParam* par1 = params[ii];
    88 
    89     _x0[ii] = par1->x0();
    90 
    91     int iOld = par1->indexOld();
    92     if (iOld < 0) {
    93       _QFlt[ii][ii] = par1->sigma0() * par1->sigma0(); // new parameter
    94     }
    95     else {
    96       _QFlt[ii][ii] = QFltOld[iOld][iOld] + par1->noise() * par1->noise();
    97       _xFlt[ii]     = xFltOld[iOld];
    98       for (unsigned jj = 0; jj < ii; jj++) {
    99         const t_pppParam* par2 = params[jj];
    100         int               jOld = par2->indexOld();
    101         if (jOld >= 0) {
    102           _QFlt[ii][jj] = QFltOld(iOld+1,jOld+1);
     413  if (h_km < 0.0) h_km = 0.0;
     414  if (h_km > 5.0) h_km = 5.0;
     415  int    ii   = int(h_km + 1);
     416  double href = ii - 1;
     417 
     418  double bCor[6];
     419  bCor[0] = 1.156;
     420  bCor[1] = 1.006;
     421  bCor[2] = 0.874;
     422  bCor[3] = 0.757;
     423  bCor[4] = 0.654;
     424  bCor[5] = 0.563;
     425 
     426  double BB = bCor[ii-1] + (bCor[ii]-bCor[ii-1]) * (h_km - href);
     427 
     428  double zen  = M_PI/2.0 - Ele;
     429
     430  return (0.002277/cos(zen)) * (pp + ((1255.0/TT)+0.05)*ee - BB*(tan(zen)*tan(zen)));
     431}
     432
     433// Prediction Step of the Filter
     434////////////////////////////////////////////////////////////////////////////
     435void t_pppFilter::predict(int iPhase, t_epoData* epoData) {
     436
     437  Tracer tracer("t_pppFilter::predict");
     438
     439  if (iPhase == 0) {
     440
     441    const double maxSolGap = 60.0;
     442
     443    bool firstCrd = false;
     444    if (!_lastTimeOK.valid() || (maxSolGap > 0.0 && _time - _lastTimeOK > maxSolGap)) {
     445      firstCrd = true;
     446      _startTime = epoData->tt;
     447      reset();
     448    }
     449   
     450    // Use different white noise for Quick-Start mode
     451    // ----------------------------------------------
     452    double sigCrdP_used = OPT->_noiseCrd(1);
     453    if ( OPT->_seedingTime > 0.0 && OPT->_seedingTime > (epoData->tt - _startTime) ) {
     454      sigCrdP_used   = 0.0;
     455    }
     456
     457    // Predict Parameter values, add white noise
     458    // -----------------------------------------
     459    for (int iPar = 1; iPar <= _params.size(); iPar++) {
     460      t_pppParam* pp = _params[iPar-1];
     461   
     462      // Coordinates
     463      // -----------
     464      if      (pp->type == t_pppParam::CRD_X) {
     465        if (firstCrd) {
     466          if (OPT->xyzAprRoverSet()) {
     467            pp->xx = OPT->_xyzAprRover[0];
     468          }
     469          else {
     470            pp->xx = _xcBanc(1);
     471          }
    103472        }
    104       }
    105     }
    106   }
    107 
    108   predictCovCrdPart(QFltOld);
    109 
    110   // Process Satellite Systems separately
    111   // ------------------------------------
    112   for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
    113     char system = OPT->systems()[iSys];
    114     vector<t_pppSatObs*> obsVector;
    115     for (unsigned jj = 0; jj < allObs.size(); jj++) {
    116       if (allObs[jj]->prn().system() == system) {
    117         obsVector.push_back(allObs[jj]);
    118       }
    119     }
    120     if ( processSystem(OPT->LCs(system), obsVector) != success ) {
    121       return failure;
    122     }
    123   }
    124    
    125   cmpDOP(allObs);
    126 
    127   _parlist->printResult(_epoTime, _QFlt, _xFlt);
    128 
    129   return success;
    130 }
    131 
    132 // Process Selected LCs
    133 ////////////////////////////////////////////////////////////////////////////
    134 t_irc t_pppFilter::processSystem(const vector<t_lc::type>& LCs,
    135                                  const vector<t_pppSatObs*>& obsVector) {
    136 
    137   LOG.setf(ios::fixed);
    138 
    139   // Detect Cycle Slips
    140   // ------------------
    141   if (detectCycleSlips(LCs, obsVector) != success) {
    142     return failure;
    143   }
    144 
    145   ColumnVector            xSav       = _xFlt;
    146   SymmetricMatrix         QSav       = _QFlt;
    147   string                  epoTimeStr = string(_epoTime);
    148   const vector<t_pppParam*>& params     = _parlist->params();
    149   unsigned                maxObs     = obsVector.size() * LCs.size();
    150    
     473        _QQ(iPar,iPar) += sigCrdP_used * sigCrdP_used;
     474      }
     475      else if (pp->type == t_pppParam::CRD_Y) {
     476        if (firstCrd) {
     477          if (OPT->xyzAprRoverSet()) {
     478            pp->xx = OPT->_xyzAprRover[1];
     479          }
     480          else {
     481            pp->xx = _xcBanc(2);
     482          }
     483        }
     484        _QQ(iPar,iPar) += sigCrdP_used * sigCrdP_used;
     485      }
     486      else if (pp->type == t_pppParam::CRD_Z) {
     487        if (firstCrd) {
     488          if (OPT->xyzAprRoverSet()) {
     489            pp->xx = OPT->_xyzAprRover[2];
     490          }
     491          else {
     492            pp->xx = _xcBanc(3);
     493          }
     494        }
     495        _QQ(iPar,iPar) += sigCrdP_used * sigCrdP_used;
     496      }   
     497   
     498      // Receiver Clocks
     499      // ---------------
     500      else if (pp->type == t_pppParam::RECCLK) {
     501        pp->xx = _xcBanc(4);
     502        for (int jj = 1; jj <= _params.size(); jj++) {
     503          _QQ(iPar, jj) = 0.0;
     504        }
     505        _QQ(iPar,iPar) = OPT->_noiseClk * OPT->_noiseClk;
     506      }
     507   
     508      // Tropospheric Delay
     509      // ------------------
     510      else if (pp->type == t_pppParam::TROPO) {
     511        _QQ(iPar,iPar) += OPT->_noiseTrp * OPT->_noiseTrp;
     512      }
     513   
     514      // Glonass Offset
     515      // --------------
     516      else if (pp->type == t_pppParam::GLONASS_OFFSET) {
     517        pp->xx = 0.0;
     518        for (int jj = 1; jj <= _params.size(); jj++) {
     519          _QQ(iPar, jj) = 0.0;
     520        }
     521        _QQ(iPar,iPar) = 1000.0 * 1000.0;
     522      }
     523
     524      // Galileo Offset
     525      // --------------
     526      else if (pp->type == t_pppParam::GALILEO_OFFSET) {
     527        _QQ(iPar,iPar) += 0.1 * 0.1;
     528      }
     529
     530      // BDS Offset
     531      // ----------
     532      else if (pp->type == t_pppParam::BDS_OFFSET) {
     533        _QQ(iPar,iPar) = 1000.0 * 1000.0;    //TODO: TEST
     534      }
     535    }
     536  }
     537
     538  // Add New Ambiguities if necessary
     539  // --------------------------------
     540  if (OPT->ambLCs('G').size() || OPT->ambLCs('R').size() ||
     541      OPT->ambLCs('E').size() || OPT->ambLCs('C').size()) {
     542
     543    // Make a copy of QQ and xx, set parameter indices
     544    // -----------------------------------------------
     545    SymmetricMatrix QQ_old = _QQ;
     546   
     547    for (int iPar = 1; iPar <= _params.size(); iPar++) {
     548      _params[iPar-1]->index_old = _params[iPar-1]->index;
     549      _params[iPar-1]->index     = 0;
     550    }
     551   
     552    // Remove Ambiguity Parameters without observations
     553    // ------------------------------------------------
     554    int iPar = 0;
     555    QMutableVectorIterator<t_pppParam*> im(_params);
     556    while (im.hasNext()) {
     557      t_pppParam* par = im.next();
     558      bool removed = false;
     559      if (par->type == t_pppParam::AMB_L3) {
     560        if (epoData->satData.find(par->prn) == epoData->satData.end()) {
     561          removed = true;
     562          delete par;
     563          im.remove();
     564        }
     565      }
     566      if (! removed) {
     567        ++iPar;
     568        par->index = iPar;
     569      }
     570    }
     571   
     572    // Add new ambiguity parameters
     573    // ----------------------------
     574    QMapIterator<QString, t_satData*> it(epoData->satData);
     575    while (it.hasNext()) {
     576      it.next();
     577      t_satData* satData = it.value();
     578      addAmb(satData);
     579    }
     580   
     581    int nPar = _params.size();
     582    _QQ.ReSize(nPar); _QQ = 0.0;
     583    for (int i1 = 1; i1 <= nPar; i1++) {
     584      t_pppParam* p1 = _params[i1-1];
     585      if (p1->index_old != 0) {
     586        _QQ(p1->index, p1->index) = QQ_old(p1->index_old, p1->index_old);
     587        for (int i2 = 1; i2 <= nPar; i2++) {
     588          t_pppParam* p2 = _params[i2-1];
     589          if (p2->index_old != 0) {
     590            _QQ(p1->index, p2->index) = QQ_old(p1->index_old, p2->index_old);
     591          }
     592        }
     593      }
     594    }
     595   
     596    for (int ii = 1; ii <= nPar; ii++) {
     597      t_pppParam* par = _params[ii-1];
     598      if (par->index_old == 0) {
     599        _QQ(par->index, par->index) = OPT->_aprSigAmb * OPT->_aprSigAmb;
     600      }
     601      par->index_old = par->index;
     602    }
     603  }
     604}
     605
     606// Update Step of the Filter (currently just a single-epoch solution)
     607////////////////////////////////////////////////////////////////////////////
     608t_irc t_pppFilter::update(t_epoData* epoData) {
     609
     610  Tracer tracer("t_pppFilter::update");
     611
     612  _time = epoData->tt; // current epoch time
     613
     614  if (OPT->useOrbClkCorr()) {
     615    LOG << "Precise Point Positioning of Epoch " << _time.datestr() <<  "_" << _time.timestr(1)
     616        << "\n---------------------------------------------------------------\n";
     617  }
     618  else {
     619    LOG << "Single Point Positioning of Epoch " << _time.datestr() <<  "_" << _time.timestr(1)
     620        << "\n--------------------------------------------------------------\n";
     621  }
     622
    151623  // Outlier Detection Loop
    152624  // ----------------------
    153   for (unsigned iOutlier = 0; iOutlier < maxObs; iOutlier++) {
    154 
    155     if (iOutlier > 0) {
    156       _xFlt = xSav;
    157       _QFlt = QSav;
    158     }
    159 
    160     // First-Design Matrix, Terms Observed-Computed, Weight Matrix
    161     // -----------------------------------------------------------
    162     Matrix                AA(maxObs, _parlist->nPar());
    163     ColumnVector          ll(maxObs);
    164     DiagonalMatrix        PP(maxObs); PP = 0.0;
    165    
    166     int iObs = -1;
    167     vector<t_pppSatObs*> usedObs;
    168     vector<t_lc::type>   usedTypes;
    169     for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    170       t_pppSatObs* obs = obsVector[ii];
    171       if (!obs->outlier()) {
    172         for (unsigned jj = 0; jj < LCs.size(); jj++) {
    173           const t_lc::type tLC = LCs[jj];
    174           ++iObs;
    175           usedObs.push_back(obs);
    176           usedTypes.push_back(tLC);
    177           for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    178             const t_pppParam* par = params[iPar];
    179             AA[iObs][iPar] = par->partial(_epoTime, obs, tLC);
     625  if (update_p(epoData) != success) {
     626    return failure;
     627  }
     628
     629  // Set Solution Vector
     630  // -------------------
     631  LOG.setf(ios::fixed);
     632  QVectorIterator<t_pppParam*> itPar(_params);
     633  while (itPar.hasNext()) {
     634    t_pppParam* par = itPar.next();
     635    if      (par->type == t_pppParam::RECCLK) {
     636      LOG << "\n    clk     = " << setw(10) << setprecision(3) << par->xx
     637          << " +- " << setw(6) << setprecision(3)
     638          << sqrt(_QQ(par->index,par->index));
     639    }
     640    else if (par->type == t_pppParam::AMB_L3) {
     641      ++par->numEpo;
     642      LOG << "\n    amb " << par->prn.mid(0,3).toAscii().data() << " = "
     643          << setw(10) << setprecision(3) << par->xx
     644          << " +- " << setw(6) << setprecision(3)
     645          << sqrt(_QQ(par->index,par->index))
     646          << "   nEpo = " << par->numEpo;
     647    }
     648    else if (par->type == t_pppParam::TROPO) {
     649      double aprTrp = delay_saast(M_PI/2.0);
     650      LOG << "\n    trp     = " << par->prn.mid(0,3).toAscii().data()
     651          << setw(7) << setprecision(3) << aprTrp << " "
     652          << setw(6) << setprecision(3) << showpos << par->xx << noshowpos
     653          << " +- " << setw(6) << setprecision(3)
     654          << sqrt(_QQ(par->index,par->index));
     655    }
     656    else if (par->type == t_pppParam::GLONASS_OFFSET) {
     657      LOG << "\n    offGlo  = " << setw(10) << setprecision(3) << par->xx
     658          << " +- " << setw(6) << setprecision(3)
     659          << sqrt(_QQ(par->index,par->index));
     660    }
     661    else if (par->type == t_pppParam::GALILEO_OFFSET) {
     662      LOG << "\n    offGal  = " << setw(10) << setprecision(3) << par->xx
     663          << " +- " << setw(6) << setprecision(3)
     664          << sqrt(_QQ(par->index,par->index));
     665    }
     666    else if (par->type == t_pppParam::BDS_OFFSET) {
     667      LOG << "\n    offBds  = " << setw(10) << setprecision(3) << par->xx
     668          << " +- " << setw(6) << setprecision(3)
     669          << sqrt(_QQ(par->index,par->index));
     670    }
     671  }
     672
     673  LOG << endl << endl;
     674
     675  // Compute dilution of precision
     676  // -----------------------------
     677  cmpDOP(epoData);
     678
     679  // Final Message (both log file and screen)
     680  // ----------------------------------------
     681  LOG << OPT->_roverName << "  PPP "
     682      << epoData->tt.timestr(1) << " " << epoData->sizeAll() << " "
     683      << setw(14) << setprecision(3) << x()                  << " +- "
     684      << setw(6)  << setprecision(3) << sqrt(_QQ(1,1))       << " "
     685      << setw(14) << setprecision(3) << y()                  << " +- "
     686      << setw(6)  << setprecision(3) << sqrt(_QQ(2,2))       << " "
     687      << setw(14) << setprecision(3) << z()                  << " +- "
     688      << setw(6)  << setprecision(3) << sqrt(_QQ(3,3));
     689
     690  // NEU Output
     691  // ----------
     692  if (OPT->xyzAprRoverSet()) {
     693    double xyz[3];
     694    xyz[0] = x() - OPT->_xyzAprRover[0];
     695    xyz[1] = y() - OPT->_xyzAprRover[1];
     696    xyz[2] = z() - OPT->_xyzAprRover[2];
     697
     698    double ellRef[3];
     699    xyz2ell(OPT->_xyzAprRover.data(), ellRef);
     700    xyz2neu(ellRef, xyz, _neu.data());
     701
     702    LOG << "  NEU "
     703        << setw(8) << setprecision(3) << _neu[0] << " "
     704        << setw(8) << setprecision(3) << _neu[1] << " "
     705        << setw(8) << setprecision(3) << _neu[2] << endl << endl;
     706  }
     707  else {
     708    LOG << endl << endl;
     709  }
     710
     711  _lastTimeOK = _time; // remember time of last successful update
     712  return success;
     713}
     714
     715// Outlier Detection
     716////////////////////////////////////////////////////////////////////////////
     717QString t_pppFilter::outlierDetection(int iPhase, const ColumnVector& vv,
     718                                   QMap<QString, t_satData*>& satData) {
     719
     720  Tracer tracer("t_pppFilter::outlierDetection");
     721
     722  QString prnGPS;
     723  QString prnGlo;
     724  double  maxResGPS = 0.0; // all the other systems except GLONASS
     725  double  maxResGlo = 0.0; // GLONASS
     726  findMaxRes(vv, satData, prnGPS, prnGlo, maxResGPS, maxResGlo);
     727
     728  if      (iPhase == 1) {
     729    if      (maxResGlo > 2.98 * OPT->_maxResL1) {
     730      LOG << "Outlier Phase " << prnGlo.mid(0,3).toAscii().data() << ' ' << maxResGlo << endl;
     731      return prnGlo;
     732    }
     733    else if (maxResGPS > 2.98 * OPT->_maxResL1) {
     734      LOG << "Outlier Phase " << prnGPS.mid(0,3).toAscii().data() << ' ' << maxResGPS << endl;
     735      return prnGPS;
     736    }
     737  }
     738  else if (iPhase == 0 && maxResGPS > 2.98 * OPT->_maxResC1) {
     739    LOG << "Outlier Code  " << prnGPS.mid(0,3).toAscii().data() << ' ' << maxResGPS << endl;
     740    return prnGPS;
     741  }
     742
     743  return QString();
     744}
     745
     746// Phase Wind-Up Correction
     747///////////////////////////////////////////////////////////////////////////
     748double t_pppFilter::windUp(const QString& prn, const ColumnVector& rSat,
     749                        const ColumnVector& rRec) {
     750
     751  Tracer tracer("t_pppFilter::windUp");
     752
     753  double Mjd = _time.mjd() + _time.daysec() / 86400.0;
     754
     755  // First time - initialize to zero
     756  // -------------------------------
     757  if (!_windUpTime.contains(prn)) {
     758    _windUpSum[prn]  = 0.0;
     759  }
     760
     761  // Compute the correction for new time
     762  // -----------------------------------
     763  if (!_windUpTime.contains(prn) || _windUpTime[prn] != Mjd) {
     764    _windUpTime[prn] = Mjd;
     765
     766    // Unit Vector GPS Satellite --> Receiver
     767    // --------------------------------------
     768    ColumnVector rho = rRec - rSat;
     769    rho /= rho.norm_Frobenius();
     770   
     771    // GPS Satellite unit Vectors sz, sy, sx
     772    // -------------------------------------
     773    ColumnVector sz = -rSat / rSat.norm_Frobenius();
     774
     775    ColumnVector xSun = t_astro::Sun(Mjd);
     776    xSun /= xSun.norm_Frobenius();
     777
     778    ColumnVector sy = crossproduct(sz, xSun);
     779    ColumnVector sx = crossproduct(sy, sz);
     780
     781    // Effective Dipole of the GPS Satellite Antenna
     782    // ---------------------------------------------
     783    ColumnVector dipSat = sx - rho * DotProduct(rho,sx)
     784                                                - crossproduct(rho, sy);
     785   
     786    // Receiver unit Vectors rx, ry
     787    // ----------------------------
     788    ColumnVector rx(3);
     789    ColumnVector ry(3);
     790
     791    double recEll[3]; xyz2ell(rRec.data(), recEll) ;
     792    double neu[3];
     793   
     794    neu[0] = 1.0;
     795    neu[1] = 0.0;
     796    neu[2] = 0.0;
     797    neu2xyz(recEll, neu, rx.data());
     798   
     799    neu[0] =  0.0;
     800    neu[1] = -1.0;
     801    neu[2] =  0.0;
     802    neu2xyz(recEll, neu, ry.data());
     803   
     804    // Effective Dipole of the Receiver Antenna
     805    // ----------------------------------------
     806    ColumnVector dipRec = rx - rho * DotProduct(rho,rx)
     807                                                   + crossproduct(rho, ry);
     808   
     809    // Resulting Effect
     810    // ----------------
     811    double alpha = DotProduct(dipSat,dipRec) /
     812                      (dipSat.norm_Frobenius() * dipRec.norm_Frobenius());
     813   
     814    if (alpha >  1.0) alpha =  1.0;
     815    if (alpha < -1.0) alpha = -1.0;
     816   
     817    double dphi = acos(alpha) / 2.0 / M_PI;  // in cycles
     818   
     819    if ( DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0 ) {
     820      dphi = -dphi;
     821    }
     822
     823    _windUpSum[prn] = floor(_windUpSum[prn] - dphi + 0.5) + dphi;
     824  }
     825
     826  return _windUpSum[prn]; 
     827}
     828
     829//
     830///////////////////////////////////////////////////////////////////////////
     831void t_pppFilter::cmpEle(t_satData* satData) {
     832  Tracer tracer("t_pppFilter::cmpEle");
     833  ColumnVector rr = satData->xx - _xcBanc.Rows(1,3);
     834  double       rho = rr.norm_Frobenius();
     835
     836  double neu[3];
     837  xyz2neu(_ellBanc.data(), rr.data(), neu);
     838
     839  satData->eleSat = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / rho );
     840  if (neu[2] < 0) {
     841    satData->eleSat *= -1.0;
     842  }
     843  satData->azSat  = atan2(neu[1], neu[0]);
     844}
     845
     846//
     847///////////////////////////////////////////////////////////////////////////
     848void t_pppFilter::addAmb(t_satData* satData) {
     849  Tracer tracer("t_pppFilter::addAmb");
     850  if (!OPT->ambLCs(satData->system()).size()){
     851    return;
     852  }
     853  bool    found = false;
     854  for (int iPar = 1; iPar <= _params.size(); iPar++) {
     855    if (_params[iPar-1]->type == t_pppParam::AMB_L3 &&
     856        _params[iPar-1]->prn == satData->prn) {
     857      found = true;
     858      break;
     859    }
     860  }
     861  if (!found) {
     862    t_pppParam* par = new t_pppParam(t_pppParam::AMB_L3,
     863                                 _params.size()+1, satData->prn);
     864    _params.push_back(par);
     865    par->xx = satData->L3 - cmpValue(satData, true);
     866  }
     867}
     868
     869//
     870///////////////////////////////////////////////////////////////////////////
     871void t_pppFilter::addObs(int iPhase, unsigned& iObs, t_satData* satData,
     872                      Matrix& AA, ColumnVector& ll, DiagonalMatrix& PP) {
     873
     874  Tracer tracer("t_pppFilter::addObs");
     875
     876  const double ELEWGHT = 20.0;
     877  double ellWgtCoef = 1.0;
     878  double eleD = satData->eleSat * 180.0 / M_PI;
     879  if (eleD < ELEWGHT) {
     880    ellWgtCoef = 1.5 - 0.5 / (ELEWGHT - 10.0) * (eleD - 10.0);
     881  }
     882
     883  // Remember Observation Index
     884  // --------------------------
     885  ++iObs;
     886  satData->obsIndex = iObs;
     887
     888  // Phase Observations
     889  // ------------------
     890
     891  if (iPhase == 1) {
     892    ll(iObs)      = satData->L3 - cmpValue(satData, true);
     893    double sigL3 = 2.98 * OPT->_sigmaL1;
     894    if (satData->system() == 'R') {
     895      sigL3 *= GLONASS_WEIGHT_FACTOR;
     896    }
     897    if  (satData->system() == 'C') {
     898      sigL3 *= BDS_WEIGHT_FACTOR;
     899    }
     900    PP(iObs,iObs) = 1.0 / (sigL3 * sigL3) / (ellWgtCoef * ellWgtCoef);
     901    for (int iPar = 1; iPar <= _params.size(); iPar++) {
     902      if (_params[iPar-1]->type == t_pppParam::AMB_L3 &&
     903          _params[iPar-1]->prn  == satData->prn) {
     904        ll(iObs) -= _params[iPar-1]->xx;
     905      }
     906      AA(iObs, iPar) = _params[iPar-1]->partial(satData, true);
     907    }
     908  }
     909
     910  // Code Observations
     911  // -----------------
     912  else {
     913    double sigP3 = 2.98 * OPT->_sigmaC1;
     914    if  (satData->system() == 'C') {
     915      sigP3 *= BDS_WEIGHT_FACTOR;
     916    }
     917    ll(iObs)      = satData->P3 - cmpValue(satData, false);
     918    PP(iObs,iObs) = 1.0 / (sigP3 * sigP3) / (ellWgtCoef * ellWgtCoef);
     919    for (int iPar = 1; iPar <= _params.size(); iPar++) {
     920      AA(iObs, iPar) = _params[iPar-1]->partial(satData, false);
     921    }
     922  }
     923}
     924
     925//
     926///////////////////////////////////////////////////////////////////////////
     927QByteArray t_pppFilter::printRes(int iPhase, const ColumnVector& vv,
     928                              const QMap<QString, t_satData*>& satDataMap) {
     929
     930  Tracer tracer("t_pppFilter::printRes");
     931
     932  ostringstream str;
     933  str.setf(ios::fixed);
     934  bool useObs;
     935  QMapIterator<QString, t_satData*> it(satDataMap);
     936  while (it.hasNext()) {
     937    it.next();
     938    t_satData* satData = it.value();
     939    (iPhase == 0) ? useObs = OPT->codeLCs(satData->system()).size() :
     940                    useObs = OPT->ambLCs(satData->system()).size();
     941    if (satData->obsIndex != 0 && useObs) {
     942      str << _time.timestr(1)
     943          << " RES " << satData->prn.mid(0,3).toAscii().data()
     944          << (iPhase ? "   L3 " : "   P3 ")
     945          << setw(9) << setprecision(4) << vv(satData->obsIndex) << endl;
     946    }
     947  }
     948
     949  return QByteArray(str.str().c_str());
     950}
     951
     952//
     953///////////////////////////////////////////////////////////////////////////
     954void t_pppFilter::findMaxRes(const ColumnVector& vv,
     955                          const QMap<QString, t_satData*>& satData,
     956                          QString& prnGPS, QString& prnGlo,
     957                          double& maxResGPS, double& maxResGlo) {
     958
     959  Tracer tracer("t_pppFilter::findMaxRes");
     960
     961  maxResGPS  = 0.0;
     962  maxResGlo  = 0.0;
     963
     964  QMapIterator<QString, t_satData*> it(satData);
     965  while (it.hasNext()) {
     966    it.next();
     967    t_satData* satData = it.value();
     968    if (satData->obsIndex != 0) {
     969      QString prn = satData->prn;
     970      if (prn[0] == 'R') {
     971        if (fabs(vv(satData->obsIndex)) > maxResGlo) {
     972          maxResGlo = fabs(vv(satData->obsIndex));
     973          prnGlo    = prn;
     974        }
     975      }
     976      else {
     977        if (fabs(vv(satData->obsIndex)) > maxResGPS) {
     978          maxResGPS = fabs(vv(satData->obsIndex));
     979          prnGPS    = prn;
     980        }
     981      }
     982    }
     983  }
     984}
     985 
     986// Update Step (private - loop over outliers)
     987////////////////////////////////////////////////////////////////////////////
     988t_irc t_pppFilter::update_p(t_epoData* epoData) {
     989
     990  Tracer tracer("t_pppFilter::update_p");
     991
     992  // Save Variance-Covariance Matrix, and Status Vector
     993  // --------------------------------------------------
     994  rememberState(epoData);
     995
     996  QString lastOutlierPrn;
     997
     998  // Try with all satellites, then with all minus one, etc.
     999  // ------------------------------------------------------
     1000  while (selectSatellites(lastOutlierPrn, epoData->satData) == success) {
     1001
     1002    QByteArray strResCode;
     1003    QByteArray strResPhase;
     1004
     1005    // Bancroft Solution
     1006    // -----------------
     1007    if (cmpBancroft(epoData) != success) {
     1008      break;
     1009    }
     1010
     1011    // First update using code observations, then phase observations
     1012    // -------------------------------------------------------------     
     1013    bool usePhase = OPT->ambLCs('G').size() || OPT->ambLCs('R').size() ||
     1014                    OPT->ambLCs('E').size() || OPT->ambLCs('C').size() ;
     1015
     1016    for (int iPhase = 0; iPhase <= (usePhase ? 1 : 0); iPhase++) {
     1017
     1018      // Status Prediction
     1019      // -----------------
     1020      predict(iPhase, epoData);
     1021     
     1022      // Create First-Design Matrix
     1023      // --------------------------
     1024      unsigned nPar = _params.size();
     1025      unsigned nObs = 0;
     1026      nObs = epoData->sizeAll();
     1027      bool useObs = false;
     1028      char sys[] ={'G', 'R', 'E', 'C'};
     1029      for (unsigned ii = 0; ii < sizeof(sys); ii++) {
     1030        const char s = sys[ii];
     1031        (iPhase == 0) ? useObs = OPT->codeLCs(s).size() : useObs = OPT->ambLCs(s).size();
     1032        if (!useObs) {
     1033          nObs -= epoData->sizeSys(s);
     1034        }
     1035      }
     1036     
     1037      // Prepare first-design Matrix, vector observed-computed
     1038      // -----------------------------------------------------
     1039      Matrix          AA(nObs, nPar);  // first design matrix
     1040      ColumnVector    ll(nObs);        // tems observed-computed
     1041      DiagonalMatrix  PP(nObs); PP = 0.0;
     1042     
     1043      unsigned iObs = 0;
     1044      QMapIterator<QString, t_satData*> it(epoData->satData);
     1045      while (it.hasNext()) {
     1046        it.next();
     1047        t_satData* satData = it.value();
     1048        QString prn = satData->prn;
     1049        (iPhase == 0) ? useObs = OPT->codeLCs(satData->system()).size() :
     1050                        useObs = OPT->ambLCs(satData->system()).size();
     1051        if (useObs) {
     1052          addObs(iPhase, iObs, satData, AA, ll, PP);
     1053        } else {
     1054          satData->obsIndex = 0;
     1055        }
     1056      }
     1057
     1058      // Compute Filter Update
     1059      // ---------------------
     1060      ColumnVector dx(nPar); dx = 0.0;
     1061      kalman(AA, ll, PP, _QQ, dx);
     1062      ColumnVector vv = ll - AA * dx;
     1063     
     1064      // Print Residuals
     1065      // ---------------
     1066      if      (iPhase == 0) {
     1067        strResCode  = printRes(iPhase, vv, epoData->satData);
     1068      }
     1069      else {
     1070        strResPhase = printRes(iPhase, vv, epoData->satData);
     1071      }
     1072
     1073      // Check the residuals
     1074      // -------------------
     1075      lastOutlierPrn = outlierDetection(iPhase, vv, epoData->satData);
     1076
     1077      // No Outlier Detected
     1078      // -------------------
     1079      if (lastOutlierPrn.isEmpty()) {
     1080
     1081        QVectorIterator<t_pppParam*> itPar(_params);
     1082        while (itPar.hasNext()) {
     1083          t_pppParam* par = itPar.next();
     1084          par->xx += dx(par->index);
     1085        }
     1086
     1087        if (!usePhase || iPhase == 1) {
     1088          if (_outlierGPS.size() > 0 || _outlierGlo.size() > 0) {
     1089            LOG << "Neglected PRNs: ";
     1090            if (!_outlierGPS.isEmpty()) {
     1091              LOG << _outlierGPS.last().mid(0,3).toAscii().data() << ' ';
     1092            }
     1093            QStringListIterator itGlo(_outlierGlo);
     1094            while (itGlo.hasNext()) {
     1095              QString prn = itGlo.next();
     1096              LOG << prn.mid(0,3).toAscii().data() << ' ';
     1097            }
    1801098          }
    181           ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1));
    182           PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC));
     1099          LOG << strResCode.data() << strResPhase.data();
     1100
     1101          return success;
    1831102        }
    1841103      }
    185     }
    186 
    187     // Check number of observations, truncate matrices
    188     // -----------------------------------------------
    189     if (iObs == -1) {
    190       return failure;
    191     }
    192     AA = AA.Rows(1, iObs+1);
    193     ll = ll.Rows(1, iObs+1);
    194     PP = PP.SymSubMatrix(1, iObs+1);
    195 
    196     // Kalman update step
    197     // ------------------
    198     kalman(AA, ll, PP, _QFlt, _xFlt);
    199 
    200     // Check Residuals
    201     // ---------------
    202     ColumnVector vv = AA * _xFlt - ll;
    203     double     maxOutlier      = 0.0;
    204     int        maxOutlierIndex = -1;
    205     t_lc::type maxOutlierLC    = t_lc::dummy;
    206     for (unsigned ii = 0; ii < usedObs.size(); ii++) {
    207       const t_lc::type tLC = usedTypes[ii];
    208       double res = fabs(vv[ii]);
    209       if (res > usedObs[ii]->maxRes(tLC)) {
    210         if (res > fabs(maxOutlier)) {
    211           maxOutlier      = vv[ii];
    212           maxOutlierIndex = ii;
    213           maxOutlierLC    = tLC;
    214         }
    215       }
    216     }
    217 
    218     // Mark outlier or break outlier detection loop
    219     // --------------------------------------------
    220     if (maxOutlierIndex > -1) {
    221       t_pppSatObs* obs = usedObs[maxOutlierIndex];
    222       t_pppParam* par = 0;
    223       LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' '
    224           << obs->prn().toString()                        << ' '
    225           << setw(8) << setprecision(4) << maxOutlier << endl;
    226       for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    227         t_pppParam* hlp = params[iPar];
    228         if (hlp->type() == t_pppParam::amb && hlp->prn()  == obs->prn() &&
    229             hlp->tLC()  == usedTypes[maxOutlierIndex]) {
    230           par = hlp;
    231         }
    232       }
    233       if (par) {
    234         if (par->ambResetCandidate()) {
    235           resetAmb(par->prn(), obsVector, &QSav, &xSav);
    236         }
    237         else {
    238           par->setAmbResetCandidate();
    239           obs->setOutlier();
    240         }
    241       }
     1104
     1105      // Outlier Found
     1106      // -------------
    2421107      else {
    243         obs->setOutlier();
    244       }
    245     }
    246 
    247     // Print Residuals
    248     // ---------------
     1108        restoreState(epoData);
     1109        break;
     1110      }
     1111
     1112    } // for iPhase
     1113
     1114  } // while selectSatellites
     1115
     1116  restoreState(epoData);
     1117  return failure;
     1118}
     1119
     1120// Remeber Original State Vector and Variance-Covariance Matrix
     1121////////////////////////////////////////////////////////////////////////////
     1122void t_pppFilter::rememberState(t_epoData* epoData) {
     1123
     1124  _QQ_sav = _QQ;
     1125
     1126  QVectorIterator<t_pppParam*> itSav(_params_sav);
     1127  while (itSav.hasNext()) {
     1128    t_pppParam* par = itSav.next();
     1129    delete par;
     1130  }
     1131  _params_sav.clear();
     1132
     1133  QVectorIterator<t_pppParam*> it(_params);
     1134  while (it.hasNext()) {
     1135    t_pppParam* par = it.next();
     1136    _params_sav.push_back(new t_pppParam(*par));
     1137  }
     1138
     1139  _epoData_sav->deepCopy(epoData);
     1140}
     1141
     1142// Restore Original State Vector and Variance-Covariance Matrix
     1143////////////////////////////////////////////////////////////////////////////
     1144void t_pppFilter::restoreState(t_epoData* epoData) {
     1145
     1146  _QQ = _QQ_sav;
     1147
     1148  QVectorIterator<t_pppParam*> it(_params);
     1149  while (it.hasNext()) {
     1150    t_pppParam* par = it.next();
     1151    delete par;
     1152  }
     1153  _params.clear();
     1154
     1155  QVectorIterator<t_pppParam*> itSav(_params_sav);
     1156  while (itSav.hasNext()) {
     1157    t_pppParam* par = itSav.next();
     1158    _params.push_back(new t_pppParam(*par));
     1159  }
     1160
     1161  epoData->deepCopy(_epoData_sav);
     1162}
     1163
     1164//
     1165////////////////////////////////////////////////////////////////////////////
     1166t_irc t_pppFilter::selectSatellites(const QString& lastOutlierPrn,
     1167                                 QMap<QString, t_satData*>& satData) {
     1168
     1169  // First Call
     1170  // ----------
     1171  if (lastOutlierPrn.isEmpty()) {
     1172    _outlierGPS.clear();
     1173    _outlierGlo.clear();
     1174    return success;
     1175  }
     1176
     1177  // Second and next trials
     1178  // ----------------------
     1179  else {
     1180
     1181    if (lastOutlierPrn[0] == 'R') {
     1182      _outlierGlo << lastOutlierPrn;
     1183    }
     1184
     1185    // Remove all Glonass Outliers
     1186    // ---------------------------
     1187    QStringListIterator it(_outlierGlo);
     1188    while (it.hasNext()) {
     1189      QString prn = it.next();
     1190      if (satData.contains(prn)) {
     1191        delete satData.take(prn);
     1192      }
     1193    }
     1194
     1195    if (lastOutlierPrn[0] == 'R') {
     1196      _outlierGPS.clear();
     1197      return success;
     1198    }
     1199
     1200    // GPS Outlier appeared for the first time - try to delete it
     1201    // ----------------------------------------------------------
     1202    if (_outlierGPS.indexOf(lastOutlierPrn) == -1) {
     1203      _outlierGPS << lastOutlierPrn;
     1204      if (satData.contains(lastOutlierPrn)) {
     1205        delete satData.take(lastOutlierPrn);
     1206      }
     1207      return success;
     1208    }
     1209
     1210  }
     1211
     1212  return failure;
     1213}
     1214
     1215//
     1216////////////////////////////////////////////////////////////////////////////
     1217double lorentz(const ColumnVector& aa, const ColumnVector& bb) {
     1218  return aa(1)*bb(1) +  aa(2)*bb(2) +  aa(3)*bb(3) -  aa(4)*bb(4);
     1219}
     1220
     1221//
     1222////////////////////////////////////////////////////////////////////////////
     1223void t_pppFilter::bancroft(const Matrix& BBpass, ColumnVector& pos) {
     1224
     1225  if (pos.Nrows() != 4) {
     1226    pos.ReSize(4);
     1227  }
     1228  pos = 0.0;
     1229
     1230  for (int iter = 1; iter <= 2; iter++) {
     1231    Matrix BB = BBpass;
     1232    int mm = BB.Nrows();
     1233    for (int ii = 1; ii <= mm; ii++) {
     1234      double xx = BB(ii,1);
     1235      double yy = BB(ii,2);
     1236      double traveltime = 0.072;
     1237      if (iter > 1) {
     1238        double zz  = BB(ii,3);
     1239        double rho = sqrt( (xx-pos(1)) * (xx-pos(1)) +
     1240                           (yy-pos(2)) * (yy-pos(2)) +
     1241                           (zz-pos(3)) * (zz-pos(3)) );
     1242        traveltime = rho / t_CST::c;
     1243      }
     1244      double angle = traveltime * t_CST::omega;
     1245      double cosa  = cos(angle);
     1246      double sina  = sin(angle);
     1247      BB(ii,1) =  cosa * xx + sina * yy;
     1248      BB(ii,2) = -sina * xx + cosa * yy;
     1249    }
     1250   
     1251    Matrix BBB;
     1252    if (mm > 4) {
     1253      SymmetricMatrix hlp; hlp << BB.t() * BB;
     1254      BBB = hlp.i() * BB.t();
     1255    }
    2491256    else {
    250       for (unsigned jj = 0; jj < LCs.size(); jj++) {
    251         for (unsigned ii = 0; ii < usedObs.size(); ii++) {
    252           const t_lc::type tLC = usedTypes[ii];
    253           t_pppSatObs*        obs = usedObs[ii];
    254           if (tLC == LCs[jj]) {
    255             obs->setRes(tLC, vv[ii]);
    256             LOG << epoTimeStr << " RES "
    257                 << left << setw(3) << t_lc::toString(tLC) << right << ' '
    258                 << obs->prn().toString().substr(0,3) << ' '
    259                 << setw(8) << setprecision(4) << vv[ii] << endl;
    260           }
    261         }
    262       }
    263       break;
    264     }
    265   }
    266 
    267   return success;
    268 }
    269 
    270 // Cycle-Slip Detection
    271 ////////////////////////////////////////////////////////////////////////////
    272 t_irc t_pppFilter::detectCycleSlips(const vector<t_lc::type>& LCs,
    273                                     const vector<t_pppSatObs*>& obsVector) {
    274 
    275   const double            SLIP       = 20.0;  // slip threshold
    276   string                  epoTimeStr = string(_epoTime);
    277   const vector<t_pppParam*>& params     = _parlist->params();
    278 
    279   for (unsigned ii = 0; ii < LCs.size(); ii++) {
    280     const t_lc::type& tLC = LCs[ii];
    281     if (t_lc::includesPhase(tLC)) {
    282       for (unsigned iObs = 0; iObs < obsVector.size(); iObs++) {
    283         const t_pppSatObs* obs = obsVector[iObs];
    284 
    285         // Check set Slips and Jump Counters
    286         // ---------------------------------
    287         bool slip = false;
    288 
    289         if (obs->slip()) {
    290           LOG << "cycle slip set (obs)" << endl;;
    291           slip = true;
    292         }
    293 
    294         if (_slips[obs->prn()]._obsSlipCounter != -1 &&
    295             _slips[obs->prn()]._obsSlipCounter != obs->slipCounter()) {
    296           LOG << "cycle slip set (obsSlipCounter)" << endl;
    297           slip = true;
    298         }
    299         _slips[obs->prn()]._obsSlipCounter = obs->slipCounter();
    300 
    301         if (_slips[obs->prn()]._biasJumpCounter != -1 &&
    302             _slips[obs->prn()]._biasJumpCounter != obs->biasJumpCounter()) {
    303           LOG << "cycle slip set (biasJumpCounter)" << endl;
    304           slip = true;
    305         }
    306         _slips[obs->prn()]._biasJumpCounter = obs->biasJumpCounter();
    307 
    308         // Slip Set
    309         // -------- 
    310         if (slip) {
    311           resetAmb(obs->prn(), obsVector);
    312         }
    313  
    314         // Check Pre-Fit Residuals
    315         // -----------------------
    316         else {
    317           ColumnVector AA(params.size());
    318           for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    319             const t_pppParam* par = params[iPar];
    320             AA[iPar] = par->partial(_epoTime, obs, tLC);
    321           }
    322          
    323           double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA);
    324           double vv = DotProduct(AA, _xFlt) - ll;
    325          
    326           if (fabs(vv) > SLIP) {
    327             LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' '
    328                 << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl;
    329             resetAmb(obs->prn(), obsVector);
    330           }
    331         }
    332       }
    333     }
    334   }
    335 
    336   return success;
    337 }
    338 
    339 // Reset Ambiguity Parameter (cycle slip)
    340 ////////////////////////////////////////////////////////////////////////////
    341 t_irc t_pppFilter::resetAmb(t_prn prn, const vector<t_pppSatObs*>& obsVector,
    342                             SymmetricMatrix* QSav, ColumnVector* xSav) {
    343   t_irc irc = failure;
    344   vector<t_pppParam*>& params = _parlist->params();
    345   for (unsigned iPar = 0; iPar < params.size(); iPar++) {
    346     t_pppParam* par = params[iPar];
    347     if (par->type() == t_pppParam::amb && par->prn() == prn) {
    348       int ind = par->indexNew();
    349       t_lc::type tLC = par->tLC();
    350       LOG << string(_epoTime) << " RESET " << par->toString() << endl;
    351       delete par; par = new t_pppParam(t_pppParam::amb, prn, tLC, &obsVector);
    352       par->setIndex(ind);
    353       params[iPar] = par;
    354       for (unsigned ii = 1; ii <= params.size(); ii++) {
    355         _QFlt(ii, ind+1) = 0.0;
    356         if (QSav) {
    357           (*QSav)(ii, ind+1) = 0.0;
    358         }
    359       }
    360       _QFlt(ind+1,ind+1) = par->sigma0() * par->sigma0();
    361       if (QSav) {
    362         (*QSav)(ind+1,ind+1) = _QFlt(ind+1,ind+1);
    363       }
    364       _xFlt[ind] = 0.0;
    365       if (xSav) {
    366         (*xSav)[ind] = _xFlt[ind];
    367       }
    368       _x0[ind] = par->x0();
    369       irc = success;
    370     }
    371   }
    372 
    373   return irc;
    374 }
    375 
    376 // Compute various DOP Values
    377 ////////////////////////////////////////////////////////////////////////////
    378 void t_pppFilter::cmpDOP(const vector<t_pppSatObs*>& obsVector) {
    379 
    380   _dop.reset();
    381 
    382   try {
    383     const unsigned numPar = 4;
    384     Matrix AA(obsVector.size(), numPar);
    385     _numSat = 0;
    386     for (unsigned ii = 0; ii < obsVector.size(); ii++) {
    387       t_pppSatObs* obs = obsVector[ii];
    388       if (obs->isValid() && !obs->outlier()) {
    389         ++_numSat;
    390         for (unsigned iPar = 0; iPar < numPar; iPar++) {
    391           const t_pppParam* par = _parlist->params()[iPar];
    392           AA[_numSat-1][iPar] = par->partial(_epoTime, obs, t_lc::c1);
    393         }
    394       }
    395     }
    396     if (_numSat < 4) {
    397       return;
    398     }
    399     AA = AA.Rows(1, _numSat);
    400     SymmetricMatrix NN; NN << AA.t() * AA; 
    401     SymmetricMatrix QQ = NN.i();
    402    
    403     _dop.P = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3));
    404     _dop.T = sqrt(QQ(4,4));
    405     _dop.G = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3) + QQ(4,4));
    406   }
    407   catch (...) {
    408   }
    409 }
    410 
    411 // Compute various DOP Values
    412 ////////////////////////////////////////////////////////////////////////////
    413 void t_pppFilter::predictCovCrdPart(const SymmetricMatrix& QFltOld) {
    414 
    415   const vector<t_pppParam*>& params = _parlist->params();
    416   if (params.size() < 3) {
     1257      BBB = BB.i();
     1258    }
     1259    ColumnVector ee(mm); ee = 1.0;
     1260    ColumnVector alpha(mm); alpha = 0.0;
     1261    for (int ii = 1; ii <= mm; ii++) {
     1262      alpha(ii) = lorentz(BB.Row(ii).t(),BB.Row(ii).t())/2.0;
     1263    }
     1264    ColumnVector BBBe     = BBB * ee;
     1265    ColumnVector BBBalpha = BBB * alpha;
     1266    double aa = lorentz(BBBe, BBBe);
     1267    double bb = lorentz(BBBe, BBBalpha)-1;
     1268    double cc = lorentz(BBBalpha, BBBalpha);
     1269    double root = sqrt(bb*bb-aa*cc);
     1270
     1271    Matrix hlpPos(4,2);
     1272    hlpPos.Column(1) = (-bb-root)/aa * BBBe + BBBalpha;
     1273    hlpPos.Column(2) = (-bb+root)/aa * BBBe + BBBalpha;
     1274
     1275    ColumnVector omc(2);
     1276    for (int pp = 1; pp <= 2; pp++) {
     1277      hlpPos(4,pp)      = -hlpPos(4,pp);
     1278      omc(pp) = BB(1,4) -
     1279                sqrt( (BB(1,1)-hlpPos(1,pp)) * (BB(1,1)-hlpPos(1,pp)) +
     1280                      (BB(1,2)-hlpPos(2,pp)) * (BB(1,2)-hlpPos(2,pp)) +
     1281                      (BB(1,3)-hlpPos(3,pp)) * (BB(1,3)-hlpPos(3,pp)) ) -
     1282                hlpPos(4,pp);
     1283    }
     1284    if ( fabs(omc(1)) > fabs(omc(2)) ) {
     1285      pos = hlpPos.Column(2);
     1286    }
     1287    else {
     1288      pos = hlpPos.Column(1);
     1289    }
     1290  }
     1291}
     1292
     1293//
     1294////////////////////////////////////////////////////////////////////////////
     1295void t_pppFilter::cmpDOP(t_epoData* epoData) {
     1296
     1297  Tracer tracer("t_pppFilter::cmpDOP");
     1298
     1299  _numSat = 0;
     1300  _pDop   = 0.0;
     1301
     1302  if (_params.size() < 4) {
    4171303    return;
    4181304  }
    4191305
    420   bool first = (params[0]->indexOld() < 0);
    421 
    422   SymmetricMatrix Qneu(3); Qneu = 0.0;
    423   for (unsigned ii = 0; ii < 3; ii++) {
    424     const t_pppParam* par = params[ii];
    425     if (first) {
    426       Qneu[ii][ii] = par->sigma0() * par->sigma0();
    427     }
    428     else {
    429       Qneu[ii][ii] = par->noise() * par->noise();
    430     }
    431   }
    432 
    433   const t_pppStation* sta = PPP_CLIENT->staRover();
    434   SymmetricMatrix Qxyz(3);
    435   covariNEU_XYZ(Qneu, sta->ellApr().data(), Qxyz);
    436 
    437   if (first) {
    438     _QFlt.SymSubMatrix(1,3) = Qxyz;
    439   }
    440   else {
    441     double dt = _epoTime - _firstEpoTime;
    442     if (dt < OPT->_seedingTime) {
    443       _QFlt.SymSubMatrix(1,3) = QFltOld.SymSubMatrix(1,3);
    444     }
    445     else {
    446       _QFlt.SymSubMatrix(1,3) = QFltOld.SymSubMatrix(1,3) + Qxyz;
    447     }
    448   }
    449 }
     1306  const unsigned numPar = 4;
     1307  Matrix AA(epoData->sizeAll(), numPar);
     1308  QMapIterator<QString, t_satData*> it(epoData->satData);
     1309  while (it.hasNext()) {
     1310    it.next();
     1311    t_satData* satData = it.value();
     1312    _numSat += 1;
     1313    for (unsigned iPar = 0; iPar < numPar; iPar++) {
     1314      AA[_numSat-1][iPar] = _params[iPar]->partial(satData, false);
     1315    }
     1316  }
     1317  if (_numSat < 4) {
     1318    return;
     1319  }
     1320  AA = AA.Rows(1, _numSat);
     1321  SymmetricMatrix NN; NN << AA.t() * AA; 
     1322  SymmetricMatrix QQ = NN.i();
     1323   
     1324  _pDop = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3));
     1325}
  • trunk/BNC/src/PPP/pppFilter.h

    r6654 r7203  
    1 #ifndef FILTER_H
    2 #define FILTER_H
    3 
    4 #include <vector>
     1// Part of BNC, a utility for retrieving decoding and
     2// converting GNSS data streams from NTRIP broadcasters.
     3//
     4// Copyright (C) 2007
     5// German Federal Agency for Cartography and Geodesy (BKG)
     6// http://www.bkg.bund.de
     7// Czech Technical University Prague, Department of Geodesy
     8// http://www.fsv.cvut.cz
     9//
     10// Email: euref-ip@bkg.bund.de
     11//
     12// This program is free software; you can redistribute it and/or
     13// modify it under the terms of the GNU General Public License
     14// as published by the Free Software Foundation, version 2.
     15//
     16// This program is distributed in the hope that it will be useful,
     17// but WITHOUT ANY WARRANTY; without even the implied warranty of
     18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     19// GNU General Public License for more details.
     20//
     21// You should have received a copy of the GNU General Public License
     22// along with this program; if not, write to the Free Software
     23// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
     24
     25#ifndef PPPFILTER_H
     26#define PPPFILTER_H
     27
     28#include <QtCore>
     29#include <QtNetwork>
    530#include <newmat.h>
    6 #include "pppInclude.h"
    7 #include "pppParlist.h"
     31
     32#include "bncconst.h"
    833#include "bnctime.h"
    9 #include "t_prn.h"
     34
     35class bncAntex;
    1036
    1137namespace BNC_PPP {
    1238
    13 class t_pppParlist;
    14 class t_pppObsPool;
    15 class t_pppSatObs;
     39class t_pppClient;
     40class t_pppOptions;
     41class t_epoData;
     42class t_satData;
     43class t_tides;
     44
     45class t_satData {
     46 public:
     47  t_satData() {
     48    obsIndex = 0;
     49    P1      = 0.0;
     50    P2      = 0.0;
     51    P5      = 0.0;
     52    P3      = 0.0;
     53    L1      = 0.0;
     54    L2      = 0.0;
     55    L5      = 0.0;
     56    L3      = 0.0;
     57    lkA     = 0.0;
     58    lkB     = 0.0;
     59  }
     60  ~t_satData() {}
     61  bncTime      tt;
     62  QString      prn;
     63  double       P1;
     64  double       P2;
     65  double       P5;
     66  double       P7;
     67  double       P3;
     68  double       L1;
     69  double       L2;
     70  double       L5;
     71  double       L7;
     72  double       L3;
     73  ColumnVector xx;
     74  ColumnVector vv;
     75  double       clk;
     76  double       eleSat;
     77  double       azSat;
     78  double       rho;
     79  bool         slipFlag;
     80  double       lambda3;
     81  double       lkA;
     82  double       lkB;
     83  unsigned     obsIndex;
     84  char system() const {return prn.toAscii()[0];}
     85};
     86
     87class t_epoData {
     88 public:
     89  t_epoData() {}
     90
     91  ~t_epoData() {
     92    clear();
     93  }
     94
     95  void clear() {
     96    QMapIterator<QString, t_satData*> it(satData);
     97    while (it.hasNext()) {
     98      it.next();
     99      delete it.value();
     100    }
     101    satData.clear();
     102    tt.reset();
     103  }
     104
     105  void deepCopy(const t_epoData* from) {
     106    clear();
     107    tt = from->tt;
     108    QMapIterator<QString, t_satData*> it(from->satData);
     109    while (it.hasNext()) {
     110      it.next();
     111      satData[it.key()] = new t_satData(*it.value());
     112    }
     113  }
     114
     115  unsigned sizeSys(char system) const {
     116    unsigned ans = 0;
     117    QMapIterator<QString, t_satData*> it(satData);
     118    while (it.hasNext()) {
     119      it.next();
     120      if (it.value()->system() == system) {
     121        ++ans;
     122      }
     123    }
     124    return ans;
     125  }
     126  unsigned sizeAll() const {return satData.size();}
     127
     128  bncTime                   tt;
     129  QMap<QString, t_satData*> satData;
     130};
     131
     132class t_pppParam {
     133 public:
     134  enum parType {CRD_X, CRD_Y, CRD_Z, RECCLK, TROPO, AMB_L3,
     135                GLONASS_OFFSET, GALILEO_OFFSET, BDS_OFFSET};
     136  t_pppParam(parType typeIn, int indexIn, const QString& prn);
     137  ~t_pppParam();
     138  double partial(t_satData* satData, bool phase);
     139  bool isCrd() const {
     140    return (type == CRD_X || type == CRD_Y || type == CRD_Z);
     141  }
     142  parType  type;
     143  double   xx;
     144  int      index;
     145  int      index_old;
     146  int      numEpo;
     147  QString  prn;
     148};
    16149
    17150class t_pppFilter {
    18151 public:
    19   t_pppFilter();
     152  t_pppFilter(t_pppClient* pppClient);
    20153  ~t_pppFilter();
    21 
    22   t_irc processEpoch(t_pppObsPool* obsPool);
    23 
    24   const ColumnVector&    x() const {return _xFlt;}
    25   const SymmetricMatrix& Q() const {return _QFlt;}
    26 
     154  t_irc update(t_epoData* epoData);
     155  bncTime time()  const {return _time;}
     156  const SymmetricMatrix& Q() const {return _QQ;}
     157  const ColumnVector& neu() const {return _neu;}
    27158  int    numSat() const {return _numSat;}
    28   double PDOP() const {return _dop.P;}
    29   double GDOP() const {return _dop.G;}
     159  double PDOP()   const {return _pDop;}
     160  double x()      const {return _params[0]->xx;}
     161  double y()      const {return _params[1]->xx;}
     162  double z()      const {return _params[2]->xx;}
     163  double clk()    const {return _params[3]->xx;}
     164  double trp0()   {return delay_saast(M_PI/2.0);}
    30165  double trp() const {
    31     const std::vector<t_pppParam*>& par = _parlist->params();
    32     for (unsigned ii = 0; ii < par.size(); ++ii) {
    33       if (par[ii]->type() == t_pppParam::trp) {
    34         return x()[ii];
    35       }
    36     }
    37     return 0.0;
    38   };
     166    for (int ii = 0; ii < _params.size(); ++ii) {
     167      t_pppParam* pp = _params[ii];
     168      if (pp->type == t_pppParam::TROPO) {
     169        return pp->xx;
     170      }
     171    }
     172    return 0.0;
     173  }
    39174  double trpStdev() const {
    40     const std::vector<t_pppParam*>& par = _parlist->params();
    41     for (unsigned ii = 0; ii < par.size(); ++ii) {
    42       if (par[ii]->type() == t_pppParam::trp) {
     175    for (int ii = 0; ii < _params.size(); ++ii) {
     176      t_pppParam* pp = _params[ii];
     177      if (pp->type == t_pppParam::TROPO) {
    43178        return sqrt(Q()[ii][ii]);
    44179      }
    45180    }
    46181    return 0.0;
    47   };
    48 
     182  }
     183  double Glonass_offset() const {
     184    for (int ii = 0; ii < _params.size(); ++ii) {
     185      t_pppParam* pp = _params[ii];
     186      if (pp->type == t_pppParam::GLONASS_OFFSET) {
     187        return pp->xx;
     188      }
     189    }
     190    return 0.0;
     191  }
     192  double Galileo_offset() const {
     193    for (int ii = 0; ii < _params.size(); ++ii) {
     194      t_pppParam* pp = _params[ii];
     195      if (pp->type == t_pppParam::GALILEO_OFFSET) {
     196        return pp->xx;
     197      }
     198    }
     199    return 0.0;
     200  }
     201  double Bds_offset() const {
     202    for (int ii = 0; ii < _params.size(); ++ii) {
     203      t_pppParam* pp = _params[ii];
     204      if (pp->type == t_pppParam::BDS_OFFSET) {
     205        return pp->xx;
     206      }
     207    }
     208    return 0.0;
     209  }
    49210 private:
    50   class t_slip {
    51    public:
    52     t_slip() {
    53       _slip            = false;
    54       _obsSlipCounter  = -1;
    55       _biasJumpCounter = -1;
    56     }
    57     bool _slip;
    58     int  _obsSlipCounter;
    59     int  _biasJumpCounter;
    60   };
    61 
    62   class t_dop {
    63    public:
    64     t_dop() {reset();}
    65     void reset() {P = T = G = 0.0;}
    66     double P;
    67     double T;
    68     double G;
    69   };
    70 
    71   t_irc processSystem(const std::vector<t_lc::type>& LCs,
    72                       const std::vector<t_pppSatObs*>& obsVector);
    73 
    74   t_irc detectCycleSlips(const std::vector<t_lc::type>& LCs,
    75                          const std::vector<t_pppSatObs*>& obsVector);
    76 
    77   t_irc resetAmb(t_prn prn, const std::vector<t_pppSatObs*>& obsVector,
    78                  SymmetricMatrix* QSav = 0, ColumnVector* xSav = 0);
    79 
    80   void cmpDOP(const std::vector<t_pppSatObs*>& obsVector);
    81 
    82   void predictCovCrdPart(const SymmetricMatrix& QFltOld);
    83 
    84   bncTime         _epoTime;
    85   t_pppParlist*   _parlist;
    86   SymmetricMatrix _QFlt;
    87   ColumnVector    _xFlt;
    88   ColumnVector    _x0;
    89   t_slip          _slips[t_prn::MAXPRN+1];
    90   int             _numSat;
    91   t_dop           _dop;
    92   bncTime         _firstEpoTime;
     211  void   reset();
     212  t_irc  cmpBancroft(t_epoData* epoData);
     213  void   cmpEle(t_satData* satData);
     214  void   addAmb(t_satData* satData);
     215  void   addObs(int iPhase, unsigned& iObs, t_satData* satData,
     216                Matrix& AA, ColumnVector& ll, DiagonalMatrix& PP);
     217  QByteArray printRes(int iPhase, const ColumnVector& vv,
     218                      const QMap<QString, t_satData*>& satDataMap);
     219  void   findMaxRes(const ColumnVector& vv,
     220                    const QMap<QString, t_satData*>& satData,
     221                    QString& prnGPS, QString& prnGlo, 
     222                    double& maxResGPS, double& maxResGlo);
     223  double cmpValue(t_satData* satData, bool phase);
     224  double delay_saast(double Ele);
     225  void   predict(int iPhase, t_epoData* epoData);
     226  t_irc  update_p(t_epoData* epoData);
     227  QString outlierDetection(int iPhase, const ColumnVector& vv,
     228                           QMap<QString, t_satData*>& satData);
     229
     230  double windUp(const QString& prn, const ColumnVector& rSat,
     231                const ColumnVector& rRec);
     232
     233  bncTime  _startTime;
     234
     235  void rememberState(t_epoData* epoData);
     236  void restoreState(t_epoData* epoData);
     237 
     238  t_irc selectSatellites(const QString& lastOutlierPrn,
     239                         QMap<QString, t_satData*>& satData);
     240
     241  void bancroft(const Matrix& BBpass, ColumnVector& pos);
     242
     243  void cmpDOP(t_epoData* epoData);
     244
     245  t_pppClient*          _pppClient;
     246  bncTime               _time;
     247  bncTime               _lastTimeOK;
     248  QVector<t_pppParam*>  _params;
     249  SymmetricMatrix       _QQ;
     250  QVector<t_pppParam*>  _params_sav;
     251  SymmetricMatrix       _QQ_sav;
     252  t_epoData*            _epoData_sav;
     253  ColumnVector          _xcBanc;
     254  ColumnVector          _ellBanc;
     255  QMap<QString, double> _windUpTime;
     256  QMap<QString, double> _windUpSum;
     257  QStringList           _outlierGPS;
     258  QStringList           _outlierGlo;
     259  bncAntex*             _antex;
     260  t_tides*              _tides;
     261  ColumnVector          _neu;
     262  int                   _numSat;
     263  double                _pDop;
    93264};
    94265
  • trunk/BNC/src/src.pri

    r6875 r7203  
    117117}
    118118
    119 exists(PPP) {
    120   INCLUDEPATH += PPP
    121   HEADERS += PPP/pppClient.h    PPP/pppObsPool.h   PPP/pppEphPool.h   \
    122              PPP/pppStation.h   PPP/pppFilter.h    PPP/pppParlist.h   \
    123              PPP/pppSatObs.h           
    124   SOURCES += PPP/pppClient.cpp  PPP/pppObsPool.cpp PPP/pppEphPool.cpp \
    125              PPP/pppStation.cpp PPP/pppFilter.cpp  PPP/pppParlist.cpp \
    126              PPP/pppSatObs.cpp     
     119exists(PPP_AR) {
     120  INCLUDEPATH += PPP_AR
     121  DEFINES += USE_PPP_AR
     122  HEADERS += PPP_AR/pppClient.h    PPP_AR/pppObsPool.h   PPP_AR/pppEphPool.h   \
     123             PPP_AR/pppStation.h   PPP_AR/pppFilter.h    PPP_AR/pppParlist.h   \
     124             PPP_AR/pppSatObs.h           
     125  SOURCES += PPP_AR/pppClient.cpp  PPP_AR/pppObsPool.cpp PPP_AR/pppEphPool.cpp \
     126             PPP_AR/pppStation.cpp PPP_AR/pppFilter.cpp  PPP_AR/pppParlist.cpp \
     127             PPP_AR/pppSatObs.cpp     
    127128}
    128129else {
    129   INCLUDEPATH += PPP_SSR_I
    130   DEFINES += USE_PPP_SSR_I
    131   HEADERS += PPP_SSR_I/pppClient.h   PPP_SSR_I/pppFilter.h
    132   SOURCES += PPP_SSR_I/pppClient.cpp PPP_SSR_I/pppFilter.cpp
     130  INCLUDEPATH += PPP
     131  DEFINES += USE_PPP
     132  HEADERS += PPP/pppClient.h   PPP/pppFilter.h
     133  SOURCES += PPP/pppClient.cpp PPP/pppFilter.cpp
    133134}
    134135
Note: See TracChangeset for help on using the changeset viewer.