tecmaps.cpp

根据GPS观测数据

            Stations[nfile].getNext = true;
         }
         else Stations[nfile].getNext = false;
      }

         // compute map(s) and output
      nepochs++;
      if(ngood > 0) {
         oflog << ngood << " data at epoch "
            << EarliestTime.printf("%Y/%m/%d %H:%M:%6.3f=%F/%10.3g")
            << ", file #" << nepochs << "." << endl;
            // compute the map(s)
         if(doVTECmap) {
            vtecmap.ComputeMap(EarliestTime,AllObs);
            OutputMapToFile(vtecmap,BaseName,EarliestTime,nepochs);
         }
         if(doMUFmap) {
            mufmap.ComputeMap(EarliestTime,AllObs);
            OutputMapToFile(mufmap,BaseName+string(".MUF"),EarliestTime,nepochs);
         }
         if(doF0F2map) {
            f0f2map.ComputeMap(EarliestTime,AllObs);
            OutputMapToFile(f0f2map,BaseName+string(".F0F2"),EarliestTime,nepochs);
         }
      }
      else oflog << "0 data at epoch "
            << EarliestTime.printf("%Y/%m/%d %H:%M:%6.3f=%F/%10.3g")
            << ", file #" << nepochs << "." << endl;

   } while(1);

      // finished...close all files
   for(nfile=0; nfile<Stations.size(); nfile++) instream[nfile].close();

   if(verbose)
      oflog << endl << "Processed " << Stations.size() << " stations\n";
}
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); }
}

//------------------------------------------------------------------------------------
// output the grid to a file
void OutputGridToFile(VTECMap& vmap, string filename) throw(Exception)
{
try {
   ofstream ofs(filename.c_str());
   if(!ofs) {
      cerr << "Failed to open grid output file " << filename << endl;
      oflog << "Failed to open grid output file " << filename << endl;
   }
   else {
      vmap.OutputGrid(ofs);
      ofs.close();
   }
}
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); }
}

//------------------------------------------------------------------------------------
// output map
void OutputMapToFile(VTECMap& vtmap, string filename, DayTime t, int n) throw(Exception)
{
try {
      // make this a function, pass it the name MUF etc, map and time
   ostringstream oss;
   oss << filename << "." << setw(4) << setfill('0') << n;
   string fn = oss.str();
   ofstream ofs(fn.c_str());
   if(!ofs) {
      cerr << "Failed to open map output file " << fn << " at epoch "
         << t.printf("%Y/%m/%d %H:%M:%6.3f=%F/%10.3g") << endl;
      oflog << "Failed to open map output file " << fn << " at epoch "
         << t.printf("%Y/%m/%d %H:%M:%6.3f=%F/%10.3g") << endl;
   }
   else {
      oflog << "Output map at epoch "
         << t.printf("%Y/%m/%d %H:%M:%6.3f=%F/%10.3g")
         << " to file " << fn << endl;
      vtmap.OutputMap(ofs,GnuplotFormat);
      ofs.close();
   }
}
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); }
}

//------------------------------------------------------------------------------------
void AddStation(string& name) throw(Exception)
{
try {
   Station s;
   s.filename = name;
   for(int i=1; i<33; i++) {
      RinexSatID p(i,SatID::systemGPS);
      s.InitTime[p] = DayTime::BEGINNING_OF_TIME;
   }
   Stations.push_back(s);
}
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,
int ProcessHeader(Station& S) throw(Exception)
{
try {
      // input header
   try {
      instream[S.nfile] >> S.header;
   }
   catch(gpstk::FFStreamError& e) {
      cerr << "Caught an FFStreamError while reading header for file "
         << S.filename << ":\n" << e.getText(0) << endl;
      oflog << "Caught an FFStreamError while reading header for file "
         << S.filename << ":\n" << e.getText(0) << endl;
      return -3;
   }
   catch(gpstk::Exception& e) {
      cerr << "Caught a gpstk exception while reading header for file "
         << S.filename << ":\n" << e.getText(0) << endl;
      oflog << "Caught a gpstk exception while reading header for file "
         << S.filename << ":\n" << e.getText(0) << endl;
      return -4;
   }
   S.nread = 0;

      // convert to geocentric LLH
   S.xyz.setECEF(S.header.antennaPosition);
   S.llr = S.xyz;
   S.llr.transformTo(Position::Geocentric);

      // save station info
   //StationName = head.markerName;
   //TotalSpan = head.lastObs.MJD()-head.firstObs.MJD();

      // dump header information
   if(verbose) {
      int i;
      oflog << "File name: " << S.filename << "  ";
      oflog << "Marker name: " << S.header.markerName << "\n";
      oflog << "Antenna position:    " << setprecision(3) << S.header.antennaPosition
         << endl;
      oflog << "ECEF Position:       " << S.xyz << endl;
      oflog << "Geocentric Position: " << S.llr << endl;
      oflog << "Observation types (" << S.header.obsTypeList.size() << ") :";
      for(i=0; i<S.header.obsTypeList.size(); i++)
         oflog << " " << RinexObsHeader::convertObsType(S.header.obsTypeList[i]);
      oflog << endl;
      oflog << "Time of first obs "
         << S.header.firstObs.printf("%04Y/%02m/%02d %02H:%02M:%010.7f")
         << " " << (S.header.firstSystem.system==RinexSatID::systemGlonass?"GLO":
                   (S.header.firstSystem.system==RinexSatID::systemGalileo?"GAL":"GPS"))
         << endl;
      oflog << "Time of  last obs "
         << S.header.lastObs.printf("%04Y/%02m/%02d %02H:%02M:%010.7f")
         << " " << (S.header.lastSystem.system==RinexSatID::systemGlonass?"GLO":
                   (S.header.lastSystem.system==RinexSatID::systemGalileo?"GAL":"GPS"))
         << 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
//        1 ok but no data read
//       -1 EOF or non-existant
//       -3 FFStream exception,
//       -4 gpstk exception,
//       -6 read error
int ReadNextObs(Station& S) throw(Exception)
{
try {
   if(S.nread == -1) return -1;
   if(!S.getNext) return 1;

   do {
      try {
         instream[S.nfile] >> S.robs;
      }
      catch(gpstk::FFStreamError& e) {
         cerr << "Caught an FFStreamError while reading obs for file "
            << S.filename << ":\n" << e.getText(0) << endl;
         oflog << "Caught an FFStreamError while reading obs for file "
            << S.filename << ":\n" << e.getText(0) << endl;
         return -3;
      }
      catch(gpstk::Exception& e) {
         cerr << "Caught a gpstk exception while reading obs for file "
            << S.filename << ":\n" << e.getText(0) << endl;
         oflog << "Caught a gpstk exception while reading obs for file "
            << S.filename << ":\n" << e.getText(0) << endl;
         return -4;
      }
      if(instream[S.nfile].eof()) {
         oflog << "End of file: " << S.filename << endl;
         return -1;
      }
      if(!instream[S.nfile].good()) {
         oflog << "Read error on file: " << S.filename << endl;
         return -6;
      }

   } while(S.robs.epochFlag != 0 && S.robs.epochFlag != 1);

   S.nread++;
   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 if data is good and was accepted
int ProcessObs(Station& S, vector<ObsData>& obsvec) throw(Exception)
{
try {
   int i,k,n;
   double EL,AZ,LA,LO,SR,VR,TP,bias,obliq;
   double la,lo;     // TEMP
   RinexSatID sat;
   RinexObsData::RinexSatMap::const_iterator it;
   RinexObsData::RinexObsTypeMap::const_iterator jt;
   Position IPP;

   //S.robs.dump(oflog);
      // loop over sat=it->first, ObsTypeMap=it->second
   for(n=0,it=S.robs.obs.begin(); it!=S.robs.obs.end(); ++it) {
      ObsData od;
      sat = it->first;
      if(sat.system != SatID::systemGPS) continue; // ignore non-GPS satellites
      k = -1;
      for(i=0; 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;
   
         // save first time
      if(S.InitTime[sat] == DayTime::BEGINNING_OF_TIME) {
         S.InitTime[sat] = S.robs.time;
      }
      
         // process this sat
      if((jt=it->second.find(ELot)) != it->second.end())
         EL = jt->second.data;
      else continue;
      if(EL < vtecmap.MinElevation) continue;   // here or inside class?
  
      if((jt=it->second.find(AZot)) != it->second.end())
         AZ = jt->second.data;
      else continue;
  
      //TEMP
      if((jt=it->second.find(LAot)) != it->second.end())
         la = jt->second.data;
      else lo = -999.0;
 
      //TEMP
      if((jt=it->second.find(LOot)) != it->second.end())
         lo = jt->second.data;
      else lo = -999.0;

      if((jt=it->second.find(SRot)) != it->second.end())
         SR = jt->second.data;
      else SR = -1.0;

      if((jt=it->second.find(VRot)) != it->second.end())
         VR = jt->second.data;
      else VR = -1.0;

      if(SR == -1.0 && VR == -1.0) continue;

      if((jt=it->second.find(TPot)) != it->second.end())
         TP = jt->second.data;
      else TP = -1.0;

         // compute the pierce point
      if(la == -999.0 || lo == -999.0) {
         IPP = S.llr.getIonosphericPiercePoint(EL,AZ,IonoHt*1000);
         LA = IPP.geocentricLatitude();
         LO = IPP.longitude();
      }
      else { LA=la; LO=lo; }

      od.elevation = EL;
      od.azimuth = AZ;
      od.latitude = LA;
      od.longitude = LO;

      if(TP != -1.0) od.AcqTime = TP;
      else od.AcqTime = S.robs.time - S.InitTime[sat];

      if(od.AcqTime < MinAcqTime) continue;

         // get the bias
      map<string,map<RinexSatID,double> >::const_iterator jt;
      jt = BiasMap.find(S.header.markerName);
         // skip sat+rx for which there are no biases
      if(jt == BiasMap.end()) continue;
      map<RinexSatID,double>::const_iterator kt;
      kt = jt->second.find(sat);
      if(kt == jt->second.end()) continue;
      bias = kt->second;
      if(debug) oflog << "Apply bias for station " << S.header.markerName
         << " and sat " << sat << " = " << fixed << setw(12) << setprecision(6)
         << bias << endl;

         // compute the obliquity
      obliq = WGS84.a()*cos(EL*DEG_TO_RAD)/(WGS84.a()+IonoHt*1000);
      obliq = SQRT(1.0-obliq*obliq);

      if(VR != -1.0) {
         od.VTEC = VR - bias*obliq;
      }
      else {
         od.VTEC = (SR - bias)*obliq;
      }

         // compute the error and save results
      od.VTECerror = vtecmap.VTECError(od.AcqTime, EL, od.VTEC);
      n++;
      obsvec.push_back(od);

         // write out
      /* */
      oflog <<        setw(4) << S.robs.time.GPSfullweek();
      oflog << " " << setw(8) << setprecision(1) << S.robs.time.GPSsow();
      oflog << " " << setw(2) << n;
      oflog << " " << setw(9) << setprecision(5) << LA; // latitude
      oflog << " " << setw(10) << setprecision(5)<< LO; // longitude
      oflog << " " << setw(4) << setprecision(2) << obliq; // obliquity
      oflog << " " << setw(8) << setprecision(3) << od.VTEC; // vertical TEC
      oflog << " " << setw(8) << setprecision(3) << od.AcqTime; // acquisition time
      oflog << " " << setw(2) << sat.id; // PRN
      oflog << " " << setw(3) << S.nfile+1; // file number
      oflog << endl;
      /* */

   }  // end for loop over sats

   if(n>0) return n;
   else return -1;
}
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); }
}

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------