vtecmap.cpp

根据GPS观测数据

#pragma ident "$Id: VTECMap.cpp 107 2006-09-08 20:11:43Z 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 VTECMap.cpp
 * Program VTECMap will TD...
 *   a simple ionospheric model using least squares and slant TEC values
 *   from multiple stations.
 */

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

#include "VTECMap.hpp"
#include "RinexUtilities.hpp"

using namespace gpstk;
using namespace std;

//------------------------------------------------------------------------------------
const double VTECMap::VTECErrorMultipath = 4.0;
const double VTECMap::VTECErrorSat = 0.9;
const double VTECMap::VTECErrorCFC[4] = {-0.000045, 0.0096, -0.6755, 15.84};
const double VTECMap::ObliqCoef[4] = {-4.316e-06, 0.001043, -0.08771, 3.57};

//------------------------------------------------------------------------------------
void VTECMap::CopyInputData(VTECMap &right)
{
   Decorrelation = right.Decorrelation;
   MinElevation = right.MinElevation;
   IonoHeight = right.IonoHeight;
   gridtype = right.gridtype;
   fittype = right.fittype;
   BeginLat = right.BeginLat;
   DeltaLat = right.DeltaLat;
   NumLat = right.NumLat;
   BeginLon = right.BeginLon;
   DeltaLon = right.DeltaLon;
   NumLon = right.NumLon;
   RefStation = right.RefStation;
}

//------------------------------------------------------------------------------------
void VTECMap::SetDefaults()
{
   Decorrelation = 3.0;    // TECU/1000km
   MinElevation = 10.0;
   gridtype = UniformLatLon;
   fittype = Constant;
   BeginLat = 21.0;
   BeginLon = 230.0;
   DeltaLat = 0.25;
   DeltaLon = 1.0;
   NumLat = NumLon = 40;
   IonoHeight = 350. * 1000.0;       // 350km in meters
}

//------------------------------------------------------------------------------------
void VTECMap::reallyMakeGrid(Station& refS, int factor)
   throw(Exception)
{
   RefStation = refS;

   grid = new GridData[NumLat*NumLon];
   if(!grid) throw Exception("VTECMap::reallyMakeGrid failed to allocate grid");

try {
   int i,j,ishift,jshift,k;

   if(gridtype == UniformSpace) {
      int ii;
      Position ptx1,ptx2,pty1,pty2,pt3,DeltaX,DeltaY;

         // azimuth = 0 : x1
      ptx1 = refS.llr.getIonosphericPiercePoint(MinElevation,  0.0,IonoHeight);
      ptx1.transformTo(Position::Cartesian);

         // azimuth = 180 : x2
      ptx2 = refS.llr.getIonosphericPiercePoint(MinElevation,180.0,IonoHeight);
      ptx2.transformTo(Position::Cartesian);

         // azimuth = 90 : y1
      pty1 = refS.llr.getIonosphericPiercePoint(MinElevation, 90.0,IonoHeight);
      pty1.transformTo(Position::Cartesian);

         // azimuth = 270 : y2
      pty2 = refS.llr.getIonosphericPiercePoint(MinElevation,270.0,IonoHeight);
      pty2.transformTo(Position::Cartesian);

         // find the center of the sheet
      pt3 = (ptx1 + ptx2)*0.5;

         // get orthogonal vectors in the plane, and compute step size info
      DeltaX = (ptx1 - ptx2)*(1.0/double(NumLon-1));
      DeltaY = (pty1 - pty2)*(1.0/double(NumLat-1));

         // create the grid
      for(i=0; i<NumLon; i++) {           // i == x == lon
         ishift = i - (NumLon/2);
         for(j=0; j<NumLat; j++) {        // j == y == lat
            k = i * NumLat + j;           // k is the index in grid array
            jshift = j - (NumLat/2);
            grid[k].XYZ = pt3 + (DeltaX*ishift + DeltaY*jshift)*factor;
            grid[k].LLR = grid[k].XYZ;
            grid[k].LLR.transformTo(Position::Geocentric);
         }  // end j loop over lon
      }  // end i loop over lat
   }

   if(gridtype == UniformLatLon) {
      double LatCenter = BeginLat + NumLat * DeltaLat/2.0;
      double LonCenter = BeginLon + NumLon * DeltaLon/2.0;
      double rad;
      {
            // this is a trick to get the radius of the ionosphere
         Position IPP = refS.llr.getIonosphericPiercePoint(90,0,IonoHeight);
         rad = IPP.radius();
      }
         // create the grid
      for(i=0; i<NumLon; i++) {           // i == x == lon
         ishift = i - (NumLon/2);
         for(j=0; j<NumLat; j++) {        // j == y == lat
            jshift = j - (NumLat/2);
            k = i * NumLat + j;           // k is the index in grid array
            grid[k].LLR.setGeocentric(
               LatCenter + factor * jshift * DeltaLat, // lat (deg)
               LonCenter + factor * ishift * DeltaLon, // lon (deg)
               rad);                                   // radius (m)
            grid[k].XYZ = grid[k].LLR;
            grid[k].XYZ.transformTo(Position::Cartesian);
         }  // end j loop over lon
      }  // end i loop over lat
   }
}
catch(gpstk::Exception& e) {
      cerr << "VTECMap:reallyMakeGrid caught an exception\n" << e;
      GPSTK_RETHROW(e);
}
catch (...) {
      cerr << "VTECMap:reallyMakeGrid caught an unknown exception\n";
}
}

//------------------------------------------------------------------------------------
void VTECMap::OutputGrid(ostream& os)
{
   int i,j,k;
   os << fixed << setprecision(3);
   for(j=0; j<NumLat; j++) {
      for(i=0; i<NumLon; i++) {
         k = i * NumLat + j;
         os << grid[k].LLR.printf(" %.8a %.8l %.3r");
         os << grid[k].XYZ.printf(" %.3x %.3y %.3z");
         os << " " << i << " " << j << endl;
      }
   }
}

//------------------------------------------------------------------------------------
void VTECMap::ComputeMap(DayTime& epoch, vector<ObsData>& data)
{
   int i,j,k,n;
      // first compute the average value
   n = 0;
   ave = 0.0;
   for(k=0; k<data.size(); k++) {
      n++;
      ave *= double(n-1)/double(n);
      ave += data[k].VTEC/double(n);
   }
      // now compute the value at each grid point
   for(i=0; i<NumLon; i++) {
      for(j=0; j<NumLat; j++) {
         k = i * NumLat + j;
         ComputeGridValue(grid[k],data);
      }
   }
}

//------------------------------------------------------------------------------------
// Compute the grid values. Called by ComputeMap.
void VTECMap::ComputeGridValue(GridData& gridpt, vector<ObsData>& data)
{
   double gridLat = gridpt.LLR.getGeocentricLatitude() * DEG_TO_RAD;
   double gridLon = gridpt.LLR.longitude();
   if(gridLon > 180.0) gridLon -= 360.0;
   gridLon *= DEG_TO_RAD;

   double destLat, destLon,dLat,dLon;
   double sg,cg,sd,dist,range,bear,d;
   vector<double> vtec,xtmp,ytmp,sigma;

      // loop over all data
   for(int k=0; k<data.size(); k++) {
      //if(data[k].elevation < MinElevation) continue;   // here?
      destLat = data[k].latitude * DEG_TO_RAD;
      destLon = data[k].longitude * DEG_TO_RAD;
      // compute distance in the plane from grid to dest(data)
      dLon = destLon - gridLon;
      sg = sin(gridLat);
      cg = cos(gridLat);
      sd = sin(destLat);
      d = sg*sd + cg*cos(destLat)*cos(dLon);
      dist = acos(d);
      // TD what is range? where does the 1.852 come from?
      // TD what are the units of range? assume km, then Decorrelation = TECU/1000km
      // decor error = range * Decorrelation
      range = 1.852 * 60 * dist * RAD_TO_DEG;   // 1.852 * distance in min of arc
      if(ABS(dist) < 0.01) dist = 0.01;
      d = (sd - sg*cos(dist)) / sin(dist)*cg;   // d = cos(bearing)
      if(ABS(d) > 1) {
         if(d > 0) d = 1.0;
         else d = -1.0;
      }
      bear = acos(d);
      if(dLon > 0) bear = TWO_PI - bear;