rescor.cpp

GPS数据预处理软件

#pragma ident "$Id: ResCor.cpp 700 2007-07-24 14:20:15Z 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 ResCor.cpp
 * 'Residuals and Corrections'
 * Open and read a single Rinex observation file, apply editing commands
 * using the RinexEditor package, compute any of several residuals and corrections
 * and register extended Rinex observation types for them, and then write
 * the edited data, along with the new extended observation types,
 * to an output Rinex observation file. Input is all on the command line.
 * ResCor is implemented by deriving a special class from class RinexEditor and
 * using its virtual functions to implement all the changes necessary to define
 * and compute the residuals and corrections.
 */

//------------------------------------------------------------------------------------
// ToDo
// catch exceptions -- elsewhere and on reading header and obs
// allow user to specify trop model, both for RAIM and for TR output
//

#include "StringUtils.hpp"
#include "DayTime.hpp"
#include "RinexSatID.hpp"
#include "CommandOptionParser.hpp"
#include "CommandOption.hpp"
#include "CommandOptionWithTimeArg.hpp"
#include "RinexObsBase.hpp"
#include "RinexObsData.hpp"
#include "RinexObsHeader.hpp"
#include "RinexObsStream.hpp"
#include "RinexNavStream.hpp"
#include "RinexNavData.hpp"
#include "SP3Stream.hpp"
#include "SP3EphemerisStore.hpp"
#include "BCEphemerisStore.hpp"
#include "EphemerisRange.hpp"
#include "TropModel.hpp"
#include "PRSolution.hpp"
#include "WGS84Geoid.hpp"           // for obliquity
#include "Stats.hpp"
#include "geometry.hpp"             // DEG_TO_RAD
#include "icd_200_constants.hpp"    // PI,C_GPS_M,OSC_FREQ,L1_MULT,L2_MULT

#include "RinexEditor.hpp"
#include "RinexUtilities.hpp"
#include "Position.hpp"

#include <time.h>
#include <string>
#include <vector>
#include <map>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <sstream>

//------------------------------------------------------------------------------------
using namespace std;
using namespace gpstk;
using namespace StringUtils;

//------------------------------------------------------------------------------------
   // prgm data
string PrgmName("ResCor");
string PrgmVers("3.7 1/22/07");

// data used in program
const double CFF=C_GPS_M/OSC_FREQ;
const double F1=L1_MULT;   // 154.0;
const double F2=L2_MULT;   // 120.0;
const double f12=F1*F1;
const double f22=F2*F2;
const double wl1=CFF/F1;
const double wl2=CFF/F2;
const double wl1r=F1/(F1+F2);
const double wl2r=F2/(F1+F2);
const double wl1p=wl1*F1/(F1-F2);
const double wl2p=-wl2*F2/(F1-F2);
const double if1r=f12/(f12-f22);
const double if2r=-f22/(f12-f22);
const double if1p=wl1*f12/(f12-f22);
const double if2p=-wl2*f22/(f12-f22);
const double gf1r=-1;
const double gf2r=1;
const double gf1p=wl1;
const double gf2p=-wl2;
const double alpha=f12/f22 - 1.0;
const double FL1=F1*10.23e6;                 // Hz
const double TECUperM=FL1*FL1*1.e-16/40.28;  // 6.1617 TECU/m (0.16229 m/TECU)

clock_t totaltime;
string Title;
   // input flags and data
bool Debug,Verbose,Callow,Cforce;
double IonoHt;
RinexSatID SVonly;
string LogFile;
ofstream logof,oferr;         // don't call it oflog - RinexEditor has that
   // Rinex headers, input and output, saved
RinexObsHeader rhead, rheadout;
   // ephemeris
string NavDir;
vector<string> NavFiles;
SP3EphemerisStore SP3EphList;
BCEphemerisStore BCEphList;
SimpleTropModel ggtm;
   // for use with current position and in RefPosMap (RAIM and/or RefPosFile)
typedef struct ReferencePositionFileData {
   Position RxPos;              // XYZT
   bool valid;
   int NPRN;
   double clk,PDOP,GDOP,RMS;
} RefPosData;
RefPosData CurrRef;        // current reference position
   // reference and RAIM solution
string RefPosFile,KnownPos;
bool doRAIM,editRAIM,outRef,headRAIM,HaveRAIM;
bool RefPosInput,KnownPosInput,KnownLLH,RefPosFlat;
double minElev;
vector<SatID> Sats;
vector<double> PRange;
//RAIMSolution RAIMSol;
PRSolution prsol;
Stats<double> ARSX,ARSY,ARSZ;      // average solution, for header output
   // computation
int inC1,inP1,inP2,inL1,inL2;      // indexes in rhead of C1, C1/P1, P2, L1 and L2
int inEP,inPS;                     // flags for input of ephemeris, Rx position
int inD1,inD2,inS1,inS2;
DayTime CurrentTime, PrgmEpoch;
// these 3 vectors parallel
vector<string> OTstrings;          // list of OTs (strings) to be computed
vector<RinexObsHeader::RinexObsType> OTList;
vector<int> OTindex;
int otC1,otP1,otP2,otL1,otL2;      // indexes in rheadout of C1, C1/P1, P2, L1 and L2
int otD1,otD2,otS1,otS2;
bool DoSVX;
WGS84Geoid WGS84;
// compute non-dispersive range, ionospheric delay, multipath (L1 and L2)
bool DoXR;
double XRM0[4],XRM1[4],XRM2[4],XRM3[4];
double *XRM[4]={XRM0,XRM1,XRM2,XRM3};
double XRdat[4],XRsol[4];
   // structure for holding raw range and phase data during computation
typedef struct range_and_phase_data {
   double L1,L2,P1,P2;
   int LL1,LL2;
} RCData;
   // map of <sat,RCData>
RCData DataStore;
map<RinexSatID,RCData> DataStoreMap;
   // debiasing output data
map<RinexObsHeader::RinexObsType,map<RinexSatID,double> > AllBiases; // (OT,SV)
   // reference position as function of time (from input)
map<DayTime,RefPosData> RefPosMap;
double RefPosMapDT;

string RxhelpString=
"\n --RxFlat <fn> : fn is a file with reference receiver positions and times:\n"
"  The first line in the file (other than comments, marked by # in column 1)\n"
"  is the format for each line of the file, using the specifications in\n"
"  DayTime::setToString() and Position::setToString().\n"
"  The second line is a pattern made up of characters T, P and X indicating the\n"
"  content of both the lines in the file and the format: (white-space-delimited)\n"
"  words on each line are either part of the time(T) or position(P) specification,\n"
"  or are to be ignored(X). For example, the file begins with these six lines:\n"
"  # format:\n"
"  t= %F %g p= %x %y %z\n"
"  # pattern:\n"
"  XTTXPPP\n"
"  # data:\n"
"  t= 1281 259200    p=   -2701232.4        6123085.7        1419837.5";

//------------------------------------------------------------------------------------
// inherit RinexEditor so that callback routines can be defined by Prgm ResCor
class RCRinexEditor : public RinexEditor
{
   public:
         /// Constructor.
      RCRinexEditor() throw() {};

         /// destructor
      virtual ~RCRinexEditor() {}
   
         /// after reading input header and before calling
         /// RinexEditor::EditHeader (pass input header)
      virtual int BeforeEditHeader(const RinexObsHeader& rhin) throw(Exception);

         /// after calling RinexEditor::EditHeader (pass output header)
      virtual int AfterEditHeader(const RinexObsHeader& rhout) throw(Exception);

         /// after reading input obs and before calling
         /// RinexEditor::EditObs (pass input obs)
      virtual int BeforeEditObs(const RinexObsData& roin) throw(Exception);

         /// before writing out header (pass output header)
      virtual int BeforeWritingHeader(RinexObsHeader& rhout) throw(Exception);

         /// before writing out filled header
      virtual int BeforeWritingFilledHeader(RinexObsHeader& rhout) throw(Exception);

         /// just before writing output obs (pass output obs)
      virtual int BeforeWritingObs(RinexObsData& roout) throw(Exception);

}; // end class RCRinexEditor

//------------------------------------------------------------------------------------
// prototypes
int GetCommandLine(int argc, char **argv, RCRinexEditor& rc) throw(Exception);
int PrepareInput(void) throw(Exception);
int LoopOverObs(void) throw(Exception);
void SaveData(const RinexObsData& rod, const RinexObsHeader& rh,
   int xL1, int xL2, int xP1, int xP2) throw(Exception);
int UpdateRxPosition(void) throw(Exception);
void ComputeNewOTs(RinexObsData& rod) throw(Exception);
void CloseOutputFile(void) throw(Exception);
void PreProcessArgs(const char *arg, vector<string>& Args) throw(Exception);
int setBiasLimit(RinexObsHeader::RinexObsType& ot, double lim) throw(Exception);
double removeBias(const RinexObsHeader::RinexObsType& ot, const RinexSatID& sat,
   bool& reset, DayTime& tt, double delta) throw(Exception);

//------------------------------------------------------------------------------------
int main(int argc, char **argv)
{
try {
   totaltime = clock();
   int iret;
   // NB. Do not instantiate editor outside main(), b/c DayTime::END_OF_TIME is a
   // static const that can produce static intialization order problems under some OS.
   RCRinexEditor REC;
   CurrentTime = DayTime::BEGINNING_OF_TIME; // for same reason, init here...

      // Title and description
   Title = PrgmName + ", part of the GPS ToolKit, Ver. " + PrgmVers
      + " (editor " + REC.getRinexEditVersion() + string("), Run ");
   time_t timer;
   struct tm *tblock;
   timer = time(NULL);
   tblock = localtime(&timer);
   PrgmEpoch.setYMDHMS(1900+tblock->tm_year,1+tblock->tm_mon,
               tblock->tm_mday,tblock->tm_hour,tblock->tm_min,tblock->tm_sec);
   Title += PrgmEpoch.printf("%04Y/%02m/%02d %02H:%02M:%02S");
   Title += "\n";
   cout << Title;

      // define extended types
   iret = RegisterARLUTExtendedTypes();
   if(iret) goto quit;

   // Set defaults, define command line and parse it.
   // Send REdit cmds to REC. Check validity of input.
   iret = GetCommandLine(argc, argv, REC);
   if(iret) goto quit;

   // Initialize, read ephemerides, set flags and prepare for processing
   iret = PrepareInput();
   if(iret) goto quit;

   // Edit the file, including callbacks
   iret = REC.EditFile();
   if(Debug) logof << "EditFile returned " << iret << endl;
   if(iret) goto quit;

   quit:
   // compute run time
   totaltime = clock()-totaltime;
   logof << "ResCor timing: " << setprecision(3)
      << double(totaltime)/double(CLOCKS_PER_SEC) << " seconds.\n";

   logof.close();
   cout << "End ResCor" << endl;
   return iret;
}
catch(gpstk::FFStreamError& e) { cerr << e; }
catch(gpstk::Exception& e) { cerr << e; }
catch(exception& e) { cerr << e.what(); }
catch (...) { cerr << "Unknown error.  Abort." << endl; }
   return 1;
}   // end main()

//------------------------------------------------------------------------------------
// Set defaults, define command line and parse it. Send REdit cmds to REC.
// Check validity of input
int GetCommandLine(int argc, char **argv, RCRinexEditor& REC) throw(Exception)
{
   bool help=false,Rxhelp=false,REChelp=false,ROThelp=false;
   int i,j,iret;
try {
      // defaults
   Debug = Verbose = false;

   doRAIM = false;
   KnownPosInput = RefPosInput = false;
   outRef = true;
   editRAIM = true;
   headRAIM = false;
   minElev = 0.0;
   
   IonoHt = 400.0;      // km

   Callow = true;
   Cforce = false;

   LogFile = string("ResCor.log");

      // -------------------------------------------------
      // required options

      // optional options

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

   // ephemeris
   CommandOption dashn(CommandOption::hasArgument, CommandOption::stdType,