Index: /trunk/BNC/src/PPP/pppClient.cpp =================================================================== --- /trunk/BNC/src/PPP/pppClient.cpp (revision 7995) +++ /trunk/BNC/src/PPP/pppClient.cpp (revision 7996) @@ -63,5 +63,5 @@ _antex = 0; } -/* + if (!_opt->_blqFileName.empty()) { _loading = new t_loading(_opt->_blqFileName.c_str()); @@ -70,5 +70,5 @@ _loading = 0; } -*/ + CLIENTS.setLocalData(this); // CLIENTS takes ownership over "this" } @@ -87,9 +87,7 @@ delete _filter; delete _tides; - /* if (_loading) { delete _loading; } - */ clearObs(); } Index: /trunk/BNC/src/PPP/pppClient.h =================================================================== --- /trunk/BNC/src/PPP/pppClient.h (revision 7995) +++ /trunk/BNC/src/PPP/pppClient.h (revision 7996) @@ -71,5 +71,5 @@ t_pppOptions* _opt; t_tides* _tides; -// t_loading* _loading; + t_loading* _loading; }; Index: /trunk/BNC/src/pppModel.cpp =================================================================== --- /trunk/BNC/src/pppModel.cpp (revision 7995) +++ /trunk/BNC/src/pppModel.cpp (revision 7996) @@ -36,5 +36,5 @@ * Created: 29-Jul-2014 * - * Changes: + * Changes: * * -----------------------------------------------------------------------*/ @@ -90,6 +90,6 @@ double Mjd_0 = floor(Mjd_UT1); double UT1 = Secs*(Mjd_UT1-Mjd_0); - double T_0 = (Mjd_0 -MJD_J2000)/36525.0; - double T = (Mjd_UT1-MJD_J2000)/36525.0; + double T_0 = (Mjd_0 -MJD_J2000)/36525.0; + double T = (Mjd_UT1-MJD_J2000)/36525.0; double gmst = 24110.54841 + 8640184.812866*T_0 + 1.002737909350795*UT1 @@ -126,5 +126,5 @@ const double T = (Mjd_1-MJD_J2000)/36525.0; const double dT = (Mjd_2-Mjd_1)/36525.0; - + double zeta = ( (2306.2181+(1.39656-0.000139*T)*T)+ ((0.30188-0.000344*T)+0.017998*dT)*dT )*dT/RHO_SEC; @@ -134,5 +134,5 @@ return rotZ(-z) * rotY(theta) * rotZ(-zeta); -} +} // Sun's position @@ -144,13 +144,13 @@ double M = 2.0*M_PI * Frac ( 0.9931267 + 99.9973583*T); - double L = 2.0*M_PI * Frac ( 0.7859444 + M/2.0/M_PI + + double L = 2.0*M_PI * Frac ( 0.7859444 + M/2.0/M_PI + (6892.0*sin(M)+72.0*sin(2.0*M)) / 1296.0e3); double r = 149.619e9 - 2.499e9*cos(M) - 0.021e9*cos(2*M); - - ColumnVector r_Sun(3); + + ColumnVector r_Sun(3); r_Sun << r*cos(L) << r*sin(L) << 0.0; r_Sun = rotX(-eps) * r_Sun; return rotZ(GMST(Mjd_TT)) - * NutMatrix(Mjd_TT) + * NutMatrix(Mjd_TT) * PrecMatrix(MJD_J2000, Mjd_TT) * r_Sun; @@ -169,6 +169,6 @@ double D = 2.0*M_PI*Frac ( 0.827361 + 1236.853086*T ); double F = 2.0*M_PI*Frac ( 0.259086 + 1342.227825*T ); - - double dL = +22640*sin(l) - 4586*sin(l-2*D) + 2370*sin(2*D) + 769*sin(2*l) + + double dL = +22640*sin(l) - 4586*sin(l-2*D) + 2370*sin(2*D) + 769*sin(2*l) -668*sin(lp) - 412*sin(2*F) - 212*sin(2*l-2*D)- 206*sin(l+lp-2*D) +192*sin(l+2*D) - 165*sin(lp-2*D) - 125*sin(D) - 110*sin(l+lp) @@ -177,28 +177,28 @@ double L = 2.0*M_PI * Frac( L_0 + dL/1296.0e3 ); - double S = F + (dL+412*sin(2*F)+541*sin(lp)) / RHO_SEC; + double S = F + (dL+412*sin(2*F)+541*sin(lp)) / RHO_SEC; double h = F-2*D; - double N = -526*sin(h) + 44*sin(l+h) - 31*sin(-l+h) - 23*sin(lp+h) + double N = -526*sin(h) + 44*sin(l+h) - 31*sin(-l+h) - 23*sin(lp+h) +11*sin(-lp+h) - 25*sin(-2*l+F) + 21*sin(-l+F); double B = ( 18520.0*sin(S) + N ) / RHO_SEC; - + double cosB = cos(B); double R = 385000e3 - 20905e3*cos(l) - 3699e3*cos(2*D-l) - 2956e3*cos(2*D) - -570e3*cos(2*l) + 246e3*cos(2*l-2*D) - 205e3*cos(lp-2*D) - -171e3*cos(l+2*D) - 152e3*cos(l+lp-2*D); - - ColumnVector r_Moon(3); + -570e3*cos(2*l) + 246e3*cos(2*l-2*D) - 205e3*cos(lp-2*D) + -171e3*cos(l+2*D) - 152e3*cos(l+lp-2*D); + + ColumnVector r_Moon(3); r_Moon << R*cos(L)*cosB << R*sin(L)*cosB << R*sin(B); r_Moon = rotX(-eps) * r_Moon; - - return rotZ(GMST(Mjd_TT)) - * NutMatrix(Mjd_TT) + + return rotZ(GMST(Mjd_TT)) + * NutMatrix(Mjd_TT) * PrecMatrix(MJD_J2000, Mjd_TT) * r_Moon; } -// Tidal Correction +// Tidal Correction //////////////////////////////////////////////////////////////////////////// ColumnVector t_tides::displacement(const bncTime& time, const ColumnVector& xyz) { @@ -239,19 +239,120 @@ double x2Sun = 3.0 * L2 * scSun; double x2Moon = 3.0 * L2 * scMoon; - + const double gmWGS = 398.6005e12; const double gms = 1.3271250e20; const double gmm = 4.9027890e12; - double facSun = gms / gmWGS * + double facSun = gms / gmWGS * (rRec * rRec * rRec * rRec) / (_rSun * _rSun * _rSun); - double facMoon = gmm / gmWGS * + double facMoon = gmm / gmWGS * (rRec * rRec * rRec * rRec) / (_rMoon * _rMoon * _rMoon); - ColumnVector dX = facSun * (x2Sun * _xSun + p2Sun * xyzUnit) + + ColumnVector dX = facSun * (x2Sun * _xSun + p2Sun * xyzUnit) + facMoon * (x2Moon * _xMoon + p2Moon * xyzUnit); return dX; +} + +// Constructor +/////////////////////////////////////////////////////////////////////////// +t_loading::t_loading(const QString& fileName) { + + newBlqData = 0; + readFile(fileName); + printAll(); +} + +t_loading::~t_loading() { + + if (newBlqData) { + delete newBlqData; + } + + QMapIterator it(blqMap); + while (it.hasNext()) { + it.next(); + delete it.value(); + } +} + +t_irc t_loading::readFile(const QString& fileName) { + QFile inFile(fileName); + inFile.open(QIODevice::ReadOnly | QIODevice::Text); + QTextStream in(&inFile); + int row = 0; + QString site = QString(); + + while ( !in.atEnd() ) { + + QString line = in.readLine(); + + // skip empty lines and comments + if (line.indexOf("$$") != -1 || line.size() == 0) { + continue; + } + + line = line.trimmed(); + QTextStream inLine(line.toAscii(), QIODevice::ReadOnly); + + switch (row) { + case 0: + site = line; + site = site.toUpper(); + if (!newBlqData) { + newBlqData = new t_blqData; + newBlqData->amplitudes.ReSize(3,11); + newBlqData->phases.ReSize(3,11); + } + break; + case 1: case 2: case 3: + for (int ii = 0; ii < 11; ii++) { + inLine >> newBlqData->amplitudes[row-1][ii]; + } + break; + case 4: case 5: case 6: + for (int ii = 0; ii < 11; ii++) { + inLine >> newBlqData->phases[row-4][ii]; + } + break; + case 7: + if (newBlqData && !site.isEmpty()) { + blqMap[site] = newBlqData; + site = QString(); + row = -1; + newBlqData = 0; + } + break; + } + row++; + } + inFile.close(); + return success; +} + +// Print +//////////////////////////////////////////////////////////////////////////// +void t_loading::printAll() const { + + QMapIterator it(blqMap); + while (it.hasNext()) { + it.next(); + t_blqData* blq = it.value(); + QString site = it.key(); + cout << site.toStdString().c_str() << "\n"; + for (int rr = 0; rr < 3; rr++) { + for (int cc = 0; cc < 11; cc++) { + cout << blq->amplitudes[rr][cc] << " "; + } + cout << endl; + } + for (int rr = 0; rr < 3; rr++) { + for (int cc = 0; cc < 11; cc++) { + cout << blq->phases[rr][cc] << " "; + } + cout << endl; + } + } } @@ -267,5 +368,5 @@ // Phase Wind-Up Correction /////////////////////////////////////////////////////////////////////////// -double t_windUp::value(const bncTime& etime, const ColumnVector& rRec, +double t_windUp::value(const bncTime& etime, const ColumnVector& rRec, t_prn prn, const ColumnVector& rSat) { @@ -276,5 +377,5 @@ ColumnVector rho = rRec - rSat; rho /= rho.norm_Frobenius(); - + // GPS Satellite unit Vectors sz, sy, sx // ------------------------------------- @@ -289,7 +390,7 @@ // Effective Dipole of the GPS Satellite Antenna // --------------------------------------------- - ColumnVector dipSat = sx - rho * DotProduct(rho,sx) + ColumnVector dipSat = sx - rho * DotProduct(rho,sx) - crossproduct(rho, sy); - + // Receiver unit Vectors rx, ry // ---------------------------- @@ -299,30 +400,30 @@ double recEll[3]; xyz2ell(rRec.data(), recEll) ; double neu[3]; - + neu[0] = 1.0; neu[1] = 0.0; neu[2] = 0.0; neu2xyz(recEll, neu, rx.data()); - + neu[0] = 0.0; neu[1] = -1.0; neu[2] = 0.0; neu2xyz(recEll, neu, ry.data()); - + // Effective Dipole of the Receiver Antenna // ---------------------------------------- - ColumnVector dipRec = rx - rho * DotProduct(rho,rx) + ColumnVector dipRec = rx - rho * DotProduct(rho,rx) + crossproduct(rho, ry); - + // Resulting Effect // ---------------- - double alpha = DotProduct(dipSat,dipRec) / + double alpha = DotProduct(dipSat,dipRec) / (dipSat.norm_Frobenius() * dipRec.norm_Frobenius()); - + if (alpha > 1.0) alpha = 1.0; if (alpha < -1.0) alpha = -1.0; - + double dphi = acos(alpha) / 2.0 / M_PI; // in cycles - + if ( DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0 ) { dphi = -dphi; @@ -339,5 +440,5 @@ } - return sumWind[prn.toInt()]; + return sumWind[prn.toInt()]; } @@ -352,5 +453,5 @@ } - double ell[3]; + double ell[3]; xyz2ell(xyz.data(), ell); double height = ell[2]; @@ -362,11 +463,11 @@ double h_km = height / 1000.0; - + if (h_km < 0.0) h_km = 0.0; if (h_km > 5.0) h_km = 5.0; int ii = int(h_km + 1); if (ii > 5) ii = 5; double href = ii - 1; - - double bCor[6]; + + double bCor[6]; bCor[0] = 1.156; bCor[1] = 1.006; @@ -375,7 +476,7 @@ bCor[4] = 0.654; bCor[5] = 0.563; - + double BB = bCor[ii-1] + (bCor[ii]-bCor[ii-1]) * (h_km - href); - + double zen = M_PI/2.0 - Ele; Index: /trunk/BNC/src/pppModel.h =================================================================== --- /trunk/BNC/src/pppModel.h (revision 7995) +++ /trunk/BNC/src/pppModel.h (revision 7996) @@ -46,4 +46,23 @@ }; +class t_loading { + public: + t_loading(const QString& fileName); + ~t_loading(); + t_irc readFile(const QString& fileName); + void printAll() const; + + private: + class t_blqData { + public: + t_blqData() {} + Matrix amplitudes; + Matrix phases; + }; + t_blqData* newBlqData; + QMap blqMap; + +}; + class t_windUp { public: