ionobias.cpp

根据GPS观测数据

   }
   else if(found_log_file || (arg[0]=='-' && arg[1]=='l')) {
      LogFile = string(arg);
      if(!found_log_file) LogFile.erase(0,2); else found_log_file = false;
   }
   else if(string(arg) == "--log")
      found_log_file = true;
   else if((arg[0]=='-' && arg[1]=='d') || string(arg)==string("--debug"))
      debug = true;
   else if((arg[0]=='-' && arg[1]=='v') || string(arg)==string("--verbose"))
      verbose = true;
   else if(string(arg) == "--file" || string(arg) == "-f")
      found_cfg_file = true;
   else Args.push_back(arg);
}
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); }
}


//------------------------------------------------------------------------------------
int Initialize(void) throw(Exception)
{
try {
   static SP3EphemerisStore SP3EphList;
   static BCEphemerisStore BCEphList;

      // open nav files and read EphemerisStore
   if(!NavDir.empty())
      for(int i=0; i<NavFiles.size(); i++)
         NavFiles[i] = NavDir + "/" + NavFiles[i];
   FillEphemerisStore(NavFiles, SP3EphList, BCEphList);
   if(SP3EphList.size()) {
      if(verbose) SP3EphList.dump(0,oflog);
   }
   else if(verbose) oflog << "SP3 Ephemeris list is empty\n";

   if(BCEphList.size()) {
      BCEphList.SearchNear();
      if(verbose) BCEphList.dump(0,oflog);
   }
   else if(verbose) oflog << "BC Ephemeris list is empty\n";

   if(SP3EphList.size()) pEph = &SP3EphList;
   else if(BCEphList.size()) pEph = &BCEphList;
   else {
      cerr << "IonoBias abort -- no ephemeris\n";
      oflog << "IonoBias abort -- no ephemeris\n";
      return -1;
   }

      // create the obs types for later use
   ELot = RinexObsHeader::convertObsType("EL");
   LAot = RinexObsHeader::convertObsType("LA");
   LOot = RinexObsHeader::convertObsType("LO");
   SRot = RinexObsHeader::convertObsType("SR");
   SSot = RinexObsHeader::convertObsType("SS");

      // initialize AT header data
   int i;
   NgoodStations = 0;
   for(i=0; i<=MAXPRN; i++) BoolVec.push_back(false);
   for(i=0; i<Filenames.size(); i++) EstimationFlag.push_back(BoolVec);

      // open output file and write zero-filled header
   fout.open(ATFileName.c_str(),ios_base::out);
   if(!fout) {
      cerr << "IonoBias abort -- failed to open AT file " << ATFileName
         << " for output." << endl;
      oflog << "IonoBias abort -- failed to open AT file " << ATFileName
         << " for output." << endl;
      return -2;
   }
   WriteATHeader();

   FoundMinLat = 90;
   FoundMinLon = 360;
   FoundMaxLat = -90;
   FoundMaxLon = 0;

   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 0 ok,
//       -2 could not open a file,
//       -3 FFStream exception,
//       -4 gpstk exception,
//       -5 no sunrise
int Process(void) throw(Exception)
{
try {
   int i,iret;
   string fname;
   //RinexObsStream instream;
   RinexObsHeader header;

      // loop over input file names
   if(verbose) oflog << "\nProcess " << Filenames.size() << " input files:\n";
   for(ndata=0,nfile=0; nfile<Filenames.size(); nfile++) {
      if(verbose) oflog << endl;
      fname = Filenames[nfile];
      //instream.clear();
      //instream.open(fname.c_str(),ios_base::in);
      RinexObsStream instream(fname.c_str(),ios_base::in);
      if(!instream) {
         oflog << " Rinex file " << fname << " could not be opened -- abort.\n";
         return -2;
      }
      instream.exceptions(ios::failbit);
      if(verbose) oflog << "Opened input file #" << nfile+1 << ": " << fname << endl;
      
      iret = ProcessHeader(instream,fname,header);
      if(iret != 0) return iret;

      if(nfile==0) {
         MJDNorm = header.firstObs.MJD();
         LonNorm = StationPosition[1]; //.getLongitude();
      }

      iret = ProcessObs(instream,fname,header);
      if(iret != 0) return iret;

      instream.close();

   }  // end loop over file names

   if(verbose) {
      oflog << endl << "Processed " << Filenames.size()
      << " files; " << NgoodStations << " of them had good data.\n";
      oflog << "Total number of data points = " << ndata << endl;
      oflog << "Found " << fixed << setprecision(2)
         << FoundMinLat << " <= raw Lat <= " << FoundMaxLat << " and "
         << FoundMinLon << " <= raw Lon <= " << FoundMaxLon << 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 0 ok, -3 FFStream exception, -4 gpstk exception, -5 no sunrise
int ProcessHeader(RinexObsStream& ins, string& filename, RinexObsHeader& head)
   throw(Exception)
{
try {
      // input header
   try {
      ins >> head;
   }
   catch(gpstk::FFStreamError& e) {
      cerr << "Caught an FFStreamError while reading header for file "
         << filename << ":\n" << e.getText(0) << endl;
      oflog << "Caught an FFStreamError while reading header for file "
         << filename << ":\n" << e.getText(0) << endl;
      return -3;
   }
   catch(gpstk::Exception& e) {
      cerr << "Caught a gpstk exception while reading header for file "
         << filename << ":\n" << e.getText(0) << endl;
      oflog << "Caught a gpstk exception while reading header for file "
         << filename << ":\n" << e.getText(0) << endl;
      return -4;
   }

      // convert to LLH
   Position xyz;
   xyz.setECEF(head.antennaPosition);
   StationPosition = xyz;
   StationPosition.transformTo(Position::Geocentric);

      // compute begin and end times
   TimeLimits(StationPosition, head.firstObs.DOY(), TimeSector, begintime, endtime);
   if(begintime == -999 || endtime == -999) return -5;

      // save station info
   StationName = lowerCase(subString(head.markerName,0,4));
   TotalSpan = head.lastObs.MJD()-head.firstObs.MJD();

      // dump header information
   if(verbose) {
      int i;
      oflog << "File name: " << filename << endl;
      oflog << "Marker name: " << head.markerName << "\n";
      oflog << "Position (XYZ,m) : " << fixed
         << setprecision(3) << head.antennaPosition << "\n";
      oflog << "Position (LLH  ) : ("
         << setprecision(8) << StationPosition[0] << ", "
         << setprecision(8) << StationPosition[1] << ", "
         << setprecision(4)
         << StationPosition[2] - StationPosition.radiusEarth()
         << ")\n";
      oflog << "Observation types (" << head.obsTypeList.size() << ") :";
      for(i=0; i<head.obsTypeList.size(); i++)
         oflog << " " << RinexObsHeader::convertObsType(head.obsTypeList[i]);
      oflog << endl;
      oflog << "Time of first obs "
         << head.firstObs.printf("%04Y/%02m/%02d %02H:%02M:%010.7f")
         << " " << (head.firstSystem.system==RinexSatID::systemGlonass?"GLO":
                   (head.firstSystem.system==RinexSatID::systemGalileo?"GAL":"GPS"))
         << endl;
      oflog << "Time of  last obs "
         << head.lastObs.printf("%04Y/%02m/%02d %02H:%02M:%010.7f")
         << " " << (head.lastSystem.system==RinexSatID::systemGlonass?"GLO":
                   (head.lastSystem.system==RinexSatID::systemGalileo?"GAL":"GPS"))
         << endl;
      oflog << "DOY = " << head.firstObs.DOY() << endl;
      oflog << "Sunrise = " << setprecision(2) << sunrise;
      oflog << "  Sunset  = " << setprecision(2) << sunset << endl;
      oflog << "Begin time = " << setprecision(2) << begintime;
      oflog << "  End time = " << setprecision(2) << endtime << 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); }
}

//------------------------------------------------------------------------------------
// Given a position (LLH), day of year, and sector flag (TimeSector),
// compute the begin and end times (hours of the day) of our data window,
// which will = sunrise + TermOffset and sunset - TermOffset.
void TimeLimits(Position llr, int doy, string& sector, double& begin, double& end)
   throw(Exception)
{
try {
   begin = 0;
   end = 24.;

   Sunrise(llr[0], llr[1], IonoHt*1000.0, doy, sunrise, sunset);
   if(sector == string("both")) return;
   if(sector == string("day")) {
      begin = sunrise + TermOffset/60.0;
      end   = sunset  - TermOffset/60.0;
   }
   if(sector == string("night")) {
      begin = sunset  + TermOffset/60.0;
      end   = sunrise - TermOffset/60.0;
   }

   while(begin <  0) begin += 24;
   while(begin >= 24) begin -= 24;
   while(end <  0) end += 24;
   while(end >= 24) end -= 24;
}
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); }
}

//------------------------------------------------------------------------------------
// Compute the position (latitude and longitude, in degrees) of the sun
// given the day of year and the hour of the day.
// Adapted from sunpos by D. Coco 12/15/94
//#include "icd_200_constants.hpp     // for TWO_PI
//#include "geometry.hpp"             // for DEG_TO_RAD and RAD_TO_DEG
void SolarPosition(int doy, double hod, double& lat, double& lon) throw(Exception)
{
try {
   lat = sin(23.5*DEG_TO_RAD)*sin(TWO_PI*double(doy-83)/365.25);
   lat = lat / ::sqrt(1.0-lat*lat);
   lat = RAD_TO_DEG*atan(lat);
   lon = 180.0 - hod*15.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); }
}

//------------------------------------------------------------------------------------
// Compute the time of day of sunrise and sunset (set to -999 if they do not exist),
// given a geographic position and day of year.
// Adapted from D. Coco 2/7/96 from equation in Supplement to the Astromonical Almanac
void Sunrise(double lat, double lon, double ht, int doy, double& rise, double& set)
   throw(Exception)
{
try {
   const double DEG_TO_HRS=(24.0/360.0); // should this be sidereal day?
   double sunlat,sunlon,hod;

   // find the position of the sun
   hod = 0;
   SolarPosition(doy, hod, sunlat, sunlon);

   double x = -1 * tan(sunlat*DEG_TO_RAD) * tan(lat*DEG_TO_RAD);
	if(x <= -1.0 || x >= 1.0) {
      rise = set = -999;
      return;
   }
   x = acos(x) * RAD_TO_DEG;

   rise = DEG_TO_HRS * (sunlon - lon - x);
   set  = DEG_TO_HRS * (sunlon - lon + x);

   // adjust for height above sea level
   double dht=0,radius=Position::radiusEarth(lat,WGS84.a(),WGS84.eccSquared());
   dht = 24.0*acos(radius/(radius+ht))/TWO_PI;
   rise -= dht;
   set += dht;

   while(rise <  0) rise += 24;
   while(rise >= 24) rise -= 24;
   while(set <  0) set += 24;
   while(set >= 24) set -= 24;
}
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 0 ok, -3 FFStream exception, -4 gpstk exception, -6 stream not good
int ProcessObs(RinexObsStream& ins, string& filename, RinexObsHeader& head)
   throw(Exception)
{
try {
   int i,j,k,npts[MAXPRN+1];
   double EL,LA,LO,SR,hours,cr,ob;
   Position LLI;
   DayTime begin[MAXPRN+1],end[MAXPRN+1];
   RinexObsData robs;
   RinexSatID sat;
   //RinexObsData::RinexObsTypeMap otmap;
   RinexObsData::RinexSatMap::const_iterator it;
   RinexObsData::RinexObsTypeMap::const_iterator jt;

   if(!ins.good()) return -6;

      // initialize for this station
   fout.seekp(0,ios_base::end);           // go to EOF
   current_header_pos = fout.tellp();     // save for later
   NgoodPoints = 0;
   WriteStationHeader(NgoodPoints,StationName,StationPosition);   // dummy

      // loop over epochs
   for(i=1; i<=MAXPRN; i++) npts[i]=0;
   do {
      try {
         ins >> robs;
      }
      catch(gpstk::FFStreamError& e) {
         cerr << "Caught an FFStreamError while reading obs for file "
            << filename << ":\n" << e.getText(0) << endl;
         oflog << "Caught an FFStreamError while reading obs for file "
            << filename << ":\n" << e.getText(0) << endl;
         return -3;
      }
      catch(gpstk::Exception& e) {
         cerr << "Caught a gpstk exception while reading obs for file "
            << filename << ":\n" << e.getText(0) << endl;
         oflog << "Caught a gpstk exception while reading obs for file "
            << filename << ":\n" << e.getText(0) << endl;
         return -4;
      }
      if(ins.eof()) break;
      if(!ins.good()) { return -6; }

      //if(verbose) oflog << " Read file " << filename
      //   << " epoch " << robs.time.printf("%Y/%m/%d_%H:%M:%6.3f=%F/%10.3g") << endl;

      hours = robs.time.secOfDay()/3600.0;          // hours of the day
         // compute co-rotating longitude CL = LO + cr
      cr = (robs.time.MJD()-MJDNorm) * 360.0;
      cr -= LonNorm + TotalSpan * 180.0;

         // loop over sat=it->first, ObsTypeMap=it->second
      for(it=robs.obs.begin(); it != robs.obs.end(); ++it) {
         sat = it->first;
         if(sat.system != SatID::systemGPS) continue; // ignore non-GPS satellites
         if(sat.id <= 0 || sat.id > MAXPRN) continue; // just in case...
         for(i=0,k=-1; i<ExSV.size(); i++) {   // Is this satellite excluded ?
            if( ExSV[i] == sat ||                                 // sat is excluded
               (ExSV[i].id==-1 && ExSV[i].system==sat.system) ) {// system excluded
               k=i;
               break;
            }
         }
         if(k != -1) continue;
   
            // process this sat
         if( (jt=it->second.find(ELot)) != it->second.end()) {
            EL = jt->second.data;
            if(EL < MinElevation) continue;
         }
         //else ...

         if( (jt=it->second.find(LAot)) != it->second.end()) {
            LA = jt->second.data;
            if(LA < FoundMinLat) FoundMinLat = LA;
            if(LA > FoundMaxLat) FoundMaxLat = LA;
            if(LA < MinLatitude || LA > MaxLatitude) continue;
         }

         if( (jt=it->second.find(LOot)) != it->second.end()) {
            LO = jt->second.data;
            while(LO < 0.0) LO+=360.0;
            if(LO < FoundMinLon) FoundMinLon = LO;
            if(LO > FoundMaxLon) FoundMaxLon = LO;
            if(LO < MinLongitude || LO > MaxLongitude) continue;
         }

         if( (jt=it->second.find(SRot)) != it->second.end()) {
            SR = jt->second.data;
            //if(jt->second.ssi == 1) continue;    // reject if ssi==1 ?? TD
         }
         else if( (jt=it->second.find(SSot)) != it->second.end()) {
            SR = jt->second.data;