Changeset 7481 in ntrip


Ignore:
Timestamp:
Sep 30, 2015, 9:37:44 PM (9 years ago)
Author:
stuerze
Message:

distinction of GEO/MEO satellites during BDS velocity determination

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/BNC/src/ephemeris.cpp

    r7278 r7481  
    3333}
    3434
    35 // 
     35//
    3636////////////////////////////////////////////////////////////////////////////
    3737void t_eph::setOrbCorr(const t_orbCorr* orbCorr) {
    38   delete _orbCorr; 
     38  delete _orbCorr;
    3939  _orbCorr = new t_orbCorr(*orbCorr);
    4040}
    4141
    42 // 
     42//
    4343////////////////////////////////////////////////////////////////////////////
    4444void t_eph::setClkCorr(const t_clkCorr* clkCorr) {
    45   delete _clkCorr; 
     45  delete _clkCorr;
    4646  _clkCorr = new t_clkCorr(*clkCorr);
    4747}
    4848
    49 // 
     49//
    5050////////////////////////////////////////////////////////////////////////////
    5151t_irc t_eph::getCrd(const bncTime& tt, ColumnVector& xc, ColumnVector& vv, bool useCorr) const {
    … …  
    319319  double xp     = r*cos(u);
    320320  double yp     = r*sin(u);
    321   double OM     = _OMEGA0 + (_OMEGADOT - omegaEarth)*tk - 
     321  double OM     = _OMEGA0 + (_OMEGADOT - omegaEarth)*tk -
    322322                   omegaEarth*_TOEsec;
    323323
    … …  
    328328  xc[0] = xp*cosom - yp*cosi*sinom;
    329329  xc[1] = xp*sinom + yp*cosi*cosom;
    330   xc[2] = yp*sini;               
    331  
     330  xc[2] = yp*sini;
     331
    332332  double tc = tt - _TOC;
    333333  xc[3] = _clock_bias + _clock_drift*tc + _clock_driftrate*tc*tc;
    … …  
    359359  // Relativistic Correction
    360360  // -----------------------
     361  // correspondent to IGS convention and GPS ICD (and SSR standard)
    361362  xc[3] -= 2.0 * (xc[0]*vv[0] + xc[1]*vv[1] + xc[2]*vv[2]) / t_CST::c / t_CST::c;
    362363
    … …  
    894895  // Relativistic Correction
    895896  // -----------------------
    896   //  xc(4) -= 4.442807633e-10 * _e * sqrt(a0) *sin(E);
    897   xc[3] -= 2.0 * (xc[0]*vv[0] + xc[1]*vv[1] + xc[2]*vv[2]) / t_CST::c / t_CST::c;
     897  // correspondent to Galileo ICD and to SSR standard
     898  xc[3] -= 4.442807633e-10 * _e * sqrt(a0) *sin(E);
     899  // correspondent to IGS convention
     900  //xc[3] -= 2.0 * (xc[0]*vv[0] + xc[1]*vv[1] + xc[2]*vv[2]) / t_CST::c / t_CST::c;
    898901
    899902  return success;
    … …  
    14641467  double cosi  = 0;
    14651468
    1466   const double iMaxGEO = 10.0 / 180.0 * M_PI;
    1467 
    1468   // MEO/IGSO satellite
    1469   // ------------------
    1470   if (_i0 > iMaxGEO) {
    1471     double OM = _OMEGA0 + (_OMEGADOT - omegaBDS)*tk - omegaBDS*toesec;
    1472 
    1473     sinom = sin(OM);
    1474     cosom = cos(OM);
    1475     sini  = sin(i);
    1476     cosi  = cos(i);
    1477 
    1478     xc[0] = xp*cosom - yp*cosi*sinom;
    1479     xc[1] = xp*sinom + yp*cosi*cosom;
    1480     xc[2] = yp*sini;
    1481   }
    1482 
    1483   // GEO satellite
    1484   // -------------
    1485   else {
    1486     double OM    = _OMEGA0 + _OMEGADOT*tk - omegaBDS*toesec;
    1487     double ll    = omegaBDS*tk;
    1488 
    1489     sinom = sin(OM);
    1490     cosom = cos(OM);
    1491     sini  = sin(i);
    1492     cosi  = cos(i);
    1493 
    1494     double xx = xp*cosom - yp*cosi*sinom;
    1495     double yy = xp*sinom + yp*cosi*cosom;
    1496     double zz = yp*sini;
    1497 
    1498     Matrix R1 = BNC_PPP::t_astro::rotX(-5.0 / 180.0 * M_PI);
    1499     Matrix R2 = BNC_PPP::t_astro::rotZ(ll);
    1500 
    1501     ColumnVector X1(3); X1 << xx << yy << zz;
    1502     ColumnVector X2 = R2*R1*X1;
    1503 
    1504     xc[0] = X2(1);
    1505     xc[1] = X2(2);
    1506     xc[2] = X2(3);
    1507   }
    1508 
    1509   double tc = tt - _TOC;
    1510   xc[3] = _clock_bias + _clock_drift*tc + _clock_driftrate*tc*tc
    1511           - 4.442807633e-10 * _e * sqrt(a0) *sin(E);
    1512 
    15131469  // Velocity
    15141470  // --------
    … …  
    15181474                 / (1 + tanv2*tanv2) * dEdM * n;
    15191475  double dotu  = dotv + (-_Cuc*sin2u0 + _Cus*cos2u0)*2*dotv;
    1520   double dotom = _OMEGADOT - t_CST::omega;
    15211476  double doti  = _IDOT + (-_Cic*sin2u0 + _Cis*cos2u0)*2*dotv;
    15221477  double dotr  = a0 * _e*sin(E) * dEdM * n
    15231478                + (-_Crc*sin2u0 + _Crs*cos2u0)*2*dotv;
     1479
    15241480  double dotx  = dotr*cos(u) - r*sin(u)*dotu;
    15251481  double doty  = dotr*sin(u) + r*cos(u)*dotu;
    15261482
    1527   vv[0]  = cosom  *dotx  - cosi*sinom   *doty   // dX / dr
    1528         - xp*sinom*dotom - yp*cosi*cosom*dotom   // dX / dOMEGA
    1529                          + yp*sini*sinom*doti;   // dX / di
    1530 
    1531   vv[1]  = sinom  *dotx  + cosi*cosom   *doty
    1532         + xp*cosom*dotom - yp*cosi*sinom*dotom
    1533                          - yp*sini*cosom*doti;
    1534 
    1535   vv[2]  = sini   *doty  + yp*cosi      *doti;
     1483
     1484  const double iMaxGEO = 10.0 / 180.0 * M_PI;
     1485
     1486  // MEO/IGSO satellite
     1487  // ------------------
     1488  if (_i0 > iMaxGEO) {
     1489    double OM = _OMEGA0 + (_OMEGADOT - omegaBDS)*tk - omegaBDS*toesec;
     1490
     1491    sinom = sin(OM);
     1492    cosom = cos(OM);
     1493    sini  = sin(i);
     1494    cosi  = cos(i);
     1495
     1496    xc[0] = xp*cosom - yp*cosi*sinom;
     1497    xc[1] = xp*sinom + yp*cosi*cosom;
     1498    xc[2] = yp*sini;
     1499
     1500    // Velocity
     1501    // --------
     1502
     1503    double dotom = _OMEGADOT - t_CST::omega;
     1504
     1505    vv[0]  = cosom  *dotx   - cosi*sinom   *doty    // dX / dr
     1506           - xp*sinom*dotom - yp*cosi*cosom*dotom   // dX / dOMEGA
     1507                            + yp*sini*sinom*doti;   // dX / di
     1508
     1509    vv[1]  = sinom  *dotx   + cosi*cosom   *doty
     1510           + xp*cosom*dotom - yp*cosi*sinom*dotom
     1511                            - yp*sini*cosom*doti;
     1512
     1513    vv[2]  = sini   *doty   + yp*cosi      *doti;
     1514
     1515  }
     1516
     1517  // GEO satellite
     1518  // -------------
     1519  else {
     1520    double OM    = _OMEGA0 + _OMEGADOT*tk - omegaBDS*toesec;
     1521    double ll    = omegaBDS*tk;
     1522
     1523    sinom = sin(OM);
     1524    cosom = cos(OM);
     1525    sini  = sin(i);
     1526    cosi  = cos(i);
     1527
     1528    double xx = xp*cosom - yp*cosi*sinom;
     1529    double yy = xp*sinom + yp*cosi*cosom;
     1530    double zz = yp*sini;
     1531
     1532    Matrix R1 = BNC_PPP::t_astro::rotX(-5.0 / 180.0 * M_PI);
     1533    Matrix R2 = BNC_PPP::t_astro::rotZ(ll);
     1534
     1535    ColumnVector X1(3); X1 << xx << yy << zz;
     1536    ColumnVector X2 = R2*R1*X1;
     1537
     1538    xc[0] = X2(1);
     1539    xc[1] = X2(2);
     1540    xc[2] = X2(3);
     1541
     1542    double dotom = _OMEGADOT;
     1543
     1544    double vx  = cosom   *dotx  - cosi*sinom   *doty
     1545               - xp*sinom*dotom - yp*cosi*cosom*dotom
     1546                                + yp*sini*sinom*doti;
     1547
     1548    double vy  = sinom   *dotx  + cosi*cosom   *doty
     1549               + xp*cosom*dotom - yp*cosi*sinom*dotom
     1550                                - yp*sini*cosom*doti;
     1551
     1552    //double vz  = sini    *doty  + yp*cosi      *doti;
     1553
     1554    ColumnVector V(3);
     1555
     1556    V[0]   = cosom  *vx     - cosi*sinom   *vy      // dX / dr
     1557           - xp*sinom*dotom - yp*cosi*cosom*dotom   // dX / dOMEGA
     1558                            + yp*sini*sinom*doti;   // dX / di
     1559
     1560    V[1]   = sinom  *vx     + cosi*cosom   *vy
     1561           + xp*cosom*dotom - yp*cosi*sinom*dotom
     1562                            - yp*sini*cosom*doti;
     1563
     1564    V[2]   = sini   *vy     + yp*cosi      *doti;
     1565
     1566    Matrix R3(3,3);
     1567    double C = cos(ll);
     1568    double S = sin(ll);
     1569    Matrix UU(3,3);
     1570    UU[0][0] =  -S;  UU[0][1] =  +C;  UU[0][2] = 0.0;
     1571    UU[1][0] =  -C;  UU[1][1] =  -S;  UU[1][2] = 0.0;
     1572    UU[2][0] = 0.0;  UU[2][1] = 0.0;  UU[2][2] = 0.0;
     1573    R3 = omegaBDS * UU;
     1574
     1575    ColumnVector VV(3);
     1576
     1577    VV = R2*R1*V + R3*R1*X1;
     1578
     1579    vv[0] = VV(1);
     1580    vv[1] = VV(2);
     1581    vv[2] = VV(3);
     1582  }
     1583
     1584  double tc = tt - _TOC;
     1585  xc[3] = _clock_bias + _clock_drift*tc + _clock_driftrate*tc*tc;
    15361586
    15371587  // dotC  = _clock_drift + _clock_driftrate*tc
    15381588  //       - 4.442807633e-10*_e*sqrt(a0)*cos(E) * dEdM * n;
     1589
     1590  // Relativistic Correction
     1591  // -----------------------
     1592  // correspondent to BDS ICD and to SSR standard
     1593    xc[3] -= 4.442807633e-10 * _e * sqrt(a0) *sin(E);
     1594  // correspondent to IGS convention
     1595  // xc[3] -= 2.0 * (xc[0]*vv[0] + xc[1]*vv[1] + xc[2]*vv[2]) / t_CST::c / t_CST::c;
    15391596
    15401597  return success;
Note: See TracChangeset for help on using the changeset viewer.