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