discfix.cpp

linux的gps应用

#pragma ident "$Id: DiscFix.cpp 327 2006-11-30 18:37:30Z ehagen $"

//============================================================================
//
//  This file is part of GPSTk, the GPS Toolkit.
//
//  The GPSTk is free software; you can redistribute it and/or modify
//  it under the terms of the GNU Lesser General Public License as published
//  by the Free Software Foundation; either version 2.1 of the License, or
//  any later version.
//
//  The GPSTk is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with GPSTk; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//  
//  Copyright 2004, The University of Texas at Austin
//
//============================================================================

//------------------------------------------------------------------------------------
// DiscFix.cpp Read a RINEX observation file containing dual frequency
//    pseudorange and phase, separate the data into satellite passes, and then
//    find and estimate discontinuities in the phase (using the GPSTk Discontinuity
//    Corrector (GDC) in DiscCorr.hpp).
//    The corrected data can be written out to another RINEX file, plus there is the
//    option to smooth the pseudorange and/or debias the phase (SatPass::smooth()).
//
//    This program is useful as a way to process RINEX data by satellite pass.
//    It reads an entire RINEX obs file, breaking it into satellite passes (SatPass)
//    and processing it (ProcessSatPass()), then writes it out again from the
//    SatPass data. It was designed so that all the input data gets into one SatPass
//    and is altered only by the routine(s) called in ProcessSatPass(). Thus by
//    modifying just that routine, this program could be used to do something else
//    to the satellite passes. Also note that there is a choice of when to write out
//    the data, either as soon as possible, or only at the end: cf. bool WriteASAP.
//------------------------------------------------------------------------------------

/**
 * @file DiscFix.cpp
 * Correct phase discontinuities (cycle slips) in dual frequency data in a RINEX
 * observation file, plus optionally smooth the pseudoranges and/or debias the phases.
 */

#include "MathBase.hpp"
#include "RinexObsBase.hpp"
#include "RinexObsData.hpp"
#include "RinexObsHeader.hpp"
#include "RinexObsStream.hpp"
#include "DayTime.hpp"
#include "StringUtils.hpp"

#include "SatPass.hpp"
#include "DiscCorr.hpp"

#include <time.h>
#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include <algorithm>

using namespace std;
using namespace gpstk;

//------------------------------------------------------------------------------------
// prgm data
string PrgmName("DiscFix");
string PrgmVers("4.0 2/28/06");

typedef struct configuration {
      // input
   string Directory;
   vector<string> InputObsName;
      // data flow
   double ith;
   DayTime begTime, endTime;
   double MaxGap;
   //int MinPts;
      // processing
   double dt;
   bool UseCA;
   vector<RinexSatID> ExSV;
   RinexSatID SVonly;
      // output files
   string LogFile,OutFile;
   ofstream oflog,ofout;
   string format;
      // output
   string OutRinexObs;
   string HDPrgm;         // header of output RINEX file
   string HDRunby;
   string HDObs;
   string HDAgency;
   string HDMarker;
   string HDNumber;
   int NrecOut;
   DayTime FirstEpoch,LastEpoch;
   bool smoothPR,smoothPH,smooth;
   //bool CAOut;
   //bool DopOut;
   bool verbose;
      // estimate dt from data
   double estdt[9];
   int ndt[9];
} DFConfig;

//------------------------------------------------------------------------------------
// data input from command line
DFConfig config;                 // for DiscFix
GDCconfiguration GDConfig;       // the discontinuity corrector configuration

// if true, write to RINEX only after ALL data has been processed (by ProcessSatPass)
// (I can't see why DiscFix would want to set this false, so it is not in the input.)
bool WriteASAP=true;

// data used in program
clock_t totaltime;
string Title;
DayTime CurrEpoch(DayTime::BEGINNING_OF_TIME), PrgmEpoch;

RinexObsStream irfstr, orfstr;      // input and output RINEX files
RinexObsHeader rhead;
int inC1,inP1,inP2,inL1,inL2;      // indexes in rhead of C1, C1/P1, P2, L1 and L2
bool UsingCA;

// Data for an entire pass is stored in SatPass object
// this contains all the SatPass's defined so far
// the parallel vector holds an iterator for use in writing out the data
vector<SatPass> SPList;
vector<unsigned int> SPIndexList;
// this is a map relating a satellite to the index in SVPList of the current pass
map<RinexSatID,int> SatToCurrentIndexMap;

//------------------------------------------------------------------------------------
// prototypes
int ReadFile(int nfile);
int ProcessOneEntireEpoch(RinexObsData& ro);
int ProcessOneSatOneEpoch(RinexSatID, DayTime, SatPassData& );

void ProcessSatPass(int index);
int AfterReadingFiles(void);
void WriteToRINEXfile(void);
void WriteRINEXheader(void);
void WriteRINEXdata(DayTime& WriteEpoch, const DayTime targetTime);

void PrintSPList(ostream&, string, const vector<SatPass>&, bool printTime);
int GetCommandLine(int argc, char **argv);
void PreProcessArgs(const char *arg, vector<string>& Args, bool& verbose);

//------------------------------------------------------------------------------------
int main(int argc, char **argv)
{
   try {
      totaltime = clock();
      int iret;

         // Title and description
      //cout << "Name " << string(argv[0]) << endl;
      Title = PrgmName + ", part of the GPS ToolKit, Ver " + PrgmVers + ", Run ";
      PrgmEpoch.setLocalTime();
      Title += PrgmEpoch.printf("%04Y/%02m/%02d %02H:%02M:%02S\n");
      cout << Title;

         // set fill char in RinexSatID
      config.SVonly.setfill('0');
      config.FirstEpoch = DayTime::BEGINNING_OF_TIME;
      config.LastEpoch = DayTime::BEGINNING_OF_TIME;

         // get command line
      iret = GetCommandLine(argc, argv);
      if(iret) return iret;

         // configure SatPass
      {
         SatPass dummy(config.SVonly,config.dt);
         dummy.setMaxGap(config.MaxGap);
         dummy.setOutputFormat(config.format);
      }

         // open output files
         // output for editing commands - write to this in ProcessSatPass()
      config.ofout.open(config.OutFile.c_str());
      if(config.oflog.fail()) {
         config.oflog << "Error: " << PrgmName << " failed to open output file "
            << config.OutFile << endl;
      }
      else
         cout << PrgmName << " is writing to output file " << config.OutFile << endl;
         // RINEX output
      orfstr.open(config.OutRinexObs.c_str(), ios::out);
      if(!config.OutRinexObs.empty()) {
         if(orfstr.fail()) {
            config.oflog << "Failed to open output file " << config.OutRinexObs
               << ". Abort.\n";
            cout << "Failed to open output file " << config.OutRinexObs
               << ". Abort.\n";
            irfstr.close();
            return 1;
         }
         else cout << PrgmName << " is writing to RINEX file "
            << config.OutRinexObs << endl;
         orfstr.exceptions(ios::failbit);
      }

         // loop over input files
      for(int nfile=0; nfile<config.InputObsName.size(); nfile++) {
         iret = ReadFile(nfile);
         if(iret < 0) break;
      }   // end loop over input files

      iret = AfterReadingFiles();

         // clean up
      orfstr.close();
      config.ofout.close();
      SatToCurrentIndexMap.clear();
      SPList.clear();
      SPIndexList.clear();

      totaltime = clock()-totaltime;
      config.oflog << PrgmName << " timing: " << fixed << setprecision(3)
         << double(totaltime)/double(CLOCKS_PER_SEC) << " seconds.\n";
      cout << PrgmName << " timing: " << fixed << setprecision(3)
         << double(totaltime)/double(CLOCKS_PER_SEC) << " seconds.\n";

      config.oflog.close();

      return iret;
   }
   catch(gpstk::Exception& e) {
      config.oflog << e;
   }
   catch (...) {
      config.oflog << PrgmName << ": Unknown error.  Abort." << endl;
      cout << PrgmName << ": Unknown error.  Abort." << endl;
   }

   return -1;

}   // end main()

//------------------------------------------------------------------------------------
// open the file, read header and check for data; then loop over the epochs
// Return 0 ok, <0 fatal error, >0 non-fatal error (ie skip this file)
// 0 ok, 1 couldn't open file, 2 file doesn't have required data
int ReadFile(int nfile)
{
   try {
      string name;
         // open input file
      name = config.Directory + string("/") + config.InputObsName[nfile];
      irfstr.open(name.c_str(),ios::in);
      if(irfstr.fail()) {
         config.oflog << "Failed to open input file " << name << ". Abort.\n";
         cout << "Failed to open input file " << name << ". Abort.\n";
         return 1;
      }
      else if(config.verbose)
         config.oflog << "Opened input file " << name << endl;
      irfstr.exceptions(ios::failbit);

         // read the header
      irfstr >> rhead;
      if(config.verbose) {
         config.oflog << "\nHere is the input header for file " << name << endl;
         rhead.dump(config.oflog);
         config.oflog << endl;
      }

         // check that file contains C1/P1,P2,L1,L2
      inC1 = inP1 = inP2 = inL1 = inL2 = -1;
      for(int j=0; j<rhead.obsTypeList.size(); j++) {
         if(rhead.obsTypeList[j] == RinexObsHeader::convertObsType("C1")) inC1=j;
         if(rhead.obsTypeList[j] == RinexObsHeader::convertObsType("L1")) inL1=j;
         if(rhead.obsTypeList[j] == RinexObsHeader::convertObsType("L2")) inL2=j;
         if(rhead.obsTypeList[j] == RinexObsHeader::convertObsType("P1")) inP1=j;
         if(rhead.obsTypeList[j] == RinexObsHeader::convertObsType("P2")) inP2=j;
      }

      if((inC1 == -1 && config.UseCA) ||             // no C1, but user wants C1
         (inP1 == -1 && inC1 == -1) ||               // no C1 and no P1
          inP2 == -1 || inL1 == -1 || inL2 == -1)
      {
         config.oflog << "Error: file " << name << " does not contain";
         if(inC1 == -1) config.oflog
            << " C1 (--useCA was" << (config.UseCA ? "" : " not") << " found)";
         if(inL1 == -1) config.oflog << " L1";
         if(inL2 == -1) config.oflog << " L2";
         if(inP1 == -1) config.oflog << " P1";
         if(inP2 == -1) config.oflog << " P2";
         config.oflog << " .. abort." << endl;
         irfstr.clear();
         irfstr.close();
         return 2;
      }
      else if(inP1==-1) {
         inP1 = inC1;
      }

      if(config.UseCA) inP1 = inC1;
      if(inP1 == inC1) UsingCA = true; else UsingCA = false;

         // loop over epochs in the file
      bool first=true;
      int iret;
      RinexObsData rodata;
      while(1) {
         irfstr >> rodata;
         if(irfstr.eof()) break;
         if(irfstr.bad()) {
            config.oflog << "input RINEX stream is bad" << endl;
            break;
         }
         iret = ProcessOneEntireEpoch(rodata);
         if(iret < -1) break;
         if(iret == -1) { iret=0; break; }  // end of file
      }

      irfstr.clear();
      irfstr.close();

      return iret;
   }
   catch(Exception& e) { GPSTK_RETHROW(e); }
   catch(exception& e)
      { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
   catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
}

//------------------------------------------------------------------------------------
// Return :
//      <-1 fatal error,
//       -1 end of file (or past end time limit),
//        0 ok,
//        1 skip this epoch : before begin time
//        2 skip this epoch : comment block,
//        3 skip this epoch : decimated
int ProcessOneEntireEpoch(RinexObsData& roe)
{
   try {
      bool ok;
      int i,j,k;
      double dt;
      string str,datastr;
      RinexSatID sat;
      SatPassData spd;
      RinexObsData::RinexObsTypeMap otmap;
      RinexObsData::RinexSatMap::iterator it;
      RinexObsData::RinexObsTypeMap::const_iterator jt;

         // stay within time limits
      if(roe.time < config.begTime) return 1;
      if(roe.time > config.endTime) return -1;

         // ignore comment blocks ...
      if(roe.epochFlag != 0 && roe.epochFlag != 1) return 2;

         // decimate data
         // if begTime is still undefined, set it to begin of week
      if(config.ith > 0.0) {
         if(fabs(config.begTime-DayTime(DayTime::BEGINNING_OF_TIME)) < 1.e-8)
            config.begTime =
               config.begTime.setGPSfullweek(roe.time.GPSfullweek(),0.0);
         double dt=fabs(roe.time - config.begTime);
         dt -= config.ith*long(0.5+dt/config.ith);
         if(fabs(dt) > 0.25) return 3;            // TD set tolerance? clock bias?
      }

         // save current time
      CurrEpoch = roe.time;