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 }
 184
 185 static inline bool
 186 maybe_copy_value (FGEnvironment * env, const SGPropertyNode * node,
 187                   const char * name, void (FGEnvironment::*setter)(double))
 188 {
 189     const SGPropertyNode * child = node->getNode(name);
 190                                 // fragile: depends on not being typed
 191                                 // as a number
 192     if (child != 0 && child->hasValue() &&
 193         child->getStringValue()[0] != '\0') {
 194         (env->*setter)(child->getDoubleValue());
 195         return true;
 196     } else {
 197         return false;
 198     }
 199 }
 200
 201 void
 202 FGEnvironment::read (const SGPropertyNode * node)
 203 {
 204     bool live_update = set_live_update( false );
 205     maybe_copy_value(this, node, "visibility-m",
 206                      &FGEnvironment::set_visibility_m);
 207
 208     maybe_copy_value(this, node, "elevation-ft",
 209                      &FGEnvironment::set_elevation_ft);
 210
 211     if (!maybe_copy_value(this, node, "temperature-sea-level-degc",
 212                           &FGEnvironment::set_temperature_sea_level_degc)) {
 213         if( maybe_copy_value(this, node, "temperature-degc",
 214                          &FGEnvironment::set_temperature_degc)) {
 215             _recalc_sl_temperature();
 216         }
 217     }
 218
 219     if (!maybe_copy_value(this, node, "dewpoint-sea-level-degc",
 220                           &FGEnvironment::set_dewpoint_sea_level_degc)) {
 221         if( maybe_copy_value(this, node, "dewpoint-degc",
 222                          &FGEnvironment::set_dewpoint_degc)) {
 223             _recalc_sl_dewpoint();
 224         }
 225     }
 226
 227     if (!maybe_copy_value(this, node, "pressure-sea-level-inhg",
 228                           &FGEnvironment::set_pressure_sea_level_inhg)) {
 229         if( maybe_copy_value(this, node, "pressure-inhg",
 230                          &FGEnvironment::set_pressure_inhg)) {
 231             _recalc_sl_pressure();
 232         }
 233    }
 234
 235     maybe_copy_value(this, node, "wind-from-heading-deg",
 236                      &FGEnvironment::set_wind_from_heading_deg);
 237
 238     maybe_copy_value(this, node, "wind-speed-kt",
 239                      &FGEnvironment::set_wind_speed_kt);
 240
 241     maybe_copy_value(this, node, "turbulence/magnitude-norm",
 242                      &FGEnvironment::set_turbulence_magnitude_norm);
 243
 244     maybe_copy_value(this, node, "turbulence/rate-hz",
 245                      &FGEnvironment::set_turbulence_rate_hz);
 246
 247     // calculate derived properties here to avoid duplicate expensive computations
 248     _recalc_ne();
 249     _recalc_alt_pt();
 250     _recalc_alt_dewpoint();
 251     _recalc_density();
 252     _recalc_relative_humidity();
 253
 254     set_live_update(live_update);
 255 }
 256
 257 void FGEnvironment::Tie( SGPropertyNode_ptr base, bool archivable )
 258 {
 259   _tiedProperties.setRoot( base );
 260
 261   _tiedProperties.Tie( "visibility-m", this,
 262       &FGEnvironment::get_visibility_m,
 263       &FGEnvironment::set_visibility_m)
 264       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 265
 266   _tiedProperties.Tie("temperature-sea-level-degc", this,
 267       &FGEnvironment::get_temperature_sea_level_degc,
 268       &FGEnvironment::set_temperature_sea_level_degc)
 269       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 270
 271   _tiedProperties.Tie("temperature-degc", this,
 272       &FGEnvironment::get_temperature_degc,
 273       &FGEnvironment::set_temperature_degc)
 274       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 275
 276   _tiedProperties.Tie("temperature-degf", this,
 277       &FGEnvironment::get_temperature_degf);
 278
 279   _tiedProperties.Tie("dewpoint-sea-level-degc", this,
 280       &FGEnvironment::get_dewpoint_sea_level_degc,
 281       &FGEnvironment::set_dewpoint_sea_level_degc)
 282       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 283
 284   _tiedProperties.Tie("dewpoint-degc", this,
 285       &FGEnvironment::get_dewpoint_degc,
 286       &FGEnvironment::set_dewpoint_degc)
 287       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 288
 289   _tiedProperties.Tie("pressure-sea-level-inhg", this,
 290       &FGEnvironment::get_pressure_sea_level_inhg,
 291       &FGEnvironment::set_pressure_sea_level_inhg)
 292       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 293
 294   _tiedProperties.Tie("pressure-inhg", this,
 295       &FGEnvironment::get_pressure_inhg,
 296       &FGEnvironment::set_pressure_inhg)
 297       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 298
 299   _tiedProperties.Tie("density-slugft3", this,
 300       &FGEnvironment::get_density_slugft3); // read-only
 301
 302   _tiedProperties.Tie("relative-humidity", this,
 303       &FGEnvironment::get_relative_humidity); //ro
 304
 305   _tiedProperties.Tie("atmosphere/density-tropo-avg", this,
 306       &FGEnvironment::get_density_tropo_avg_kgm3); //ro
 307
 308   _tiedProperties.Tie("atmosphere/altitude-half-to-sun", this,
 309       &FGEnvironment::get_altitude_half_to_sun_m,
 310       &FGEnvironment::set_altitude_half_to_sun_m)
 311       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 312
 313   _tiedProperties.Tie("atmosphere/altitude-troposphere-top", this,
 314       &FGEnvironment::get_altitude_tropo_top_m,
 315       &FGEnvironment::set_altitude_tropo_top_m)
 316       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 317
 318   _tiedProperties.Tie("wind-from-heading-deg", this,
 319       &FGEnvironment::get_wind_from_heading_deg,
 320       &FGEnvironment::set_wind_from_heading_deg)
 321       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 322
 323   _tiedProperties.Tie("wind-speed-kt", this,
 324       &FGEnvironment::get_wind_speed_kt,
 325       &FGEnvironment::set_wind_speed_kt)
 326       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 327
 328   _tiedProperties.Tie("wind-from-north-fps", this,
 329       &FGEnvironment::get_wind_from_north_fps,
 330       &FGEnvironment::set_wind_from_north_fps)
 331       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 332
 333   _tiedProperties.Tie("wind-from-east-fps", this,
 334       &FGEnvironment::get_wind_from_east_fps,
 335       &FGEnvironment::set_wind_from_east_fps)
 336       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 337
 338   _tiedProperties.Tie("wind-from-down-fps", this,
 339       &FGEnvironment::get_wind_from_down_fps,
 340       &FGEnvironment::set_wind_from_down_fps)
 341       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 342
 343   _tiedProperties.Tie("turbulence/magnitude-norm", this,
 344       &FGEnvironment::get_turbulence_magnitude_norm,
 345       &FGEnvironment::set_turbulence_magnitude_norm)
 346       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 347
 348   _tiedProperties.Tie("turbulence/rate-hz", this,
 349       &FGEnvironment::get_turbulence_rate_hz,
 350       &FGEnvironment::set_turbulence_rate_hz)
 351       ->setAttribute( SGPropertyNode::ARCHIVE, archivable );
 352 }
 353
 354 void FGEnvironment::Untie()
 355 {
 356     _tiedProperties.Untie();
 357 }
 358
 359 double
 360 FGEnvironment::get_visibility_m () const
 361 {
 362   return visibility_m;
 363 }
 364
 365 double
 366 FGEnvironment::get_temperature_sea_level_degc () const
 367 {
 368   return temperature_sea_level_degc;
 369 }
 370
 371 double
 372 FGEnvironment::get_temperature_degc () const
 373 {
 374   return temperature_degc;
 375 }
 376
 377 double
 378 FGEnvironment::get_temperature_degf () const
 379 {
 380   return (temperature_degc * 9.0 / 5.0) + 32.0;
 381 }
 382
 383 double
 384 FGEnvironment::get_dewpoint_sea_level_degc () const
 385 {
 386   return dewpoint_sea_level_degc;
 387 }
 388
 389 double
 390 FGEnvironment::get_dewpoint_degc () const
 391 {
 392   return dewpoint_degc;
 393 }
 394
 395 double
 396 FGEnvironment::get_pressure_sea_level_inhg () const
 397 {
 398   return pressure_sea_level_inhg;
 399 }
 400
 401 double
 402 FGEnvironment::get_pressure_inhg () const
 403 {
 404   return pressure_inhg;
 405 }
 406
 407 double
 408 FGEnvironment::get_density_slugft3 () const
 409 {
 410   return density_slugft3;
 411 }
 412
 413 double
 414 FGEnvironment::get_relative_humidity () const
 415 {
 416   return relative_humidity;
 417 }
 418
 419 double
 420 FGEnvironment::get_density_tropo_avg_kgm3 () const
 421 {
 422   return density_tropo_avg_kgm3;
 423 }
 424
 425 double
 426 FGEnvironment::get_altitude_half_to_sun_m () const
 427 {
 428   return altitude_half_to_sun_m;
 429 }
 430
 431 double
 432 FGEnvironment::get_altitude_tropo_top_m () const
 433 {
 434   return altitude_tropo_top_m;
 435 }
 436
 437 double
 438 FGEnvironment::get_wind_from_heading_deg () const
 439 {
 440   return wind_from_heading_deg;
 441 }
 442
 443 double
 444 FGEnvironment::get_wind_speed_kt () const
 445 {
 446   return wind_speed_kt;
 447 }
 448
 449 double
 450 FGEnvironment::get_wind_from_north_fps () const
 451 {
 452   return wind_from_north_fps;
 453 }
 454
 455 double
 456 FGEnvironment::get_wind_from_east_fps () const
 457 {
 458   return wind_from_east_fps;
 459 }
 460
 461 double
 462 FGEnvironment::get_wind_from_down_fps () const
 463 {
 464   return wind_from_down_fps;
 465 }
 466
 467 double
 468 FGEnvironment::get_turbulence_magnitude_norm () const
 469 {
 470   return turbulence_magnitude_norm;
 471 }
 472
 473 double
 474 FGEnvironment::get_turbulence_rate_hz () const
 475 {
 476   return turbulence_rate_hz;
 477 }
 478
 479 double
 480 FGEnvironment::get_elevation_ft () const
 481 {
 482   return elevation_ft;
 483 }
 484
 485 void
 486 FGEnvironment::set_visibility_m (double v)
 487 {
 488   visibility_m = v;
 489 }
 490
 491 void
 492 FGEnvironment::set_temperature_sea_level_degc (double t)
 493 {
 494   temperature_sea_level_degc = t;
 495   if (dewpoint_sea_level_degc > t)
 496       dewpoint_sea_level_degc = t;
 497   if( live_update ) {
 498     _recalc_alt_pt();
 499     _recalc_density();
 500   }
 501 }
 502
 503 void
 504 FGEnvironment::set_temperature_degc (double t)
 505 {
 506   temperature_degc = t;
 507   if( live_update ) {
 508     _recalc_sl_temperature();
 509     _recalc_sl_pressure();
 510     _recalc_alt_pt();
 511     _recalc_density();
 512     _recalc_relative_humidity();
 513   }
 514 }
 515
 516 void
 517 FGEnvironment::set_dewpoint_sea_level_degc (double t)
 518 {
 519   dewpoint_sea_level_degc = t;
 520   if (temperature_sea_level_degc < t)
 521       temperature_sea_level_degc = t;
 522   if( live_update ) {
 523     _recalc_alt_dewpoint();
 524     _recalc_density();
 525   }
 526 }
 527
 528 void
 529 FGEnvironment::set_dewpoint_degc (double t)
 530 {
 531   dewpoint_degc = t;
 532   if( live_update ) {
 533     _recalc_sl_dewpoint();
 534     _recalc_density();
 535     _recalc_relative_humidity();
 536   }
 537 }
 538
 539 void
 540 FGEnvironment::set_pressure_sea_level_inhg (double p)
 541 {
 542   pressure_sea_level_inhg = p;
 543   if( live_update ) {
 544     _recalc_alt_pt();
 545     _recalc_density();
 546   }
 547 }
 548
 549 void
 550 FGEnvironment::set_pressure_inhg (double p)
 551 {
 552   pressure_inhg = p;
 553   if( live_update ) {
 554     _recalc_sl_pressure();
 555     _recalc_density();
 556   }
 557 }
 558
 559 void
 560 FGEnvironment::set_wind_from_heading_deg (double h)
 561 {
 562   wind_from_heading_deg = h;
 563   if( live_update ) {
 564     _recalc_ne();
 565   }
 566 }
 567
 568 void
 569 FGEnvironment::set_wind_speed_kt (double s)
 570 {
 571   wind_speed_kt = s;
 572   if( live_update ) {
 573     _recalc_ne();
 574   }
 575 }
 576
 577 void
 578 FGEnvironment::set_wind_from_north_fps (double n)
 579 {
 580   wind_from_north_fps = n;
 581   if( live_update ) {
 582     _recalc_hdgspd();
 583   }
 584 }
 585
 586 void
 587 FGEnvironment::set_wind_from_east_fps (double e)
 588 {
 589   wind_from_east_fps = e;
 590   if( live_update ) {
 591     _recalc_hdgspd();
 592   }
 593 }
 594
 595 void
 596 FGEnvironment::set_wind_from_down_fps (double d)
 597 {
 598   wind_from_down_fps = d;
 599   if( live_update ) {
 600     _recalc_hdgspd();
 601   }
 602 }
 603
 604 void
 605 FGEnvironment::set_turbulence_magnitude_norm (double t)
 606 {
 607   turbulence_magnitude_norm = t;
 608 }
 609
 610 void
 611 FGEnvironment::set_turbulence_rate_hz (double r)
 612 {
 613   turbulence_rate_hz = r;
 614 }
 615
 616 void
 617 FGEnvironment::set_elevation_ft (double e)
 618 {
 619   elevation_ft = e;
 620   if( live_update ) {
 621     _recalc_alt_pt();
 622     _recalc_alt_dewpoint();
 623     _recalc_density();
 624     _recalc_relative_humidity();
 625   }
 626 }
 627
 628 void
 629 FGEnvironment::set_altitude_half_to_sun_m (double alt)
 630 {
 631   altitude_half_to_sun_m = alt;
 632   if( live_update ) {
 633     _recalc_density_tropo_avg_kgm3();
 634   }
 635 }
 636
 637 void
 638 FGEnvironment::set_altitude_tropo_top_m (double alt)
 639 {
 640   altitude_tropo_top_m = alt;
 641   if( live_update ) {
 642     _recalc_density_tropo_avg_kgm3();
 643   }
 644 }
 645
 646
 647 void
 648 FGEnvironment::_recalc_hdgspd ()
 649 {
 650   wind_from_heading_deg =
 651     atan2(wind_from_east_fps, wind_from_north_fps) * SGD_RADIANS_TO_DEGREES;
 652
 653   if( wind_from_heading_deg < 0 )
 654     wind_from_heading_deg += 360.0;
 655
 656   wind_speed_kt = sqrt(wind_from_north_fps * wind_from_north_fps +
 657                        wind_from_east_fps * wind_from_east_fps)
 658                   * SG_METER_TO_NM * SG_FEET_TO_METER * 3600;
 659 }
 660
 661 void
 662 FGEnvironment::_recalc_ne ()
 663 {
 664   double speed_fps =
 665     wind_speed_kt * SG_NM_TO_METER * SG_METER_TO_FEET * (1.0/3600);
 666
 667   wind_from_north_fps = speed_fps *
 668     cos(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
 669   wind_from_east_fps = speed_fps *
 670     sin(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
 671 }
 672
 673 // Intended to help with the interpretation of METAR data,
 674 // not for random in-flight outside-air temperatures.
 675 void
 676 FGEnvironment::_recalc_sl_temperature ()
 677 {
 678
 679 #if 0
 680   {
 681     SG_LOG(SG_ENVIRONMENT, SG_DEBUG, "recalc_sl_temperature: using "
 682       << temperature_degc << " @ " << elevation_ft << " :: " << this);
 683   }
 684 #endif
 685
 686   if (elevation_ft * atmodel::foot >= ISA_def[1].height) {
 687     SG_LOG(SG_ENVIRONMENT, SG_ALERT, "recalc_sl_temperature: "
 688         << "valid only in troposphere, not " << elevation_ft);
 689     return;
 690   }
 691
 692 // Clamp: temperature of the stratosphere, in degrees C:
 693   double t_strato = ISA_def[1].temp - atmodel::freezing;
 694   if (temperature_degc < t_strato) temperature_sea_level_degc = t_strato;
 695   else temperature_sea_level_degc =
 696       temperature_degc + elevation_ft * atmodel::foot * ISA_def[0].lapse;
 697
 698 // Alternative implemenation:
 699 //  else temperature_sea_level_inhg = T_layer(0., elevation_ft * foot,
 700 //      pressure_inhg * inHg, temperature_degc + freezing, ISA_def[0].lapse) - freezing;
 701 }
 702
 703 void
 704 FGEnvironment::_recalc_sl_dewpoint ()
 705 {
 706                                 // 0.2degC/1000ft
 707                                 // FIXME: this will work only for low
 708                                 // elevations
 709   dewpoint_sea_level_degc = dewpoint_degc + (elevation_ft * .0002);
 710   if (dewpoint_sea_level_degc > temperature_sea_level_degc)
 711     dewpoint_sea_level_degc = temperature_sea_level_degc;
 712 }
 713
 714 void
 715 FGEnvironment::_recalc_alt_dewpoint ()
 716 {
 717                                 // 0.2degC/1000ft
 718                                 // FIXME: this will work only for low
 719                                 // elevations
 720   dewpoint_degc = dewpoint_sea_level_degc + (elevation_ft * .0002);
 721   if (dewpoint_degc > temperature_degc)
 722     dewpoint_degc = temperature_degc;
 723 }
 724
 725 void
 726 FGEnvironment::_recalc_sl_pressure ()
 727 {
 728   using namespace atmodel;
 729 #if 0
 730   {
 731     SG_LOG(SG_ENVIRONMENT, SG_ALERT, "recalc_sl_pressure: using "
 732       << pressure_inhg << " and "
 733       << temperature_degc << " @ " << elevation_ft << " :: " << this);
 734   }
 735 #endif
 736   pressure_sea_level_inhg = P_layer(0., elevation_ft * foot,
 737       pressure_inhg * inHg, temperature_degc + freezing, ISA_def[0].lapse) / inHg;
 738 }
 739
 740 // This gets called at frame rate, to account for the aircraft's
 741 // changing altitude.
 742 // Called by set_elevation_ft() which is called by FGEnvironmentMgr::update
 743
 744 void
 745 FGEnvironment::_recalc_alt_pt ()
 746 {
 747   using namespace atmodel;
 748 #if 0
 749   {
 750     static int count(0);
 751     if (++count % 1000 == 0) {
 752       SG_LOG(SG_ENVIRONMENT, SG_ALERT,
 753            "recalc_alt_pt for: " << elevation_ft
 754         << "  using "  << pressure_sea_level_inhg
 755         << "  and "  << temperature_sea_level_degc
 756         << " :: " << this
 757         << "  # " << count);
 758     }
 759   }
 760 #endif
 761   double press = pressure_inhg * inHg;
 762   double temp = temperature_degc + freezing;
 763   boost::tie(press, temp) = PT_vs_hpt(elevation_ft * foot,
 764         pressure_sea_level_inhg * inHg, temperature_sea_level_degc + freezing);
 765   temperature_degc = temp - freezing;
 766   pressure_inhg = press / inHg;
 767 }
 768
 769 void
 770 FGEnvironment::_recalc_density ()
 771 {
 772   double pressure_psf = pressure_inhg * 70.7487;
 773
 774   // adjust for humidity
 775   // calculations taken from USA Today (oops!) at
 776   // http://www.usatoday.com/weather/basics/density-calculations.htm
 777   double temperature_degk = temperature_degc + 273.15;
 778   double pressure_mb = pressure_inhg * 33.86;
 779   double vapor_pressure_mb =
 780     6.11 * pow(10.0, 7.5 * dewpoint_degc / (237.7 + dewpoint_degc));
 781   double virtual_temperature_degk = temperature_degk / (1 - (vapor_pressure_mb / pressure_mb) * (1.0 - 0.622));
 782   double virtual_temperature_degr = virtual_temperature_degk * 1.8;
 783
 784   density_slugft3 = pressure_psf / (virtual_temperature_degr * 1718);
 785   _recalc_density_tropo_avg_kgm3();
 786 }
 787
 788 // This is used to calculate the average density on the path
 789 // of sunlight to the observer for calculating sun-color
 790 void
 791 FGEnvironment::_recalc_density_tropo_avg_kgm3 ()
 792 {
 793   double pressure_mb = pressure_inhg * 33.86;
 794   double vaporpressure = 6.11 * pow(10.0, ((7.5 * dewpoint_degc) / (237.7 + dewpoint_degc)));
 795
 796   double virtual_temp = (temperature_degc + 273.15) / (1 - 0.379 * (vaporpressure/pressure_mb));
 797
 798   double density_half = (100 * pressure_mb * exp(-altitude_half_to_sun_m / 8000))
 799       / (287.05 * virtual_temp);
 800   double density_tropo = (100 * pressure_mb * exp((-1 * altitude_tropo_top_m) / 8000))
 801       / ( 287.05 * virtual_temp);
 802
 803   density_tropo_avg_kgm3 = ((density_slugft3 * 515.379) + density_half + density_tropo) / 3;
 804 }
 805
 806 void
 807 FGEnvironment::_recalc_relative_humidity ()
 808 {
 809 /*
 810   double vaporpressure = 6.11 * pow(10.0, ((7.5 * dewpoint_degc) / ( 237.7 + dewpoint_degc)));
 811   double sat_vaporpressure = 6.11 * pow(10.0, ((7.5 * temperature_degc)
 812       / ( 237.7 + temperature_degc)) );
 813   relative_humidity = 100 * vaporpressure / sat_vaporpressure ;
 814
 815   with a little algebra, this gets the same result and spares two multiplications and one pow()
 816 */
 817   double a = (7.5 * dewpoint_degc)    / ( 237.7 + dewpoint_degc);
 818   double b = (7.5 * temperature_degc) / ( 237.7 + temperature_degc);
 819   relative_humidity = 100 * pow(10.0,a-b);
 820 }
 821
 822 bool
 823 FGEnvironment::set_live_update( bool _live_update )
 824 {
 825   bool b = live_update;
 826   live_update = _live_update;
 827   return b;
 828 }
 829
 830
 831 ////////////////////////////////////////////////////////////////////////
 832 // Functions.
 833 ////////////////////////////////////////////////////////////////////////
 834
 835 static inline double
 836 do_interp (double a, double b, double fraction)
 837 {
 838     double retval = (a + ((b - a) * fraction));
 839     return retval;
 840 }
 841
 842 static inline double
 843 do_interp_deg (double a, double b, double fraction)
 844 {
 845     a = fmod(a, 360);
 846     b = fmod(b, 360);
 847     if (fabs(b-a) > 180) {
 848         if (a < b)
 849             a += 360;
 850         else
 851             b += 360;
 852     }
 853     return fmod(do_interp(a, b, fraction), 360);
 854 }
 855
 856 FGEnvironment &
 857 FGEnvironment::interpolate( const FGEnvironment & env2,
 858              double fraction, FGEnvironment * result) const
 859 {
 860     // don't calculate each internal property every time we set a single value
 861     // we trigger that at the end of the interpolation process
 862     bool live_update = result->set_live_update( false );
 863
 864     result->set_visibility_m
 865         (do_interp(get_visibility_m(),
 866                    env2.get_visibility_m(),
 867                    fraction));
 868
 869     result->set_temperature_sea_level_degc
 870         (do_interp(get_temperature_sea_level_degc(),
 871                    env2.get_temperature_sea_level_degc(),
 872                    fraction));
 873
 874     result->set_dewpoint_sea_level_degc
 875         (do_interp(get_dewpoint_sea_level_degc(),
 876                    env2.get_dewpoint_sea_level_degc(),
 877                    fraction));
 878
 879     result->set_pressure_sea_level_inhg
 880         (do_interp(get_pressure_sea_level_inhg(),
 881                    env2.get_pressure_sea_level_inhg(),
 882                    fraction));
 883
 884     result->set_wind_from_heading_deg
 885         (do_interp_deg(get_wind_from_heading_deg(),
 886                        env2.get_wind_from_heading_deg(),
 887                        fraction));
 888
 889     result->set_wind_speed_kt
 890         (do_interp(get_wind_speed_kt(),
 891                    env2.get_wind_speed_kt(),
 892                    fraction));
 893
 894     result->set_elevation_ft
 895         (do_interp(get_elevation_ft(),
 896                    env2.get_elevation_ft(),
 897                    fraction));
 898
 899     result->set_turbulence_magnitude_norm
 900         (do_interp(get_turbulence_magnitude_norm(),
 901                    env2.get_turbulence_magnitude_norm(),
 902                    fraction));
 903
 904     result->set_turbulence_rate_hz
 905         (do_interp(get_turbulence_rate_hz(),
 906                    env2.get_turbulence_rate_hz(),
 907                    fraction));
 908
 909     // calculate derived properties here to avoid duplicate expensive computations
 910     result->_recalc_ne();
 911     result->_recalc_alt_pt();
 912     result->_recalc_alt_dewpoint();
 913     result->_recalc_density();
 914     result->_recalc_relative_humidity();
 915
 916     result->set_live_update(live_update);
 917
 918     return *result;
 919 }
 920
 921 // end of environment.cxx