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