tropmodel.cpp

GPS定位中对流层模型改正

   {
      if(!valid)
         GPSTK_THROW(InvalidTropModel("Invalid model"));

      return (Cwetdelay * GGwetscale);
   }  // end GGTropModel::wet_zenith_delay()

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

      if(elevation < 0.0) return 0.0;

      GPSGeoid geoid;
      double ce=std::cos(elevation*DEG_TO_RAD), se=std::sin(elevation*DEG_TO_RAD);
      double ad = -se/Cdrymap;
      double bd = -ce*ce/(2.0*geoid.a()*Cdrymap);
      double Rd = SQRT((geoid.a()+Cdrymap)*(geoid.a()+Cdrymap)
                - geoid.a()*geoid.a()*ce*ce) - geoid.a()*se;

      double Ad[9], ad2=ad*ad, bd2=bd*bd;
      Ad[0] = 1.0;
      Ad[1] = 4.0*ad;
      Ad[2] = 6.0*ad2 + 4.0*bd;
      Ad[3] = 4.0*ad*(ad2+3.0*bd);
      Ad[4] = ad2*ad2 + 12.0*ad2*bd + 6.0*bd2;
      Ad[5] = 4.0*ad*bd*(ad2+3.0*bd);
      Ad[6] = bd2*(6.0*ad2+4.0*bd);
      Ad[7] = 4.0*ad*bd*bd2;
      Ad[8] = bd2*bd2;

         // compute dry component of the mapping function
      double sumd=0.0;
      for(int j=9; j>=1; j--) {
         sumd += Ad[j-1]/double(j);
         sumd *= Rd;
      }
      return sumd/GGdryscale;

   }  // end GGTropModel::dry_mapping_function(elev)

      // compute wet component of the mapping function
   double GGTropModel::wet_mapping_function(double elevation) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid)
         GPSTK_THROW(InvalidTropModel("Invalid model"));

      if(elevation < 0.0) return 0.0;

      GPSGeoid geoid;
      double ce = std::cos(elevation*DEG_TO_RAD), se = std::sin(elevation*DEG_TO_RAD);
      double aw = -se/Cwetmap;
      double bw = -ce*ce/(2.0*geoid.a()*Cwetmap);
      double Rw = SQRT((geoid.a()+Cwetmap)*(geoid.a()+Cwetmap)
                - geoid.a()*geoid.a()*ce*ce) - geoid.a()*se;

      double Aw[9], aw2=aw*aw, bw2=bw*bw;
      Aw[0] = 1.0;
      Aw[1] = 4.0*aw;
      Aw[2] = 6.0*aw2 + 4.0*bw;
      Aw[3] = 4.0*aw*(aw2+3.0*bw);
      Aw[4] = aw2*aw2 + 12.0*aw2*bw + 6.0*bw2;
      Aw[5] = 4.0*aw*bw*(aw2+3.0*bw);
      Aw[6] = bw2*(6.0*aw2+4.0*bw);
      Aw[7] = 4.0*aw*bw*bw2;
      Aw[8] = bw2*bw2;

      double sumw=0.0;
      for(int j=9; j>=1; j--) {
         sumw += Aw[j-1]/double(j);
         sumw *= Rw;
      }
      return sumw/GGwetscale;

   }  // end GGTropModel::wet_mapping_function(elev)

   void GGTropModel::setWeather(const double& T,
                                const double& P,
                                const double& H)
      throw(InvalidParameter)
   {
      TropModel::setWeather(T,P,H);
      double th=300./temp;
         // water vapor partial pressure (mb)
         // this comes from Leick and is not good.
         // double wvpp=6.108*(RHum*0.01)*exp((17.15*Tk-4684.0)/(Tk-38.45));
      double wvpp=2.409e9*humid*th*th*th*th*std::exp(-22.64*th);
      Cdrydelay = 7.7624e-5*press/temp;
      Cwetdelay = 1.0e-6*(-12.92+3.719e+05/temp)*(wvpp/temp);
      Cdrymap = (5.0*0.002277*press)/Cdrydelay;
      Cwetmap = (5.0*0.002277/Cwetdelay)*(1255.0/temp+0.5)*wvpp;
      valid = true;
   }  // end GGTropModel::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 GGTropModel::setWeather(const WxObservation& wx)
      throw(InvalidParameter)
   {
      TropModel::setWeather(wx);
   }

   // -------------------------------------------------------------------------------
   // Tropospheric model with heights 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.
   // (Not the same as GGTropModel because this has height dependence,
   // and the computation of this model does not break cleanly into
   // wet and dry components.)

      // Default constructor
   GGHeightTropModel::GGHeightTropModel(void)
   {
      validWeather = false; //setWeather(20.0,980.0,50.0);
      validHeights = false; //setHeights(0.0,0.0,0.0);
      validRxHeight = false;
   }

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

      // 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
   GGHeightTropModel::GGHeightTropModel(const double& T,
                                        const double& P,
                                        const double& H)
      throw(InvalidParameter)
   {
      validRxHeight = validHeights = false;
      setWeather(T,P,H);
   }

      // Create a valid model from explicit input.
      // @param T temperature in degrees Celsius
      // @param P atmospheric pressure in millibars
      // @param H relative humidity in percent
      // @param hT height at which temperature applies in meters.
      // @param hP height at which atmospheric pressure applies in meters.
      // @param hH height at which relative humidity applies in meters.
   GGHeightTropModel::GGHeightTropModel(const double& T,
                                        const double& P,
                                        const double& H,
                                        const double hT,
                                        const double hP,
                                        const double hH)
      throw(InvalidParameter)
   {
      validRxHeight = false;
      setWeather(T,P,H);
      setHeights(hT,hP,hH);
   }

      // re-define this to get the throws
   double GGHeightTropModel::correction(double elevation) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid)
      {
         if(!validWeather)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
         if(!validHeights)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
         if(!validRxHeight)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
      }
      if(elevation < 0.0) return 0.0;

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

   }  // end GGHeightTropModel::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 GGHeightTropModel::correction(const Position& RX,
                                        const Position& SV,
                                        const DayTime& tt)
      throw(TropModel::InvalidTropModel)
   {
      if(!valid)
      {
         if(!validWeather)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
         if(!validHeights)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
         if(!validRxHeight)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
      }

      // compute height from RX
      setReceiverHeight(RX.getHeight());

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

   }  // end GGHeightTropModel::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 GGHeightTropModel::correction(const Xvt& RX,
                                        const Xvt& SV,
                                        const DayTime& tt)
      throw(TropModel::InvalidTropModel)
   {
      Position R(RX),S(SV);
      return GGHeightTropModel::correction(R,S,tt);
   }

      // Compute and return the zenith delay for dry component of the troposphere
   double GGHeightTropModel::dry_zenith_delay(void) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
         if(!validHeights)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
         if(!validRxHeight)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
      }
      double hrate=6.5e-3;
      double Ts=temp+hrate*height;
      double em=978.77/(2.8704e4*hrate);
      double Tp=Ts-hrate*hpress;
      double ps=press*std::pow(Ts/Tp,em)/1000.0;
      double rs=77.624e-3/Ts;
      double ho=11.385/rs;
      rs *= ps;
      double zen=(ho-height)/ho;
      zen = rs*zen*zen*zen*zen;
         // normalize
      zen *= (ho-height)/5;
      return zen;

   }  // end GGHeightTropModel::dry_zenith_delay()

      // Compute and return the zenith delay for wet component of the troposphere
   double GGHeightTropModel::wet_zenith_delay(void) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
         if(!validHeights)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
         if(!validRxHeight)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
      }
      double hrate=6.5e-3; //   deg K / m
      double Th=temp-273.15-hrate*(hhumid-htemp);
      double Ta=7.5*Th/(237.3+Th);
         // water vapor partial pressure
      double e0=6.11e-5*humid*std::pow(10.0,Ta);
      double Ts=temp+hrate*htemp;
      double em=978.77/(2.8704e4*hrate);
      double Tk=Ts-hrate*hhumid;
      double es=e0*std::pow(Ts/Tk,4.0*em);
      double rs=(371900.0e-3/Ts-12.92e-3)/Ts;
      double ho=11.385*(1255/Ts+0.05)/rs;
      double zen=(ho-height)/ho;
      zen = rs*es*zen*zen*zen*zen;
         //normalize
      zen *= (ho-height)/5;
      return zen;
      
   }  // end GGHeightTropModel::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 GGHeightTropModel::dry_mapping_function(double elevation) const
      throw(TropModel::InvalidTropModel)
   {
      if(!valid) {
         if(!validWeather)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
         if(!validHeights)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
         if(!validRxHeight)
            GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
      }
      if(elevation < 0.0) return 0.0;

      double hrate=6.5e-3;
      double Ts=temp+hrate*htemp;
      double ho=(11.385/77.624e-3)*Ts;
      double se=std::sin(elevation*DEG_TO_RAD);
      if(se < 0.0) se=0.0;

      GPSGeoid geoid;
      double rt,a,b,rn[8],al[8],er=geoid.a();
      rt = (er+ho)/(er+height);
      rt = rt*rt - (1.0-se*se);
      if(rt < 0) rt=0.0;
      rt = (er+height)*(SQRT(rt)-se);
      a = -se/(ho-height);
      b = -(1.0-se*se)/(2.0*er*(ho-height));
      rn[0] = rt*rt;
      for(int j=1; j<8; j++) rn[j]=rn[j-1]*rt;
      al[0] = 2*a;
      al[1] = 2*a*a+4*b/3;
      al[2] = a*(a*a+3*b);
      al[3] = a*a*a*a/5+2.4*a*a*b+1.2*b*b;
      al[4] = 2*a*b*(a*a+3*b)/3;
      al[5] = b*b*(6*a*a+4*b)*0.1428571;
      if(b*b > 1.0e-35) {
         al[6] = a*b*b*b/2;
         al[7] = b*b*b*b/9;
      } else {
         al[6] = 0.0;
         al[7] = 0.0;
      }
      double map=rt;
      for(int k=0; k<8; k++) map += al[k]*rn[k];