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