rinsum.cpp

GPS数据预处理软件

               }
               // dump this data
               if(debug) *pout << " " << RinexObsHeader::convertObsType(jt->first)
                  << " " << setw(13) << setprecision(3) << jt->second.data << " "
                  << jt->second.lli << " " << jt->second.ssi;
            }  // end loop over obs types
            if(debug) *pout << endl;

            // test for millisecond clock adjusts -
            // sometimes they are applied to range but not phase or vice-versa
            if(prev != DayTime::BEGINNING_OF_TIME && L1 != 0 && ptab->prevL1 != 0) {
               int nms;
               double test;
               nsats++;
               if(P1 != 0 && ptab->prevP1 != 0)
                  test = P1-L1_WAVELENGTH*L1
                     - (ptab->prevP1-L1_WAVELENGTH*ptab->prevL1);
               else if(C1 != 0 && ptab->prevC1 != 0)
                  test = C1-L1_WAVELENGTH*L1
                     - (ptab->prevC1-L1_WAVELENGTH*ptab->prevL1);
               else
                  test = 0.0;
               if(fabs(test) > 0.5) {      // test must be > 150 km =~ 1/2 millisecond
                  // is it nearly an even multiple of 1 millisecond?
                  //test *= 1000.0/C_GPS_M;
                  if(debug) *pout << "possible clock jump: test = "
                                 << setprecision(9) << test;
                  nms = long(test + (test > 0 ? 0.5 : -0.5));
                  if(fabs(test - double(nms)) < 0.001) {
                     if(debug) *pout << " -> " << setprecision(9)
                              << fabs(test - double(nms));
                     // keep clkjumpave = sequential average nms, clkjumpvar=variance
                     if(test < 0) nms *= -1;
                     nclkjumps++;
                     clkjumpave += (double(nms)-clkjumpave)/double(nclkjumps);
                     if(nclkjumps > 1)
                        clkjumpvar = (clkjumpvar*(nclkjumps-2)
                         + nclkjumps*(double(nms)-clkjumpave)*(double(nms)-clkjumpave)
                            /(nclkjumps-1))/(nclkjumps-1);
                  }
                  else if(debug) *pout << " - failed.";
                  if(debug && L1lli != 0) { *pout << " LLI is set"; }
                  if(debug) *pout << " " << RinexSatID(it->first)
                     << " " << last.printf("%4F %.3g") << endl;
               }
            }
            // save C1,L1,P1 for this sat for next time
            ptab->prevC1 = C1;
            ptab->prevL1 = L1;
            ptab->prevP1 = P1;

         }  // end loop over sats

         //out << robs;

         // if more than half the sats saw a clk jump, call it
         if(nclkjumps > nsats/2) {
            if(debug) *pout << "test nclkjumps is " << nclkjumps
               << " and nsats is " << nsats
               << ", ave is " << fixed << setprecision(3) << clkjumpave
               << " and stddev is " << setprecision(3) << sqrt(clkjumpvar)
               << endl;
            clkjumpTimes.push_back(last);
            clkjumpMillsecs.push_back(clkjumpave);
            clkjumpAgree.push_back(nsats-nclkjumps);
            clkjumpUncertainty.push_back(sqrt(clkjumpvar));
         }

         if(prev != DayTime::BEGINNING_OF_TIME) {
            dt = last-prev;
            if(dt > 0.0) {
               for(i=0; i<ndtmax; i++) {
                  if(ndt[i] <= 0) { bestdt[i]=dt; ndt[i]=1; break; }
                  if(fabs(dt-bestdt[i]) < 0.0001) { ndt[i]++; break; }
                  if(i == ndtmax-1) {
                     k = 0;
                     int nleast=ndt[k];
                     for(j=1; j<ndtmax; j++) if(ndt[j] <= nleast) {
                        k=j; nleast=ndt[j];
                     }
                     ndt[k]=1; bestdt[k]=dt;
                  }
               }
            }
            else {
               cerr << " WARNING time tags out of order: "
                  //<< " prev > curr : "
                  << prev.printf("%F/%.0g = %04Y/%02m/%02d %02H:%02M:%02S")
                  << " > "
                  << last.printf("%F/%.0g = %04Y/%02m/%02d %02H:%02M:%02S")
                  << endl;
            }
         }
         prev = last;

      }  // end loop over epochs in the file
      InStream.close();

         // check that we found some data
      if(nepochs <= 0) {
         *pout << "File " << filename << " : no data found. Are time limits wrong?\n";
         continue;
      }

         // compute interval
      for(i=1,j=0; i<ndtmax; i++) if(ndt[i]>ndt[j]) j=i;
      dt = bestdt[j];

         // summary info
      *pout << "Computed interval "
         << fixed << setw(5) << setprecision(2) << dt << " seconds." << endl;
      *pout << "Computed first epoch: " << ftime.printf("%4F %14.7g") << " = "
            << ftime.printf("%04Y/%02m/%02d %02H:%02M:%010.7f") << endl;
      *pout << "Computed last  epoch: " << last.printf("%4F %14.7g") << " = "
            << last.printf("%04Y/%02m/%02d %02H:%02M:%010.7f") << endl;

      *pout << "Computed time span:";
      double secs=last-ftime;
      int iday = int(secs/86400.0);
      if(iday > 0) *pout << " " << iday << "d";
      DayTime delta;
      delta.setSecOfDay(secs - iday*86400);
      *pout << " " << delta.hour() << "h "
         << delta.minute() << "m "
         << delta.second() << "s = "
         << secs << " seconds." << endl;

      i = 1+int(0.5+(last-ftime)/dt);
      if(!brief) *pout << "There were " << nepochs << " epochs ("
         << setprecision(2) << double(nepochs*100)/i
         << "% of " << i << " possible epochs in this timespan) and "
         << ncommentblocks << " inline header blocks.\n";

         // sort table
      sort(table.begin(),table.end(),TableSATLessThan());
      if(TimeSortTable) sort(table.begin(),table.end(),TableBegLessThan());

         // output table
         // header
      vector<TableData>::iterator tit;
      if(table.size() > 0) table.begin()->sat.setfill('0');
      if(!brief) {
         *pout << "\n          Summary of data available in this file: "
            << "(Totals are based on times and interval)\n";
         *pout << "Sat  OT:";
         for(k=0; k<n; k++)
            *pout << setw(7) << rheader.obsTypeList[k].type;
         *pout << "  Total             Begin time - End time\n";
            // loop
         for(tit=table.begin(); tit!=table.end(); ++tit) {
            *pout << "Sat " << tit->sat << " ";
            for(k=0; k<n; k++) *pout << setw(7) << tit->nobs[k];
            // compute total based on times
            *pout << setw(7) << 1+int(0.5+(tit->end-tit->begin)/dt);
            if(GPSTimeOutput) {
               *pout << "  " << tit->begin.printf("%4F %10.3g")
                  << " - " << tit->end.printf("%4F %10.3g") << endl;
            }
            else {
               *pout
                  << "  " << tit->begin.printf("%04Y/%02m/%02d %02H:%02M:%04.1f")
                  << " - " << tit->end.printf("%04Y/%02m/%02d %02H:%02M:%04.1f")
                  << endl;
            }
         }
         *pout << "TOTAL   "; for(k=0; k<n; k++) *pout << setw(7) << totals[k];
         *pout << endl;
      }
      else {
         *pout << "SATs(" << table.size() << "):";
         for(tit=table.begin(); tit!=table.end(); ++tit)
            *pout << " " << tit->sat;
         *pout << endl;

         *pout << "Obs types(" << rheader.obsTypeList.size() << "): ";
         for(i=0; i<rheader.obsTypeList.size(); i++)
            *pout << " " << rheader.obsTypeList[i].type;
         *pout << endl;
      }

         // warnings
      if((rheader.valid & RinexObsHeader::intervalValid)
            && fabs(dt-rheader.interval) > 1.e-3)
         *pout << " WARNING: Computed interval is " << setprecision(2)
            << dt << " sec, while input header has " << setprecision(2)
            << rheader.interval << " sec.\n";
      if(fabs(ftime-rheader.firstObs) > 1.e-8)
         *pout << " WARNING: Computed first time does not agree with header\n";
      if((rheader.valid & RinexObsHeader::lastTimeValid)
            && fabs(last-rheader.lastObs) > 1.e-8)
         *pout << " WARNING: Computed last time does not agree with header\n";

      if(clkjumpTimes.size() > 0) {
         *pout << " WARNING: millisecond clock adjusts at these times:\n";
         for(i=0; i<clkjumpTimes.size(); i++) {
            *pout << "   "
             << clkjumpTimes[i].printf("%4F %10.3g = %04Y/%02m/%02d %02H:%02M:%06.3f")
             << " " << setprecision(2) << clkjumpMillsecs[i] << " ms_clock_adjust";
             if(clkjumpAgree[i] > 0 || clkjumpUncertainty[i] > 0.01)
               *pout << " (low quality determination; data may be irredeemable)";
            *pout << endl;
         }
      }
         // look for 'empty' obs types
      for(k=0; k<n; k++) {
         if(totals[k] <= 0) *pout << " WARNING: ObsType "
            << rheader.obsTypeList[k].type
            << " should be deleted from header.\n";
      }

      if(ReplaceHeader) {
            // modify the header
         rheader.version = 2.1; rheader.valid |= RinexObsHeader::versionValid;
         rheader.interval = dt; rheader.valid |= RinexObsHeader::intervalValid;
         rheader.lastObs = last; rheader.valid |= RinexObsHeader::lastTimeValid;
            // now the table
         rheader.numSVs = table.size(); rheader.valid |= RinexObsHeader::numSatsValid;
         rheader.numObsForSat.clear();
         for(tit=table.begin(); tit!=table.end(); ++tit) {      // tit defined above
            rheader.numObsForSat.insert(
               map<SatID, vector<int> >::value_type(tit->sat,tit->nobs) );
         }
         rheader.valid |= RinexObsHeader::prnObsValid;
         //*pout << "\nNew header\n";
         //rheader.dump(*pout);

            // now re-open the file and replace the header
#ifdef _MSC_VER
         char newname[L_tmpnam];
         if(!tmpnam(newname)) {
            cerr << "Could not create temporary file name - abort\n";
            return -1;
         }
#else
         char newname[]="RinSumTemp.XXXXXX";
         if(mkstemp(newname)==-1) {
            cerr << "Could not create temporary file name - abort\n";
            return -1;
         }
#endif
         remove(newname);

         RinexObsHeader rhjunk;
         RinexObsStream ROutStr(newname, ios::out);
         RinexObsStream InAgain(filename.c_str());
         InAgain.exceptions(ios::failbit);

         InAgain >> rhjunk;
         ROutStr << rheader;
         while(InAgain >> robs) {
            last = robs.time;
            if(last < BegTime) continue;
            if(last > EndTime) break;
            ROutStr << robs;
         }
         InAgain.close();
         ROutStr.close();
            // delete original file and rename the temporary
         iret = remove(filename.c_str());
         if(iret) *pout << "RinSum: Error: Could not remove existing file: "
            << filename << endl;
         else {
            iret = rename(newname,filename.c_str());
            if(iret) *pout << "RinSum: Error: Could not rename new file " << newname
               << " using old name " << filename << endl;
            else *pout << "\nRinSum: Replaced original header with complete one,"
               << " using temporary file name "
               << newname << endl;
         }
      }

      if(!brief) *pout << "\n+++++++++++++ End of RinSum summary of " << filename
         << " +++++++++++++\n";
   }

   if(pout != &cout) {
      ((ofstream *)pout)->close();
      delete pout;
   }

   return 0;
}
catch(gpstk::FFStreamError& e) { cerr << "FFStreamError: " << e; }
catch(gpstk::Exception& e) { cerr << "Exception: " << e; }
catch (...) { cerr << "Unknown exception.  Abort." << endl; }