00001 #pragma ident "$Id: TropModel.cpp 1431 2008-10-31 16:43:03Z btolman $"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00045 #include "TropModel.hpp"
00046 #include "EphemerisRange.hpp"
00047 #include "MathBase.hpp"
00048 #include "geometry.hpp"
00049 #include "GPSGeoid.hpp"
00050 #include "icd_200_constants.hpp"
00051 #include "Geodetic.hpp"
00052 #include "ECEF.hpp"
00053
00054 namespace gpstk
00055 {
00056
00057 static const double CELSIUS_TO_KELVIN = 273.15;
00058
00059
00060
00061
00062 double TropModel::correction(double elevation) const
00063 throw(TropModel::InvalidTropModel)
00064 {
00065 if(!valid)
00066 GPSTK_THROW(InvalidTropModel("Invalid model"));
00067
00068 if(elevation < 0.0)
00069 return 0.0;
00070
00071 return (dry_zenith_delay() * dry_mapping_function(elevation)
00072 + wet_zenith_delay() * wet_mapping_function(elevation));
00073
00074 }
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085 double TropModel::correction(const Position& RX,
00086 const Position& SV,
00087 const DayTime& tt)
00088 throw(TropModel::InvalidTropModel)
00089 {
00090 if(!valid)
00091 GPSTK_THROW(InvalidTropModel("Invalid model"));
00092
00093 double c;
00094 try
00095 {
00096 c = correction(RX.elevation(SV));
00097 }
00098 catch(InvalidTropModel& e)
00099 {
00100 GPSTK_RETHROW(e);
00101 }
00102 return c;
00103 }
00104
00105
00106
00107
00108
00109
00110 void TropModel::setWeather(const double& T,
00111 const double& P,
00112 const double& H)
00113 throw(InvalidParameter)
00114 {
00115 temp = T + CELSIUS_TO_KELVIN;
00116 press = P;
00117 humid = H;
00118 if (temp < 0.0)
00119 {
00120 valid = false;
00121 GPSTK_THROW(InvalidParameter("Invalid temperature parameter."));
00122 }
00123 if (press < 0.0)
00124 {
00125 valid = false;
00126 GPSTK_THROW(InvalidParameter("Invalid pressure parameter."));
00127 }
00128 if (humid < 0.0 || humid > 100.0)
00129 {
00130 valid = false;
00131 GPSTK_THROW(InvalidParameter("Invalid humidity parameter."));
00132 }
00133 }
00134
00135
00136
00137
00138 void TropModel::setWeather(const WxObservation& wx)
00139 throw(InvalidParameter)
00140 {
00141 if (wx.isAllValid())
00142 {
00143 try
00144 {
00145 setWeather(wx.temperature, wx.pressure, wx.humidity);
00146 valid = true;
00147 }
00148 catch(InvalidParameter& e)
00149 {
00150 valid = false;
00151 GPSTK_RETHROW(e);
00152 }
00153 }
00154 else
00155 {
00156 valid = false;
00157 GPSTK_THROW(InvalidParameter("Invalid weather data"));
00158 }
00159 }
00160
00161
00162
00163
00164
00165 SimpleTropModel::SimpleTropModel(void)
00166 {
00167 setWeather(20.0, 980.0, 50.0);
00168 Cwetdelay = 0.122382715318184;
00169 Cdrydelay = 2.235486646978727;
00170 Cwetmap = 1.000282213715744;
00171 Cdrymap = 1.001012704615527;
00172 valid = true;
00173 }
00174
00175
00176
00177 SimpleTropModel::SimpleTropModel(const WxObservation& wx)
00178 throw(InvalidParameter)
00179 {
00180 setWeather(wx);
00181 valid = true;
00182 }
00183
00184
00185
00186
00187
00188 SimpleTropModel::SimpleTropModel(const double& T,
00189 const double& P,
00190 const double& H)
00191 throw(InvalidParameter)
00192 {
00193 setWeather(T,P,H);
00194 valid = true;
00195 }
00196
00197
00198
00199
00200
00201
00202 void SimpleTropModel::setWeather(const double& T,
00203 const double& P,
00204 const double& H)
00205 throw(InvalidParameter)
00206 {
00207 TropModel::setWeather(T,P,H);
00208 GPSGeoid geoid;
00209 Cdrydelay = 2.343*(press/1013.25)*(temp-3.96)/temp;
00210 double tks = temp * temp;
00211 Cwetdelay = 8.952/tks*humid*std::exp(-37.2465+0.213166*temp-(0.256908e-3)*tks);
00212 Cdrymap =1.0+(0.15)*148.98*(temp-3.96)/geoid.a();
00213 Cwetmap =1.0+(0.15)*12000.0/geoid.a();
00214 valid = true;
00215 }
00216
00217
00218
00219
00220 void SimpleTropModel::setWeather(const WxObservation& wx)
00221 throw(InvalidParameter)
00222 {
00223 TropModel::setWeather(wx);
00224 }
00225
00226
00227 double SimpleTropModel::dry_zenith_delay(void) const
00228 throw(TropModel::InvalidTropModel)
00229 {
00230 if(!valid)
00231 GPSTK_THROW(InvalidTropModel("Invalid model"));
00232
00233 return Cdrydelay;
00234 }
00235
00236
00237 double SimpleTropModel::wet_zenith_delay(void) const
00238 throw(TropModel::InvalidTropModel)
00239 {
00240 if(!valid)
00241 GPSTK_THROW(InvalidTropModel("Invalid model"));
00242
00243 return Cwetdelay;
00244 }
00245
00246
00247
00248
00249
00250 double SimpleTropModel::dry_mapping_function(double elevation) const
00251 throw(TropModel::InvalidTropModel)
00252 {
00253 if(!valid)
00254 GPSTK_THROW(InvalidTropModel("Invalid model"));
00255
00256 if(elevation < 0.0) return 0.0;
00257
00258 double d = std::cos(elevation*DEG_TO_RAD);
00259 d /= Cdrymap;
00260 return (1.0/SQRT(1.0-d*d));
00261 }
00262
00263
00264
00265
00266
00267 double SimpleTropModel::wet_mapping_function(double elevation) const
00268 throw(TropModel::InvalidTropModel)
00269 {
00270 if(!valid)
00271 GPSTK_THROW(InvalidTropModel("Invalid model"));
00272
00273 if(elevation < 0.0) return 0.0;
00274
00275 double d = std::cos(elevation*DEG_TO_RAD);
00276 d /= Cwetmap;
00277 return (1.0/SQRT(1.0-d*d));
00278 }
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289 static const double GGdryscale = 8594.777388436570600;
00290 static const double GGwetscale = 2540.042008403690900;
00291
00292
00293 GGTropModel::GGTropModel(void)
00294 {
00295 TropModel::setWeather(20.0, 980.0, 50.0);
00296 Cdrydelay = 2.59629761092150147e-4;
00297 Cwetdelay = 4.9982784999977412e-5;
00298 Cdrymap = 42973.886942182834900;
00299 Cwetmap = 12700.210042018454260;
00300 valid = true;
00301 }
00302
00303
00304
00305 GGTropModel::GGTropModel(const WxObservation& wx)
00306 throw(InvalidParameter)
00307 {
00308 setWeather(wx);
00309 valid = true;
00310 }
00311
00312
00313
00314
00315
00316 GGTropModel::GGTropModel(const double& T,
00317 const double& P,
00318 const double& H)
00319 throw(InvalidParameter)
00320 {
00321 setWeather(T,P,H);
00322 valid = true;
00323 }
00324
00325 double GGTropModel::dry_zenith_delay(void) const
00326 throw(TropModel::InvalidTropModel)
00327 {
00328 if(!valid)
00329 GPSTK_THROW(InvalidTropModel("Invalid model"));
00330
00331 return (Cdrydelay * GGdryscale);
00332 }
00333
00334 double GGTropModel::wet_zenith_delay(void) const
00335 throw(TropModel::InvalidTropModel)
00336 {
00337 if(!valid)
00338 GPSTK_THROW(InvalidTropModel("Invalid model"));
00339
00340 return (Cwetdelay * GGwetscale);
00341 }
00342
00343 double GGTropModel::dry_mapping_function(double elevation) const
00344 throw(TropModel::InvalidTropModel)
00345 {
00346 if(!valid)
00347 GPSTK_THROW(InvalidTropModel("Invalid model"));
00348
00349 if(elevation < 0.0) return 0.0;
00350
00351 GPSGeoid geoid;
00352 double ce=std::cos(elevation*DEG_TO_RAD), se=std::sin(elevation*DEG_TO_RAD);
00353 double ad = -se/Cdrymap;
00354 double bd = -ce*ce/(2.0*geoid.a()*Cdrymap);
00355 double Rd = SQRT((geoid.a()+Cdrymap)*(geoid.a()+Cdrymap)
00356 - geoid.a()*geoid.a()*ce*ce) - geoid.a()*se;
00357
00358 double Ad[9], ad2=ad*ad, bd2=bd*bd;
00359 Ad[0] = 1.0;
00360 Ad[1] = 4.0*ad;
00361 Ad[2] = 6.0*ad2 + 4.0*bd;
00362 Ad[3] = 4.0*ad*(ad2+3.0*bd);
00363 Ad[4] = ad2*ad2 + 12.0*ad2*bd + 6.0*bd2;
00364 Ad[5] = 4.0*ad*bd*(ad2+3.0*bd);
00365 Ad[6] = bd2*(6.0*ad2+4.0*bd);
00366 Ad[7] = 4.0*ad*bd*bd2;
00367 Ad[8] = bd2*bd2;
00368
00369
00370 double sumd=0.0;
00371 for(int j=9; j>=1; j--) {
00372 sumd += Ad[j-1]/double(j);
00373 sumd *= Rd;
00374 }
00375 return sumd/GGdryscale;
00376
00377 }
00378
00379
00380 double GGTropModel::wet_mapping_function(double elevation) const
00381 throw(TropModel::InvalidTropModel)
00382 {
00383 if(!valid)
00384 GPSTK_THROW(InvalidTropModel("Invalid model"));
00385
00386 if(elevation < 0.0) return 0.0;
00387
00388 GPSGeoid geoid;
00389 double ce = std::cos(elevation*DEG_TO_RAD), se = std::sin(elevation*DEG_TO_RAD);
00390 double aw = -se/Cwetmap;
00391 double bw = -ce*ce/(2.0*geoid.a()*Cwetmap);
00392 double Rw = SQRT((geoid.a()+Cwetmap)*(geoid.a()+Cwetmap)
00393 - geoid.a()*geoid.a()*ce*ce) - geoid.a()*se;
00394
00395 double Aw[9], aw2=aw*aw, bw2=bw*bw;
00396 Aw[0] = 1.0;
00397 Aw[1] = 4.0*aw;
00398 Aw[2] = 6.0*aw2 + 4.0*bw;
00399 Aw[3] = 4.0*aw*(aw2+3.0*bw);
00400 Aw[4] = aw2*aw2 + 12.0*aw2*bw + 6.0*bw2;
00401 Aw[5] = 4.0*aw*bw*(aw2+3.0*bw);
00402 Aw[6] = bw2*(6.0*aw2+4.0*bw);
00403 Aw[7] = 4.0*aw*bw*bw2;
00404 Aw[8] = bw2*bw2;
00405
00406 double sumw=0.0;
00407 for(int j=9; j>=1; j--) {
00408 sumw += Aw[j-1]/double(j);
00409 sumw *= Rw;
00410 }
00411 return sumw/GGwetscale;
00412
00413 }
00414
00415 void GGTropModel::setWeather(const double& T,
00416 const double& P,
00417 const double& H)
00418 throw(InvalidParameter)
00419 {
00420 TropModel::setWeather(T,P,H);
00421 double th=300./temp;
00422
00423
00424
00425 double wvpp=2.409e9*humid*th*th*th*th*std::exp(-22.64*th);
00426 Cdrydelay = 7.7624e-5*press/temp;
00427 Cwetdelay = 1.0e-6*(-12.92+3.719e+05/temp)*(wvpp/temp);
00428 Cdrymap = (5.0*0.002277*press)/Cdrydelay;
00429 Cwetmap = (5.0*0.002277/Cwetdelay)*(1255.0/temp+0.5)*wvpp;
00430 valid = true;
00431 }
00432
00433
00434
00435
00436 void GGTropModel::setWeather(const WxObservation& wx)
00437 throw(InvalidParameter)
00438 {
00439 TropModel::setWeather(wx);
00440 }
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452 GGHeightTropModel::GGHeightTropModel(void)
00453 {
00454 validWeather = false;
00455 validHeights = false;
00456 validRxHeight = false;
00457 }
00458
00459
00460
00461 GGHeightTropModel::GGHeightTropModel(const WxObservation& wx)
00462 throw(InvalidParameter)
00463 {
00464 valid = validRxHeight = validHeights = false;
00465 setWeather(wx);
00466 }
00467
00468
00469
00470
00471
00472 GGHeightTropModel::GGHeightTropModel(const double& T,
00473 const double& P,
00474 const double& H)
00475 throw(InvalidParameter)
00476 {
00477 validRxHeight = validHeights = false;
00478 setWeather(T,P,H);
00479 }
00480
00481
00482
00483
00484
00485
00486
00487
00488 GGHeightTropModel::GGHeightTropModel(const double& T,
00489 const double& P,
00490 const double& H,
00491 const double hT,
00492 const double hP,
00493 const double hH)
00494 throw(InvalidParameter)
00495 {
00496 validRxHeight = false;
00497 setWeather(T,P,H);
00498 setHeights(hT,hP,hH);
00499 }
00500
00501
00502 double GGHeightTropModel::correction(double elevation) const
00503 throw(TropModel::InvalidTropModel)
00504 {
00505 if(!valid)
00506 {
00507 if(!validWeather)
00508 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
00509 if(!validHeights)
00510 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
00511 if(!validRxHeight)
00512 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
00513 }
00514 if(elevation < 0.0) return 0.0;
00515
00516 return (dry_zenith_delay() * dry_mapping_function(elevation)
00517 + wet_zenith_delay() * wet_mapping_function(elevation));
00518
00519 }
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530 double GGHeightTropModel::correction(const Position& RX,
00531 const Position& SV,
00532 const DayTime& tt)
00533 throw(TropModel::InvalidTropModel)
00534 {
00535 if(!valid)
00536 {
00537 if(!validWeather)
00538 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
00539 if(!validHeights)
00540 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
00541 if(!validRxHeight)
00542 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
00543 }
00544
00545
00546 setReceiverHeight(RX.getHeight());
00547
00548 return TropModel::correction(RX.elevation(SV));
00549
00550 }
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562 double GGHeightTropModel::correction(const Xvt& RX,
00563 const Xvt& SV,
00564 const DayTime& tt)
00565 throw(TropModel::InvalidTropModel)
00566 {
00567 Position R(RX),S(SV);
00568 return GGHeightTropModel::correction(R,S,tt);
00569 }
00570
00571
00572 double GGHeightTropModel::dry_zenith_delay(void) const
00573 throw(TropModel::InvalidTropModel)
00574 {
00575 if(!valid) {
00576 if(!validWeather)
00577 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
00578 if(!validHeights)
00579 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
00580 if(!validRxHeight)
00581 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
00582 }
00583 double hrate=6.5e-3;
00584 double Ts=temp+hrate*height;
00585 double em=978.77/(2.8704e4*hrate);
00586 double Tp=Ts-hrate*hpress;
00587 double ps=press*std::pow(Ts/Tp,em)/1000.0;
00588 double rs=77.624e-3/Ts;
00589 double ho=11.385/rs;
00590 rs *= ps;
00591 double zen=(ho-height)/ho;
00592 zen = rs*zen*zen*zen*zen;
00593
00594 zen *= (ho-height)/5;
00595 return zen;
00596
00597 }
00598
00599
00600 double GGHeightTropModel::wet_zenith_delay(void) const
00601 throw(TropModel::InvalidTropModel)
00602 {
00603 if(!valid) {
00604 if(!validWeather)
00605 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
00606 if(!validHeights)
00607 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
00608 if(!validRxHeight)
00609 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
00610 }
00611 double hrate=6.5e-3;
00612 double Th=temp-273.15-hrate*(hhumid-htemp);
00613 double Ta=7.5*Th/(237.3+Th);
00614
00615 double e0=6.11e-5*humid*std::pow(10.0,Ta);
00616 double Ts=temp+hrate*htemp;
00617 double em=978.77/(2.8704e4*hrate);
00618 double Tk=Ts-hrate*hhumid;
00619 double es=e0*std::pow(Ts/Tk,4.0*em);
00620 double rs=(371900.0e-3/Ts-12.92e-3)/Ts;
00621 double ho=11.385*(1255/Ts+0.05)/rs;
00622 double zen=(ho-height)/ho;
00623 zen = rs*es*zen*zen*zen*zen;
00624
00625 zen *= (ho-height)/5;
00626 return zen;
00627
00628 }
00629
00630
00631
00632
00633
00634 double GGHeightTropModel::dry_mapping_function(double elevation) const
00635 throw(TropModel::InvalidTropModel)
00636 {
00637 if(!valid) {
00638 if(!validWeather)
00639 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
00640 if(!validHeights)
00641 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
00642 if(!validRxHeight)
00643 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
00644 }
00645 if(elevation < 0.0) return 0.0;
00646
00647 double hrate=6.5e-3;
00648 double Ts=temp+hrate*htemp;
00649 double ho=(11.385/77.624e-3)*Ts;
00650 double se=std::sin(elevation*DEG_TO_RAD);
00651 if(se < 0.0) se=0.0;
00652
00653 GPSGeoid geoid;
00654 double rt,a,b,rn[8],al[8],er=geoid.a();
00655 rt = (er+ho)/(er+height);
00656 rt = rt*rt - (1.0-se*se);
00657 if(rt < 0) rt=0.0;
00658 rt = (er+height)*(SQRT(rt)-se);
00659 a = -se/(ho-height);
00660 b = -(1.0-se*se)/(2.0*er*(ho-height));
00661 rn[0] = rt*rt;
00662 for(int j=1; j<8; j++) rn[j]=rn[j-1]*rt;
00663 al[0] = 2*a;
00664 al[1] = 2*a*a+4*b/3;
00665 al[2] = a*(a*a+3*b);
00666 al[3] = a*a*a*a/5+2.4*a*a*b+1.2*b*b;
00667 al[4] = 2*a*b*(a*a+3*b)/3;
00668 al[5] = b*b*(6*a*a+4*b)*0.1428571;
00669 if(b*b > 1.0e-35) {
00670 al[6] = a*b*b*b/2;
00671 al[7] = b*b*b*b/9;
00672 } else {
00673 al[6] = 0.0;
00674 al[7] = 0.0;
00675 }
00676 double map=rt;
00677 for(int k=0; k<8; k++) map += al[k]*rn[k];
00678
00679 double norm=(ho-height)/5;
00680 return map/norm;
00681
00682 }
00683
00684
00685
00686
00687
00688 double GGHeightTropModel::wet_mapping_function(double elevation) const
00689 throw(TropModel::InvalidTropModel)
00690 {
00691 if(!valid) {
00692 if(!validWeather)
00693 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Weather"));
00694 if(!validHeights)
00695 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Heights"));
00696 if(!validRxHeight)
00697 GPSTK_THROW(InvalidTropModel("Invalid GGH trop model: Rx Height"));
00698 }
00699 if(elevation < 0.0) return 0.0;
00700
00701 double hrate=6.5e-3;
00702 double Ts=temp+hrate*htemp;
00703 double rs=(371900.0e-3/Ts-12.92e-3)/Ts;
00704 double ho=11.385*(1255/Ts+0.05)/rs;
00705 double se=std::sin(elevation*DEG_TO_RAD);
00706 if(se < 0.0) se=0.0;
00707
00708 GPSGeoid geoid;
00709 double rt,a,b,rn[8],al[8],er=geoid.a();
00710 rt = (er+ho)/(er+height);
00711 rt = rt*rt - (1.0-se*se);
00712 if(rt < 0) rt=0.0;
00713 rt = (er+height)*(SQRT(rt)-se);
00714 a = -se/(ho-height);
00715 b = -(1.0-se*se)/(2.0*er*(ho-height));
00716 rn[0] = rt*rt;
00717 for(int i=1; i<8; i++) rn[i]=rn[i-1]*rt;
00718 al[0] = 2*a;
00719 al[1] = 2*a*a+4*b/3;
00720 al[2] = a*(a*a+3*b);
00721 al[3] = a*a*a*a/5+2.4*a*a*b+1.2*b*b;
00722 al[4] = 2*a*b*(a*a+3*b)/3;
00723 al[5] = b*b*(6*a*a+4*b)*0.1428571;
00724 if(b*b > 1.0e-35) {
00725 al[6] = a*b*b*b/2;
00726 al[7] = b*b*b*b/9;
00727 } else {
00728 al[6] = 0.0;
00729 al[7] = 0.0;
00730 }
00731 double map=rt;
00732 for(int j=0; j<8; j++) map += al[j]*rn[j];
00733
00734 double norm=(ho-height)/5;
00735 return map/norm;
00736
00737 }
00738
00739
00740
00741
00742
00743
00744 void GGHeightTropModel::setWeather(const double& T,
00745 const double& P,
00746 const double& H)
00747 throw(InvalidParameter)
00748 {
00749 try
00750 {
00751 TropModel::setWeather(T,P,H);
00752 validWeather = true;
00753 valid = validWeather && validHeights && validRxHeight;
00754 }
00755 catch(InvalidParameter& e)
00756 {
00757 valid = validWeather = false;
00758 GPSTK_RETHROW(e);
00759 }
00760 }
00761
00762
00763
00764
00765 void GGHeightTropModel::setWeather(const WxObservation& wx)
00766 throw(InvalidParameter)
00767 {
00768 try
00769 {
00770 TropModel::setWeather(wx);
00771 validWeather = true;
00772 valid = validWeather && validHeights && validRxHeight;
00773 }
00774 catch(InvalidParameter& e)
00775 {
00776 valid = validWeather = false;
00777 GPSTK_RETHROW(e);
00778 }
00779 }
00780
00781
00782
00783
00784
00785
00786
00787 void GGHeightTropModel::setHeights(const double& hT,
00788 const double& hP,
00789 const double& hH)
00790 {
00791 htemp = hT;
00792 hpress = hP;
00793 hhumid = hH;
00794 validHeights = true;
00795 valid = validWeather && validHeights && validRxHeight;
00796 }
00797
00798
00799
00800 void GGHeightTropModel::setReceiverHeight(const double& ht)
00801 {
00802 height = ht;
00803 validRxHeight = true;
00804 if(!validHeights) {
00805 htemp = hpress = hhumid = ht;
00806 validHeights = true;
00807 }
00808 valid = validWeather && validHeights && validRxHeight;
00809 }
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828 static const double NBRd=287.054;
00829 static const double NBg=9.80665;
00830 static const double NBk1=77.604;
00831 static const double NBk3p=382000;
00832
00833 static const double NBLat[5]={ 15.0, 30.0, 45.0, 60.0, 75.0};
00834
00835
00836 static const double NBZP0[5]={1013.25,1017.25,1015.75,1011.75,1013.00};
00837 static const double NBZT0[5]={ 299.65, 294.15, 283.15, 272.15, 263.65};
00838 static const double NBZW0[5]={ 26.31, 21.79, 11.66, 6.78, 4.11};
00839 static const double NBZB0[5]={6.30e-3,6.05e-3,5.58e-3,5.39e-3,4.53e-3};
00840 static const double NBZL0[5]={ 2.77, 3.15, 2.57, 1.81, 1.55};
00841
00842
00843 static const double NBZPa[5]={ 0.0, -3.75, -2.25, -1.75, -0.50};
00844 static const double NBZTa[5]={ 0.0, 7.0, 11.0, 15.0, 14.5};
00845 static const double NBZWa[5]={ 0.0, 8.85, 7.24, 5.36, 3.39};
00846 static const double NBZBa[5]={ 0.0,0.25e-3,0.32e-3,0.81e-3,0.62e-3};
00847 static const double NBZLa[5]={ 0.0, 0.33, 0.46, 0.74, 0.30};
00848
00849
00850 static const double NBMad[5]={1.2769934e-3,1.2683230e-3,1.2465397e-3,1.2196049e-3,
00851 1.2045996e-3};
00852 static const double NBMbd[5]={2.9153695e-3,2.9152299e-3,2.9288445e-3,2.9022565e-3,
00853 2.9024912e-3};
00854 static const double NBMcd[5]={62.610505e-3,62.837393e-3,63.721774e-3,63.824265e-3,
00855 64.258455e-3};
00856
00857
00858 static const double NBMaa[5]={0.0,1.2709626e-5,2.6523662e-5,3.4000452e-5,
00859 4.1202191e-5};
00860 static const double NBMba[5]={0.0,2.1414979e-5,3.0160779e-5,7.2562722e-5,
00861 11.723375e-5};
00862 static const double NBMca[5]={0.0,9.0128400e-5,4.3497037e-5,84.795348e-5,
00863 170.37206e-5};
00864
00865
00866 static const double NBMaw[5]={5.8021897e-4,5.6794847e-4,5.8118019e-4,5.9727542e-4,
00867 6.1641693e-4};
00868 static const double NBMbw[5]={1.4275268e-3,1.5138625e-3,1.4572752e-3,1.5007428e-3,
00869 1.7599082e-3};
00870 static const double NBMcw[5]={4.3472961e-2,4.6729510e-2,4.3908931e-2,4.4626982e-2,
00871 5.4736038e-2};
00872
00873
00874 enum TableEntry { ZP=1, ZT, ZW, ZB, ZL, Mad, Mbd, Mcd, Maw, Mbw, Mcw };
00875
00876
00877 static double NB_Interpolate(double lat, int doy, TableEntry entry)
00878 {
00879 const double *pave, *pamp;
00880 double ret, day=double(doy);
00881
00882
00883 switch(entry) {
00884 case ZP: pave=&NBZP0[0]; pamp=&NBZPa[0]; break;
00885 case ZT: pave=&NBZT0[0]; pamp=&NBZTa[0]; break;
00886 case ZW: pave=&NBZW0[0]; pamp=&NBZWa[0]; break;
00887 case ZB: pave=&NBZB0[0]; pamp=&NBZBa[0]; break;
00888 case ZL: pave=&NBZL0[0]; pamp=&NBZLa[0]; break;
00889 case Mad: pave=&NBMad[0]; pamp=&NBMaa[0]; break;
00890 case Mbd: pave=&NBMbd[0]; pamp=&NBMba[0]; break;
00891 case Mcd: pave=&NBMcd[0]; pamp=&NBMca[0]; break;
00892 case Maw: pave=&NBMaw[0]; break;
00893 case Mbw: pave=&NBMbw[0]; break;
00894 case Mcw: pave=&NBMcw[0]; break;
00895 }
00896
00897
00898 int i = int(ABS(lat)/15.0)-1;
00899
00900 if(i>=0 && i<=3) {
00901 double m=(ABS(lat)-NBLat[i])/(NBLat[i+1]-NBLat[i]);
00902 ret = pave[i]+m*(pave[i+1]-pave[i]);
00903 if(entry < Maw)
00904 ret -= (pamp[i]+m*(pamp[i+1]-pamp[i]))
00905 * std::cos(TWO_PI*(day-28.0)/365.25);
00906 }
00907 else {
00908 if(i<0) i=0; else i=4;
00909 ret = pave[i];
00910 if(entry < Maw)
00911 ret -= pamp[i]*std::cos(TWO_PI*(day-28.0)/365.25);
00912 }
00913
00914 return ret;
00915
00916 }
00917
00918
00919 NBTropModel::NBTropModel(void)
00920 {
00921 validWeather = false;
00922 validRxLatitude = false;
00923 validDOY = false;
00924 validRxHeight = false;
00925 }
00926
00927
00928
00929
00930
00931 NBTropModel::NBTropModel(const double& lat,
00932 const int& day)
00933 {
00934 validRxHeight = false;
00935 setReceiverLatitude(lat);
00936 setDayOfYear(day);
00937 setWeather();
00938 }
00939
00940
00941
00942
00943
00944 NBTropModel::NBTropModel(const double& lat,
00945 const int& day,
00946 const WxObservation& wx)
00947 throw(InvalidParameter)
00948 {
00949 validRxHeight = false;
00950 setReceiverLatitude(lat);
00951 setDayOfYear(day);
00952 setWeather(wx);
00953 }
00954
00955
00956
00957
00958
00959
00960
00961 NBTropModel::NBTropModel(const double& lat,
00962 const int& day,
00963 const double& T,
00964 const double& P,
00965 const double& H)
00966 throw(InvalidParameter)
00967 {
00968 validRxHeight = false;
00969 setReceiverLatitude(lat);
00970 setDayOfYear(day);
00971 setWeather(T,P,H);
00972 }
00973
00974
00975
00976
00977
00978
00979 NBTropModel::NBTropModel(const double& ht,
00980 const double& lat,
00981 const int& day)
00982 {
00983 setReceiverHeight(ht);
00984 setReceiverLatitude(lat);
00985 setDayOfYear(day);
00986 setWeather();
00987 }
00988
00989
00990 double NBTropModel::correction(double elevation) const
00991 throw(TropModel::InvalidTropModel)
00992 {
00993 if(!valid) {
00994 if(!validWeather)
00995 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
00996 if(!validRxLatitude)
00997 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
00998 if(!validRxHeight)
00999 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
01000 if(!validDOY)
01001 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
01002 }
01003
01004 if(elevation < 0.0) return 0.0;
01005
01006 return (dry_zenith_delay() * dry_mapping_function(elevation)
01007 + wet_zenith_delay() * wet_mapping_function(elevation));
01008 }
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019 double NBTropModel::correction(const Position& RX,
01020 const Position& SV,
01021 const DayTime& tt)
01022 throw(TropModel::InvalidTropModel)
01023 {
01024 if(!valid) {
01025 if(!validWeather)
01026
01027 if(!validRxLatitude)
01028 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
01029 if(!validRxHeight)
01030 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
01031 if(!validDOY)
01032 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
01033 }
01034
01035
01036 setReceiverHeight(RX.getHeight());
01037 setReceiverLatitude(RX.getGeodeticLatitude());
01038
01039
01040 setDayOfYear(int(tt.DOYday()));
01041
01042 return TropModel::correction(RX.elevation(SV));
01043
01044 }
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056 double NBTropModel::correction(const Xvt& RX,
01057 const Xvt& SV,
01058 const DayTime& tt)
01059 throw(TropModel::InvalidTropModel)
01060 {
01061 Position R(RX),S(SV);
01062 return NBTropModel::correction(R,S,tt);
01063 }
01064
01065
01066 double NBTropModel::dry_zenith_delay(void) const
01067 throw(TropModel::InvalidTropModel)
01068 {
01069 if(!valid) {
01070 if(!validWeather)
01071 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
01072 if(!validRxLatitude)
01073 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
01074 if(!validRxHeight)
01075 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
01076 if(!validDOY)
01077 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
01078 }
01079 double beta = NB_Interpolate(latitude,doy,ZB);
01080 double gm = 9.784*(1.0-2.66e-3*std::cos(2.0*latitude*DEG_TO_RAD)-2.8e-7*height);
01081
01082
01083
01084 double kd=1, base=std::log(1.0-beta*height/temp);
01085 if(interpolateWeather)
01086 kd = std::exp(base*NBg/(NBRd*beta));
01087
01088
01089 return ((1.0e-6*NBk1*NBRd/gm) * kd * press);
01090
01091 }
01092
01093
01094 double NBTropModel::wet_zenith_delay(void) const
01095 throw(TropModel::InvalidTropModel)
01096 {
01097 if(!valid) {
01098 if(!validWeather)
01099 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
01100 if(!validRxLatitude)
01101 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
01102 if(!validRxHeight)
01103 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
01104 if(!validDOY)
01105 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
01106 }
01107 double beta = NB_Interpolate(latitude,doy,ZB);
01108 double lam = NB_Interpolate(latitude,doy,ZL);
01109 double gm = 9.784*(1.0-2.66e-3*std::cos(2.0*latitude*DEG_TO_RAD)-2.8e-7*height);
01110
01111
01112
01113 double kw=1, base=std::log(1.0-beta*height/temp);
01114 if(interpolateWeather)
01115 kw = std::exp(base*(-1.0+(lam+1)*NBg/(NBRd*beta)));
01116
01117
01118 return ((1.0e-6*NBk3p*NBRd/(gm*(lam+1)-beta*NBRd)) * kw * humid/temp);
01119
01120 }
01121
01122
01123
01124
01125
01126 double NBTropModel::dry_mapping_function(double elevation) const
01127 throw(TropModel::InvalidTropModel)
01128 {
01129 if(!valid) {
01130 if(!validWeather)
01131 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
01132 if(!validRxLatitude)
01133 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
01134 if(!validRxHeight)
01135 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
01136 if(!validDOY)
01137 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
01138 }
01139 if(elevation < 0.0) return 0.0;
01140
01141 double a,b,c,se,map;
01142 se = std::sin(elevation*DEG_TO_RAD);
01143
01144 a = NB_Interpolate(latitude,doy,Mad);
01145 b = NB_Interpolate(latitude,doy,Mbd);
01146 c = NB_Interpolate(latitude,doy,Mcd);
01147 map = (1.0+a/(1.0+b/(1.0+c))) / (se+a/(se+b/(se+c)));
01148
01149 a = 2.53e-5;
01150 b = 5.49e-3;
01151 c = 1.14e-3;
01152 if(ABS(elevation)<=0.001) se=0.001;
01153 map += ((1.0/se)-(1.0+a/(1.0+b/(1.0+c)))/(se+a/(se+b/(se+c))))*height/1000.0;
01154
01155 return map;
01156
01157 }
01158
01159
01160
01161
01162
01163 double NBTropModel::wet_mapping_function(double elevation) const
01164 throw(TropModel::InvalidTropModel)
01165 {
01166 if(!valid) {
01167 if(!validWeather)
01168 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: weather"));
01169 if(!validRxLatitude)
01170 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Latitude"));
01171 if(!validRxHeight)
01172 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: Rx Height"));
01173 if(!validDOY)
01174 GPSTK_THROW(InvalidTropModel("Invalid NB trop model: day of year"));
01175 }
01176 if(elevation < 0.0) return 0.0;
01177
01178 double a,b,c,se;
01179 se = std::sin(elevation*DEG_TO_RAD);
01180 a = NB_Interpolate(latitude,doy,Maw);
01181 b = NB_Interpolate(latitude,doy,Mbw);
01182 c = NB_Interpolate(latitude,doy,Mcw);
01183
01184 return ( (1.0+a/(1.0+b/(1.0+c))) / (se+a/(se+b/(se+c))) );
01185
01186 }
01187
01188
01189
01190
01191
01192
01193 void NBTropModel::setWeather(const double& T,
01194 const double& P,
01195 const double& H)
01196 throw(InvalidParameter)
01197 {
01198 interpolateWeather=false;
01199 TropModel::setWeather(T,P,H);
01200
01201 double th=300./temp;
01202 humid = 2.409e9*H*th*th*th*th*std::exp(-22.64*th);
01203 validWeather = true;
01204 valid = validWeather && validRxHeight && validRxLatitude && validDOY;
01205
01206 }
01207
01208
01209
01210
01211 void NBTropModel::setWeather(const WxObservation& wx)
01212 throw(InvalidParameter)
01213 {
01214 interpolateWeather = false;
01215 try
01216 {
01217 TropModel::setWeather(wx);
01218
01219 double th=300./temp;
01220 humid = 2.409e9*humid*th*th*th*th*std::exp(-22.64*th);
01221 validWeather = true;
01222 valid = validWeather && validRxHeight && validRxLatitude && validDOY;
01223 }
01224 catch(InvalidParameter& e)
01225 {
01226 valid = validWeather = false;
01227 GPSTK_RETHROW(e);
01228 }
01229 }
01230
01231
01232
01233 void NBTropModel::setWeather()
01234 throw(TropModel::InvalidTropModel)
01235 {
01236 interpolateWeather = true;
01237 if(!validRxLatitude)
01238 {
01239 valid = validWeather = false;
01240 GPSTK_THROW(InvalidTropModel(
01241 "NBTropModel must have Rx latitude before interpolating weather"));
01242 }
01243 if(!validDOY)
01244 {
01245 valid = validWeather = false;
01246 GPSTK_THROW(InvalidTropModel(
01247 "NBTropModel must have day of year before interpolating weather"));
01248 }
01249 temp = NB_Interpolate(latitude,doy,ZT);
01250 press = NB_Interpolate(latitude,doy,ZP);
01251 humid = NB_Interpolate(latitude,doy,ZW);
01252 validWeather = true;
01253 valid = validWeather && validRxHeight && validRxLatitude && validDOY;
01254 }
01255
01256
01257
01258 void NBTropModel::setReceiverHeight(const double& ht)
01259 {
01260 height = ht;
01261 validRxHeight = true;
01262 valid = validWeather && validRxHeight && validRxLatitude && validDOY;
01263 if(!validWeather && validRxLatitude && validDOY)
01264 setWeather();
01265 }
01266
01267
01268
01269 void NBTropModel::setReceiverLatitude(const double& lat)
01270 {
01271 latitude = lat;
01272 validRxLatitude = true;
01273 valid = validWeather && validRxHeight && validRxLatitude && validDOY;
01274 if(!validWeather && validRxLatitude && validDOY)
01275 setWeather();
01276 }
01277
01278
01279
01280 void NBTropModel::setDayOfYear(const int& d)
01281 {
01282 doy = d;
01283 if(doy > 0 && doy < 367) validDOY=true; else validDOY = false;
01284 valid = validWeather && validRxHeight && validRxLatitude && validDOY;
01285 if(!validWeather && validRxLatitude && validDOY)
01286 setWeather();
01287 }
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320 static const double SaasWetA[5]=
01321 { 0.00058021897, 0.00056794847, 0.00058118019, 0.00059727542, 0.00061641693 };
01322 static const double SaasWetB[5]=
01323 { 0.0014275268, 0.0015138625, 0.0014572752, 0.0015007428, 0.0017599082 };
01324 static const double SaasWetC[5]=
01325 { 0.043472961, 0.046729510, 0.043908931, 0.044626982, 0.054736038 };
01326
01327
01328 static const double SaasDryA[5]=
01329 { 0.0012769934, 0.0012683230, 0.0012465397, 0.0012196049, 0.0012045996 };
01330 static const double SaasDryB[5]=
01331 { 0.0029153695, 0.0029152299, 0.0029288445, 0.0029022565, 0.0029024912 };
01332 static const double SaasDryC[5]=
01333 { 0.062610505, 0.062837393, 0.063721774, 0.063824265, 0.064258455 };
01334
01335 static const double SaasDryA1[5]=
01336 { 0.0, 0.000012709626, 0.000026523662, 0.000034000452, 0.000041202191 };
01337 static const double SaasDryB1[5]=
01338 { 0.0, 0.000021414979, 0.000030160779, 0.000072562722, 0.00011723375 };
01339 static const double SaasDryC1[5]=
01340 { 0.0, 0.000090128400, 0.000043497037, 0.00084795348, 0.0017037206 };
01341
01342
01343 SaasTropModel::SaasTropModel(void)
01344 {
01345 validWeather = false;
01346 validRxLatitude = false;
01347 validDOY = false;
01348 validRxHeight = false;
01349 }
01350
01351
01352
01353
01354
01355 SaasTropModel::SaasTropModel(const double& lat,
01356 const int& day)
01357 {
01358 validWeather = false;
01359 validRxHeight = false;
01360 SaasTropModel::setReceiverLatitude(lat);
01361 SaasTropModel::setDayOfYear(day);
01362 }
01363
01364
01365
01366
01367
01368 SaasTropModel::SaasTropModel(const double& lat,
01369 const int& day,
01370 const WxObservation& wx)
01371 throw(InvalidParameter)
01372 {
01373 validRxHeight = false;
01374 SaasTropModel::setReceiverLatitude(lat);
01375 SaasTropModel::setDayOfYear(day);
01376 SaasTropModel::setWeather(wx);
01377 }
01378
01379
01380
01381
01382
01383
01384
01385 SaasTropModel::SaasTropModel(const double& lat,
01386 const int& day,
01387 const double& T,
01388 const double& P,
01389 const double& H)
01390 throw(InvalidParameter)
01391 {
01392 validRxHeight = false;
01393 SaasTropModel::setReceiverLatitude(lat);
01394 SaasTropModel::setDayOfYear(day);
01395 SaasTropModel::setWeather(T,P,H);
01396 }
01397
01398
01399 double SaasTropModel::correction(double elevation) const
01400 throw(TropModel::InvalidTropModel)
01401 {
01402 if(!valid) {
01403 if(!validWeather) GPSTK_THROW(
01404 InvalidTropModel("Invalid Saastamoinen trop model: weather"));
01405 if(!validRxLatitude) GPSTK_THROW(
01406 InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
01407 if(!validRxHeight) GPSTK_THROW(
01408 InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
01409 if(!validDOY) GPSTK_THROW(
01410 InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
01411 GPSTK_THROW(
01412 InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
01413 }
01414
01415 if(elevation < 0.0) return 0.0;
01416
01417 double corr=0.0;
01418 try {
01419 corr = (dry_zenith_delay() * dry_mapping_function(elevation)
01420 + wet_zenith_delay() * wet_mapping_function(elevation));
01421 }
01422 catch(Exception& e) { GPSTK_RETHROW(e); }
01423
01424 return corr;
01425
01426 }
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437 double SaasTropModel::correction(const Position& RX,
01438 const Position& SV,
01439 const DayTime& tt)
01440 throw(TropModel::InvalidTropModel)
01441 {
01442 SaasTropModel::setReceiverHeight(RX.getHeight());
01443 SaasTropModel::setReceiverLatitude(RX.getGeodeticLatitude());
01444 SaasTropModel::setDayOfYear(int(tt.DOYday()));
01445
01446 if(!valid) {
01447 if(!validWeather) GPSTK_THROW(
01448 InvalidTropModel("Invalid Saastamoinen trop model: weather"));
01449 if(!validRxLatitude) GPSTK_THROW(
01450 InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
01451 if(!validRxHeight) GPSTK_THROW(
01452 InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
01453 if(!validDOY) GPSTK_THROW(
01454 InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
01455 valid = true;
01456 }
01457
01458 double corr=0.0;
01459 try {
01460 corr = SaasTropModel::correction(RX.elevation(SV));
01461 }
01462 catch(Exception& e) { GPSTK_RETHROW(e); }
01463
01464 return corr;
01465
01466 }
01467
01468 double SaasTropModel::correction(const Xvt& RX,
01469 const Xvt& SV,
01470 const DayTime& tt)
01471 throw(TropModel::InvalidTropModel)
01472 {
01473 Position R(RX),S(SV);
01474 return SaasTropModel::correction(R,S,tt);
01475 }
01476
01477
01478 double SaasTropModel::dry_zenith_delay(void) const
01479 throw(TropModel::InvalidTropModel)
01480 {
01481 if(!valid) {
01482 if(!validWeather) GPSTK_THROW(
01483 InvalidTropModel("Invalid Saastamoinen trop model: weather"));
01484 if(!validRxLatitude) GPSTK_THROW(
01485 InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
01486 if(!validRxHeight) GPSTK_THROW(
01487 InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
01488 if(!validDOY) GPSTK_THROW(
01489 InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
01490 GPSTK_THROW(
01491 InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
01492 }
01493
01494
01495
01496
01497
01498 double delay = 0.0022768 * press
01499 / (1 - 0.00266 * ::cos(2*latitude*DEG_TO_RAD) - 0.00028 * height/1000.);
01500
01501 return delay;
01502
01503 }
01504
01505
01506 double SaasTropModel::wet_zenith_delay(void) const
01507 throw(TropModel::InvalidTropModel)
01508 {
01509 if(!valid) {
01510 if(!validWeather) GPSTK_THROW(
01511 InvalidTropModel("Invalid Saastamoinen trop model: weather"));
01512 if(!validRxLatitude) GPSTK_THROW(
01513 InvalidTropModel("Invalid Saastamoinen trop model: Rx Latitude"));
01514 if(!validRxHeight) GPSTK_THROW(
01515 InvalidTropModel("Invalid Saastamoinen trop model: Rx Height"));
01516 if(!validDOY) GPSTK_THROW(
01517 InvalidTropModel("Invalid Saastamoinen trop model: day of year"));
01518 GPSTK_THROW(
01519 InvalidTropModel("Valid flag corrupted in Saastamoinen trop model"));
01520 }
01521
01522
01523 double delay = 0.0022768 * humid * 1255/((temp+CELSIUS_TO_KELVIN) + 0.05)
01524 / (1 - 0.00266 * ::cos(2*latitude*DEG_TO_RAD) - 0.00028 * height/1000.);
01525
01526 return delay;
01527
01528 }
01529
01530
01531
01532 double SaasTropModel::dry_mapping_function(