ionobias.cpp

根据GPS观测数据

#pragma ident "$Id: IonoBias.cpp 584 2007-06-06 14:16:05Z btolman $"

//============================================================================
//
//  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
//
//============================================================================

/**
 * @file IonoBias.cpp
 * Program IonoBias will estimate satellite and receiver biases and compute
 *   a simple ionospheric model using least squares and slant TEC values
 *   from multiple stations.
 */

//------------------------------------------------------------------------------------
#include "StringUtils.hpp"
#include "DayTime.hpp"
#include "RinexSatID.hpp"
#include "CommandOption.hpp"
#include "CommandOptionWithTimeArg.hpp"
#include "CommandOptionParser.hpp"

#include "RinexObsData.hpp"
#include "RinexObsHeader.hpp"
#include "RinexObsStream.hpp"

#include "Vector.hpp"
#include "Matrix.hpp"

#include "Position.hpp"
#include "WGS84Geoid.hpp"
#include "icd_200_constants.hpp"     // for TWO_PI
#include "geometry.hpp"              // for DEG_TO_RAD and RAD_TO_DEG

#include "RinexUtilities.hpp"

#include <iostream>
#include <time.h>
#include <string>
#include <vector>
#include <map>
#include <utility>      // for pair

using namespace std;
using namespace gpstk;
using namespace StringUtils;

//------------------------------------------------------------------------------------
// TD. Robustness!
// Check for duplicate file names or site names - this makes it singular
// Check for unreasonable TEC values - this also can make it singular
// What is the purpose of max and min longitude?
// look at the german objections
// I don't understand co-rotating longitude - perhaps I can see making the sun stand
// still, but why rotate so middle of each pass is aligned?
//------------------------------------------------------------------------------------
// Max PRN 
const int MAXPRN=32;
// data input from command line
bool verbose,debug;
   // log file
string LogFile;
ofstream oflog;
string Title;
   // output file
string ATFileName,BiasFileName;
ofstream fout;
ios::pos_type current_header_pos;
   // input path
string InputPath;
vector<string> Filenames;
   // excluded satellites
vector<RinexSatID> ExSV;
   // ephemeris
string NavDir;
vector<string> NavFiles;
EphemerisStore *pEph;
   // obs types needed
RinexObsHeader::RinexObsType ELot,LAot,LOot,SRot,SSot;
   // geoid
WGS84Geoid WGS84;
   // Start and stop times
DayTime BegTime,EndTime;

   // processing
int MinPoints;
double MinTimeSpan;  // TD not implemented
double MinElevation;
double MinLatitude,MaxLatitude;
double MinLongitude,MaxLongitude;
double FoundMinLat,FoundMinLon,FoundMaxLat,FoundMaxLon;
string TimeSector;
double TermOffset;
double IonoHt;
   // 
double sunrise,sunset;     // times in hours of day
double begintime,endtime;  // "
   // normalizations
double MJDNorm,LonNorm;
   // data that goes into output file headers
long NgoodStations;
vector<vector<bool> > EstimationFlag;
vector<bool> BoolVec;
   // data per station that goes into AT output file
int nfile;     // current file number (0..Filenames.size()-1)
long NgoodPoints;
double TotalSpan;       // time in days covered by the file
string StationName;
Position StationPosition;    // station position in geographic lat,lon,radius
   // least squares
bool ComputeSatBiases,DoEstimation,SkipPreproc;
string Model("linear");
int NIonoParam,NBiasParam,NTotalParam;
Vector<double> Sol,InfData;
Matrix<double> Cov;
int ndata;
double MaxLat,MinLat,MaxCRLon,MinCRLon,PM[10];
vector<pair<string,int> > ComponentIDs;
map<string,int> mapN;
map<string,string> mapFilename;

//------------------------------------------------------------------------------------
// prototypes
void ConfigureAndDefaults(void) throw(Exception);
int GetCommandLine(int argc, char **argv) throw(Exception);
void PreProcessArgs(const char *arg, vector<string>& Args) throw(Exception);
int Initialize(void) throw(Exception);
int Process(void) throw(Exception);
int ProcessHeader(RinexObsStream& ins, string& filename, RinexObsHeader& head)
   throw(Exception);
void TimeLimits(Position llr, int doy, string& sector, double& begin, double& end)
   throw(Exception);
void SolarPosition(int doy, double hod, double& lat, double& lon) throw(Exception);
void Sunrise(double lat, double lon, double ht, int doy, double& rise, double& set)
   throw(Exception);
int ProcessObs(RinexObsStream& ins, string& filename, RinexObsHeader& head)
   throw(Exception);
void WriteATHeader(void) throw(Exception);
void WriteStationHeader(int npts, string sta_name, Position llr) throw(Exception);
void ParseLine(string& str, vector<string>& wds) throw(Exception);
int ReadATandCompute(void) throw(Exception);
double obliquity(double elevation) throw(Exception);
//void PartialsMatrix(Matrix<double>& P,int index,double lat,double lon,double obq);

//------------------------------------------------------------------------------------
// utility routines
//------------------------------------------------------------------------------------
// find the index of first occurance of item t (of type T) in vector<T> v;
// i.e. v[index]=t  Return -1 if t is not found.
template<class T> int index(const std::vector<T> v, const T& t) 
{
   for(int i=v.size()-1; i>=0; i--) {
      if(v[i] == t) return i;
   }
   return -1;
}

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

   BegTime = DayTime::BEGINNING_OF_TIME;
   EndTime = DayTime::END_OF_TIME;

      // Title description and run time
   CurrEpoch.setLocalTime();
   Title = "IonoBias, built on the GPSTK ToolKit, Ver 1.0 6/25/04, Run ";
   Title += CurrEpoch.printf("%04Y/%02m/%02d %02H:%02M:%02S\n");
   cout << Title;

      // set configuration and default values
   ConfigureAndDefaults();

      // define extended types
   iret = RegisterARLUTExtendedTypes();
   if(iret) goto quit;
   iret = RegisterExtendedRinexObsType("SS","Slant TEC (Phase smoothed)","TECU",0x1E);
   if(iret) goto quit;

      // get command line arguments
   iret = GetCommandLine(argc, argv);
   if(iret) goto quit;

   if(!SkipPreproc) {
         // initialize
      iret = Initialize();
      if(iret) goto quit;

         // process the data
      iret = Process();
      if(iret) goto quit;

         // write the revised header
      WriteATHeader();
      fout.close();
   }

   if(DoEstimation) {
      // read the AT file and compute biases and model
      iret = ReadATandCompute();
   }

quit:
      // compute run time
   totaltime = clock()-totaltime;
   cout << "IonoBias timing: " << fixed << setprecision(3)
      << double(totaltime)/double(CLOCKS_PER_SEC) << " seconds.\n";
   oflog << "\nIonoBias timing: " << fixed << setprecision(3)
      << double(totaltime)/double(CLOCKS_PER_SEC) << " seconds.\n";
   oflog.close();

   return iret;
}
catch(gpstk::Exception& e) {
   cerr << e;
   return 1;
}
catch (...) {
   cerr << "Unknown error.  Abort." << endl;
   return 1;
}
   return 0;
}   // end main()

//------------------------------------------------------------------------------------
// set defaults
void ConfigureAndDefaults(void) throw(Exception)
{
try {
   verbose = false;
   debug = false;
   LogFile = string("IonoBias.log");
   BiasFileName = string("");          // no output

   MinPoints = 0;
   MinTimeSpan = 0.0;      // minutes
   MinElevation = 0.0;
   MinLatitude = -90.0;
   MaxLatitude = 90.0;
   MinLongitude = 0.0;
   MaxLongitude = 360.0;
   TimeSector = string("night");
   TermOffset = 0.0;       // min
   IonoHt = 350.0;         // km

   DoEstimation=true;      // if false, quit after writing the AT file
   SkipPreproc =false;     // if true, assume AT file exists and don't generate it
   ComputeSatBiases=true;  // if true, compute Sat+Rx biases, else Rx biases only
}
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); }
}

//------------------------------------------------------------------------------------
// Define, parse and evaluate command line
int GetCommandLine(int argc, char **argv) throw(Exception)
{
try {
   bool help=false;
   int i,j;
   RinexSatID sat;
   sat.setfill('0');

      // required options

   RequiredOption dashin(CommandOption::hasArgument, CommandOption::stdType,
      0,"input"," --input <file>       Input Rinex obs file name(s)");
   //dashin.setMaxCount(1);

      // optional options

   // this only so it will show up in help page...
   CommandOption dashf(CommandOption::hasArgument, CommandOption::stdType,
      'f',""," -f<file>             file containing more options");

   CommandOption dashp(CommandOption::hasArgument, CommandOption::stdType,
      0,"inputdir"," --inputdir <path>    Path for input file(s)");
   dashp.setMaxCount(1);

   // ephemeris
   CommandOption dashnd(CommandOption::hasArgument, CommandOption::stdType, 0,
      "navdir"," Ephemeris input:\n --navdir <dir>       Path of navigation file(s)");
   dashnd.setMaxCount(1);

   CommandOption dashn(CommandOption::hasArgument, CommandOption::stdType,
      0,"nav"," --nav <file>         Navigation (Rinex Nav OR SP3) file(s)");

   CommandOption dashat(CommandOption::hasArgument, CommandOption::stdType,
      0,"datafile",
      " Output:\n --datafile <file>    Data (AT) file name, for output and/or input");
   dashat.setMaxCount(1);
   
   CommandOption dashl(CommandOption::hasArgument, CommandOption::stdType,
      0,"log"," --log <file>         Output log file name");
   dashl.setMaxCount(1);

   CommandOption dashout(CommandOption::hasArgument, CommandOption::stdType,
      0,"biasout"," --biasout <file>     Output satellite+receiver biases file name");
   dashout.setMaxCount(1);

   // time
   CommandOptionWithTimeArg dasheb(0,"BeginTime","%Y,%m,%d,%H,%M,%f",
      " Time limits:\n --BeginTime <arg>    Start time, arg is of the form "
      "YYYY,MM,DD,HH,Min,Sec");
   CommandOptionWithTimeArg dashgb(0,"BeginGPSTime","%F,%g",
      " --BeginGPSTime <arg> Start time, arg is of the form GPSweek,GPSsow");

   CommandOptionWithTimeArg dashee(0,"EndTime","%Y,%m,%d,%H,%M,%f",
      " --EndTime <arg>      End time, arg is of the form YYYY,MM,DD,HH,Min,Sec");
   CommandOptionWithTimeArg dashge(0,"EndGPSTime","%F,%g",
      " --EndGPSTime <arg>   End time, arg is of the form GPSweek,GPSsow");

   // allow ONLY one start time (use startmutex(true) if one is required)
   CommandOptionMutex startmutex(false);
   startmutex.addOption(&dasheb);
   startmutex.addOption(&dashgb);
   CommandOptionMutex stopmutex(false);
   stopmutex.addOption(&dashee);
   stopmutex.addOption(&dashge);

   CommandOptionNoArg dashde(0, "NoEstimation"," Processing:\n"
      " --NoEstimation       Do NOT perform the estimation (default=false).");

   CommandOptionNoArg dashwo(0, "NoPreprocess",
      " --NoPreprocess       Skip preprocessing; read (existing) AT file "
      "(false).");

   CommandOptionNoArg dashsb(0, "NoSatBiases",
      " --NoSatBiases        Compute Receiver biases ONLY (not Rx+Sat biases) "
      "(false).");

   CommandOption dashmod(CommandOption::hasArgument, CommandOption::stdType,
      0,"Model"," --Model <type>       Ionospheric model: type is linear, "
      "quadratic or cubic");

   CommandOption dashMinPoints(CommandOption::hasArgument, CommandOption::stdType,
      0,"MinPoints",
     " --MinPoints <n>      Minimum points per satellite required");
   dashMinPoints.setMaxCount(1);

   CommandOption dashMinTimeSpan(CommandOption::hasArgument, CommandOption::stdType,
      0,"MinTimeSpan",
     " --MinTimeSpan <n>    Minimum timespan per satellite required (minutes)");
   dashMinTimeSpan.setMaxCount(1);

   CommandOption dashMinElevation(CommandOption::hasArgument, CommandOption::stdType,
      0,"MinElevation", " --MinElevation <n>   Minimum elevation angle (degrees)");
   dashMinElevation.setMaxCount(1);

   CommandOption dashMinLatitude(CommandOption::hasArgument, CommandOption::stdType,
      0,"MinLatitude", " --MinLatitude <n>    Minimum latitude (degrees)");
   dashMinLatitude.setMaxCount(1);

   CommandOption dashMaxLatitude(CommandOption::hasArgument, CommandOption::stdType,
      0,"MaxLatitude", " --MaxLatitude <n>    Maximum latitude (degrees)");
   dashMaxLatitude.setMaxCount(1);

   CommandOption dashMinLongitude(CommandOption::hasArgument, CommandOption::stdType,
      0,"MinLongitude", " --MinLongitude <n>   Minimum longitude (degrees)");
   dashMinLongitude.setMaxCount(1);

   CommandOption dashMaxLongitude(CommandOption::hasArgument, CommandOption::stdType,
      0,"MaxLongitude", " --MaxLongitude <n>   Maximum longitude (degrees)");
   dashMaxLongitude.setMaxCount(1);

   CommandOption dashTimeSector(CommandOption::hasArgument, CommandOption::stdType,
      0,"TimeSector", " --TimeSector <n>     Time sector (day | night | both)");
   dashTimeSector.setMaxCount(1);

   CommandOption dashTermOffset(CommandOption::hasArgument, CommandOption::stdType,
      0,"TerminOffset", " --TerminOffset <n>   Terminator offset (minutes)");
   dashTermOffset.setMaxCount(1);

   CommandOption dashIonoHt(CommandOption::hasArgument, CommandOption::stdType,
      0,"IonoHeight", " --IonoHeight <n>     Ionosphere height (km)");
   dashIonoHt.setMaxCount(1);

   CommandOption dashXsat(CommandOption::hasArgument, CommandOption::stdType,
      '0', "XSat", " Other options:\n --XSat <sat>         Exclude this satellite "
      "(<sat> may be <system> only)");
   
   // ... other options
   CommandOptionNoArg dashv('v', "verbose",
         " Help:\n [-v|--verbose]       print extended output info.");

   CommandOptionNoArg dashd('d', "debug",
         " [-d|--debug]         print extended output info.");

   CommandOptionNoArg dashh('h', "help",
         " [-h|--help]          print syntax and quit.");

   // ... rest of the command line
   CommandOptionRest Rest("");

   CommandOptionParser Par(