rescor.cpp

GPS数据预处理软件

                  have = havefmt & havepat;
                  break;
               }
               fmtT = fmtP = lineT = lineP = string("");
               for(i=0; i<numWords(line); i++) {
                  word = words(line,i,1);
                  fword = words(format,i,1);
                  if(pattern[i] == 'X') continue;
                  else if(pattern[i] == 'T') {
                     lineT += string(" ") + word;
                     fmtT += string(" ") + fword;
                  }
                  else if(pattern[i] == 'P') {
                     lineP += string(" ") + word;
                     fmtP += string(" ") + fword;
                  }
               }
               try {
                  timetag.setToString(lineT,fmtT);
               }
               catch(Exception& dte) {
                  logof << "ERROR: reading the receiver position flat file threw"
                     << " a DayTime exception:\n"
                     << "  This is the time format: " << fmtT << endl;
                  ok = have = havefmt = false;
                  break;
               }
               try {
                  pos.setToString(lineP,fmtP);
                  pos.transformTo(Position::Cartesian);
                  CurrRef.RxPos = pos;
               }
               catch(Exception& ge) {
                  logof << "ERROR: reading the receiver position flat file threw"
                     << " a Position exception:\n"
                     << "  This is the position format: " << fmtP << endl;
                  ok = have = havefmt = havepat = false;
               }
               if(ok) {
                  if(Debug)logof << "Result: t= " << timetag << " p= " << pos << endl;
                  RefPosMap[timetag] = CurrRef;
                  CurrRef.valid = true;
               }
               break;
            }
            if(!ok) break;
            getline(inf,line);
            stripTrailing(line,'\r');
         }
         inf.close();
         if(!have) {
            logof << "ERROR in reading receiver position file: ";
            if(!havefmt) logof << "format ";
            if(!havepat) {
               if(!havefmt) logof << "and pattern ";
               else logof << "pattern ";
            }
            logof << ((havepat || havefmt) ? "was " : "were ")
               << "wrong or not found!\n";
            logof << RxhelpString << endl;
            logof << "  [The input format is " << format << "]" << endl;
            logof << "  [The input pattern is " << pattern << "]" << endl;
            return -2;
         }
         inPS = 1;
      }  // end flat file input

      // compute the nominal time spacing of the map
      {
         const int ndtmax=15;
         double dt,bestdt[ndtmax];
         int j,k,nleast,ndt[ndtmax]={-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};
         DayTime prev(DayTime::BEGINNING_OF_TIME);
         map<DayTime,RefPosData>::const_iterator it;

         if(Debug) logof << "Here is the reference position map\n";
         for(it=RefPosMap.begin(); it != RefPosMap.end(); it++) {
            if(Debug) logof << "   " << it->first << " " << fixed
                  << " " << setw(13) << setprecision(3) << it->second.RxPos.X()
                  << " " << setw(13) << setprecision(3) << it->second.RxPos.Y()
                  << " " << setw(13) << setprecision(3) << it->second.RxPos.Z()
                  << endl;
            if(prev != DayTime::BEGINNING_OF_TIME) {
               dt = it->first - prev;
               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; nleast = ndt[k];
                     for(j=1; j<ndtmax; j++) if(ndt[j] <= nleast) {
                        k=j; nleast=ndt[j];
                     }
                     ndt[k]=1; bestdt[k]=dt;
                  }
               }
            }
            prev = it->first;
         }
         for(i=1,j=0; i<ndtmax; i++) if(ndt[i] > ndt[j]) j=i;
         RefPosMapDT = bestdt[j];
      }

   }  // end non-empty RefPosFile name

   else if(doRAIM) {
      // if(Debug) prsol.Debug = true; // write to cout ...
      prsol.Algebraic = false;
      //prsol.MaxNIterations = PIC.NIter;    // TD add to command line?
      //prsol.Convergence = PIC.Conv;
      // set inPS below, when you know you can do RAIM
      logof << "Reference position will come from RAIM\n";
   }
   else if(RefPosInput) {
      logof << "Reference position will come from the input file\n";
      inPS = 1;
   }
 
      // reset average RAIM solution
   if(headRAIM) {
      ARSX.Reset();
      ARSY.Reset();
      ARSZ.Reset();
   }

      // --------------------------------------------------------------------
      // misc
      // IonoHt used in meters
   IonoHt *= 1000.0;

      // search for SX,Y,Z input and set DoSX flag, also XR,XI,X1,X2 and DoXR
   DoSVX = DoXR = false;
   for(i=0; i<OTstrings.size(); i++) {
      if(OTstrings[i]==string("SX")
            || OTstrings[i]==string("SY")
            || OTstrings[i]==string("SZ")) DoSVX = true;
      if(OTstrings[i]==string("XR") || OTstrings[i]==string("XI")
            || OTstrings[i]==string("X1") || OTstrings[i]==string("X2")) DoXR = true;
   }

   if(DoXR) {
      int j;
      // transformation matrix is constant
      XRM0[0] = alpha+1;      XRM0[1] = -1;      XRM0[2] = 0;     XRM0[3] = 0;
      XRM1[0] = 1;            XRM1[1] = -1;      XRM1[2] = 0;     XRM1[3] = 0;
      XRM2[0] = -alpha-2;     XRM2[1] = 2;       XRM2[2] = alpha; XRM2[3] = 0;
      XRM3[0] = -2*(alpha+1); XRM3[1] = alpha+2; XRM3[2] = 0;     XRM3[3] = alpha;
      for(i=0; i<4; i++) for(j=0; j<4; j++) XRM[i][j] /= alpha;
      if(Debug) {
         logof << "XRM matrix is:\n" << fixed;
         for(i=0; i<4; i++) {
            for(j=0; j<4; j++) {
               logof << " " << setw(20) << setprecision(4) << XRM[i][j];
            }
            logof << endl;
         }
      }
   }

   CurrRef.valid = false;
   if(Debug) logof << "Return from PrepareInput" << endl;

   return 0;
}
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;
}

//------------------------------------------------------------------------------------
// after reading input header and before calling REC.EditHeader (pass input header)
int RCRinexEditor::BeforeEditHeader(const RinexObsHeader& rhin) throw(Exception)
{
try {
   int i;

      // save the header for later use by SaveData and ComputeNewOTs
   rhead = rhin;

      // -----------------------------------------------------------------------
      // get indexes of input obs types, for dependence checking and fast access
   for(i=0; i<rhin.obsTypeList.size(); i++) {
      if(rhin.obsTypeList[i] == RinexObsHeader::convertObsType("C1")) inC1=i;
      if(rhin.obsTypeList[i] == RinexObsHeader::convertObsType("L1")) inL1=i;
      if(rhin.obsTypeList[i] == RinexObsHeader::convertObsType("L2")) inL2=i;
      if(rhin.obsTypeList[i] == RinexObsHeader::convertObsType("P1")) inP1=i;
      if(rhin.obsTypeList[i] == RinexObsHeader::convertObsType("P2")) inP2=i;
      if(rhin.obsTypeList[i] == RinexObsHeader::convertObsType("D1")) inD1=i;
      if(rhin.obsTypeList[i] == RinexObsHeader::convertObsType("D2")) inD2=i;
      if(rhin.obsTypeList[i] == RinexObsHeader::convertObsType("S1")) inS1=i;
      if(rhin.obsTypeList[i] == RinexObsHeader::convertObsType("S2")) inS2=i;
   }

      // redefine inP1 based on inC1, Callow and Cforce
   if(Callow && inC1 > -1 && inP1 == -1) inP1=inC1;
   if(Cforce && inC1 > -1)               inP1=inC1;

      // -----------------------------------------------------------------------
      // Check dependences of input and output data types
      // -----------------------------------------------------------------------
      // check that we can do RAIM
   if(doRAIM) {
      if(inP1>-1 && inP2>-1) inPS=1;
      else {
         ostringstream stst;
         stst << "Error: cannot compute RAIM solution: missing";
         if(inP1 == -1) stst << " P1";
         if(inP2 == -1) stst << " P2";
         if(inEP == -1) stst << " EP";
         stst << "; abort.\n";
         logof << stst.str();
         oferr << stst.str();
         return -2;
      }
   }

      // -----------------------------------------------------------------------
      // Define bit flags for input data types
   unsigned int InputData=0;
   if(Verbose) logof << "Input data:\n";
   if(inP1 > -1) {
      InputData |= 0x08;
      if(Verbose) logof << " P1(" << inP1 << ")";
   }
   if(inP2 > -1) {
      InputData |= 0x10;
      if(Verbose) logof << " P2(" << inP2 << ")";
   }
   if(inL1 > -1) {
      InputData |= 0x02;
      if(Verbose) logof << " L1(" << inL1 << ")";
   }
   if(inL2 > -1) {
      InputData |= 0x04;
      if(Verbose) logof << " L2(" << inL2 << ")";
   }
   if(inEP > -1) {
      InputData |= RinexObsHeader::RinexObsType::EPdepend;
      if(Verbose) logof << " EP";
   }
   if(inPS > -1) {
      InputData |= 0x40;
      if(Verbose) logof << " PS";
   }
   if(Verbose) logof << "(" << hex << InputData << ")" << dec << endl;

      // -----------------------------------------------------------------------
      // create list OTList of RinexObsTypes here, for use later
      // check dependencies of requested output OTs
   if(Verbose) logof << "Here is the list of added OTs:";
   for(i=0; i<OTstrings.size(); i++) {
      if(Verbose) logof << " " << OTstrings[i];
      OTList.push_back(RinexObsHeader::convertObsType(OTstrings[i]));
   }
   if(Verbose) logof << endl;
   bool ok=true;
   for(i=0; i<OTList.size(); i++) {
      if((InputData & OTList[i].depend) != OTList[i].depend) {
         ostringstream stst;
         ok = false;
         stst << "ResCor Error: Abort: Output OT " << OTstrings[i]
            << " requires missing input:";
         unsigned int test=(InputData & OTList[i].depend);
         test ^= OTList[i].depend;
         if(test & rhin.obsTypeList[inL1].depend) stst << " L1";
         if(test & rhin.obsTypeList[inL2].depend) stst << " L2";
         if(test & rhin.obsTypeList[inP1].depend) stst << " P1";
         if(test & rhin.obsTypeList[inP2].depend) stst << " P2";
         if(test & RinexObsHeader::RinexObsType::EPdepend) stst << " EP";
         if(test & RinexObsHeader::RinexObsType::PSdepend) stst << " PS";
         stst << endl;
         logof << stst.str();
         oferr << stst.str();
      }
   }
   if(!ok) return -3;

   return 0;
}
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); }
}

//------------------------------------------------------------------------------------
// after calling REC.EditHeader (pass output header)
int RCRinexEditor::AfterEditHeader(const RinexObsHeader& rhout) throw(Exception)
{
try {
   int i,j;

      // save header for later use by SaveData
   rheadout = rhout;

      // -----------------------------------------------------------------------
      // define indexes of raw data in output header
   for(i=0; i<rhout.obsTypeList.size(); i++) {
      if(rhout.obsTypeList[i] == RinexObsHeader::convertObsType("C1")) otC1=i;
      if(rhout.obsTypeList[i] == RinexObsHeader::convertObsType("L1")) otL1=i;
      if(rhout.obsTypeList[i] == RinexObsHeader::convertObsType("L2")) otL2=i;
      if(rhout.obsTypeList[i] == RinexObsHeader::convertObsType("P1")) otP1=i;
      if(rhout.obsTypeList[i] == RinexObsHeader::convertObsType("P2")) otP2=i;
      if(rhout.obsTypeList[i] == RinexObsHeader::convertObsType("D1")) otD1=i;
      if(rhout.obsTypeList[i] == RinexObsHeader::convertObsType("D2")) otD2=i;
      if(rhout.obsTypeList[i] == RinexObsHeader::convertObsType("S1")) otS1=i;
      if(rhout.obsTypeList[i] == RinexObsHeader::convertObsType("S2")) otS2=i;
   }

      // redefine otP1 based on otC1, inP1, inC1, Callow and Cforce
   if(Callow && otC1 > -1 && inC1 > -1 && inP1 == -1) otP1=otC1;
   if(Cforce && otC1 > -1)                            otP1=otC1;

      // -----------------------------------------------------------------------
      // create a list of indexes parallel to OTstrings and OTList
   for(j=0; j<OTList.size(); j++) {
      for(i=0; i<rhout.obsTypeList.size(); i++) {
         if(rhout.obsTypeList[i] == OTList[j]) OTindex.push_back(i);
      }
   }

   return 0;
}
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); }
}

//------------------------------------------------------------------------------------
// after reading input obs and before calling EditObs (pass input obs)
int RCRinexEditor::BeforeEditObs(const RinexObsData& roin) throw(Exception)
{
try {
   if(Debug) logof << "\n----------------------------- " << roin.time
      << " ------------------------" << endl;

   // -----------------------------------------------------------------------
   // in-line header info
   // note that often these have a bad (all zeros) epoch