Changeset 6348 in ntrip for trunk/BNC/src


Ignore:
Timestamp:
Nov 26, 2014, 6:06:04 PM (10 years ago)
Author:
mervart
Message:
 
Location:
trunk/BNC/src/orbComp
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/BNC/src/orbComp/sp3Comp.cpp

    r6346 r6348  
    4040
    4141#include <iostream>
     42#include <iomanip>
    4243#include "sp3Comp.h"
    4344#include "bnccore.h"
    4445#include "bncsettings.h"
    4546#include "bncutils.h"
     47#include "bncsp3.h"
    4648
    4749using namespace std;
    … …  
    175177}
    176178
    177 //// Program
    178 //////////////////////////////////////////////////////////////////////////////////
    179 //int main(int argc, char* argv[]) {
    180 //
    181 //  t_sysout* sysout = new t_sysout("", 0, false);
    182 //  rtnet_core_log  = sysout->getSysoutCore();
    183 //
    184 //  // Parse Input Options
    185 //  // -------------------
    186 //  map<string, string> OPT;
    187 //  parseCmdLine(argc, argv, OPT);
    188 //  if (OPT.find("sp3File1") == OPT.end() || OPT.find("sp3File2") == OPT.end()) {
    189 //    cerr << "usage: sp3comp --sp3File1 <fileName> --sp3File2 <fileName> " << endl;
    190 //    return 1;
    191 //  }
    192 //
    193 //  // Synchronize reading of two sp3 files
    194 //  // ------------------------------------
    195 //  vector<t_prn> prnList;
    196 //  for (unsigned prn = 1; prn <= t_prn::MAXPRN; prn++) {
    197 //    prnList.push_back(t_prn(prn));
    198 //  }
    199 //  t_sp3 in1(prnList); if (in1.assignFile(OPT["sp3File1"].c_str())) return 1;
    200 //  t_sp3 in2(prnList); if (in2.assignFile(OPT["sp3File2"].c_str())) return 1;
    201 //
    202 //  vector<t_epoch*> epochs;
    203 //  set<t_prn> clkSats;
    204 //  while (true) {
    205 //    in1.readNextEpoch();
    206 //    bncTime tt = in1.etime();
    207 //    if (!tt.valid()) {
    208 //      break;
    209 //    }
    210 //    if (in2.readEpoch(tt) == t_irc::success) {
    211 //      t_epoch* epo = new t_epoch; epo->_tt = tt;
    212 //      bool epochOK = false;
    213 //      for (unsigned iPrn = 1; iPrn <= t_prn::MAXPRN; iPrn++) {
    214 //        t_prn prn(iPrn);
    215 //        ColumnVector xyz1(3), xyz2(3);
    216 //        ColumnVector vel1(3), vel2(3);
    217 //        double       clk1, clk2;
    218 //        bool         posSet1, posSet2;
    219 //        bool         clkSet1, clkSet2;
    220 //        bool         velSet1, velSet2;
    221 //        if (in1.position(prn, xyz1[0], xyz1[1], xyz1[2], posSet1, clk1, clkSet1,
    222 //                         vel1[0], vel1[1], vel1[2], velSet1) == 0 &&
    223 //            in2.position(prn, xyz2[0], xyz2[1], xyz2[2], posSet2, clk2, clkSet2,
    224 //                         vel1[0], vel1[1], vel1[2], velSet1) == 0) {
    225 //          if (posSet1 && posSet2) {
    226 //            epochOK       = true;
    227 //            epo->_dr[prn]  = xyz1 - xyz2;
    228 //            epo->_xyz[prn] = xyz1;
    229 //            if (clkSet1 && clkSet2) {
    230 //              epo->_dc[prn] = (clk1 - clk2) * t_genConst::c;
    231 //              clkSats.insert(prn);
    232 //            }
    233 //          }
    234 //        }
    235 //      }
    236 //      if (epochOK) {
    237 //        epochs.push_back(epo);
    238 //      }
    239 //      else {
    240 //        delete epo;
    241 //      }
    242 //    }
    243 //    if (in1.finished()) {
    244 //      break;
    245 //    }
    246 //  }
    247 //
    248 //  // Transform xyz into radial, along-track, and out-of-plane
    249 //  // --------------------------------------------------------
    250 //  for (unsigned ie = 0; ie < epochs.size(); ie++) {
    251 //    t_epoch* epoch  = epochs[ie];
    252 //    t_epoch* epoch2 = 0;
    253 //    if (ie == 0 && epochs.size() > 1) {
    254 //      epoch2 = epochs[ie+1];
    255 //    }
    256 //    else {
    257 //      epoch2 = epochs[ie-1];
    258 //    }
    259 //    double dt = epoch->_tt - epoch2->_tt;
    260 //    map<t_prn, ColumnVector>& dr   = epoch->_dr;
    261 //    map<t_prn, ColumnVector>& xyz  = epoch->_xyz;
    262 //    map<t_prn, ColumnVector>& xyz2 = epoch2->_xyz;
    263 //    for (map<t_prn, ColumnVector>::const_iterator it = dr.begin(); it != dr.end(); it++) {
    264 //      const t_prn&  prn = it->first;
    265 //      if (xyz2.find(prn) != xyz2.end()) {
    266 //        const ColumnVector  dx = dr[prn];
    267 //        const ColumnVector& x1 = xyz[prn];
    268 //        const ColumnVector& x2 = xyz2[prn];
    269 //        ColumnVector vel = (x1 - x2) / dt;
    270 //        t_astro::XYZ_to_RSW(x1, vel, dx, dr[prn]);
    271 //      }
    272 //      else {
    273 //        cerr << "not found: " << prn << endl;
    274 //      }
    275 //    }
    276 //  }
    277 //
    278 //  map<string, t_stat> stat;
    279 //
    280 //  // Estimate Clock Offsets
    281 //  // ----------------------
    282 //  processClocks(clkSats, epochs, stat);
    283 //
    284 //  // Print Residuals
    285 //  // ---------------
    286 //  const string all = "ZZZ";
    287 //
    288 //  cout.setf(ios::fixed);
    289 //  cout << "!\n!  MJD       PRN  radial   along   out        clk    clkRed   iPRN"
    290 //           "\n! ----------------------------------------------------------------\n";
    291 //  for (unsigned ii = 0; ii < epochs.size(); ii++) {
    292 //    const t_epoch* epo = epochs[ii];
    293 //    const map<t_prn, ColumnVector>& dr = epochs[ii]->_dr;
    294 //    const map<t_prn, double>&       dc = epochs[ii]->_dc;
    295 //    for (map<t_prn, ColumnVector>::const_iterator it = dr.begin(); it != dr.end(); it++) {
    296 //      const t_prn&  prn = it->first;
    297 //      const ColumnVector& rao = it->second;
    298 //      cout << setprecision(6) << epo->_tt.mjddec() << ' ' << prn << ' '
    299 //           << setw(7) << setprecision(4) << rao[0] << ' '
    300 //           << setw(7) << setprecision(4) << rao[1] << ' '
    301 //           << setw(7) << setprecision(4) << rao[2] << "    ";
    302 //      stat[prn]._rao += SP(rao, rao); // Schur product
    303 //      stat[prn]._nr  += 1;
    304 //      stat[all]._rao += SP(rao, rao);
    305 //      stat[all]._nr  += 1;
    306 //      if (dc.find(prn) != dc.end()) {
    307 //        double clkRes    = dc.find(prn)->second;
    308 //        double clkResRed = clkRes - it->second[0]; // clock minus radial component
    309 //        cout << setw(7) << setprecision(4) << clkRes << ' '
    310 //             << setw(7) << setprecision(4) << clkResRed;
    311 //        stat[prn]._dc    += clkRes * clkRes;
    312 //        stat[prn]._dcRed += clkResRed * clkResRed;
    313 //        stat[prn]._nc    += 1;
    314 //        stat[all]._dc    += clkRes * clkRes;
    315 //        stat[all]._dcRed += clkResRed * clkResRed;
    316 //        stat[all]._nc    += 1;
    317 //      }
    318 //      else {
    319 //        cout << "  .       .    ";
    320 //      }
    321 //      cout << "    " << setw(2) << int(prn) << endl;
    322 //    }
    323 //    delete epo;
    324 //  }
    325 //
    326 //  // Print Summary
    327 //  // -------------
    328 //  cout << "!\n! RMS:\n";
    329 //  for (map<string, t_stat>::iterator it = stat.begin(); it != stat.end(); it++) {
    330 //    const string& prn  = it->first;
    331 //    t_stat&       stat = it->second;
    332 //    if (stat._nr > 0) {
    333 //      stat._rao[0] = sqrt(stat._rao[0] / stat._nr);
    334 //      stat._rao[1] = sqrt(stat._rao[1] / stat._nr);
    335 //      stat._rao[2] = sqrt(stat._rao[2] / stat._nr);
    336 //      if (prn == all) {
    337 //        cout << "!\n!          Total ";
    338 //      }
    339 //      else {
    340 //        cout << "!            " << prn << ' ';
    341 //      }
    342 //      cout << setw(7) << setprecision(4) << stat._rao[0] << ' '
    343 //           << setw(7) << setprecision(4) << stat._rao[1] << ' '
    344 //           << setw(7) << setprecision(4) << stat._rao[2] << "    ";
    345 //      if (stat._nc > 0) {
    346 //        stat._dc    = sqrt(stat._dc / stat._nc);
    347 //        stat._dcRed = sqrt(stat._dcRed / stat._nc);
    348 //        cout << setw(7) << setprecision(4) << stat._dc << ' '
    349 //             << setw(7) << setprecision(4) << stat._dcRed;
    350 //        if (prn != all) {
    351 //          cout << "    offset " << setw(7) << setprecision(4) << stat._offset;
    352 //        }
    353 //      }
    354 //      else {
    355 //        cout << "  .       .    ";
    356 //      }
    357 //      cout << endl;
    358 //    }
    359 //  }
    360 //
    361 //  return 0;
    362 //}
     179// Main Routine
     180////////////////////////////////////////////////////////////////////////////////
     181void t_sp3Comp::compare() const {
     182
     183  // Synchronize reading of two sp3 files
     184  // ------------------------------------
     185  bncSP3 in1(_sp3FileNames[0]);
     186  bncSP3 in2(_sp3FileNames[1]);
     187
     188  vector<t_epoch*> epochs;
     189  set<t_prn>       clkSats;
     190  while (in1.nextEpoch()) {
     191    bncTime tt = in1.currEpoch()->_tt;
     192    while (in2.nextEpoch()) {
     193      if (tt == in2.currEpoch()->_tt) {
     194        t_epoch* epo = new t_epoch; epo->_tt = tt;
     195
     196        bool epochOK = false;
     197        for (int i1 = 0; i1 < in1.currEpoch()->_sp3Sat.size(); i1++) {
     198          bncSP3::t_sp3Sat* sat1 = in1.currEpoch()->_sp3Sat[i1];
     199          for (int i2 = 0; i2 < in2.currEpoch()->_sp3Sat.size(); i2++) {
     200            bncSP3::t_sp3Sat* sat2 = in2.currEpoch()->_sp3Sat[i2];
     201            if (sat1->_prn == sat2->_prn) {
     202              epochOK        = true;
     203              epo->_dr[sat1->_prn]  = sat1->_xyz - sat2->_xyz;
     204              epo->_xyz[sat1->_prn] = sat1->_xyz;
     205              if (true) {  //// TODO: check whether clocks set in SP3 files
     206                epo->_dc[sat1->_prn] = (sat1->_clk - sat2->_clk) * t_CST::c;
     207                clkSats.insert(sat1->_prn);
     208              }
     209            }
     210          }
     211        }
     212        if (epochOK) {
     213          epochs.push_back(epo);
     214        }
     215        else {
     216          delete epo;
     217        }
     218      }
     219    }
     220  }
     221
     222  // Transform xyz into radial, along-track, and out-of-plane
     223  // --------------------------------------------------------
     224  for (unsigned ie = 0; ie < epochs.size(); ie++) {
     225    t_epoch* epoch  = epochs[ie];
     226    t_epoch* epoch2 = 0;
     227    if (ie == 0 && epochs.size() > 1) {
     228      epoch2 = epochs[ie+1];
     229    }
     230    else {
     231      epoch2 = epochs[ie-1];
     232    }
     233    double dt = epoch->_tt - epoch2->_tt;
     234    map<t_prn, ColumnVector>& dr   = epoch->_dr;
     235    map<t_prn, ColumnVector>& xyz  = epoch->_xyz;
     236    map<t_prn, ColumnVector>& xyz2 = epoch2->_xyz;
     237    for (map<t_prn, ColumnVector>::const_iterator it = dr.begin(); it != dr.end(); it++) {
     238      const t_prn&  prn = it->first;
     239      if (xyz2.find(prn) != xyz2.end()) {
     240        const ColumnVector  dx = dr[prn];
     241        const ColumnVector& x1 = xyz[prn];
     242        const ColumnVector& x2 = xyz2[prn];
     243        ColumnVector vel = (x1 - x2) / dt;
     244        XYZ_to_RSW(x1, vel, dx, dr[prn]);
     245      }
     246      else {
     247        cerr << "not found: " << prn << endl;
     248      }
     249    }
     250  }
     251
     252  map<string, t_stat> stat;
     253
     254  // Estimate Clock Offsets
     255  // ----------------------
     256  processClocks(clkSats, epochs, stat);
     257
     258  // Print Residuals
     259  // ---------------
     260  const string all = "ZZZ";
     261
     262  cout.setf(ios::fixed);
     263  cout << "!\n!  MJD       PRN  radial   along   out        clk    clkRed   iPRN"
     264           "\n! ----------------------------------------------------------------\n";
     265  for (unsigned ii = 0; ii < epochs.size(); ii++) {
     266    const t_epoch* epo = epochs[ii];
     267    const map<t_prn, ColumnVector>& dr = epochs[ii]->_dr;
     268    const map<t_prn, double>&       dc = epochs[ii]->_dc;
     269    for (map<t_prn, ColumnVector>::const_iterator it = dr.begin(); it != dr.end(); it++) {
     270      const t_prn&  prn = it->first;
     271      const ColumnVector& rao = it->second;
     272      cout << setprecision(6) << epo->_tt.mjddec() << ' ' << prn << ' '
     273           << setw(7) << setprecision(4) << rao[0] << ' '
     274           << setw(7) << setprecision(4) << rao[1] << ' '
     275           << setw(7) << setprecision(4) << rao[2] << "    ";
     276      stat[prn.toString()]._rao += SP(rao, rao); // Schur product
     277      stat[prn.toString()]._nr  += 1;
     278      stat[all]._rao            += SP(rao, rao);
     279      stat[all]._nr             += 1;
     280      if (dc.find(prn) != dc.end()) {
     281        double clkRes    = dc.find(prn)->second;
     282        double clkResRed = clkRes - it->second[0]; // clock minus radial component
     283        cout << setw(7) << setprecision(4) << clkRes << ' '
     284             << setw(7) << setprecision(4) << clkResRed;
     285        stat[prn.toString()]._dc    += clkRes * clkRes;
     286        stat[prn.toString()]._dcRed += clkResRed * clkResRed;
     287        stat[prn.toString()]._nc    += 1;
     288        stat[all]._dc               += clkRes * clkRes;
     289        stat[all]._dcRed            += clkResRed * clkResRed;
     290        stat[all]._nc               += 1;
     291      }
     292      else {
     293        cout << "  .       .    ";
     294      }
     295      cout << "    " << setw(2) << int(prn) << endl;
     296    }
     297    delete epo;
     298  }
     299
     300  // Print Summary
     301  // -------------
     302  cout << "!\n! RMS:\n";
     303  for (map<string, t_stat>::iterator it = stat.begin(); it != stat.end(); it++) {
     304    const string& prn  = it->first;
     305    t_stat&       stat = it->second;
     306    if (stat._nr > 0) {
     307      stat._rao[0] = sqrt(stat._rao[0] / stat._nr);
     308      stat._rao[1] = sqrt(stat._rao[1] / stat._nr);
     309      stat._rao[2] = sqrt(stat._rao[2] / stat._nr);
     310      if (prn == all) {
     311        cout << "!\n!          Total ";
     312      }
     313      else {
     314        cout << "!            " << prn << ' ';
     315      }
     316      cout << setw(7) << setprecision(4) << stat._rao[0] << ' '
     317           << setw(7) << setprecision(4) << stat._rao[1] << ' '
     318           << setw(7) << setprecision(4) << stat._rao[2] << "    ";
     319      if (stat._nc > 0) {
     320        stat._dc    = sqrt(stat._dc / stat._nc);
     321        stat._dcRed = sqrt(stat._dcRed / stat._nc);
     322        cout << setw(7) << setprecision(4) << stat._dc << ' '
     323             << setw(7) << setprecision(4) << stat._dcRed;
     324        if (prn != all) {
     325          cout << "    offset " << setw(7) << setprecision(4) << stat._offset;
     326        }
     327      }
     328      else {
     329        cout << "  .       .    ";
     330      }
     331      cout << endl;
     332    }
     333  }
     334}

  • trunk/BNC/src/orbComp/sp3Comp.h

    r6345 r6348  
    8383  void processClocks(const std::set<t_prn>& clkSats, std::vector<t_epoch*>& epochs,
    8484                     std::map<std::string, t_stat>& stat) const;
     85  void compare() const;
    8586
    8687  QStringList  _sp3FileNames;
Note: See TracChangeset for help on using the changeset viewer.