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., 675 Mass Ave, Cambridge, MA 02139, USA.
  20 //
  21 // $Id$
  22
  23
  24 #ifdef HAVE_CONFIG_H
  25 #  include <config.h>
  26 #endif
  27
  28 #ifdef HAVE_WINDOWS_H
  29 #  include <windows.h>
  30 #endif
  31
  32 #include <math.h>
  33
  34 #include <plib/sg.h>
  35
  36 #include <simgear/constants.h>
  37 #include <simgear/debug/logstream.hxx>
  38 #include <simgear/math/interpolater.hxx>
  39
  40 #include <Main/fg_props.hxx>
  41 #include "environment.hxx"
  42
  43
  44 FGEnvironment::FGEnvironment()
  45   : _temperature_degc_table(new SGInterpTable),
  46     _pressure_inhg_table(new SGInterpTable),
  47     visibility_m(32000),
  48     temperature_sea_level_degc(15),
  49     dewpoint_sea_level_degc(5), // guess
  50     pressure_sea_level_inhg(29.92),
  51     wind_from_heading_deg(0),
  52     wind_speed_kt(0),
  53     wind_from_north_fps(0),
  54     wind_from_east_fps(0),
  55     wind_from_down_fps(0)
  56 {
  57   _setup_tables();
  58 }
  59
  60 FGEnvironment::FGEnvironment (const FGEnvironment &env)
  61   : _temperature_degc_table(new SGInterpTable),
  62     _pressure_inhg_table(new SGInterpTable),
  63     elevation_ft(env.elevation_ft),
  64     visibility_m(env.visibility_m),
  65     temperature_sea_level_degc(env.temperature_sea_level_degc),
  66     dewpoint_sea_level_degc(env.dewpoint_sea_level_degc),
  67     pressure_sea_level_inhg(env.pressure_sea_level_inhg),
  68     wind_from_heading_deg(env.wind_from_heading_deg),
  69     wind_speed_kt(env.wind_speed_kt),
  70     wind_from_north_fps(env.wind_from_north_fps),
  71     wind_from_east_fps(env.wind_from_east_fps),
  72     wind_from_down_fps(env.wind_from_down_fps)
  73 {
  74   _setup_tables();
  75 }
  76
  77 FGEnvironment::~FGEnvironment()
  78 {
  79   delete _temperature_degc_table;
  80   delete _pressure_inhg_table;
  81 }
  82
  83
  84 double
  85 FGEnvironment::get_visibility_m () const
  86 {
  87   return visibility_m;
  88 }
  89
  90 double
  91 FGEnvironment::get_temperature_sea_level_degc () const
  92 {
  93   return temperature_sea_level_degc;
  94 }
  95
  96 double
  97 FGEnvironment::get_temperature_degc () const
  98 {
  99   return temperature_degc;
 100 }
 101
 102 double
 103 FGEnvironment::get_dewpoint_sea_level_degc () const
 104 {
 105   return dewpoint_sea_level_degc;
 106 }
 107
 108 double
 109 FGEnvironment::get_dewpoint_degc () const
 110 {
 111   return dewpoint_degc;
 112 }
 113
 114 double
 115 FGEnvironment::get_pressure_sea_level_inhg () const
 116 {
 117   return pressure_sea_level_inhg;
 118 }
 119
 120 double
 121 FGEnvironment::get_pressure_inhg () const
 122 {
 123   return pressure_inhg;
 124 }
 125
 126 double
 127 FGEnvironment::get_density_slugft3 () const
 128 {
 129   return density_slugft3;
 130 }
 131
 132 double
 133 FGEnvironment::get_wind_from_heading_deg () const
 134 {
 135   return wind_from_heading_deg;
 136 }
 137
 138 double
 139 FGEnvironment::get_wind_speed_kt () const
 140 {
 141   return wind_speed_kt;
 142 }
 143
 144 double
 145 FGEnvironment::get_wind_from_north_fps () const
 146 {
 147   return wind_from_north_fps;
 148 }
 149
 150 double
 151 FGEnvironment::get_wind_from_east_fps () const
 152 {
 153   return wind_from_east_fps;
 154 }
 155
 156 double
 157 FGEnvironment::get_wind_from_down_fps () const
 158 {
 159   return wind_from_down_fps;
 160 }
 161
 162 double
 163 FGEnvironment::get_elevation_ft () const
 164 {
 165   return elevation_ft;
 166 }
 167
 168 void
 169 FGEnvironment::set_visibility_m (double v)
 170 {
 171   visibility_m = v;
 172 }
 173
 174 void
 175 FGEnvironment::set_temperature_sea_level_degc (double t)
 176 {
 177   temperature_sea_level_degc = t;
 178   _recalc_alt_temperature();
 179   _recalc_density();
 180 }
 181
 182 void
 183 FGEnvironment::set_temperature_degc (double t)
 184 {
 185   temperature_degc = t;
 186   _recalc_sl_temperature();
 187   _recalc_density();
 188 }
 189
 190 void
 191 FGEnvironment::set_dewpoint_sea_level_degc (double t)
 192 {
 193   dewpoint_sea_level_degc = t;
 194   _recalc_alt_dewpoint();
 195   _recalc_density();
 196 }
 197
 198 void
 199 FGEnvironment::set_dewpoint_degc (double t)
 200 {
 201   dewpoint_degc = t;
 202   _recalc_sl_dewpoint();
 203   _recalc_density();
 204 }
 205
 206 void
 207 FGEnvironment::set_pressure_sea_level_inhg (double p)
 208 {
 209   pressure_sea_level_inhg = p;
 210   _recalc_alt_pressure();
 211   _recalc_density();
 212 }
 213
 214 void
 215 FGEnvironment::set_pressure_inhg (double p)
 216 {
 217   pressure_inhg = p;
 218   _recalc_sl_pressure();
 219   _recalc_density();
 220 }
 221
 222 void
 223 FGEnvironment::set_wind_from_heading_deg (double h)
 224 {
 225   wind_from_heading_deg = h;
 226   _recalc_ne();
 227 }
 228
 229 void
 230 FGEnvironment::set_wind_speed_kt (double s)
 231 {
 232   wind_speed_kt = s;
 233   _recalc_ne();
 234 }
 235
 236 void
 237 FGEnvironment::set_wind_from_north_fps (double n)
 238 {
 239   wind_from_north_fps = n;
 240   _recalc_hdgspd();
 241 }
 242
 243 void
 244 FGEnvironment::set_wind_from_east_fps (double e)
 245 {
 246   wind_from_east_fps = e;
 247   _recalc_hdgspd();
 248 }
 249
 250 void
 251 FGEnvironment::set_wind_from_down_fps (double d)
 252 {
 253   wind_from_down_fps = d;
 254   _recalc_hdgspd();
 255 }
 256
 257 void
 258 FGEnvironment::set_elevation_ft (double e)
 259 {
 260   elevation_ft = e;
 261   _recalc_alt_temperature();
 262   _recalc_alt_dewpoint();
 263   _recalc_alt_pressure();
 264   _recalc_density();
 265 }
 266
 267 // Atmosphere model.
 268
 269 // Copied from YASim Atmosphere.cxx, with m converted to ft, degK
 270 // converted to degC, Pa converted to inHG, and kg/m^3 converted to
 271 // slug/ft^3; they were then converted to deltas from the sea-level
 272 // defaults (approx. 15degC, 29.92inHG, and 0.00237slugs/ft^3).
 273
 274 // Original comment from YASim:
 275
 276 // Copied from McCormick, who got it from "The ARDC Model Atmosphere"
 277 // Note that there's an error in the text in the first entry,
 278 // McCormick lists 299.16/101325/1.22500, but those don't agree with
 279 // R=287.  I chose to correct the temperature to 288.20, since 79F is
 280 // pretty hot for a "standard" atmosphere.
 281
 282 // Elevation (ft), temperature factor (degK), pressure factor (inHG)
 283 static double atmosphere_data[][3] = {
 284   0.00, 1.00, 1.000,
 285   2952.76, 0.98, 0.898,
 286   5905.51, 0.96, 0.804,
 287   8858.27, 0.94, 0.719,
 288   11811.02, 0.92, 0.641,
 289   14763.78, 0.90, 0.570,
 290   17716.54, 0.88, 0.506,
 291   20669.29, 0.86, 0.447,
 292   23622.05, 0.84, 0.394,
 293   26574.80, 0.82, 0.347,
 294   29527.56, 0.80, 0.304,
 295   32480.31, 0.78, 0.266,
 296   35433.07, 0.76, 0.231,
 297   38385.83, 0.75, 0.201,
 298   41338.58, 0.75, 0.174,
 299   44291.34, 0.75, 0.151,
 300   47244.09, 0.75, 0.131,
 301   50196.85, 0.75, 0.114,
 302   53149.61, 0.75, 0.099,
 303   56102.36, 0.75, 0.086,
 304   59055.12, 0.75, 0.075,
 305   62007.87, 0.75, 0.065,
 306   -1, -1, -1
 307 };
 308
 309 void
 310 FGEnvironment::_setup_tables ()
 311 {
 312   for (int i = 0; atmosphere_data[i][0] != -1; i++) {
 313     _temperature_degc_table->addEntry(atmosphere_data[i][0],
 314                                       atmosphere_data[i][1]);
 315     _pressure_inhg_table->addEntry(atmosphere_data[i][0],
 316                                    atmosphere_data[i][2]);
 317   }
 318 }
 319
 320 void
 321 FGEnvironment::_recalc_hdgspd ()
 322 {
 323   double angle_rad;
 324
 325   if (wind_from_east_fps == 0) {
 326     angle_rad = (wind_from_north_fps >= 0 ? SGD_PI/2 : -SGD_PI/2);
 327   } else {
 328     angle_rad = atan(wind_from_north_fps/wind_from_east_fps);
 329   }
 330   wind_from_heading_deg = angle_rad * SGD_RADIANS_TO_DEGREES;
 331   if (wind_from_east_fps >= 0)
 332     wind_from_heading_deg = 90 - wind_from_heading_deg;
 333   else
 334     wind_from_heading_deg = 270 - wind_from_heading_deg;
 335   if (angle_rad == 0)
 336     wind_speed_kt = fabs(wind_from_east_fps
 337                          * SG_METER_TO_NM * SG_FEET_TO_METER * 3600);
 338   else
 339     wind_speed_kt = (wind_from_north_fps / sin(angle_rad))
 340       * SG_METER_TO_NM * SG_FEET_TO_METER * 3600;
 341 }
 342
 343 void
 344 FGEnvironment::_recalc_ne ()
 345 {
 346   double speed_fps =
 347     wind_speed_kt * SG_NM_TO_METER * SG_METER_TO_FEET * (1.0/3600);
 348
 349   wind_from_north_fps = speed_fps *
 350     cos(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
 351   wind_from_east_fps = speed_fps *
 352     sin(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
 353 }
 354
 355 void
 356 FGEnvironment::_recalc_sl_temperature ()
 357 {
 358   // If we're in the stratosphere, leave sea-level temp alone
 359   if (elevation_ft < 38000) {
 360     temperature_sea_level_degc =
 361       (temperature_degc + 273.15)
 362       /_temperature_degc_table->interpolate(elevation_ft)
 363       - 273.15;
 364   }
 365 }
 366
 367 void
 368 FGEnvironment::_recalc_alt_temperature ()
 369 {
 370   if (elevation_ft < 38000) {
 371     temperature_degc =
 372       (temperature_sea_level_degc + 273.15) *
 373       _temperature_degc_table->interpolate(elevation_ft) - 273.15;
 374   } else {
 375     temperature_degc = -56.49;  // Stratosphere is constant
 376   }
 377 }
 378
 379 void
 380 FGEnvironment::_recalc_sl_dewpoint ()
 381 {
 382                                 // 0.2degC/1000ft
 383                                 // FIXME: this will work only for low
 384                                 // elevations
 385   dewpoint_sea_level_degc = dewpoint_degc + (elevation_ft * .0002);
 386   if (dewpoint_sea_level_degc > temperature_sea_level_degc)
 387     dewpoint_sea_level_degc = temperature_sea_level_degc;
 388 }
 389
 390 void
 391 FGEnvironment::_recalc_alt_dewpoint ()
 392 {
 393                                 // 0.2degC/1000ft
 394                                 // FIXME: this will work only for low
 395                                 // elevations
 396   dewpoint_degc = dewpoint_sea_level_degc + (elevation_ft * .0002);
 397   if (dewpoint_degc > temperature_degc)
 398     dewpoint_degc = temperature_degc;
 399 }
 400
 401 void
 402 FGEnvironment::_recalc_sl_pressure ()
 403 {
 404   pressure_sea_level_inhg =
 405     pressure_inhg / _pressure_inhg_table->interpolate(elevation_ft);
 406 }
 407
 408 void
 409 FGEnvironment::_recalc_alt_pressure ()
 410 {
 411   pressure_inhg =
 412     pressure_sea_level_inhg * _pressure_inhg_table->interpolate(elevation_ft);
 413 }
 414
 415 void
 416 FGEnvironment::_recalc_density ()
 417 {
 418   double pressure_psf = pressure_inhg * 70.7487;
 419
 420   // adjust for humidity
 421   // calculations taken from USA Today (oops!) at
 422   // http://www.usatoday.com/weather/basics/density-calculations.htm
 423   double temperature_degk = temperature_degc + 273.15;
 424   double pressure_mb = pressure_inhg * 33.86;
 425   double vapor_pressure_mb =
 426     6.11 * pow(10.0, 7.5 * dewpoint_degc / (237.7 + dewpoint_degc));
 427   double virtual_temperature_degk = temperature_degk / (1 - (vapor_pressure_mb / pressure_mb) * (1.0 - 0.622));
 428   double virtual_temperature_degr = virtual_temperature_degk * 1.8;
 429
 430   density_slugft3 = pressure_psf / (virtual_temperature_degr * 1718);
 431 }
 432
 433 // end of environment.cxx