1  #ifndef FILTER_H


2  #define FILTER_H


3 


4  #include <vector>


5  #include <newmat.h>


6  #include "pppInclude.h"


7  #include "pppParlist.h"


8  #include "bnctime.h"


9  #include "t_prn.h"


10  #include "pppClient.h"


11 


12  namespace BNC_PPP {


13 


14  class t_pppParlist;


15  class t_pppObsPool;


16  class t_pppSatObs;


17 


18  class t_pppFilter {


19  public:


20  t_pppFilter(t_pppObsPool* obsPool);


21  ~t_pppFilter();


22 


23  t_irc processEpoch();


24 


25  t_irc datumTransformation();


26  void initDatumTransformation(const std::vector<t_pppSatObs*>& allObs, bool pseudoObsIono);


27  unsigned setTrafoObs();


28 


29  const ColumnVector& x() const {return _xFlt;}


30  const SymmetricMatrix& Q() const {return _QFlt;}


31 


32  int numSat() const {return _numSat;}


33  double HDOP() const {return _dop.H;}


34  double HDOV() const {return _dop.V;}


35  double PDOP() const {return _dop.P;}


36  double GDOP() const {return _dop.G;}


37  double trp() const {


38  const std::vector<t_pppParam*>& par = _parlist>params();


39  for (unsigned ii = 0; ii < par.size(); ++ii) {


40  if (par[ii]>type() == t_pppParam::trp) {


41  return x()[ii];


42  }


43  }


44  return 0.0;


45  };


46  double trpStdev() const {


47  const std::vector<t_pppParam*>& par = _parlist>params();


48  for (unsigned ii = 0; ii < par.size(); ++ii) {


49  if (par[ii]>type() == t_pppParam::trp) {


50  return sqrt(Q()[ii][ii]);


51  }


52  }


53  return 0.0;


54  };


55 


56  private:


57  class t_slip {


58  public:


59  t_slip() {


60  _slip = false;


61  _obsSlipCounter = 1;


62  _biasJumpCounter = 1;


63  }


64  bool _slip;


65  int _obsSlipCounter;


66  int _biasJumpCounter;


67  };


68 


69  class t_dop {


70  public:


71  t_dop() {reset();}


72  void reset() {H = V = P = T = G = 0.0;}


73  double H;


74  double V;


75  double P;


76  double T;


77  double G;


78  };


79 


80  class t_datumTrafo {


81  public:


82  t_datumTrafo () {


83  initIndices();


84  }


85  ~t_datumTrafo () {


86 


87  }


88 


89  void initIndices() {


90  _firstRow = 1; _lastRow = 0;


91  }


92 


93  void setFirstSystem(char firstSys) {


94  _firstSys = firstSys;


95  }


96 


97  bool firstSystem(char sys) {


98  if (_firstSys == sys) {


99  return true;


100  }


101  return false;


102  }


103 


104  void updateIndices(char sys, int maxObsSys) {


105  if (firstSystem(sys)) {


106  initIndices();


107  }


108  else {


109  _firstRow = _lastRow + 1;


110  }


111  _lastRow += maxObsSys;//LOG << "updateIndices: lastRow: " << _lastRow << endl;


112  };


113 


114  void setNumObs(int maxObs) {_maxObs = maxObs;}


115  void setNumPar(int numPar) { _numPar = numPar;}


116  void setLastEpoParlist(t_pppParlist* parlist) {_parlist = parlist;}


117  t_pppParlist* lastEpoParlist() {return _parlist;}


118  int numPar() {return _numPar;}


119  int numObs() {return _maxObs;}


120  void updateNumObs() {//LOG << "updateObsNum _maxObs " << _maxObs;


121  _maxObs = _lastRow;//LOG << "=> _maxObs " << _maxObs << " _numPar: " << _numPar << endl;


122  _AA1 = _AA1.SubMatrix(1, _lastRow, 1, _numPar);


123  _AA2 = _AA2.SubMatrix(1, _lastRow, 1, _numPar);


124  }


125 


126  const Matrix& AA1() {return _AA1;}


127  const Matrix& AA2() {return _AA2;}


128  const Matrix& D21() {return _D21;}


129 


130  void initAA() {//LOG << "initAA: _maxObs: " << _maxObs << " _numPar: " << _numPar << endl;


131  _AA1.ReSize(_maxObs, _numPar); _AA1 = 0.0;


132  _AA2.ReSize(_maxObs, _numPar); _AA2 = 0.0;


133  _D21.ReSize(_numPar, _numPar); _D21 = 0.0;


134  }


135 


136  void prepareAA(const Matrix& AA, int ind) {


137  Matrix* Prep = &_AA2;


138  if (ind == 1) {


139  Prep = &_AA1;


140  }


141  Prep>SubMatrix(_firstRow, _lastRow, 1, _numPar) << AA;


142  }


143 


144  void switchAA() {


145  _AA1 = _AA2;


146  }


147 


148  t_irc computeTrafoMatrix() {


149  if (((_AA2.t() * _AA2)).Determinant() == 0.0) {


150  LOG << "(_AA2.t() * _AA2).inv() is singular" << endl;


151  return failure;


152  }


153  _D21 = ((_AA2.t() * _AA2).i()) * _AA2.t() * _AA1;


154  return success;


155  }


156 


157  void printMatrix(const Matrix& X, int nRow, int nCol) {


158  for (int rr = 0; rr < nRow; rr++) {


159  for (int cc = 0; cc < nCol; cc++) {


160  LOG << setw(6) << setprecision(3) << X[rr][cc] << " ;";


161  }


162  LOG << endl;


163  }


164  LOG << endl;


165  }


166  private:


167  int _firstRow;


168  int _lastRow;


169  Matrix _AA1;


170  Matrix _AA2;


171  Matrix _D21;


172  char _firstSys;


173  int _maxObs;


174  int _numPar;


175  QMap<char, t_prn> _refSatMapPseudoObs;


176  t_pppParlist* _parlist;


177  };


178 


179  t_irc processSystem(const std::vector<t_lc::type>& LCs,


180  const std::vector<t_pppSatObs*>& obsVector,


181  const t_prn& refPrn,


182  bool pseudoObsIonoAvailable,


183  bool preProcessing);


184 


185  t_irc detectCycleSlips(const std::vector<t_lc::type>& LCs,


186  const std::vector<t_pppSatObs*>& obsVector,


187  const t_prn& refPrn,


188  bool preProcessing);


189 


190  t_irc resetAmb(t_prn prn, const std::vector<t_pppSatObs*>& obsVector,


191  SymmetricMatrix* QSav = 0, ColumnVector* xSav = 0);


192 


193  void cmpDOP(const std::vector<t_pppSatObs*>& obsVector);


194 


195  void predictCovCrdPart(const SymmetricMatrix& QFltOld);


196 


197  t_irc addNoiseToIono(char sys);


198 


199  bool resetRefSatellitesLastEpoch(std::vector<t_pppSatObs*>& obsVector);


200 


201  bncTime _epoTime;


202  t_pppParlist* _parlist;


203  t_pppObsPool* _obsPool;


204  t_datumTrafo* _datumTrafo;


205  SymmetricMatrix _QFlt;


206  ColumnVector _xFlt;


207  ColumnVector _x0;


208  t_slip _slips[t_prn::MAXPRN+1];


209  int _numSat;


210  t_dop _dop;


211  bncTime _firstEpoTime;


212  bncTime _lastEpoTimeOK;


213  t_prn _refPrn;


214  };


215 


216  }


217 


218  #endif

