tropmodel.cpp

GPS定位中对流层模型改正

                                  const Position& SV,
                                  const DayTime& tt)
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather)
            
         if(!validRxLatitude)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
         if(!validRxHeight)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
         if(!validDOY)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
      }

         // compute height and latitude from RX
      setReceiverHeight(RX.getHeight());
      setReceiverLatitude(RX.getGeodeticLatitude());

         // compute day of year from tt
      setDayOfYear(int(tt.DOYday()));

      return TropModel::correction(RX.elevation(SV));

   }  // end NBTropModel::correction(RX,SV,TT)

      // 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 
      // This function is deprecated; use the Position version
   double NBTropModel::correction(const Xvt& RX,
                                  const Xvt& SV,
                                  const DayTime& tt)
      throw(TropModel::InvalidTropModel)
   {
      Position R(RX),S(SV);
      return NBTropModel::correction(R,S,tt);
   }

      // Compute and return the zenith delay for dry component of the troposphere
   double NBTropModel::dry_zenith_delay(void) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
         if(!validRxLatitude)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
         if(!validRxHeight)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
         if(!validDOY)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
      }
      double beta = NB_Interpolate(latitude,doy,ZB);
      double gm = 9.784*(1.0-2.66e-3*std::cos(2.0*latitude*DEG_TO_RAD)-2.8e-7*height);

         // scale factors for height above mean sea level
         // if weather is given, assume it's measured at ht -> kw=kd=1
      double kd=1, base=std::log(1.0-beta*height/temp);
      if(interpolateWeather)
         kd = std::exp(base*NBg/(NBRd*beta));

         // compute the zenith delay for dry component
      return ((1.0e-6*NBk1*NBRd/gm) * kd * press);

   }  // end NBTropModel::dry_zenith_delay()

      // Compute and return the zenith delay for wet component of the troposphere
   double NBTropModel::wet_zenith_delay(void) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
         if(!validRxLatitude)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
         if(!validRxHeight)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
         if(!validDOY)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
      }
      double beta = NB_Interpolate(latitude,doy,ZB);
      double lam = NB_Interpolate(latitude,doy,ZL);
      double gm = 9.784*(1.0-2.66e-3*std::cos(2.0*latitude*DEG_TO_RAD)-2.8e-7*height);

         // scale factors for height above mean sea level
         // if weather is given, assume it's measured at ht -> kw=kd=1
      double kw=1, base=std::log(1.0-beta*height/temp);
      if(interpolateWeather)
         kw = std::exp(base*(-1.0+(lam+1)*NBg/(NBRd*beta)));

         // compute the zenith delay for wet component
      return ((1.0e-6*NBk3p*NBRd/(gm*(lam+1)-beta*NBRd)) * kw * humid/temp);

   }  // end NBTropModel::wet_zenith_delay()

      // Compute and return the mapping function for dry component
      // of the troposphere
      // @param elevation Elevation of satellite as seen at receiver,
      //                  in degrees
   double NBTropModel::dry_mapping_function(double elevation) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
         if(!validRxLatitude)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
         if(!validRxHeight)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
         if(!validDOY)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
      }
      if(elevation < 0.0) return 0.0;

      double a,b,c,se,map;
      se = std::sin(elevation*DEG_TO_RAD);

      a = NB_Interpolate(latitude,doy,Mad);
      b = NB_Interpolate(latitude,doy,Mbd);
      c = NB_Interpolate(latitude,doy,Mcd);
      map = (1.0+a/(1.0+b/(1.0+c))) / (se+a/(se+b/(se+c)));

      a = 2.53e-5;
      b = 5.49e-3;
      c = 1.14e-3;
      if(ABS(elevation)<=0.001) se=0.001;
      map += ((1.0/se)-(1.0+a/(1.0+b/(1.0+c)))/(se+a/(se+b/(se+c))))*height/1000.0;

      return map;

   }  // end NBTropModel::dry_mapping_function()

      // Compute and return the mapping function for wet component
      // of the troposphere
      // @param elevation Elevation of satellite as seen at receiver,
      //                  in degrees
   double NBTropModel::wet_mapping_function(double elevation) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
         if(!validRxLatitude)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
         if(!validRxHeight)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
         if(!validDOY)
            GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
      }
      if(elevation < 0.0) return 0.0;

      double a,b,c,se;
      se = std::sin(elevation*DEG_TO_RAD);
      a = NB_Interpolate(latitude,doy,Maw);
      b = NB_Interpolate(latitude,doy,Mbw);
      c = NB_Interpolate(latitude,doy,Mcw);

      return ( (1.0+a/(1.0+b/(1.0+c))) / (se+a/(se+b/(se+c))) );

   }  // end NBTropModel::wet_mapping_function()

      // Re-define the weather data.
      // If called, typically called before any calls to correction().
      // @param T temperature in degrees Celsius
      // @param P atmospheric pressure in millibars
      // @param H relative humidity in percent
   void NBTropModel::setWeather(const double& T,
                                const double& P,
                                const double& H)
      throw(InvalidParameter)
   {
      interpolateWeather=false;
      TropModel::setWeather(T,P,H);
            // humid actually stores water vapor partial pressure
      double th=300./temp;
      humid = 2.409e9*H*th*th*th*th*std::exp(-22.64*th);
      validWeather = true;
      valid = validWeather && validRxHeight && validRxLatitude && validDOY;

   }  // end NBTropModel::setWeather()
   
      // Re-define the tropospheric model with explicit weather data.
      // Typically called just before correction().
      // @param wx the weather to use for this correction       
   void NBTropModel::setWeather(const WxObservation& wx)
      throw(InvalidParameter)
   {
      interpolateWeather = false;
      try
      {
         TropModel::setWeather(wx);
            // humid actually stores vapor partial pressure
         double th=300./temp;
         humid = 2.409e9*humid*th*th*th*th*std::exp(-22.64*th);
         validWeather = true;         
         valid = validWeather && validRxHeight && validRxLatitude && validDOY;
      }
      catch(InvalidParameter& e)
      {
         valid = validWeather = false;
         GPSTK_RETHROW(e);
      }
   }
   
      // configure the model to estimate the weather from the internal model,
      // using lat and doy
   void NBTropModel::setWeather()
      throw(TropModel::InvalidTropModel)
   {
      interpolateWeather = true;
      if(!validRxLatitude)
      {
         valid = validWeather = false;
         GPSTK_THROW(InvalidTropModel(
            "NBTropModel must have Rx latitude before interpolating weather"));
      }
      if(!validDOY)
      {
         valid = validWeather = false;
         GPSTK_THROW(InvalidTropModel(
            "NBTropModel must have day of year before interpolating weather"));
      }
      temp = NB_Interpolate(latitude,doy,ZT);
      press = NB_Interpolate(latitude,doy,ZP);
      humid = NB_Interpolate(latitude,doy,ZW);
      validWeather = true;
      valid = validWeather && validRxHeight && validRxLatitude && validDOY;
   }

      // Define the receiver height; this required before calling
      // correction() or any of the zenith_delay or mapping_function routines.
   void NBTropModel::setReceiverHeight(const double& ht)
   {
      height = ht;
      validRxHeight = true;
      valid = validWeather && validRxHeight && validRxLatitude && validDOY;
      if(!validWeather && validRxLatitude && validDOY)
         setWeather();
   }  // end NBTropModel::setReceiverHeight()

      // Define the latitude of the receiver; this is required before calling
      // correction() or any of the zenith_delay or mapping_function routines.
   void NBTropModel::setReceiverLatitude(const double& lat)
   {
      latitude = lat;
      validRxLatitude = true;
      valid = validWeather && validRxHeight && validRxLatitude && validDOY;
      if(!validWeather && validRxLatitude && validDOY)
         setWeather();
   }  // end NBTropModel::setReceiverLatitude(lat)

      // Define the day of year; this is required before calling
      // correction() or any of the zenith_delay or mapping_function routines.
   void NBTropModel::setDayOfYear(const int& d)
   {
      doy = d;
      if(doy > 0 && doy < 367) validDOY=true; else validDOY = false;
      valid = validWeather && validRxHeight && validRxLatitude && validDOY;
      if(!validWeather && validRxLatitude && validDOY)
         setWeather();
   }  // end NBTropModel::setDayOfYear(doy)

   // -------------------------------------------------------------------------------
   // Saastamoinen tropospheric model based on Saastamoinen, J., 'Atmospheric
   // Correction for the Troposphere and Stratosphere in Radio Ranging of
   // Satellites,' Geophysical Monograph 15, American Geophysical Union, 1972,
   // and Ch. 9 of McCarthy, D. and Petit, G., IERS Conventions (2003), IERS
   // Technical Note 32, IERS, 2004. The mapping functions are from
   // Neill, A.E., 1996, 'Global Mapping Functions for the Atmosphere Delay of
   // Radio Wavelengths,' J. Geophys. Res., 101, pp. 3227-3246 (also see IERS TN 32).
   //
   // This model includes a wet and dry component, and requires input of the
   // geodetic latitude, day of year and height above the ellipsoid of the receiver.
   //
   // Usually, the caller will set the latitude and day of year at the same
   // time the weather is set
   //   SaasTropModel stm;
   //   stm.setReceiverLatitude(lat);
   //   stm.setDayOfYear(doy);
   //   stm.setWeather(T,P,H);
   // Then, when the correction (and/or delay and map) is computed, receiver height
   // should be set before the call to correction(elevation):
   //   stm.setReceiverHeight(height);
   //   trop_corr = stm.correction(elevation);
   //
   // NB in this model, units of 'temp' are degrees Celcius and
   // humid actually stores water vapor partial pressure in mbars
   //

   // constants for wet mapping function
   static const double SaasWetA[5]=
     { 0.00058021897, 0.00056794847, 0.00058118019, 0.00059727542, 0.00061641693 };
   static const double SaasWetB[5]=
     { 0.0014275268, 0.0015138625, 0.0014572752, 0.0015007428, 0.0017599082 };
   static const double SaasWetC[5]=
     { 0.043472961, 0.046729510, 0.043908931, 0.044626982, 0.054736038 };

   // constants for dry mapping function
   static const double SaasDryA[5]=
     { 0.0012769934, 0.0012683230, 0.0012465397, 0.0012196049, 0.0012045996 };
   static const double SaasDryB[5]=
     { 0.0029153695, 0.0029152299, 0.0029288445, 0.0029022565, 0.0029024912 };
   static const double SaasDryC[5]=
     { 0.062610505, 0.062837393, 0.063721774, 0.063824265, 0.064258455 };

   static const double SaasDryA1[5]=
     { 0.0, 0.000012709626, 0.000026523662, 0.000034000452, 0.000041202191 };
   static const double SaasDryB1[5]=
     { 0.0, 0.000021414979, 0.000030160779, 0.000072562722, 0.00011723375 };
   static const double SaasDryC1[5]=
     { 0.0, 0.000090128400, 0.000043497037, 0.00084795348, 0.0017037206 };

      // Default constructor
   SaasTropModel::SaasTropModel(void)
   {
      validWeather = false;
      validRxLatitude = false;
      validDOY = false;
      validRxHeight = false;
   } // end SaasTropModel::SaasTropModel()

      // Create a trop model using the minimum information: latitude and doy.
      // Interpolate the weather unless setWeather (optional) is called.
      // @param lat Latitude of the receiver in degrees.
      // @param day Day of year.
   SaasTropModel::SaasTropModel(const double& lat,
                                const int& day)
   {
      validWeather = false;
      validRxHeight = false;
      SaasTropModel::setReceiverLatitude(lat);
      SaasTropModel::setDayOfYear(day);
   } // end SaasTropModel::SaasTropModel

      // Create a trop model with weather.