tropmodel.cpp

GPS定位中对流层模型改正

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

//============================================================================
//
//This software developed by Applied Research Laboratories at the University of
//Texas at Austin, under contract to an agency or agencies within the U.S. 
//Department of Defense. The U.S. Government retains all rights to use,
//duplicate, distribute, disclose, or release this software. 
//
//Pursuant to DoD Directive 523024 
//
// DISTRIBUTION STATEMENT A: This software has been approved for public 
//                           release, distribution is unlimited.
//
//=============================================================================






/**
 * @file TropModel.cpp
 * Base class for tropospheric models, plus implementations of several
 * published models
 */

#include "TropModel.hpp"
#include "EphemerisRange.hpp"             // for Elevation()
#include "MathBase.hpp"                   // SQRT
#include "geometry.hpp"                   // DEG_TO_RAD
#include "GPSGeoid.hpp"                   // geoid.a() = R earth
#include "icd_200_constants.hpp"          // TWO_PI
#include "Geodetic.hpp"
#include "ECEF.hpp"

namespace gpstk
{
      // for temperature conversion from Celcius to Kelvin
   static const double CELSIUS_TO_KELVIN = 273.15;

      // Compute and return the full tropospheric delay. Typically call
      // setWeather(T,P,H) before making this call.
      // @param elevation Elevation of satellite as seen at receiver, in degrees
   double TropModel::correction(double elevation) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid)
         GPSTK_THROW(InvalidTropModel("Invalid model"));

      if(elevation < 0.0)
         return 0.0;

      return (dry_zenith_delay() * dry_mapping_function(elevation)
            + wet_zenith_delay() * wet_mapping_function(elevation));

   }  // end TropModel::correction(elevation)

      // Compute and return the full tropospheric delay, given the positions of
      // receiver and satellite and the time tag. This version is most useful
      // within positioning algorithms, where the receiver position and timetag may
      // vary; it computes the elevation (and other receiver location information)
      // and passes them to appropriate set...() routines and the
      // correction(elevation) routine.
      // @param RX  Receiver position in ECEF cartesian coordinates (meters)
      // @param SV  Satellite position in ECEF cartesian coordinates (meters)
      // @param tt  Time tag of the signal 
   double TropModel::correction(const Position& RX,
                                const Position& SV,
                                const DayTime& tt)
      throw(TropModel::InvalidTropModel)
   {
      if(!valid)
         GPSTK_THROW(InvalidTropModel("Invalid model"));

      double c;
      try
      {
         c = correction(RX.elevation(SV));
      }
      catch(InvalidTropModel& e)
      {
         GPSTK_RETHROW(e);
      }
      return c;
   }  // end TropModel::correction(RX,SV,TT)

      // Re-define the tropospheric model with explicit weather data.
      // Typically called just before correction().
      // @param T temperature in degrees Celsius
      // @param P atmospheric pressure in millibars
      // @param H relative humidity in percent
   void TropModel::setWeather(const double& T,
                              const double& P,
                              const double& H)
      throw(InvalidParameter)
   {
      temp = T + CELSIUS_TO_KELVIN;
      press = P;
      humid = H;
      if (temp < 0.0)
      {
         valid = false;
         GPSTK_THROW(InvalidParameter("Invalid temperature parameter."));
      }
      if (press < 0.0)
      {
         valid = false;
         GPSTK_THROW(InvalidParameter("Invalid pressure parameter."));
      }
      if (humid < 0.0 || humid > 100.0)
      {
         valid = false;
         GPSTK_THROW(InvalidParameter("Invalid humidity parameter."));
      }         
   }  // end TropModel::setWeather(T,P,H)
   
      // Re-define the tropospheric model with explicit weather data.
      // Typically called just before correction().
      // @param wx the weather to use for this correction       
   void TropModel::setWeather(const WxObservation& wx)
      throw(InvalidParameter)
   {
      if (wx.isAllValid())
      {
         try
         {
            setWeather(wx.temperature, wx.pressure, wx.humidity);
            valid = true;
         }
         catch(InvalidParameter& e)
         {
            valid = false;
            GPSTK_RETHROW(e);
         }
      }
      else
      {
         valid = false;
         GPSTK_THROW(InvalidParameter("Invalid weather data"));
      }
   }
         
   // -------------------------------------------------------------------------------
   // Simple Black model. This has been used as the 'default' for many years.

      // Default constructor
   SimpleTropModel::SimpleTropModel(void)
   {
      setWeather(20.0, 980.0, 50.0);
      Cwetdelay = 0.122382715318184;
      Cdrydelay = 2.235486646978727;
      Cwetmap = 1.000282213715744;
      Cdrymap = 1.001012704615527;
      valid = true;
   }

      // Creates a trop model from a weather observation
      // @param wx the weather to use for this correction.
   SimpleTropModel::SimpleTropModel(const WxObservation& wx)
      throw(InvalidParameter)
   {
      setWeather(wx);
      valid = true;
   }

      // Create a tropospheric model from explicit weather data
      // @param T temperature in degrees Celsius
      // @param P atmospheric pressure in millibars
      // @param H relative humidity in percent
   SimpleTropModel::SimpleTropModel(const double& T,
                                    const double& P,
                                    const double& H)
      throw(InvalidParameter)
   {
      setWeather(T,P,H);
      valid = true;
   }

      // Re-define the tropospheric model with explicit weather data.
      // Typically called just before correction().
      // @param T temperature in degrees Celsius
      // @param P atmospheric pressure in millibars
      // @param H relative humidity in percent
   void SimpleTropModel::setWeather(const double& T,
                                    const double& P,
                                    const double& H)
      throw(InvalidParameter)
   {
      TropModel::setWeather(T,P,H);
      GPSGeoid geoid;
      Cdrydelay = 2.343*(press/1013.25)*(temp-3.96)/temp;
      double tks = temp * temp;
      Cwetdelay = 8.952/tks*humid*std::exp(-37.2465+0.213166*temp-(0.256908e-3)*tks);
      Cdrymap =1.0+(0.15)*148.98*(temp-3.96)/geoid.a();
      Cwetmap =1.0+(0.15)*12000.0/geoid.a();
      valid = true;
   }  // end SimpleTropModel::setWeather(T,P,H)

      // Re-define the tropospheric model with explicit weather data.
      // Typically called just before correction().
      // @param wx the weather to use for this correction       
   void SimpleTropModel::setWeather(const WxObservation& wx)
      throw(InvalidParameter)
   {
      TropModel::setWeather(wx);
   }

      // Compute and return the zenith delay for dry component of the troposphere
   double SimpleTropModel::dry_zenith_delay(void) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid)
         GPSTK_THROW(InvalidTropModel("Invalid model"));

      return Cdrydelay;
   }

      // Compute and return the zenith delay for wet component of the troposphere
   double SimpleTropModel::wet_zenith_delay(void) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid)
         GPSTK_THROW(InvalidTropModel("Invalid model"));

      return Cwetdelay;
   }

      // Compute and return the mapping function for dry component
      // of the troposphere
      // @param elevation is the Elevation of satellite as seen at receiver,
      //                  in degrees
   double SimpleTropModel::dry_mapping_function(double elevation) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid)
         GPSTK_THROW(InvalidTropModel("Invalid model"));

      if(elevation < 0.0) return 0.0;

      double d = std::cos(elevation*DEG_TO_RAD);
      d /= Cdrymap;
      return (1.0/SQRT(1.0-d*d));
   }

      // Compute and return the mapping function for wet component
      // of the troposphere
      // @param elevation is the Elevation of satellite as seen at receiver,
      //                  in degrees
   double SimpleTropModel::wet_mapping_function(double elevation) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid)
         GPSTK_THROW(InvalidTropModel("Invalid model"));

      if(elevation < 0.0) return 0.0;

      double d = std::cos(elevation*DEG_TO_RAD);
      d /= Cwetmap;
      return (1.0/SQRT(1.0-d*d));
   }

   // -------------------------------------------------------------------------------
   // Tropospheric model based on Goad and Goodman(1974),
   // "A Modified Hopfield Tropospheric Refraction Correction Model," Paper
   // presented at the Fall Annual Meeting of the American Geophysical Union,
   // San Francisco, December 1974.
   // See Leick, "GPS Satellite Surveying," Wiley, NY, 1990, Chapter 9,
   // particularly Table 9.1.
   // -------------------------------------------------------------------------------
 
   static const double GGdryscale = 8594.777388436570600;
   static const double GGwetscale = 2540.042008403690900;

      // Default constructor
   GGTropModel::GGTropModel(void)
   {
      TropModel::setWeather(20.0, 980.0, 50.0);
      Cdrydelay = 2.59629761092150147e-4;    // zenith delay, dry
      Cwetdelay = 4.9982784999977412e-5;     // zenith delay, wet
      Cdrymap = 42973.886942182834900;       // height for mapping, dry
      Cwetmap = 12700.210042018454260;       // height for mapping, wet
      valid = true;
   }  // end GGTropModel::GGTropModel()

      // Creates a trop model from a weather observation
      // @param wx the weather to use for this correction.
   GGTropModel::GGTropModel(const WxObservation& wx)
      throw(InvalidParameter)
   {
      setWeather(wx);
      valid = true;
   }  // end GGTropModel::GGTropModel(weather)

      // Create a tropospheric model from explicit weather data
      // @param T temperature in degrees Celsius
      // @param P atmospheric pressure in millibars
      // @param H relative humidity in percent
   GGTropModel::GGTropModel(const double& T,
                            const double& P,
                            const double& H)
      throw(InvalidParameter)
   {
      setWeather(T,P,H);
      valid = true;
   } // end GGTropModel::GGTropModel()

   double GGTropModel::dry_zenith_delay(void) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid)
         GPSTK_THROW(InvalidTropModel("Invalid model"));

      return (Cdrydelay * GGdryscale);
   }  // end GGTropModel::dry_zenith_delay()

   double GGTropModel::wet_zenith_delay(void) const
      throw(TropModel::InvalidTropModel)