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