ionobias.cpp

根据GPS观测数据

         }

         LLI = Position(LA,LO,IonoHt*1000.0);     // 3rd entry is actually not used.
         TimeLimits(LLI, robs.time.DOY(), TimeSector, begintime, endtime);
         if(endtime >= begintime) {
            if(hours < begintime || hours > endtime) continue;
         }
         else {
            if(hours < begintime && hours > endtime) continue;
         }

            // compute the obliquity
         ob = obliquity(EL);

            // write out
         fout <<        setw(4)                    << robs.time.GPSfullweek();
         fout << " " << setw(8) << setprecision(1) << robs.time.GPSsow();
         fout << " " << setw(9) << setprecision(5) << LA; // latitude
         fout << " " << setw(10) << setprecision(5) << LO+cr; // co-rotating longitude
         fout << " " << setw(4) << setprecision(2) << ob; // 1/obliquity
         fout << " " << setw(8) << setprecision(3) << SR; // slant TEC
         fout << " " << setw(6) << setprecision(2) << 1;  // sigma ?? TD
         fout << " " << setw(2) << sat.id; // PRN
         fout << " " << setw(3) << nfile+1; // file number
         fout << endl;

         EstimationFlag[nfile][sat.id] = true;
         NgoodPoints++;
         npts[sat.id]++; // Npts for this sat
         if(npts[sat.id]==1) begin[sat.id] = robs.time;
         end[sat.id] = robs.time;

      }  // end for loop over sats


   } while(1);

      // revised header
   WriteStationHeader(NgoodPoints,StationName,StationPosition);

      // revise estimation flags
   if(verbose) oflog << "PRN  Points  Timespan   Begin       End  (hrs)\n";
   for(i=1; i<=MAXPRN; i++) {
      if(npts[i] > 0) {
         if(verbose) oflog << "G" << setfill('0') << setw(2) << i << setfill(' ')
            << setw(6) << npts[i]
            << setw(10) << setprecision(2) << (end[i]-begin[i])/3600.0
            << setw(10) << setprecision(2) << begin[i].secOfDay()/3600.0
            << setw(10) << setprecision(2) << end[i].secOfDay()/3600.0;
         if((end[i]-begin[i] < MinTimeSpan*60.0 || npts[i] < MinPoints)){
            if(verbose) {
               oflog << " reject(";
               if(end[i]-begin[i] < MinTimeSpan*60.0) oflog << " time ";
               if(npts[i] < MinPoints) oflog << " pts ";
               oflog << ")";
            }
            EstimationFlag[nfile][i] = false;
            NgoodPoints -= npts[i];
         }
         if(verbose) oflog << endl;
      }
   }

   if(NgoodPoints > 0) {
      NgoodStations++;
      ndata += NgoodPoints;
   }

   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); }
}

//------------------------------------------------------------------------------------
void WriteATHeader(void) throw(Exception)
{
try {
   int i,j;
   fout.seekp(0);
   fout << setw(5) << Filenames.size() << " " << setw(5) << NgoodStations
      << " Number (max, good) stations in this file \n";
   for(i=0; i<Filenames.size(); i++) {
      for(j=0; j<MAXPRN+1; j++) fout << (EstimationFlag[i][j] ? '1' : '0');
      fout << "\n";
   }
   fout << fixed;
}
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 WriteStationHeader(int npts, string sta_name, Position llr) throw(Exception)
{
try {
   fout.seekp(current_header_pos);
   fout << "Npt " << setw(5) << npts;
   //fout << in_file_s;
   //fout << sta_id;
   fout << " Sta " << sta_name;
   fout << " LLH " << setw(10) << setprecision(4) << llr[0]; //gllh.getLatitude();
   fout << " " << setw(10) << setprecision(4) << llr[1]; //gllh.getLongitude();
   fout << " " << setw(10) << setprecision(4)
      << llr[2]-Position::radiusEarth(llr[0],WGS84.a(),WGS84.eccSquared());
                //gllh.getAltitude();
   fout << endl;
}
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 ParseLine(string& str, vector<string>& wds) throw(Exception)
{
try {
   istringstream iss(str);
   string wd;
   wds.clear();
   while(iss >> wd) {
      wds.push_back(wd);
   }
}
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 ReadATandCompute(void) throw(Exception)
{
try {
   ifstream ifs;
   ifs.open(ATFileName.c_str());       // output file is now the input
   if(!ifs) {
      cerr << "Failed to open AT file " << ATFileName << " for input" << endl;
      return -1;
   }
   else if(verbose) oflog << "\nOpened AT file " << ATFileName << " for input\n";

      // read the AT header
   int i,j,k,ii,jj;
   long N,n;
   string line;

   EstimationFlag.clear();
   ifs >> N >> n;
   getline(ifs,line);      // read to eol
   for(i=0; i<N; i++) {
      getline(ifs,line);
      for(j=0; j<line.size(); j++) BoolVec[j] = (line[j] == '1');
      EstimationFlag.push_back(BoolVec);
   }
   
   //if(verbose) {
   //   oflog << "Read AT (" <<N<< "," <<n<< "," << EstimationFlag.size() << ")\n";
   //   for(i=0; i<EstimationFlag.size(); i++) {
   //      for(j=0; j<MAXPRN+1; j++) oflog << (EstimationFlag[i][j] ? '1' : '0');
   //      oflog << "\n";
   //   }
   //}

   if(N != EstimationFlag.size()) { //oops
   }

      // dimension and initialize the LS problem
   if(Model == "cubic") {
      oflog << "Model is cubic\n";
      NIonoParam = 10;
   }
   else if(Model == "quadratic") {
      oflog << "Model is quadratic\n";
      NIonoParam = 6;
   }
   else {   // linear and default
      oflog << "Model is linear\n";
      NIonoParam = 3;
   }

   for(NBiasParam=0,i=0; i<N; i++) {
      if(!ComputeSatBiases) NBiasParam++;
      else for(j=0; j<MAXPRN+1; j++) if(EstimationFlag[i][j]) NBiasParam++;
   }
   NTotalParam = NIonoParam + NBiasParam;

   Sol.resize(NTotalParam,0.0);
   Cov.resize(NTotalParam,NTotalParam,0.0);
   InfData.resize(NTotalParam,0.0);

      // Read the rest of the AT file
   int wn,prn,nfile,in;
   double sow,lat,lon,obq,sr,sig,d2=0;
   string stationID;
   vector<string> words;
   pair<string,int> Comp;

      // loop over stations
   oflog << setw(2) << N << "  Number of stations (N data and filename follow).\n";
   for(ndata=0,i=0; i<N; i++) {
         // read station header
      getline(ifs,line);
      ParseLine(line,words);
      if(words[0] != string("Npt")) { //oops
      }
      n = asInt(words[1]);
      stationID = words[3];

      if(n > 0 && verbose) {
         oflog << setw(3) << i+1 << "  " << stationID << " " << setw(4) << n << " ";
         for(j=0; j<=MAXPRN; j++) oflog << (EstimationFlag[i][j] ? '1' : '0');
         oflog << " " << Filenames[i];
         oflog << endl;
         mapN[stationID] = n;
         mapFilename[stationID] = Filenames[i];
      }

         // read data
      for(j=0; j<n; j++) {
         getline(ifs,line);
         ParseLine(line,words);
         wn = asInt(words[0]);
         sow = asDouble(words[1]);
         lat = asDouble(words[2]);
         lon = asDouble(words[3]);
         obq = asDouble(words[4]);
         sr = asDouble(words[5]);
         sig = asDouble(words[6]);
         prn = asInt(words[7]);
         nfile = asInt(words[8]);

         // do not include rejected data
         if(!(EstimationFlag[i][prn])) continue;

         // if NOT computing satellite biases, lump all data together into "PRN 0"
         if(!ComputeSatBiases) { prn = 0; }

         // find min and max lat and lon
         if(ndata == 0) {
            MaxLat = MinLat = lat;
            MaxCRLon = MinCRLon = lon;
         }
         else {
            if(fabs(lat) > MaxLat) MaxLat=lat;
            if(fabs(lat) < MinLat) MinLat=lat;
            if(fabs(lon) > MaxCRLon) MaxCRLon=lon;
            if(fabs(lon) < MinCRLon) MinCRLon=lon;
         }
         ndata++;

         // add this data to the LS
         //d2 += sr*sr;
         // find the index in partials matrix for this station-satellite pair
         Comp = make_pair(stationID,prn);
         in = index(ComponentIDs,Comp);
         if(in == -1) {
            in = ComponentIDs.size();
            ComponentIDs.push_back(Comp);
         }

         //PartialsMatrix(Par,in,lat,lon,obq);
         // note that obq is 1/obliquity
         // row of partials matrix has [in] = 1 and if nb=NBiasParam
            PM[0] =       obq; // [nb+0]               (all models)
            PM[1] = lat * obq; // [nb+1]               (all models)
            PM[2] = lon * obq; // [nnb2]               (all models)
         if(NIonoParam > 3) {
            PM[3] = lat * lat * obq; // [nb+3]         (quadratic and cubic)
            PM[4] = lon * lon * obq; // [nb+4]         (quadratic and cubic)
            PM[5] = lat * lon * obq; // [nb+5]         (quadratic and cubic)
         }
         if(NIonoParam > 6) {
            PM[6] = lat * lat * lat * obq; // [nb+6]   (cubic only)
            PM[7] = lon * lon * lon * obq; // [nb+7]   (cubic only)
            PM[8] = lat * lat * lon * obq; // [nb+8]   (cubic only)
            PM[9] = lat * lon * lon * obq; // [nb+9]   (cubic only)
         }

         //LS.Add(Par,Dat,Wgt); do the sequential LS by hand for efficiency
         //
         // Inf += transpose(partials) * partials (weight = 1)
         // InfData += transpose(partials) * data
         //
         Cov(in,in) += 1.0;
         InfData(in) += sr;
         for(ii=0; ii<NIonoParam; ii++) {
            k = NBiasParam + ii;
            InfData(k) += sr * PM[ii];
            Cov(k,in) += PM[ii];
            Cov(in,k) += PM[ii];
            for(jj=0; jj<NIonoParam; jj++) {
               Cov(k,NBiasParam+jj) += PM[ii]*PM[jj];
            }
         }

      }  // end loop over points for this station

   }  // end loop over stations

   ifs.close();

   oflog << setw(9) << setprecision(2) << MinLat << "  Minimum Latitude\n";
   oflog << setw(9) << setprecision(2) << MaxLat << "  Maximum Latitude\n";
   oflog << setw(9) << setprecision(2) << MinCRLon << "  Minimum Co-rot lon\n";
   oflog << setw(9) << setprecision(2) << MaxCRLon << "  Maximum Co-rot lon\n";
   oflog << setw(5) << ndata << " data points used." << endl << endl;

   // solve the LS problem
   // Cov = inverse(information)
   // X = Cov * InfData
   try { Cov = inverse(Cov); }
   catch(Exception& e) {
      oflog << "Least squares failed because the problem is singular\n";
      return -2;
   }

   // Invert Cov via SVD - also expensive - maybe make option, see SVs and conditionN
   //SVD<double> svd;
   //svd(Cov);
   //oflog << "Singular Values range " << svd.S(0)
   //   << " to " << svd.S(NTotalParam-1) << endl;
   //for(i=1; i<NIonoParam; i++) {
   //   if(svd.S(i) < 1.e-14 * svd.S(0)) {
   //      oflog << "Edit SingularValue(" << i << ") = " << svd.S(i) << endl;
   //      svd.S(i) = 0;
   //   }
   //}
   //Vector<double> T(NIonoParam);
   //for(j=0; j<NIonoParam; j++) { // loop over columns
   //   T = 0.0;
   //   T(j) = 1.0;
   //   svd.backSub(T);
   //   for(i=0; i<NTotalParam; i++) Cov(i,j)=T(i);
   //}

   // compute solution
   Sol = Cov * InfData;
   //if(verbose) oflog << "Least squares solved successfully.\n";

   // print solution and sigma - remember lat and lon may be scaled by 1/1000
   bool biasout=false;
   if(BiasFileName.length() > 0) {
      fout.open(BiasFileName.c_str(),ios_base::out);
      if(!fout)
         cerr << "Failed to open output biases file " << BiasFileName << endl;
      else {
         biasout = true;
         fout << Title;
      }
   }
   oflog << setw(2) << NBiasParam << "  Number of SPR biases\n";
   if(biasout) fout << setw(2) << NBiasParam << "  Number of SPR biases\n";
   for(i=0; i<NBiasParam; i++) {
      ostringstream oss;
      oss << "BIAS " << setw(3) << i+1                                  // number
         << "  " << ComponentIDs[i].first                               // station id
         << " G" << setw(2) << setfill('0') << ComponentIDs[i].second   // sat G<prn>
         << setfill(' ') << fixed
         << " " << setw(4) << mapN[ComponentIDs[i].first]
         << " " << setw(12) << setprecision(6) << Sol(i)                // bias
         << scientific
         << " " << setw(10) << setprecision(3) << ::sqrt(Cov(i,i))        // sigma
         << " " << mapFilename[ComponentIDs[i].first]
         << endl;
      oflog << oss.str();
      if(biasout) fout << oss.str();
   }
   oflog << setw(2) << NTotalParam-NBiasParam << "  Number of ion model parameters\n";
   for(i=NBiasParam; i<NTotalParam; i++) {
      ostringstream oss;
      oss << setw(3) << i+1-NBiasParam << fixed                         // number
         << " " << setw(12) << setprecision(6) << Sol(i)                // solution
         << scientific
         << " " << setw(10) << setprecision(3) << ::sqrt(Cov(i,i))        // sigma
         << endl;
      oflog << oss.str();
      if(biasout) fout << oss.str();
   }

   // compute standard error estimates
   // TBD


   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); }
}

//------------------------------------------------------------------------------------
// elevation in degrees. Return the inverse of the obliquity factor.
double obliquity(double elevation) throw(Exception)
{
try {
   double ob;
   //const double coef[4]={1.02056,0.466332,3.50523,-1.84119};
   //double x2=(1-elevation/90.)*(1-elevation/90.);
   //ob = coef[3];
   //for(int i=2; i>=0; i--) ob = ob*x2 + coef[i];

   ob = WGS84.a()*cos(elevation*DEG_TO_RAD)/(WGS84.a()+IonoHt*1000);
   ob = ::sqrt(1.0-ob*ob);

   return ob;
}
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); }
}

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