earthorientation.cpp

linux的gps应用

         << setw(10) << setprecision(6) << eopp.K1
         << setw(10) << setprecision(6) << eopp.K2
         << setw(10) << setprecision(6) << eopp.K3
         << setw(10) << setprecision(6) << eopp.K4
         << "          " << endl;
      os << setw(10) << setprecision(6) << eopp.L1
         << setw(10) << setprecision(6) << eopp.L2
         << setw(10) << setprecision(6) << eopp.L3
         << setw(10) << setprecision(6) << eopp.L4
         << setw( 9) << setprecision(4) << eopp.R1
         << setw( 9) << setprecision(4) << eopp.R2
         << setw( 9) << setprecision(4) << eopp.R3
         << setw( 9) << setprecision(4) << eopp.R4
         << "    " << endl;
      os << setw(4) << eopp.TAIUTC
         << setw(5) << eopp.SerialNo
         << setw(6) << int(eopp.tv+0.5)
         << " " << eopp.Info
         << "                    "
         << "                    "
         << "      ";
      return os;
   }

   //---------------------------------------------------------------------------------
   // class EOPStore
   //---------------------------------------------------------------------------------
   // Add to the store directly -- not recommended,
   // use the form that takes EOPPrediction
   void EOPStore::addEOP(int mjd, EarthOrientation& eop)
      throw()
   {
      mapMJD_EOP[mjd] = eop;

      if(begMJD == -1 || endMJD == -1) {
         begMJD = endMJD = mjd;
      }
      else if(mjd < begMJD) {
         begMJD = mjd;
      }
      else if(mjd > endMJD) {
         endMJD = mjd;
      }
   }

   //---------------------------------------------------------------------------------
   // Add to the store by computing using an EOPPrediction,
   // this is the usual way.
   // @param MJD integer MJD at which to add EOPs
   // @return non-0 if MJD is outside range
   int EOPStore::addEOP(int mjd, EOPPrediction& eopp)
      throw(DayTime::DayTimeException)
   {
      EarthOrientation eo;
      try {
         eo = eopp.computeEOP(mjd);
      }
      catch(DayTime::DayTimeException& dte)
      {
         GPSTK_RETHROW(dte);
      }

      addEOP(mjd,eo);

      return 0;
   }

   //---------------------------------------------------------------------------------
   // Add EOPs to the store via an inpu file: either an EOPP file
   // or a flat file produced by USNO (see http://maia.usno.navy.mil/
   // and get either file 'finals.data' or finals2000A.data').
   // @param filename Name of file to read, including path.
   void EOPStore::addFile(const string& filename)
      throw(FileMissingException)
   {
      try {
         addEOPPFile(filename);
      }
      catch(FileMissingException& fme)
      {
         if(StringUtils::matches(fme.getText(),string("wrong format")).empty()) {
            GPSTK_RETHROW(fme);
         }

         // try other format
         try {
            addIERSFile(filename);
         }
         catch(FileMissingException& fme)
         {
            GPSTK_RETHROW(fme);
         }
      }
   }

   //---------------------------------------------------------------------------------
   // Add EOPs to the store via an EOPP file:
   // read the EOPPrediction from the file and then compute EOPs
   // for all days within the valid range.
   // @param filename Name of file to read, including path.
   void EOPStore::addEOPPFile(const string& filename)
      throw(FileMissingException)
   {
      // read the file into an EOPPrediction
      EOPPrediction eopp;
      try {
         eopp.loadFile(filename);
      }
      catch(FileMissingException& fme)
      {
         GPSTK_RETHROW(fme);
      }

      // pull out the beginning of the valid time range
      int mjd;
      mjd = eopp.getValidTime();
      // add all 7 days
      for(int i=0; i<7; i++) {
         EarthOrientation eo;
         eo = eopp.computeEOP(mjd);
         addEOP(mjd,eo);
         mjd++;
      }
   }

   //---------------------------------------------------------------------------------
   // see http://maia.usno.navy.mil/readme.finals
   void EOPStore::addIERSFile(const string& filename)
      throw(FileMissingException)
   {
      bool ok;
      int n,mjd;
      double fracmjd;
      string line,word;

      ifstream inpf(filename.c_str());
      if(!inpf) {
         FileMissingException fme("Could not open IERS file " + filename);
         GPSTK_THROW(fme);
      }

      ok = true;
      while(!inpf.eof() && inpf.good()) {
         getline(inpf,line);
         StringUtils::stripTrailing(line,'\r');
         if(inpf.eof()) break;
            // line length is actually 187
         if(inpf.bad() || line.size() < 70) { ok = false; break; }
         EarthOrientation eo;
         mjd = StringUtils::asInt(line.substr(7,5));
         eo.xp = StringUtils::asDouble(line.substr(18,9));       // arcseconds
         eo.yp = StringUtils::asDouble(line.substr(37,9));       // arcseconds
         eo.UT1mUTC = StringUtils::asDouble(line.substr(58,10)); // seconds

         addEOP(mjd,eo);
      };
      inpf.close();

      if(!ok) {
         FileMissingException fme("IERS File " + filename
            + " is corrupted or wrong format");
         GPSTK_THROW(fme);
      }
   }

   //---------------------------------------------------------------------------------
   // Edit the store by deleting all entries before(after)
   //  the given min(max) MJDs. If mjdmin is later than mjdmax,
   //  the two times are switched.
   //  @param mjdmin integer MJD desired earliest store time.
   //  @param mjdmax integer MJD desired latest store time.
   void EOPStore::edit(int mjdmin, int mjdmax)
      throw()
   {
      if(mjdmin > mjdmax) {
         int m=mjdmin;
         mjdmin = mjdmax;
         mjdmax = m;
      }

      map<int,EarthOrientation>::iterator it;
      it = mapMJD_EOP.lower_bound(mjdmin);
      if(it != mapMJD_EOP.begin())
         mapMJD_EOP.erase(mapMJD_EOP.begin(), it);

      it = mapMJD_EOP.upper_bound(mjdmax);
      if(it != mapMJD_EOP.end())
         mapMJD_EOP.erase(it, mapMJD_EOP.end());

      it = mapMJD_EOP.begin();
      if(it == mapMJD_EOP.end())
         begMJD = -1;
      else 
         begMJD = it->first;

      it = mapMJD_EOP.end();
      if(--it == mapMJD_EOP.end())
         endMJD = -1;
      else 
         endMJD = it->first;
   }

   //---------------------------------------------------------------------------------
   // Dump the store to cout.
   // @param detail determines how much detail to include in the output
   //   0 start and stop times (MJD), and number of EOPs.
   //   1 list of all times and EOPs.
   void EOPStore::dump(short detail, ostream& os) const
      throw()
   {
      DayTime tt;
      os << "EOPStore dump (" << mapMJD_EOP.size() << " entries):\n";
      os << " Time limits: [MJD " << begMJD << " - " << endMJD << "]";
      tt.setMJD(double(begMJD));
      os << " = [m/d/y " << tt.printf("%m/%d/%Y");
      tt.setMJD(double(endMJD));
      os << " - " << tt.printf("%m/%d/%Y") << "]" << endl;
      if(detail > 0) {
         int lastmjd=-1;
         map<int,EarthOrientation>::const_iterator it;
         for(it=mapMJD_EOP.begin(); it != mapMJD_EOP.end(); it++) {
            if(lastmjd != -1 && it->first - lastmjd > 1)
               os << " ....." << endl;
            os << " " << it->first << " " << it->second
               << "     (" << setfill('0') << setw(3)
               << EOPPrediction::getSerialNumber(it->first) << setfill(' ') << ")"
               << endl;
            lastmjd = it->first;
         }
      }
   }

   //---------------------------------------------------------------------------------
   // Get the EOP at the given epoch and return it.
   // @param t DayTime at which to compute the EOPs.
   // @return EarthOrientation EOPs at time t.
   // @throw InvalidRequest if the (int) MJDs on either side of t
   // cannot be found in the map.
   EarthOrientation EOPStore::getEOP(DayTime& t) const
      throw(InvalidRequest)
   {
         // find the EOs before and after epoch
      int loMJD = int(t.MJD());
      int hiMJD = loMJD + 1;
         // find these EOPs
      map<int,EarthOrientation>::const_iterator hit,lit;
      lit = mapMJD_EOP.find(loMJD);
      hit = mapMJD_EOP.find(hiMJD);
      if(lit == mapMJD_EOP.end() || hit == mapMJD_EOP.end()) {
         InvalidRequest ire(string("Time tag (MJD=")
         + (lit == mapMJD_EOP.end() ?
               StringUtils::asString(loMJD) : StringUtils::asString(hiMJD))
         + string(") not found within the EOP store - EOPP files are out-of-date"));
         GPSTK_THROW(ire);
      }
         // linearly interpolate to get EOP at the desired time
      EarthOrientation eo;
      double dt=t.MJD()-double(loMJD);
      eo.xp = (1.0-dt) * lit->second.xp + dt * hit->second.xp;
      eo.yp = (1.0-dt) * lit->second.yp + dt * hit->second.yp;
      eo.UT1mUTC = (1.0-dt) * lit->second.UT1mUTC + dt * hit->second.UT1mUTC;

      return eo;
   }

} // end namespace gpstk
//------------------------------------------------------------------------------------