src/Environment/environment.cxx

   1 // environment.cxx -- routines to model the natural environment
   2 //
   3 // Written by David Megginson, started February 2002.
   4 //
   5 // Copyright (C) 2002  David Megginson - david@megginson.com
   6 //
   7 // This program is free software; you can redistribute it and/or
   8 // modify it under the terms of the GNU General Public License as
   9 // published by the Free Software Foundation; either version 2 of the
  10 // License, or (at your option) any later version.
  11 //
  12 // This program is distributed in the hope that it will be useful, but
  13 // WITHOUT ANY WARRANTY; without even the implied warranty of
  14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15 // General Public License for more details.
  16 //
  17 // You should have received a copy of the GNU General Public License
  18 // along with this program; if not, write to the Free Software
  19 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  20 //
  21
  22 #ifdef HAVE_CONFIG_H
  23 #  include <config.h>
  24 #endif
  25
  26 #include <math.h>
  27
  28 #include <boost/tuple/tuple.hpp>
  29
  30 #include <simgear/props/props.hxx>
  31 #include <simgear/math/SGMath.hxx>
  32
  33 #include <Main/fg_props.hxx>
  34
  35 #include "environment.hxx"
  36 #include "atmosphere.hxx"
  37
  38 \f
  39 ////////////////////////////////////////////////////////////////////////
  40 // Atmosphere model.
  41 ////////////////////////////////////////////////////////////////////////
  42
  43 #ifdef USING_TABLES
  44
  45 // Calculated based on the ISA standard day, as found at e.g.
  46 // http://www.av8n.com/physics/altimetry.htm
  47
  48 // Each line of data has 3 elements:
  49 //   Elevation (ft),
  50 //   temperature factor (dimensionless ratio of absolute temp),
  51 //   pressure factor (dimensionless ratio)
  52 static double atmosphere_data[][3] = {
  53  {  -3000.00,   1.021,  1.1133 },
  54  {      0.00,   1.000,  1.0000 },
  55  {   2952.76,   0.980,  0.8978 },
  56  {   5905.51,   0.959,  0.8042 },
  57  {   8858.27,   0.939,  0.7187 },
  58  {  11811.02,   0.919,  0.6407 },
  59  {  14763.78,   0.898,  0.5697 },
  60  {  17716.54,   0.878,  0.5052 },
  61  {  20669.29,   0.858,  0.4468 },
  62  {  23622.05,   0.838,  0.3940 },
  63  {  26574.80,   0.817,  0.3463 },
  64  {  29527.56,   0.797,  0.3034 },
  65  {  32480.31,   0.777,  0.2649 },
  66  {  35433.07,   0.756,  0.2305 },
  67  {  38385.83,   0.752,  0.2000 },
  68  {  41338.58,   0.752,  0.1736 },
  69  {  44291.34,   0.752,  0.1506 },
  70  {  47244.09,   0.752,  0.1307 },
  71  {  50196.85,   0.752,  0.1134 },
  72  {  53149.61,   0.752,  0.0984 },
  73  {  56102.36,   0.752,  0.0854 },
  74  {  59055.12,   0.752,  0.0741 },
  75  {  62007.87,   0.752,  0.0643 },
  76  {  65000.00,   0.752,  0.0557 },
  77  {  68000.00,   0.754,  0.0482 },
  78  {  71000.00,   0.758,  0.0418 },
  79  {  74000.00,   0.761,  0.0362 },
  80  {  77000.00,   0.764,  0.0314 },
  81  {  80000.00,   0.767,  0.0273 },
  82  {  83000.00,   0.770,  0.0237 },
  83  {  86000.00,   0.773,  0.0206 },
  84  {  89000.00,   0.777,  0.0179 },
  85  {  92000.00,   0.780,  0.0156 },
  86  {  95000.00,   0.783,  0.0135 },
  87  {  98000.00,   0.786,  0.0118 },
  88  { 101000.00,   0.789,  0.0103 },
  89  { -1, -1, -1 }
  90 };
  91
  92 static SGInterpTable * _temperature_degc_table = 0;
  93 static SGInterpTable * _pressure_inhg_table = 0;
  94
  95 static void
  96 _setup_tables ()
  97 {
  98   if (_temperature_degc_table != 0)
  99       return;
 100
 101   _temperature_degc_table = new SGInterpTable;
 102   _pressure_inhg_table = new SGInterpTable;
 103
 104   for (int i = 0; atmosphere_data[i][0] != -1; i++) {
 105     _temperature_degc_table->addEntry(atmosphere_data[i][0],
 106                                       atmosphere_data[i][1]);
 107     _pressure_inhg_table->addEntry(atmosphere_data[i][0],
 108                                    atmosphere_data[i][2]);
 109   }
 110 }
 111 #endif
 112
 113 \f
 114 ////////////////////////////////////////////////////////////////////////
 115 // Implementation of FGEnvironment.
 116 ////////////////////////////////////////////////////////////////////////
 117
 118 void FGEnvironment::_init()
 119 {
 120     live_update = false;
 121     elevation_ft = 0;
 122     visibility_m = 32000;
 123     temperature_sea_level_degc = 15;
 124     temperature_degc = 15;
 125     dewpoint_sea_level_degc = 5; // guess
 126     dewpoint_degc = 5;
 127     pressure_sea_level_inhg = 29.92;
 128     pressure_inhg = 29.92;
 129     density_slugft3 = 0;
 130     turbulence_magnitude_norm = 0;
 131     turbulence_rate_hz = 1;
 132     wind_from_heading_deg = 0;
 133     wind_speed_kt = 0;
 134     wind_from_north_fps = 0;
 135     wind_from_east_fps = 0;
 136     wind_from_down_fps = 0;
 137     altitude_half_to_sun_m = 1000;
 138     altitude_tropo_top_m = 10000;
 139 #ifdef USING_TABLES
 140     _setup_tables();
 141 #endif
 142     _recalc_density();
 143     _recalc_relative_humidity();
 144     live_update = true;
 145 }
 146
 147 FGEnvironment::FGEnvironment()
 148 {
 149     _init();
 150 }
 151
 152 FGEnvironment::FGEnvironment (const FGEnvironment &env)
 153 {
 154     _init();
 155     copy(env);
 156 }
 157
 158 FGEnvironment::~FGEnvironment()
 159 {
 160     Untie();
 161 }
 162
 163 FGEnvironment & FGEnvironment::operator = ( const FGEnvironment & other )
 164 {
 165     copy( other );
 166     return *this;
 167 }
 168
 169 void
 170 FGEnvironment::copy (const FGEnvironment &env)
 171 {
 172     elevation_ft = env.elevation_ft;
 173     visibility_m = env.visibility_m;
 174     temperature_sea_level_degc = env.temperature_sea_level_degc;
 175     temperature_degc = env.temperature_degc;
 176     dewpoint_sea_level_degc = env.dewpoint_sea_level_degc;
 177     dewpoint_degc = env.dewpoint_degc;
 178     pressure_sea_level_inhg = env.pressure_sea_level_inhg;
 179     wind_from_heading_deg = env.wind_from_heading_deg;
 180     wind_speed_kt = env.wind_speed_kt;
 181     wind_from_north_fps = env.wind_from_north_fps;
 182     wind_from_east_fps = env.wind_from_east_fps;
 183     wind_from_down_fps = env.wind_from_down_fps;
 184     turbulence_magnitude_norm = env.turbulence_magnitude_norm;
 185     turbulence_rate_hz = env.turbulence_rate_hz;
 186     pressure_inhg = env.pressure_inhg;
 187     density_slugft3 = env.density_slugft3;
 188     density_tropo_avg_kgm3 = env.density_tropo_avg_kgm3;
 189     relative_humidity = env.relative_humidity;
 190     altitude_half_to_sun_m = env.altitude_half_to_sun_m;
 191     altitude_tropo_top_m = env.altitude_tropo_top_m;
 192     live_update = env.live_update;
 193 }
 194
 195 static inline bool
 196 maybe_copy_value (FGEnvironment * env, const SGPropertyNode * node,
 197                   const char * name, void (FGEnvironment::*setter)(double))
 198 {
 199     const SGPropertyNode * child = node->getNode(name);
 200                                 // fragile: depends on not being typed
 201                                 // as a number
 202     if (child != 0 && child->hasValue() &&
 203         child->getStringValue()[0] != '\0') {
 204         (env->*setter)(child->getDoubleValue());
 205         return true;
 206     } else {
 207         return false;
 208     }
 209 }
 210
 211 void
 212 FGEnvironment::read (const SGPropertyNode * node)
 213 {
 214     bool live_update = set_live_update( false );
 215     maybe_copy_value(this, node, "visibility-m",
 216                      &FGEnvironment::set_visibility_m);
 217
 218     maybe_copy_value(this, node, "elevation-ft",
 219                      &FGEnvironment::set_elevation_ft);
 220
 221     if (!maybe_copy_value(this, node, "temperature-sea-level-degc",
 222                           &FGEnvironment::set_temperature_sea_level_degc)) {
 223         if( maybe_copy_value(this, node, "temperature-degc",
 224                          &FGEnvironment::set_temperature_degc)) {
 225             _recalc_sl_temperature();
 226         }
 227     }
 228
 229     if (!maybe_copy_value(this, node, "dewpoint-sea-level-degc",
 230                           &FGEnvironment::set_dewpoint_sea_level_degc)) {
 231         if( maybe_copy_value(this, node, "dewpoint-degc",
 232                          &FGEnvironment::set_dewpoint_degc)) {
 233             _recalc_sl_dewpoint();
 234         }
 235     }
 236
 237     if (!maybe_copy_value(this, node, "pressure-sea-level-inhg",
 238                           &FGEnvironment::set_pressure_sea_level_inhg)) {
 239         if( maybe_copy_value(this, node, "pressure-inhg",
 240                          &FGEnvironment::set_pressure_inhg)) {
 241             _recalc_sl_pressure();
 242         }
 243    }
 244
 245     maybe_copy_value(this, node, "wind-from-heading-deg",
 246                      &FGEnvironment::set_wind_from_heading_deg);
 247
 248     maybe_copy_value(this, node, "wind-speed-kt",
 249                      &FGEnvironment::set_wind_speed_kt);
 250
 251     maybe_copy_value(this, node, "turbulence/magnitude-norm",
 252                      &FGEnvironment::set_turbulence_magnitude_norm);
 253
 254     maybe_copy_value(this, node, "turbulence/rate-hz",
 255                      &FGEnvironment::set_turbulence_rate_hz);
 256
 257     // calculate derived properties here to avoid duplicate expensive computations
 258     _recalc_ne();
 259     _recalc_alt_pt();
 260     _recalc_alt_dewpoint();
 261     _recalc_density();
 262     _recalc_relative_humidity();
 263
 264     set_live_update(live_update);
 265 }
 266
 267 void FGEnvironment::Tie( SGPropertyNode_ptr base, bool archivable )
 268 {
 269   _tiedProperties.setRoot( base );
 270
 271   _tiedProperties.Tie( "visibility-m", this,
 272       &FGEnvironment::get_visibility_m,
 273       &FGEnvironment::set_visibility_m);
 274
 275   _tiedProperties.Tie("temperature-sea-level-degc", this,
 276       &FGEnvironment::get_temperature_sea_level_degc,
 277       &FGEnvironment::set_temperature_sea_level_degc);
 278
 279   _tiedProperties.Tie("temperature-degc", this,
 280       &FGEnvironment::get_temperature_degc,
 281       &FGEnvironment::set_temperature_degc);
 282
 283   _tiedProperties.Tie("dewpoint-sea-level-degc", this,
 284       &FGEnvironment::get_dewpoint_sea_level_degc,
 285       &FGEnvironment::set_dewpoint_sea_level_degc);
 286
 287   _tiedProperties.Tie("dewpoint-degc", this,
 288       &FGEnvironment::get_dewpoint_degc,
 289       &FGEnvironment::set_dewpoint_degc);
 290
 291   _tiedProperties.Tie("pressure-sea-level-inhg", this,
 292       &FGEnvironment::get_pressure_sea_level_inhg,
 293       &FGEnvironment::set_pressure_sea_level_inhg);
 294
 295   _tiedProperties.Tie("pressure-inhg", this,
 296       &FGEnvironment::get_pressure_inhg,
 297       &FGEnvironment::set_pressure_inhg);
 298
 299   _tiedProperties.Tie("atmosphere/altitude-half-to-sun", this,
 300       &FGEnvironment::get_altitude_half_to_sun_m,
 301       &FGEnvironment::set_altitude_half_to_sun_m);
 302
 303   _tiedProperties.Tie("atmosphere/altitude-troposphere-top", this,
 304       &FGEnvironment::get_altitude_tropo_top_m,
 305       &FGEnvironment::set_altitude_tropo_top_m);
 306
 307   _tiedProperties.Tie("wind-from-heading-deg", this,
 308       &FGEnvironment::get_wind_from_heading_deg,
 309       &FGEnvironment::set_wind_from_heading_deg);
 310
 311   _tiedProperties.Tie("wind-speed-kt", this,
 312       &FGEnvironment::get_wind_speed_kt,
 313       &FGEnvironment::set_wind_speed_kt);
 314
 315   _tiedProperties.Tie("wind-from-north-fps", this,
 316       &FGEnvironment::get_wind_from_north_fps,
 317       &FGEnvironment::set_wind_from_north_fps);
 318
 319   _tiedProperties.Tie("wind-from-east-fps", this,
 320       &FGEnvironment::get_wind_from_east_fps,
 321       &FGEnvironment::set_wind_from_east_fps);
 322
 323   _tiedProperties.Tie("wind-from-down-fps", this,
 324       &FGEnvironment::get_wind_from_down_fps,
 325       &FGEnvironment::set_wind_from_down_fps);
 326
 327   _tiedProperties.Tie("turbulence/magnitude-norm", this,
 328       &FGEnvironment::get_turbulence_magnitude_norm,
 329       &FGEnvironment::set_turbulence_magnitude_norm);
 330
 331   _tiedProperties.Tie("turbulence/rate-hz", this,
 332       &FGEnvironment::get_turbulence_rate_hz,
 333       &FGEnvironment::set_turbulence_rate_hz);
 334
 335   _tiedProperties.setAttribute( SGPropertyNode::ARCHIVE, archivable );
 336
 337   _tiedProperties.Tie("temperature-degf", this,
 338       &FGEnvironment::get_temperature_degf);
 339
 340   _tiedProperties.Tie("density-slugft3", this,
 341       &FGEnvironment::get_density_slugft3); // read-only
 342
 343   _tiedProperties.Tie("relative-humidity", this,
 344       &FGEnvironment::get_relative_humidity); //ro
 345
 346   _tiedProperties.Tie("atmosphere/density-tropo-avg", this,
 347       &FGEnvironment::get_density_tropo_avg_kgm3); //ro
 348 }
 349
 350 void FGEnvironment::Untie()
 351 {
 352     _tiedProperties.Untie();
 353 }
 354
 355 double
 356 FGEnvironment::get_visibility_m () const
 357 {
 358   return visibility_m;
 359 }
 360
 361 double
 362 FGEnvironment::get_temperature_sea_level_degc () const
 363 {
 364   return temperature_sea_level_degc;
 365 }
 366
 367 double
 368 FGEnvironment::get_temperature_degc () const
 369 {
 370   return temperature_degc;
 371 }
 372
 373 double
 374 FGEnvironment::get_temperature_degf () const
 375 {
 376   return (temperature_degc * 9.0 / 5.0) + 32.0;
 377 }
 378
 379 double
 380 FGEnvironment::get_dewpoint_sea_level_degc () const
 381 {
 382   return dewpoint_sea_level_degc;
 383 }
 384
 385 double
 386 FGEnvironment::get_dewpoint_degc () const
 387 {
 388   return dewpoint_degc;
 389 }
 390
 391 double
 392 FGEnvironment::get_pressure_sea_level_inhg () const
 393 {
 394   return pressure_sea_level_inhg;
 395 }
 396
 397 double
 398 FGEnvironment::get_pressure_inhg () const
 399 {
 400   return pressure_inhg;
 401 }
 402
 403 double
 404 FGEnvironment::get_density_slugft3 () const
 405 {
 406   return density_slugft3;
 407 }
 408
 409 double
 410 FGEnvironment::get_relative_humidity () const
 411 {
 412   return relative_humidity;
 413 }
 414
 415 double
 416 FGEnvironment::get_density_tropo_avg_kgm3 () const
 417 {
 418   return density_tropo_avg_kgm3;
 419 }
 420
 421 double
 422 FGEnvironment::get_altitude_half_to_sun_m () const
 423 {
 424   return altitude_half_to_sun_m;
 425 }
 426
 427 double
 428 FGEnvironment::get_altitude_tropo_top_m () const
 429 {
 430   return altitude_tropo_top_m;
 431 }
 432
 433 double
 434 FGEnvironment::get_wind_from_heading_deg () const
 435 {
 436   return wind_from_heading_deg;
 437 }
 438
 439 double
 440 FGEnvironment::get_wind_speed_kt () const
 441 {
 442   return wind_speed_kt;
 443 }
 444
 445 double
 446 FGEnvironment::get_wind_from_north_fps () const
 447 {
 448   return wind_from_north_fps;
 449 }
 450
 451 double
 452 FGEnvironment::get_wind_from_east_fps () const
 453 {
 454   return wind_from_east_fps;
 455 }
 456
 457 double
 458 FGEnvironment::get_wind_from_down_fps () const
 459 {
 460   return wind_from_down_fps;
 461 }
 462
 463 double
 464 FGEnvironment::get_turbulence_magnitude_norm () const
 465 {
 466   return turbulence_magnitude_norm;
 467 }
 468
 469 double
 470 FGEnvironment::get_turbulence_rate_hz () const
 471 {
 472   return turbulence_rate_hz;
 473 }
 474
 475 double
 476 FGEnvironment::get_elevation_ft () const
 477 {
 478   return elevation_ft;
 479 }
 480
 481 void
 482 FGEnvironment::set_visibility_m (double v)
 483 {
 484   visibility_m = v;
 485 }
 486
 487 void
 488 FGEnvironment::set_temperature_sea_level_degc (double t)
 489 {
 490   temperature_sea_level_degc = t;
 491   if (dewpoint_sea_level_degc > t)
 492       dewpoint_sea_level_degc = t;
 493   if( live_update ) {
 494     _recalc_alt_pt();
 495     _recalc_density();
 496   }
 497 }
 498
 499 void
 500 FGEnvironment::set_temperature_degc (double t)
 501 {
 502   temperature_degc = t;
 503   if( live_update ) {
 504     _recalc_sl_temperature();
 505     _recalc_sl_pressure();
 506     _recalc_alt_pt();
 507     _recalc_density();
 508     _recalc_relative_humidity();
 509   }
 510 }
 511
 512 void
 513 FGEnvironment::set_dewpoint_sea_level_degc (double t)
 514 {
 515   dewpoint_sea_level_degc = t;
 516   if (temperature_sea_level_degc < t)
 517       temperature_sea_level_degc = t;
 518   if( live_update ) {
 519     _recalc_alt_dewpoint();
 520     _recalc_density();
 521   }
 522 }
 523
 524 void
 525 FGEnvironment::set_dewpoint_degc (double t)
 526 {
 527   dewpoint_degc = t;
 528   if( live_update ) {
 529     _recalc_sl_dewpoint();
 530     _recalc_density();
 531     _recalc_relative_humidity();
 532   }
 533 }
 534
 535 void
 536 FGEnvironment::set_pressure_sea_level_inhg (double p)
 537 {
 538   pressure_sea_level_inhg = p;
 539   if( live_update ) {
 540     _recalc_alt_pt();
 541     _recalc_density();
 542   }
 543 }
 544
 545 void
 546 FGEnvironment::set_pressure_inhg (double p)
 547 {
 548   pressure_inhg = p;
 549   if( live_update ) {
 550     _recalc_sl_pressure();
 551     _recalc_density();
 552   }
 553 }
 554
 555 void
 556 FGEnvironment::set_wind_from_heading_deg (double h)
 557 {
 558   wind_from_heading_deg = h;
 559   if( live_update ) {
 560     _recalc_ne();
 561   }
 562 }
 563
 564 void
 565 FGEnvironment::set_wind_speed_kt (double s)
 566 {
 567   wind_speed_kt = s;
 568   if( live_update ) {
 569     _recalc_ne();
 570   }
 571 }
 572
 573 void
 574 FGEnvironment::set_wind_from_north_fps (double n)
 575 {
 576   wind_from_north_fps = n;
 577   if( live_update ) {
 578     _recalc_hdgspd();
 579   }
 580 }
 581
 582 void
 583 FGEnvironment::set_wind_from_east_fps (double e)
 584 {
 585   wind_from_east_fps = e;
 586   if( live_update ) {
 587     _recalc_hdgspd();
 588   }
 589 }
 590
 591 void
 592 FGEnvironment::set_wind_from_down_fps (double d)
 593 {
 594   wind_from_down_fps = d;
 595   if( live_update ) {
 596     _recalc_hdgspd();
 597   }
 598 }
 599
 600 void
 601 FGEnvironment::set_turbulence_magnitude_norm (double t)
 602 {
 603   turbulence_magnitude_norm = t;
 604 }
 605
 606 void
 607 FGEnvironment::set_turbulence_rate_hz (double r)
 608 {
 609   turbulence_rate_hz = r;
 610 }
 611
 612 void
 613 FGEnvironment::set_elevation_ft (double e)
 614 {
 615   elevation_ft = e;
 616   if( live_update ) {
 617     _recalc_alt_pt();
 618     _recalc_alt_dewpoint();
 619     _recalc_density();
 620     _recalc_relative_humidity();
 621   }
 622 }
 623
 624 void
 625 FGEnvironment::set_altitude_half_to_sun_m (double alt)
 626 {
 627   altitude_half_to_sun_m = alt;
 628   if( live_update ) {
 629     _recalc_density_tropo_avg_kgm3();
 630   }
 631 }
 632
 633 void
 634 FGEnvironment::set_altitude_tropo_top_m (double alt)
 635 {
 636   altitude_tropo_top_m = alt;
 637   if( live_update ) {
 638     _recalc_density_tropo_avg_kgm3();
 639   }
 640 }
 641
 642
 643 void
 644 FGEnvironment::_recalc_hdgspd ()
 645 {
 646   wind_from_heading_deg =
 647     atan2(wind_from_east_fps, wind_from_north_fps) * SGD_RADIANS_TO_DEGREES;
 648
 649   if( wind_from_heading_deg < 0 )
 650     wind_from_heading_deg += 360.0;
 651
 652   wind_speed_kt = sqrt(wind_from_north_fps * wind_from_north_fps +
 653                        wind_from_east_fps * wind_from_east_fps)
 654                   * SG_METER_TO_NM * SG_FEET_TO_METER * 3600;
 655 }
 656
 657 void
 658 FGEnvironment::_recalc_ne ()
 659 {
 660   double speed_fps =
 661     wind_speed_kt * SG_NM_TO_METER * SG_METER_TO_FEET * (1.0/3600);
 662
 663   wind_from_north_fps = speed_fps *
 664     cos(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
 665   wind_from_east_fps = speed_fps *
 666     sin(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
 667 }
 668
 669 // Intended to help with the interpretation of METAR data,
 670 // not for random in-flight outside-air temperatures.
 671 void
 672 FGEnvironment::_recalc_sl_temperature ()
 673 {
 674
 675 #if 0
 676   {
 677     SG_LOG(SG_ENVIRONMENT, SG_DEBUG, "recalc_sl_temperature: using "
 678       << temperature_degc << " @ " << elevation_ft << " :: " << this);
 679   }
 680 #endif
 681
 682   if (elevation_ft * atmodel::foot >= ISA_def[1].height) {
 683     SG_LOG(SG_ENVIRONMENT, SG_ALERT, "recalc_sl_temperature: "
 684         << "valid only in troposphere, not " << elevation_ft);
 685     return;
 686   }
 687
 688 // Clamp: temperature of the stratosphere, in degrees C:
 689   double t_strato = ISA_def[1].temp - atmodel::freezing;
 690   if (temperature_degc < t_strato) temperature_sea_level_degc = t_strato;
 691   else temperature_sea_level_degc =
 692       temperature_degc + elevation_ft * atmodel::foot * ISA_def[0].lapse;
 693
 694 // Alternative implemenation:
 695 //  else temperature_sea_level_inhg = T_layer(0., elevation_ft * foot,
 696 //      pressure_inhg * inHg, temperature_degc + freezing, ISA_def[0].lapse) - freezing;
 697 }
 698
 699 void
 700 FGEnvironment::_recalc_sl_dewpoint ()
 701 {
 702                                 // 0.2degC/1000ft
 703                                 // FIXME: this will work only for low
 704                                 // elevations
 705   dewpoint_sea_level_degc = dewpoint_degc + (elevation_ft * .0002);
 706   if (dewpoint_sea_level_degc > temperature_sea_level_degc)
 707     dewpoint_sea_level_degc = temperature_sea_level_degc;
 708 }
 709
 710 void
 711 FGEnvironment::_recalc_alt_dewpoint ()
 712 {
 713                                 // 0.2degC/1000ft
 714                                 // FIXME: this will work only for low
 715                                 // elevations
 716   dewpoint_degc = dewpoint_sea_level_degc + (elevation_ft * .0002);
 717   if (dewpoint_degc > temperature_degc)
 718     dewpoint_degc = temperature_degc;
 719 }
 720
 721 void
 722 FGEnvironment::_recalc_sl_pressure ()
 723 {
 724   using namespace atmodel;
 725 #if 0
 726   {
 727     SG_LOG(SG_ENVIRONMENT, SG_ALERT, "recalc_sl_pressure: using "
 728       << pressure_inhg << " and "
 729       << temperature_degc << " @ " << elevation_ft << " :: " << this);
 730   }
 731 #endif
 732   pressure_sea_level_inhg = P_layer(0., elevation_ft * foot,
 733       pressure_inhg * inHg, temperature_degc + freezing, ISA_def[0].lapse) / inHg;
 734 }
 735
 736 // This gets called at frame rate, to account for the aircraft's
 737 // changing altitude.
 738 // Called by set_elevation_ft() which is called by FGEnvironmentMgr::update
 739
 740 void
 741 FGEnvironment::_recalc_alt_pt ()
 742 {
 743   using namespace atmodel;
 744 #if 0
 745   {
 746     static int count(0);
 747     if (++count % 1000 == 0) {
 748       SG_LOG(SG_ENVIRONMENT, SG_ALERT,
 749            "recalc_alt_pt for: " << elevation_ft
 750         << "  using "  << pressure_sea_level_inhg
 751         << "  and "  << temperature_sea_level_degc
 752         << " :: " << this
 753         << "  # " << count);
 754     }
 755   }
 756 #endif
 757   double press = pressure_inhg * inHg;
 758   double temp = temperature_degc + freezing;
 759   boost::tie(press, temp) = PT_vs_hpt(elevation_ft * foot,
 760         pressure_sea_level_inhg * inHg, temperature_sea_level_degc + freezing);
 761   temperature_degc = temp - freezing;
 762   pressure_inhg = press / inHg;
 763 }
 764
 765 void
 766 FGEnvironment::_recalc_density ()
 767 {
 768   double pressure_psf = pressure_inhg * 70.7487;
 769
 770   // adjust for humidity
 771   // calculations taken from USA Today (oops!) at
 772   // http://www.usatoday.com/weather/basics/density-calculations.htm
 773   double temperature_degk = temperature_degc + 273.15;
 774   double pressure_mb = pressure_inhg * 33.86;
 775   double vapor_pressure_mb =
 776     6.11 * pow(10.0, 7.5 * dewpoint_degc / (237.7 + dewpoint_degc));
 777   double virtual_temperature_degk = temperature_degk / (1 - (vapor_pressure_mb / pressure_mb) * (1.0 - 0.622));
 778   double virtual_temperature_degr = virtual_temperature_degk * 1.8;
 779
 780   density_slugft3 = pressure_psf / (virtual_temperature_degr * 1718);
 781   _recalc_density_tropo_avg_kgm3();
 782 }
 783
 784 // This is used to calculate the average density on the path
 785 // of sunlight to the observer for calculating sun-color
 786 void
 787 FGEnvironment::_recalc_density_tropo_avg_kgm3 ()
 788 {
 789   double pressure_mb = pressure_inhg * 33.86;
 790   double vaporpressure = 6.11 * pow(10.0, ((7.5 * dewpoint_degc) / (237.7 + dewpoint_degc)));
 791
 792   double virtual_temp = (temperature_degc + 273.15) / (1 - 0.379 * (vaporpressure/pressure_mb));
 793
 794   double density_half = (100 * pressure_mb * exp(-altitude_half_to_sun_m / 8000))
 795       / (287.05 * virtual_temp);
 796   double density_tropo = (100 * pressure_mb * exp((-1 * altitude_tropo_top_m) / 8000))
 797       / ( 287.05 * virtual_temp);
 798
 799   density_tropo_avg_kgm3 = ((density_slugft3 * 515.379) + density_half + density_tropo) / 3;
 800 }
 801
 802 void
 803 FGEnvironment::_recalc_relative_humidity ()
 804 {
 805 /*
 806   double vaporpressure = 6.11 * pow(10.0, ((7.5 * dewpoint_degc) / ( 237.7 + dewpoint_degc)));
 807   double sat_vaporpressure = 6.11 * pow(10.0, ((7.5 * temperature_degc)
 808       / ( 237.7 + temperature_degc)) );
 809   relative_humidity = 100 * vaporpressure / sat_vaporpressure ;
 810
 811   with a little algebra, this gets the same result and spares two multiplications and one pow()
 812 */
 813   double a = (7.5 * dewpoint_degc)    / ( 237.7 + dewpoint_degc);
 814   double b = (7.5 * temperature_degc) / ( 237.7 + temperature_degc);
 815   relative_humidity = 100 * pow(10.0,a-b);
 816 }
 817
 818 bool
 819 FGEnvironment::set_live_update( bool _live_update )
 820 {
 821   bool b = live_update;
 822   live_update = _live_update;
 823   return b;
 824 }
 825
 826
 827 ////////////////////////////////////////////////////////////////////////
 828 // Functions.
 829 ////////////////////////////////////////////////////////////////////////
 830
 831 static inline double
 832 do_interp (double a, double b, double fraction)
 833 {
 834     double retval = (a + ((b - a) * fraction));
 835     return retval;
 836 }
 837
 838 static inline double
 839 do_interp_deg (double a, double b, double fraction)
 840 {
 841     a = fmod(a, 360);
 842     b = fmod(b, 360);
 843     if (fabs(b-a) > 180) {
 844         if (a < b)
 845             a += 360;
 846         else
 847             b += 360;
 848     }
 849     return fmod(do_interp(a, b, fraction), 360);
 850 }
 851
 852 FGEnvironment &
 853 FGEnvironment::interpolate( const FGEnvironment & env2,
 854              double fraction, FGEnvironment * result) const
 855 {
 856     // don't calculate each internal property every time we set a single value
 857     // we trigger that at the end of the interpolation process
 858     bool live_update = result->set_live_update( false );
 859
 860     result->set_visibility_m
 861         (do_interp(get_visibility_m(),
 862                    env2.get_visibility_m(),
 863                    fraction));
 864
 865     result->set_temperature_sea_level_degc
 866         (do_interp(get_temperature_sea_level_degc(),
 867                    env2.get_temperature_sea_level_degc(),
 868                    fraction));
 869
 870     result->set_dewpoint_sea_level_degc
 871         (do_interp(get_dewpoint_sea_level_degc(),
 872                    env2.get_dewpoint_sea_level_degc(),
 873                    fraction));
 874
 875     result->set_pressure_sea_level_inhg
 876         (do_interp(get_pressure_sea_level_inhg(),
 877                    env2.get_pressure_sea_level_inhg(),
 878                    fraction));
 879
 880     result->set_wind_from_heading_deg
 881         (do_interp_deg(get_wind_from_heading_deg(),
 882                        env2.get_wind_from_heading_deg(),
 883                        fraction));
 884
 885     result->set_wind_speed_kt
 886         (do_interp(get_wind_speed_kt(),
 887                    env2.get_wind_speed_kt(),
 888                    fraction));
 889
 890     result->set_elevation_ft
 891         (do_interp(get_elevation_ft(),
 892                    env2.get_elevation_ft(),
 893                    fraction));
 894
 895     result->set_turbulence_magnitude_norm
 896         (do_interp(get_turbulence_magnitude_norm(),
 897                    env2.get_turbulence_magnitude_norm(),
 898                    fraction));
 899
 900     result->set_turbulence_rate_hz
 901         (do_interp(get_turbulence_rate_hz(),
 902                    env2.get_turbulence_rate_hz(),
 903                    fraction));
 904
 905     // calculate derived properties here to avoid duplicate expensive computations
 906     result->_recalc_ne();
 907     result->_recalc_alt_pt();
 908     result->_recalc_alt_dewpoint();
 909     result->_recalc_density();
 910     result->_recalc_relative_humidity();
 911
 912     result->set_live_update(live_update);
 913
 914     return *result;
 915 }
 916
 917 // end of environment.cxx