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     elevation_ft = 0;
 121     visibility_m = 32000;
 122     temperature_sea_level_degc = 15;
 123     temperature_degc = 15;
 124     dewpoint_sea_level_degc = 5; // guess
 125     dewpoint_degc = 5;
 126     pressure_sea_level_inhg = 29.92;
 127     pressure_inhg = 29.92;
 128     turbulence_magnitude_norm = 0;
 129     turbulence_rate_hz = 1;
 130     wind_from_heading_deg = 0;
 131     wind_speed_kt = 0;
 132     wind_from_north_fps = 0;
 133     wind_from_east_fps = 0;
 134     wind_from_down_fps = 0;
 135     altitude_half_to_sun_m = 1000;
 136     altitude_tropo_top_m = 10000;
 137 #ifdef USING_TABLES
 138     _setup_tables();
 139 #endif
 140     _recalc_density();
 141     _recalc_relative_humidity();
 142     live_update = true;
 143 }
 144
 145 FGEnvironment::FGEnvironment()
 146 {
 147     _init();
 148 }
 149
 150 FGEnvironment::FGEnvironment (const FGEnvironment &env)
 151 {
 152     _init();
 153     copy(env);
 154 }
 155
 156 FGEnvironment::~FGEnvironment()
 157 {
 158     Untie();
 159 }
 160
 161 FGEnvironment & FGEnvironment::operator = ( const FGEnvironment & other )
 162 {
 163     copy( other );
 164     return *this;
 165 }
 166
 167 void
 168 FGEnvironment::copy (const FGEnvironment &env)
 169 {
 170     elevation_ft = env.elevation_ft;
 171     visibility_m = env.visibility_m;
 172     temperature_sea_level_degc = env.temperature_sea_level_degc;
 173     temperature_degc = env.temperature_degc;
 174     dewpoint_sea_level_degc = env.dewpoint_sea_level_degc;
 175     dewpoint_degc = env.dewpoint_degc;
 176     pressure_sea_level_inhg = env.pressure_sea_level_inhg;
 177     wind_from_heading_deg = env.wind_from_heading_deg;
 178     wind_speed_kt = env.wind_speed_kt;
 179     wind_from_north_fps = env.wind_from_north_fps;
 180     wind_from_east_fps = env.wind_from_east_fps;
 181     wind_from_down_fps = env.wind_from_down_fps;
 182     turbulence_magnitude_norm = env.turbulence_magnitude_norm;
 183     turbulence_rate_hz = env.turbulence_rate_hz;
 184     pressure_inhg = env.pressure_inhg;
 185     density_slugft3 = env.density_slugft3;
 186     density_tropo_avg_kgm3 = env.density_tropo_avg_kgm3;
 187     relative_humidity = env.relative_humidity;
 188     altitude_half_to_sun_m = env.altitude_half_to_sun_m;
 189     altitude_tropo_top_m = env.altitude_tropo_top_m;
 190     live_update = env.live_update;
 191 }
 192
 193 static inline bool
 194 maybe_copy_value (FGEnvironment * env, const SGPropertyNode * node,
 195                   const char * name, void (FGEnvironment::*setter)(double))
 196 {
 197     const SGPropertyNode * child = node->getNode(name);
 198                                 // fragile: depends on not being typed
 199                                 // as a number
 200     if (child != 0 && child->hasValue() &&
 201         child->getStringValue()[0] != '\0') {
 202         (env->*setter)(child->getDoubleValue());
 203         return true;
 204     } else {
 205         return false;
 206     }
 207 }
 208
 209 void
 210 FGEnvironment::read (const SGPropertyNode * node)
 211 {
 212     bool live_update = set_live_update( false );
 213     maybe_copy_value(this, node, "visibility-m",
 214                      &FGEnvironment::set_visibility_m);
 215
 216     maybe_copy_value(this, node, "elevation-ft",
 217                      &FGEnvironment::set_elevation_ft);
 218
 219     if (!maybe_copy_value(this, node, "temperature-sea-level-degc",
 220                           &FGEnvironment::set_temperature_sea_level_degc)) {
 221         if( maybe_copy_value(this, node, "temperature-degc",
 222                          &FGEnvironment::set_temperature_degc)) {
 223             _recalc_sl_temperature();
 224         }
 225     }
 226
 227     if (!maybe_copy_value(this, node, "dewpoint-sea-level-degc",
 228                           &FGEnvironment::set_dewpoint_sea_level_degc)) {
 229         if( maybe_copy_value(this, node, "dewpoint-degc",
 230                          &FGEnvironment::set_dewpoint_degc)) {
 231             _recalc_sl_dewpoint();
 232         }
 233     }
 234
 235     if (!maybe_copy_value(this, node, "pressure-sea-level-inhg",
 236                           &FGEnvironment::set_pressure_sea_level_inhg)) {
 237         if( maybe_copy_value(this, node, "pressure-inhg",
 238                          &FGEnvironment::set_pressure_inhg)) {
 239             _recalc_sl_pressure();
 240         }
 241    }
 242
 243     maybe_copy_value(this, node, "wind-from-heading-deg",
 244                      &FGEnvironment::set_wind_from_heading_deg);
 245
 246     maybe_copy_value(this, node, "wind-speed-kt",
 247                      &FGEnvironment::set_wind_speed_kt);
 248
 249     maybe_copy_value(this, node, "turbulence/magnitude-norm",
 250                      &FGEnvironment::set_turbulence_magnitude_norm);
 251
 252     maybe_copy_value(this, node, "turbulence/rate-hz",
 253                      &FGEnvironment::set_turbulence_rate_hz);
 254
 255     // calculate derived properties here to avoid duplicate expensive computations
 256     _recalc_ne();
 257     _recalc_alt_pt();
 258     _recalc_alt_dewpoint();
 259     _recalc_density();
 260     _recalc_relative_humidity();
 261
 262     set_live_update(live_update);
 263 }
 264
 265 void FGEnvironment::Tie( SGPropertyNode_ptr base, bool archivable )
 266 {
 267   _tiedProperties.setRoot( base );
 268
 269   _tiedProperties.Tie( "visibility-m", this,
 270       &FGEnvironment::get_visibility_m,
 271       &FGEnvironment::set_visibility_m)
 272       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 273
 274   _tiedProperties.Tie("temperature-sea-level-degc", this,
 275       &FGEnvironment::get_temperature_sea_level_degc,
 276       &FGEnvironment::set_temperature_sea_level_degc)
 277       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 278
 279   _tiedProperties.Tie("temperature-degc", this,
 280       &FGEnvironment::get_temperature_degc,
 281       &FGEnvironment::set_temperature_degc)
 282       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 283
 284   _tiedProperties.Tie("temperature-degf", this,
 285       &FGEnvironment::get_temperature_degf);
 286
 287   _tiedProperties.Tie("dewpoint-sea-level-degc", this,
 288       &FGEnvironment::get_dewpoint_sea_level_degc,
 289       &FGEnvironment::set_dewpoint_sea_level_degc)
 290       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 291
 292   _tiedProperties.Tie("dewpoint-degc", this,
 293       &FGEnvironment::get_dewpoint_degc,
 294       &FGEnvironment::set_dewpoint_degc)
 295       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 296
 297   _tiedProperties.Tie("pressure-sea-level-inhg", this,
 298       &FGEnvironment::get_pressure_sea_level_inhg,
 299       &FGEnvironment::set_pressure_sea_level_inhg)
 300       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 301
 302   _tiedProperties.Tie("pressure-inhg", this,
 303       &FGEnvironment::get_pressure_inhg,
 304       &FGEnvironment::set_pressure_inhg)
 305       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 306
 307   _tiedProperties.Tie("density-slugft3", this,
 308       &FGEnvironment::get_density_slugft3); // read-only
 309
 310   _tiedProperties.Tie("relative-humidity", this,
 311       &FGEnvironment::get_relative_humidity); //ro
 312
 313   _tiedProperties.Tie("atmosphere/density-tropo-avg", this,
 314       &FGEnvironment::get_density_tropo_avg_kgm3); //ro
 315
 316   _tiedProperties.Tie("atmosphere/altitude-half-to-sun", this,
 317       &FGEnvironment::get_altitude_half_to_sun_m,
 318       &FGEnvironment::set_altitude_half_to_sun_m)
 319       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 320
 321   _tiedProperties.Tie("atmosphere/altitude-troposphere-top", this,
 322       &FGEnvironment::get_altitude_tropo_top_m,
 323       &FGEnvironment::set_altitude_tropo_top_m)
 324       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 325
 326   _tiedProperties.Tie("wind-from-heading-deg", this,
 327       &FGEnvironment::get_wind_from_heading_deg,
 328       &FGEnvironment::set_wind_from_heading_deg)
 329       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 330
 331   _tiedProperties.Tie("wind-speed-kt", this,
 332       &FGEnvironment::get_wind_speed_kt,
 333       &FGEnvironment::set_wind_speed_kt)
 334       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 335
 336   _tiedProperties.Tie("wind-from-north-fps", this,
 337       &FGEnvironment::get_wind_from_north_fps,
 338       &FGEnvironment::set_wind_from_north_fps)
 339       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 340
 341   _tiedProperties.Tie("wind-from-east-fps", this,
 342       &FGEnvironment::get_wind_from_east_fps,
 343       &FGEnvironment::set_wind_from_east_fps)
 344       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 345
 346   _tiedProperties.Tie("wind-from-down-fps", this,
 347       &FGEnvironment::get_wind_from_down_fps,
 348       &FGEnvironment::set_wind_from_down_fps)
 349       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 350
 351   _tiedProperties.Tie("turbulence/magnitude-norm", this,
 352       &FGEnvironment::get_turbulence_magnitude_norm,
 353       &FGEnvironment::set_turbulence_magnitude_norm)
 354       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 355
 356   _tiedProperties.Tie("turbulence/rate-hz", this,
 357       &FGEnvironment::get_turbulence_rate_hz,
 358       &FGEnvironment::set_turbulence_rate_hz)
 359       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 360 }
 361
 362 void FGEnvironment::Untie()
 363 {
 364     _tiedProperties.Untie();
 365 }
 366
 367 double
 368 FGEnvironment::get_visibility_m () const
 369 {
 370   return visibility_m;
 371 }
 372
 373 double
 374 FGEnvironment::get_temperature_sea_level_degc () const
 375 {
 376   return temperature_sea_level_degc;
 377 }
 378
 379 double
 380 FGEnvironment::get_temperature_degc () const
 381 {
 382   return temperature_degc;
 383 }
 384
 385 double
 386 FGEnvironment::get_temperature_degf () const
 387 {
 388   return (temperature_degc * 9.0 / 5.0) + 32.0;
 389 }
 390
 391 double
 392 FGEnvironment::get_dewpoint_sea_level_degc () const
 393 {
 394   return dewpoint_sea_level_degc;
 395 }
 396
 397 double
 398 FGEnvironment::get_dewpoint_degc () const
 399 {
 400   return dewpoint_degc;
 401 }
 402
 403 double
 404 FGEnvironment::get_pressure_sea_level_inhg () const
 405 {
 406   return pressure_sea_level_inhg;
 407 }
 408
 409 double
 410 FGEnvironment::get_pressure_inhg () const
 411 {
 412   return pressure_inhg;
 413 }
 414
 415 double
 416 FGEnvironment::get_density_slugft3 () const
 417 {
 418   return density_slugft3;
 419 }
 420
 421 double
 422 FGEnvironment::get_relative_humidity () const
 423 {
 424   return relative_humidity;
 425 }
 426
 427 double
 428 FGEnvironment::get_density_tropo_avg_kgm3 () const
 429 {
 430   return density_tropo_avg_kgm3;
 431 }
 432
 433 double
 434 FGEnvironment::get_altitude_half_to_sun_m () const
 435 {
 436   return altitude_half_to_sun_m;
 437 }
 438
 439 double
 440 FGEnvironment::get_altitude_tropo_top_m () const
 441 {
 442   return altitude_tropo_top_m;
 443 }
 444
 445 double
 446 FGEnvironment::get_wind_from_heading_deg () const
 447 {
 448   return wind_from_heading_deg;
 449 }
 450
 451 double
 452 FGEnvironment::get_wind_speed_kt () const
 453 {
 454   return wind_speed_kt;
 455 }
 456
 457 double
 458 FGEnvironment::get_wind_from_north_fps () const
 459 {
 460   return wind_from_north_fps;
 461 }
 462
 463 double
 464 FGEnvironment::get_wind_from_east_fps () const
 465 {
 466   return wind_from_east_fps;
 467 }
 468
 469 double
 470 FGEnvironment::get_wind_from_down_fps () const
 471 {
 472   return wind_from_down_fps;
 473 }
 474
 475 double
 476 FGEnvironment::get_turbulence_magnitude_norm () const
 477 {
 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_ENVIRONMENT, 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_ENVIRONMENT, 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_ENVIRONMENT, 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_ENVIRONMENT, 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