tropmodel.cpp

GPS定位中对流层模型改正

      // @param lat Latitude of the receiver in degrees.
      // @param day Day of year.
      // @param wx the weather to use for this correction.
   SaasTropModel::SaasTropModel(const double& lat,
                                const int& day,
                                const WxObservation& wx)
      throw(InvalidParameter)
   {
      validRxHeight = false;
      SaasTropModel::setReceiverLatitude(lat);
      SaasTropModel::setDayOfYear(day);
      SaasTropModel::setWeather(wx);
   }  // end SaasTropModel::SaasTropModel(weather)

      // Create a tropospheric model from explicit weather data
      // @param lat Latitude of the receiver in degrees.
      // @param day Day of year.
      // @param T temperature in degrees Celsius
      // @param P atmospheric pressure in millibars
      // @param H relative humidity in percent
   SaasTropModel::SaasTropModel(const double& lat,
                                const int& day,
                                const double& T,
                                const double& P,
                                const double& H)
      throw(InvalidParameter)
   {
      validRxHeight = false;
      SaasTropModel::setReceiverLatitude(lat);
      SaasTropModel::setDayOfYear(day);
      SaasTropModel::setWeather(T,P,H);
   } // end SaasTropModel::SaasTropModel()

      // re-define this to get the throws correct
   double SaasTropModel::correction(double elevation) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: weather"));
         if(!validRxLatitude) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
         if(!validRxHeight) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
         if(!validDOY) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
         GPSTK_THROW(
            InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
      }

      if(elevation < 0.0) return 0.0;

      double corr=0.0;
      try {
         corr = (dry_zenith_delay() * dry_mapping_function(elevation)
            + wet_zenith_delay() * wet_mapping_function(elevation));
      }
      catch(Exception& e) { GPSTK_RETHROW(e); }

      return corr;

   }  // end SaasTropModel::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
      // @param SV  Satellite position
      // @param tt  Time tag of the signal 
   double SaasTropModel::correction(const Position& RX,
                                    const Position& SV,
                                    const DayTime& tt)
      throw(TropModel::InvalidTropModel)
   {
      SaasTropModel::setReceiverHeight(RX.getHeight());
      SaasTropModel::setReceiverLatitude(RX.getGeodeticLatitude());
      SaasTropModel::setDayOfYear(int(tt.DOYday()));

      if(!valid) {
         if(!validWeather) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: weather"));
         if(!validRxLatitude) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
         if(!validRxHeight) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
         if(!validDOY) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
         valid = true;
      }

      double corr=0.0;
      try {
         corr = SaasTropModel::correction(RX.elevation(SV));
      }
      catch(Exception& e) { GPSTK_RETHROW(e); }

      return corr;

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

   double SaasTropModel::correction(const Xvt& RX,
                                    const Xvt& SV,
                                    const DayTime& tt)
      throw(TropModel::InvalidTropModel)
   {
      Position R(RX),S(SV);
      return SaasTropModel::correction(R,S,tt);
   }

      // Compute and return the zenith delay for dry component of the troposphere
   double SaasTropModel::dry_zenith_delay(void) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: weather"));
         if(!validRxLatitude) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
         if(!validRxHeight) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
         if(!validDOY) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
         GPSTK_THROW(   
            InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
      }

      // correct pressure for height
      double press_at_h =
         press * std::pow((temp+273.16-4.5*height/1000.0)/(temp+273.16),34.1/4.5);
      // humid is zero for the dry component
      double delay = 0.0022768 * press_at_h
            / (1 - 0.00266 * ::cos(2*latitude*DEG_TO_RAD) - 0.00028 * height/1000.);

      return delay;

   }  // end SaasTropModel::dry_zenith_delay()

      // Compute and return the zenith delay for wet component of the troposphere
   double SaasTropModel::wet_zenith_delay(void) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: weather"));
         if(!validRxLatitude) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
         if(!validRxHeight) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
         if(!validDOY) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
         GPSTK_THROW(
            InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
      }

      // press is zero for the wet component
      double delay = 0.0022768 * humid * 1255/(temp+273.20)
            / (1 - 0.00266 * ::cos(2*latitude*DEG_TO_RAD) - 0.00028 * height/1000.);

      return delay;

   }  // end SaasTropModel::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 SaasTropModel::dry_mapping_function(double elevation) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: weather"));
         if(!validRxLatitude) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
         if(!validRxHeight) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
         if(!validDOY) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
         GPSTK_THROW(
            InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
      }
      if(elevation < 0.0) return 0.0;

      double lat,t,ct;
      lat = fabs(latitude);         // degrees
      t = doy - 28.;                // mid-winter
      if(latitude < 0)              // southern hemisphere
         t += 365.25/2.;
      t *= 360.0/365.25;            // convert to degrees
      ct = ::cos(t*DEG_TO_RAD);

      double a,b,c;
      if(lat < 15.) {
         a = SaasDryA[0];
         b = SaasDryB[0];
         c = SaasDryC[0];
      }
      else if(lat < 75.) {          // coefficients are for 15,30,45,60,75 deg
         int i=int(lat/15.0)-1;
         double frac=(lat-15.*(i+1))/15.;
         a = SaasDryA[i] + frac*(SaasDryA[i+1]-SaasDryA[i]);
         b = SaasDryB[i] + frac*(SaasDryB[i+1]-SaasDryB[i]);
         c = SaasDryC[i] + frac*(SaasDryC[i+1]-SaasDryC[i]);

         a -= ct * (SaasDryA1[i] + frac*(SaasDryA1[i+1]-SaasDryA1[i]));
         b -= ct * (SaasDryB1[i] + frac*(SaasDryB1[i+1]-SaasDryB1[i]));
         c -= ct * (SaasDryC1[i] + frac*(SaasDryC1[i+1]-SaasDryC1[i]));
      }
      else {
         a = SaasDryA[4] - ct * SaasDryA1[4];
         b = SaasDryB[4] - ct * SaasDryB1[4];
         c = SaasDryC[4] - ct * SaasDryC1[4];
      }

      double se = ::sin(elevation*DEG_TO_RAD);
      double map = (1.+a/(1.+b/(1.+c)))/(se+a/(se+b/(se+c)));

      a = 0.0000253;
      b = 0.00549;
      c = 0.00114;
      map += (height/1000.0)*(1./se-(1+a/(1.+b/(1.+c)))/(se+a/(se+b/(se+c))));

      return map;

   }  // end SaasTropModel::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 SaasTropModel::wet_mapping_function(double elevation) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: weather"));
         if(!validRxLatitude) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
         if(!validRxHeight) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
         if(!validDOY) GPSTK_THROW(
            InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
         GPSTK_THROW(
            InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
      }
      if(elevation < 0.0) return 0.0;

      double a,b,c,lat;
      lat = fabs(latitude);         // degrees
      if(lat < 15.) {
         a = SaasWetA[0];
         b = SaasWetB[0];
         c = SaasWetC[0];
      }
      else if(lat < 75.) {          // coefficients are for 15,30,45,60,75 deg
         int i=int(lat/15.0)-1;
         double frac=(lat-15.*(i+1))/15.;
         a = SaasWetA[i] + frac*(SaasWetA[i+1]-SaasWetA[i]);
         b = SaasWetB[i] + frac*(SaasWetB[i+1]-SaasWetB[i]);
         c = SaasWetC[i] + frac*(SaasWetC[i+1]-SaasWetC[i]);
      }
      else {
         a = SaasWetA[4];
         b = SaasWetB[4];
         c = SaasWetC[4];
      }

      double se = ::sin(elevation*DEG_TO_RAD);
      double map = (1.+a/(1.+b/(1.+c)))/(se+a/(se+b/(se+c)));

      return map;

   }  // end SaasTropModel::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 SaasTropModel::setWeather(const double& T,
                                  const double& P,
                                  const double& H)
      throw(InvalidParameter)
   {
      temp = T;
      press = P;
         // humid actually stores water vapor partial pressure
      double exp=7.5*T/(T+237.3);
      humid = 6.11 * (H/100.) * std::pow(10.0,exp);

      validWeather = true;
      valid = (validWeather && validRxHeight && validRxLatitude && validDOY);

   }  // end SaasTropModel::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 SaasTropModel::setWeather(const WxObservation& wx)
      throw(InvalidParameter)
   {
      try
      {
         SaasTropModel::setWeather(wx.temperature,wx.pressure,wx.humidity);
      }
      catch(InvalidParameter& e)
      {
         valid = validWeather = false;
         GPSTK_RETHROW(e);
      }
   }
   
      // Define the receiver height; this required before calling
      // correction() or any of the zenith_delay or mapping_function routines.
   void SaasTropModel::setReceiverHeight(const double& ht)
   {
      height = ht;
      validRxHeight = true;
      valid = (validWeather && validRxHeight && validRxLatitude && validDOY);
   }  // end SaasTropModel::setReceiverHeight()

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

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


    //-----