Changeset 7996 in ntrip for trunk/BNC/src


Ignore:
Timestamp:
Aug 4, 2016, 9:51:01 AM (8 years ago)
Author:
stuerze
Message:

reading of ocean loading files is added

Location:
trunk/BNC/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/BNC/src/PPP/pppClient.cpp

    r7972 r7996  
    6363    _antex = 0;
    6464  }
    65 /*
     65
    6666  if (!_opt->_blqFileName.empty()) {
    6767    _loading = new t_loading(_opt->_blqFileName.c_str());
    … …  
    7070    _loading = 0;
    7171  }
    72 */
     72
    7373  CLIENTS.setLocalData(this);  // CLIENTS takes ownership over "this"
    7474}
    … …  
    8787  delete _filter;
    8888  delete _tides;
    89   /*
    9089  if (_loading) {
    9190    delete _loading;
    9291  }
    93   */
    9492  clearObs();
    9593}

  • trunk/BNC/src/PPP/pppClient.h

    r7972 r7996  
    7171  t_pppOptions*             _opt;
    7272  t_tides*                  _tides;
    73 //  t_loading*                _loading;
     73  t_loading*                _loading;
    7474};
    7575

  • trunk/BNC/src/pppModel.cpp

    r7259 r7996  
    3636 * Created:    29-Jul-2014
    3737 *
    38  * Changes:   
     38 * Changes:
    3939 *
    4040 * -----------------------------------------------------------------------*/
    … …  
    9090  double Mjd_0 = floor(Mjd_UT1);
    9191  double UT1   = Secs*(Mjd_UT1-Mjd_0);
    92   double T_0   = (Mjd_0  -MJD_J2000)/36525.0; 
    93   double T     = (Mjd_UT1-MJD_J2000)/36525.0; 
     92  double T_0   = (Mjd_0  -MJD_J2000)/36525.0;
     93  double T     = (Mjd_UT1-MJD_J2000)/36525.0;
    9494
    9595  double gmst  = 24110.54841 + 8640184.812866*T_0 + 1.002737909350795*UT1
    … …  
    126126  const double T  = (Mjd_1-MJD_J2000)/36525.0;
    127127  const double dT = (Mjd_2-Mjd_1)/36525.0;
    128  
     128
    129129  double zeta  =  ( (2306.2181+(1.39656-0.000139*T)*T)+
    130130                        ((0.30188-0.000344*T)+0.017998*dT)*dT )*dT/RHO_SEC;
    … …  
    134134
    135135  return rotZ(-z) * rotY(theta) * rotZ(-zeta);
    136 }   
     136}
    137137
    138138// Sun's position
    … …  
    144144
    145145  double M = 2.0*M_PI * Frac ( 0.9931267 + 99.9973583*T);
    146   double L = 2.0*M_PI * Frac ( 0.7859444 + M/2.0/M_PI + 
     146  double L = 2.0*M_PI * Frac ( 0.7859444 + M/2.0/M_PI +
    147147                        (6892.0*sin(M)+72.0*sin(2.0*M)) / 1296.0e3);
    148148  double r = 149.619e9 - 2.499e9*cos(M) - 0.021e9*cos(2*M);
    149  
    150   ColumnVector r_Sun(3); 
     149
     150  ColumnVector r_Sun(3);
    151151  r_Sun << r*cos(L) << r*sin(L) << 0.0; r_Sun = rotX(-eps) * r_Sun;
    152152
    153153  return    rotZ(GMST(Mjd_TT))
    154           * NutMatrix(Mjd_TT) 
     154          * NutMatrix(Mjd_TT)
    155155          * PrecMatrix(MJD_J2000, Mjd_TT)
    156156          * r_Sun;
    … …  
    169169  double D   = 2.0*M_PI*Frac ( 0.827361 + 1236.853086*T );
    170170  double F   = 2.0*M_PI*Frac ( 0.259086 + 1342.227825*T );
    171    
    172   double dL = +22640*sin(l) - 4586*sin(l-2*D) + 2370*sin(2*D) +  769*sin(2*l) 
     171
     172  double dL = +22640*sin(l) - 4586*sin(l-2*D) + 2370*sin(2*D) +  769*sin(2*l)
    173173              -668*sin(lp) - 412*sin(2*F) - 212*sin(2*l-2*D)- 206*sin(l+lp-2*D)
    174174              +192*sin(l+2*D) - 165*sin(lp-2*D) - 125*sin(D) - 110*sin(l+lp)
    … …  
    177177  double L = 2.0*M_PI * Frac( L_0 + dL/1296.0e3 );
    178178
    179   double S  = F + (dL+412*sin(2*F)+541*sin(lp)) / RHO_SEC; 
     179  double S  = F + (dL+412*sin(2*F)+541*sin(lp)) / RHO_SEC;
    180180  double h  = F-2*D;
    181   double N  = -526*sin(h) + 44*sin(l+h) - 31*sin(-l+h) - 23*sin(lp+h) 
     181  double N  = -526*sin(h) + 44*sin(l+h) - 31*sin(-l+h) - 23*sin(lp+h)
    182182              +11*sin(-lp+h) - 25*sin(-2*l+F) + 21*sin(-l+F);
    183183
    184184  double B = ( 18520.0*sin(S) + N ) / RHO_SEC;
    185    
     185
    186186  double cosB = cos(B);
    187187
    188188  double R = 385000e3 - 20905e3*cos(l) - 3699e3*cos(2*D-l) - 2956e3*cos(2*D)
    189       -570e3*cos(2*l) + 246e3*cos(2*l-2*D) - 205e3*cos(lp-2*D) 
    190       -171e3*cos(l+2*D) - 152e3*cos(l+lp-2*D);   
    191 
    192   ColumnVector r_Moon(3); 
     189      -570e3*cos(2*l) + 246e3*cos(2*l-2*D) - 205e3*cos(lp-2*D)
     190      -171e3*cos(l+2*D) - 152e3*cos(l+lp-2*D);
     191
     192  ColumnVector r_Moon(3);
    193193  r_Moon << R*cos(L)*cosB << R*sin(L)*cosB << R*sin(B);
    194194  r_Moon = rotX(-eps) * r_Moon;
    195    
    196   return    rotZ(GMST(Mjd_TT)) 
    197           * NutMatrix(Mjd_TT) 
     195
     196  return    rotZ(GMST(Mjd_TT))
     197          * NutMatrix(Mjd_TT)
    198198          * PrecMatrix(MJD_J2000, Mjd_TT)
    199199          * r_Moon;
    200200}
    201201
    202 // Tidal Correction 
     202// Tidal Correction
    203203////////////////////////////////////////////////////////////////////////////
    204204ColumnVector t_tides::displacement(const bncTime& time, const ColumnVector& xyz) {
    … …  
    239239  double x2Sun  = 3.0 * L2 * scSun;
    240240  double x2Moon = 3.0 * L2 * scMoon;
    241  
     241
    242242  const double gmWGS = 398.6005e12;
    243243  const double gms   = 1.3271250e20;
    244244  const double gmm   = 4.9027890e12;
    245245
    246   double facSun  = gms / gmWGS * 
     246  double facSun  = gms / gmWGS *
    247247                   (rRec * rRec * rRec * rRec) / (_rSun * _rSun * _rSun);
    248248
    249   double facMoon = gmm / gmWGS * 
     249  double facMoon = gmm / gmWGS *
    250250                   (rRec * rRec * rRec * rRec) / (_rMoon * _rMoon * _rMoon);
    251251
    252   ColumnVector dX = facSun  * (x2Sun  * _xSun  + p2Sun  * xyzUnit) + 
     252  ColumnVector dX = facSun  * (x2Sun  * _xSun  + p2Sun  * xyzUnit) +
    253253                    facMoon * (x2Moon * _xMoon + p2Moon * xyzUnit);
    254254
    255255  return dX;
     256}
     257
     258// Constructor
     259///////////////////////////////////////////////////////////////////////////
     260t_loading::t_loading(const QString& fileName) {
     261
     262  newBlqData = 0;
     263  readFile(fileName);
     264  printAll();
     265}
     266
     267t_loading::~t_loading() {
     268
     269  if (newBlqData) {
     270    delete newBlqData;
     271  }
     272
     273  QMapIterator<QString, t_blqData*> it(blqMap);
     274  while (it.hasNext()) {
     275    it.next();
     276    delete it.value();
     277  }
     278}
     279
     280t_irc t_loading::readFile(const QString& fileName) {
     281  QFile inFile(fileName);
     282  inFile.open(QIODevice::ReadOnly | QIODevice::Text);
     283  QTextStream in(&inFile);
     284  int row = 0;
     285  QString site = QString();
     286
     287  while ( !in.atEnd() ) {
     288
     289    QString line = in.readLine();
     290
     291    // skip empty lines and comments
     292    if (line.indexOf("$$") != -1 || line.size() == 0) {
     293      continue;
     294    }
     295
     296    line = line.trimmed();
     297    QTextStream inLine(line.toAscii(), QIODevice::ReadOnly);
     298
     299    switch (row) {
     300      case 0:
     301        site = line;
     302        site = site.toUpper();
     303        if (!newBlqData) {
     304          newBlqData = new t_blqData;
     305          newBlqData->amplitudes.ReSize(3,11);
     306          newBlqData->phases.ReSize(3,11);
     307        }
     308        break;
     309      case 1: case 2: case 3:
     310        for (int ii = 0; ii < 11; ii++) {
     311          inLine >> newBlqData->amplitudes[row-1][ii];
     312        }
     313        break;
     314      case 4: case 5: case 6:
     315        for (int ii = 0; ii < 11; ii++) {
     316          inLine >> newBlqData->phases[row-4][ii];
     317        }
     318        break;
     319      case 7:
     320        if (newBlqData && !site.isEmpty()) {
     321          blqMap[site] = newBlqData;
     322          site = QString();
     323          row = -1;
     324          newBlqData = 0;
     325        }
     326        break;
     327    }
     328    row++;
     329  }
     330  inFile.close();
     331  return success;
     332}
     333
     334// Print
     335////////////////////////////////////////////////////////////////////////////
     336void t_loading::printAll() const {
     337
     338  QMapIterator<QString, t_blqData*> it(blqMap);
     339  while (it.hasNext()) {
     340    it.next();
     341    t_blqData* blq = it.value();
     342    QString site = it.key();
     343    cout << site.toStdString().c_str() << "\n";
     344    for (int rr = 0; rr < 3; rr++) {
     345      for (int cc = 0; cc < 11; cc++) {
     346        cout << blq->amplitudes[rr][cc] << " ";
     347      }
     348      cout << endl;
     349    }
     350    for (int rr = 0; rr < 3; rr++) {
     351      for (int cc = 0; cc < 11; cc++) {
     352        cout << blq->phases[rr][cc] << " ";
     353      }
     354      cout << endl;
     355    }
     356  }
    256357}
    257358
    … …  
    267368// Phase Wind-Up Correction
    268369///////////////////////////////////////////////////////////////////////////
    269 double t_windUp::value(const bncTime& etime, const ColumnVector& rRec, 
     370double t_windUp::value(const bncTime& etime, const ColumnVector& rRec,
    270371                       t_prn prn, const ColumnVector& rSat) {
    271372
    … …  
    276377    ColumnVector rho = rRec - rSat;
    277378    rho /= rho.norm_Frobenius();
    278    
     379
    279380    // GPS Satellite unit Vectors sz, sy, sx
    280381    // -------------------------------------
    … …  
    289390    // Effective Dipole of the GPS Satellite Antenna
    290391    // ---------------------------------------------
    291     ColumnVector dipSat = sx - rho * DotProduct(rho,sx) 
     392    ColumnVector dipSat = sx - rho * DotProduct(rho,sx)
    292393                                                - crossproduct(rho, sy);
    293    
     394
    294395    // Receiver unit Vectors rx, ry
    295396    // ----------------------------
    … …  
    299400    double recEll[3]; xyz2ell(rRec.data(), recEll) ;
    300401    double neu[3];
    301    
     402
    302403    neu[0] = 1.0;
    303404    neu[1] = 0.0;
    304405    neu[2] = 0.0;
    305406    neu2xyz(recEll, neu, rx.data());
    306    
     407
    307408    neu[0] =  0.0;
    308409    neu[1] = -1.0;
    309410    neu[2] =  0.0;
    310411    neu2xyz(recEll, neu, ry.data());
    311    
     412
    312413    // Effective Dipole of the Receiver Antenna
    313414    // ----------------------------------------
    314     ColumnVector dipRec = rx - rho * DotProduct(rho,rx) 
     415    ColumnVector dipRec = rx - rho * DotProduct(rho,rx)
    315416                                                   + crossproduct(rho, ry);
    316    
     417
    317418    // Resulting Effect
    318419    // ----------------
    319     double alpha = DotProduct(dipSat,dipRec) / 
     420    double alpha = DotProduct(dipSat,dipRec) /
    320421                      (dipSat.norm_Frobenius() * dipRec.norm_Frobenius());
    321    
     422
    322423    if (alpha >  1.0) alpha =  1.0;
    323424    if (alpha < -1.0) alpha = -1.0;
    324    
     425
    325426    double dphi = acos(alpha) / 2.0 / M_PI;  // in cycles
    326    
     427
    327428    if ( DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0 ) {
    328429      dphi = -dphi;
    … …  
    339440  }
    340441
    341   return sumWind[prn.toInt()]; 
     442  return sumWind[prn.toInt()];
    342443}
    343444
    … …  
    352453  }
    353454
    354   double ell[3]; 
     455  double ell[3];
    355456  xyz2ell(xyz.data(), ell);
    356457  double height = ell[2];
    … …  
    362463
    363464  double h_km = height / 1000.0;
    364  
     465
    365466  if (h_km < 0.0) h_km = 0.0;
    366467  if (h_km > 5.0) h_km = 5.0;
    367468  int    ii   = int(h_km + 1); if (ii > 5) ii = 5;
    368469  double href = ii - 1;
    369  
    370   double bCor[6]; 
     470
     471  double bCor[6];
    371472  bCor[0] = 1.156;
    372473  bCor[1] = 1.006;
    … …  
    375476  bCor[4] = 0.654;
    376477  bCor[5] = 0.563;
    377  
     478
    378479  double BB = bCor[ii-1] + (bCor[ii]-bCor[ii-1]) * (h_km - href);
    379  
     480
    380481  double zen  = M_PI/2.0 - Ele;
    381482

  • trunk/BNC/src/pppModel.h

    r7625 r7996  
    4646};
    4747
     48class t_loading {
     49 public:
     50  t_loading(const QString& fileName);
     51  ~t_loading();
     52  t_irc   readFile(const QString& fileName);
     53  void    printAll() const;
     54
     55 private:
     56  class t_blqData {
     57   public:
     58    t_blqData() {}
     59    Matrix amplitudes;
     60    Matrix phases;
     61  };
     62  t_blqData*                 newBlqData;
     63  QMap <QString, t_blqData*> blqMap;
     64
     65};
     66
    4867class t_windUp {
    4968 public:
Note: See TracChangeset for help on using the changeset viewer.