 Timestamp:
 Aug 4, 2016, 9:51:01 AM (3 years ago)
 Location:
 trunk/BNC/src
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/BNC/src/pppModel.cpp
r7259 r7996 36 36 * Created: 29Jul2014 37 37 * 38 * Changes: 38 * Changes: 39 39 * 40 40 * */ … … 90 90 double Mjd_0 = floor(Mjd_UT1); 91 91 double UT1 = Secs*(Mjd_UT1Mjd_0); 92 double T_0 = (Mjd_0 MJD_J2000)/36525.0; 93 double T = (Mjd_UT1MJD_J2000)/36525.0; 92 double T_0 = (Mjd_0 MJD_J2000)/36525.0; 93 double T = (Mjd_UT1MJD_J2000)/36525.0; 94 94 95 95 double gmst = 24110.54841 + 8640184.812866*T_0 + 1.002737909350795*UT1 … … 126 126 const double T = (Mjd_1MJD_J2000)/36525.0; 127 127 const double dT = (Mjd_2Mjd_1)/36525.0; 128 128 129 129 double zeta = ( (2306.2181+(1.396560.000139*T)*T)+ 130 130 ((0.301880.000344*T)+0.017998*dT)*dT )*dT/RHO_SEC; … … 134 134 135 135 return rotZ(z) * rotY(theta) * rotZ(zeta); 136 } 136 } 137 137 138 138 // Sun's position … … 144 144 145 145 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 + 147 147 (6892.0*sin(M)+72.0*sin(2.0*M)) / 1296.0e3); 148 148 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); 151 151 r_Sun << r*cos(L) << r*sin(L) << 0.0; r_Sun = rotX(eps) * r_Sun; 152 152 153 153 return rotZ(GMST(Mjd_TT)) 154 * NutMatrix(Mjd_TT) 154 * NutMatrix(Mjd_TT) 155 155 * PrecMatrix(MJD_J2000, Mjd_TT) 156 156 * r_Sun; … … 169 169 double D = 2.0*M_PI*Frac ( 0.827361 + 1236.853086*T ); 170 170 double F = 2.0*M_PI*Frac ( 0.259086 + 1342.227825*T ); 171 172 double dL = +22640*sin(l)  4586*sin(l2*D) + 2370*sin(2*D) + 769*sin(2*l) 171 172 double dL = +22640*sin(l)  4586*sin(l2*D) + 2370*sin(2*D) + 769*sin(2*l) 173 173 668*sin(lp)  412*sin(2*F)  212*sin(2*l2*D) 206*sin(l+lp2*D) 174 174 +192*sin(l+2*D)  165*sin(lp2*D)  125*sin(D)  110*sin(l+lp) … … 177 177 double L = 2.0*M_PI * Frac( L_0 + dL/1296.0e3 ); 178 178 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; 180 180 double h = F2*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) 182 182 +11*sin(lp+h)  25*sin(2*l+F) + 21*sin(l+F); 183 183 184 184 double B = ( 18520.0*sin(S) + N ) / RHO_SEC; 185 185 186 186 double cosB = cos(B); 187 187 188 188 double R = 385000e3  20905e3*cos(l)  3699e3*cos(2*Dl)  2956e3*cos(2*D) 189 570e3*cos(2*l) + 246e3*cos(2*l2*D)  205e3*cos(lp2*D) 190 171e3*cos(l+2*D)  152e3*cos(l+lp2*D); 191 192 ColumnVector r_Moon(3); 189 570e3*cos(2*l) + 246e3*cos(2*l2*D)  205e3*cos(lp2*D) 190 171e3*cos(l+2*D)  152e3*cos(l+lp2*D); 191 192 ColumnVector r_Moon(3); 193 193 r_Moon << R*cos(L)*cosB << R*sin(L)*cosB << R*sin(B); 194 194 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) 198 198 * PrecMatrix(MJD_J2000, Mjd_TT) 199 199 * r_Moon; 200 200 } 201 201 202 // Tidal Correction 202 // Tidal Correction 203 203 //////////////////////////////////////////////////////////////////////////// 204 204 ColumnVector t_tides::displacement(const bncTime& time, const ColumnVector& xyz) { … … 239 239 double x2Sun = 3.0 * L2 * scSun; 240 240 double x2Moon = 3.0 * L2 * scMoon; 241 241 242 242 const double gmWGS = 398.6005e12; 243 243 const double gms = 1.3271250e20; 244 244 const double gmm = 4.9027890e12; 245 245 246 double facSun = gms / gmWGS * 246 double facSun = gms / gmWGS * 247 247 (rRec * rRec * rRec * rRec) / (_rSun * _rSun * _rSun); 248 248 249 double facMoon = gmm / gmWGS * 249 double facMoon = gmm / gmWGS * 250 250 (rRec * rRec * rRec * rRec) / (_rMoon * _rMoon * _rMoon); 251 251 252 ColumnVector dX = facSun * (x2Sun * _xSun + p2Sun * xyzUnit) + 252 ColumnVector dX = facSun * (x2Sun * _xSun + p2Sun * xyzUnit) + 253 253 facMoon * (x2Moon * _xMoon + p2Moon * xyzUnit); 254 254 255 255 return dX; 256 } 257 258 // Constructor 259 /////////////////////////////////////////////////////////////////////////// 260 t_loading::t_loading(const QString& fileName) { 261 262 newBlqData = 0; 263 readFile(fileName); 264 printAll(); 265 } 266 267 t_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 280 t_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[row1][ii]; 312 } 313 break; 314 case 4: case 5: case 6: 315 for (int ii = 0; ii < 11; ii++) { 316 inLine >> newBlqData>phases[row4][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 //////////////////////////////////////////////////////////////////////////// 336 void 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 } 256 357 } 257 358 … … 267 368 // Phase WindUp Correction 268 369 /////////////////////////////////////////////////////////////////////////// 269 double t_windUp::value(const bncTime& etime, const ColumnVector& rRec, 370 double t_windUp::value(const bncTime& etime, const ColumnVector& rRec, 270 371 t_prn prn, const ColumnVector& rSat) { 271 372 … … 276 377 ColumnVector rho = rRec  rSat; 277 378 rho /= rho.norm_Frobenius(); 278 379 279 380 // GPS Satellite unit Vectors sz, sy, sx 280 381 //  … … 289 390 // Effective Dipole of the GPS Satellite Antenna 290 391 //  291 ColumnVector dipSat = sx  rho * DotProduct(rho,sx) 392 ColumnVector dipSat = sx  rho * DotProduct(rho,sx) 292 393  crossproduct(rho, sy); 293 394 294 395 // Receiver unit Vectors rx, ry 295 396 //  … … 299 400 double recEll[3]; xyz2ell(rRec.data(), recEll) ; 300 401 double neu[3]; 301 402 302 403 neu[0] = 1.0; 303 404 neu[1] = 0.0; 304 405 neu[2] = 0.0; 305 406 neu2xyz(recEll, neu, rx.data()); 306 407 307 408 neu[0] = 0.0; 308 409 neu[1] = 1.0; 309 410 neu[2] = 0.0; 310 411 neu2xyz(recEll, neu, ry.data()); 311 412 312 413 // Effective Dipole of the Receiver Antenna 313 414 //  314 ColumnVector dipRec = rx  rho * DotProduct(rho,rx) 415 ColumnVector dipRec = rx  rho * DotProduct(rho,rx) 315 416 + crossproduct(rho, ry); 316 417 317 418 // Resulting Effect 318 419 //  319 double alpha = DotProduct(dipSat,dipRec) / 420 double alpha = DotProduct(dipSat,dipRec) / 320 421 (dipSat.norm_Frobenius() * dipRec.norm_Frobenius()); 321 422 322 423 if (alpha > 1.0) alpha = 1.0; 323 424 if (alpha < 1.0) alpha = 1.0; 324 425 325 426 double dphi = acos(alpha) / 2.0 / M_PI; // in cycles 326 427 327 428 if ( DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0 ) { 328 429 dphi = dphi; … … 339 440 } 340 441 341 return sumWind[prn.toInt()]; 442 return sumWind[prn.toInt()]; 342 443 } 343 444 … … 352 453 } 353 454 354 double ell[3]; 455 double ell[3]; 355 456 xyz2ell(xyz.data(), ell); 356 457 double height = ell[2]; … … 362 463 363 464 double h_km = height / 1000.0; 364 465 365 466 if (h_km < 0.0) h_km = 0.0; 366 467 if (h_km > 5.0) h_km = 5.0; 367 468 int ii = int(h_km + 1); if (ii > 5) ii = 5; 368 469 double href = ii  1; 369 370 double bCor[6]; 470 471 double bCor[6]; 371 472 bCor[0] = 1.156; 372 473 bCor[1] = 1.006; … … 375 476 bCor[4] = 0.654; 376 477 bCor[5] = 0.563; 377 478 378 479 double BB = bCor[ii1] + (bCor[ii]bCor[ii1]) * (h_km  href); 379 480 380 481 double zen = M_PI/2.0  Ele; 381 482 
trunk/BNC/src/pppModel.h
r7625 r7996 46 46 }; 47 47 48 class 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 48 67 class t_windUp { 49 68 public:
Note: See TracChangeset
for help on using the changeset viewer.