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   if (dewpoint_sea_level_degc > t)
 183       dewpoint_sea_level_degc = t;
 184   _recalc_alt_temperature();
 185   _recalc_density();
 186 }
 187
 188 void
 189 FGEnvironment::set_temperature_degc (double t)
 190 {
 191   temperature_degc = t;
 192   _recalc_sl_temperature();
 193   _recalc_density();
 194 }
 195
 196 void
 197 FGEnvironment::set_dewpoint_sea_level_degc (double t)
 198 {
 199   dewpoint_sea_level_degc = t;
 200   if (temperature_sea_level_degc < t)
 201       temperature_sea_level_degc = t;
 202   _recalc_alt_dewpoint();
 203   _recalc_density();
 204 }
 205
 206 void
 207 FGEnvironment::set_dewpoint_degc (double t)
 208 {
 209   dewpoint_degc = t;
 210   _recalc_sl_dewpoint();
 211   _recalc_density();
 212 }
 213
 214 void
 215 FGEnvironment::set_pressure_sea_level_inhg (double p)
 216 {
 217   pressure_sea_level_inhg = p;
 218   _recalc_alt_pressure();
 219   _recalc_density();
 220 }
 221
 222 void
 223 FGEnvironment::set_pressure_inhg (double p)
 224 {
 225   pressure_inhg = p;
 226   _recalc_sl_pressure();
 227   _recalc_density();
 228 }
 229
 230 void
 231 FGEnvironment::set_wind_from_heading_deg (double h)
 232 {
 233   wind_from_heading_deg = h;
 234   _recalc_ne();
 235 }
 236
 237 void
 238 FGEnvironment::set_wind_speed_kt (double s)
 239 {
 240   wind_speed_kt = s;
 241   _recalc_ne();
 242 }
 243
 244 void
 245 FGEnvironment::set_wind_from_north_fps (double n)
 246 {
 247   wind_from_north_fps = n;
 248   _recalc_hdgspd();
 249 }
 250
 251 void
 252 FGEnvironment::set_wind_from_east_fps (double e)
 253 {
 254   wind_from_east_fps = e;
 255   _recalc_hdgspd();
 256 }
 257
 258 void
 259 FGEnvironment::set_wind_from_down_fps (double d)
 260 {
 261   wind_from_down_fps = d;
 262   _recalc_hdgspd();
 263 }
 264
 265 void
 266 FGEnvironment::set_elevation_ft (double e)
 267 {
 268   elevation_ft = e;
 269   _recalc_alt_temperature();
 270   _recalc_alt_dewpoint();
 271   _recalc_alt_pressure();
 272   _recalc_density();
 273 }
 274
 275 // Atmosphere model.
 276
 277 // Copied from YASim Atmosphere.cxx, with m converted to ft, degK
 278 // converted to degC, Pa converted to inHG, and kg/m^3 converted to
 279 // slug/ft^3; they were then converted to deltas from the sea-level
 280 // defaults (approx. 15degC, 29.92inHG, and 0.00237slugs/ft^3).
 281
 282 // Original comment from YASim:
 283
 284 // Copied from McCormick, who got it from "The ARDC Model Atmosphere"
 285 // Note that there's an error in the text in the first entry,
 286 // McCormick lists 299.16/101325/1.22500, but those don't agree with
 287 // R=287.  I chose to correct the temperature to 288.20, since 79F is
 288 // pretty hot for a "standard" atmosphere.
 289
 290 // Elevation (ft), temperature factor (degK), pressure factor (inHG)
 291 static double atmosphere_data[][3] = {
 292  { 0.00, 1.00, 1.000 },
 293  { 2952.76, 0.98, 0.898 },
 294  { 5905.51, 0.96, 0.804 },
 295  { 8858.27, 0.94, 0.719 },
 296  { 11811.02, 0.92, 0.641 },
 297  { 14763.78, 0.90, 0.570 },
 298  { 17716.54, 0.88, 0.506 },
 299  { 20669.29, 0.86, 0.447 },
 300  { 23622.05, 0.84, 0.394 },
 301  { 26574.80, 0.82, 0.347 },
 302  { 29527.56, 0.80, 0.304 },
 303  { 32480.31, 0.78, 0.266 },
 304  { 35433.07, 0.76, 0.231 },
 305  { 38385.83, 0.75, 0.201 },
 306  { 41338.58, 0.75, 0.174 },
 307  { 44291.34, 0.75, 0.151 },
 308  { 47244.09, 0.75, 0.131 },
 309  { 50196.85, 0.75, 0.114 },
 310  { 53149.61, 0.75, 0.099 },
 311  { 56102.36, 0.75, 0.086 },
 312  { 59055.12, 0.75, 0.075 },
 313  { 62007.87, 0.75, 0.065 },
 314  { -1, -1, -1 }
 315 };
 316
 317 void
 318 FGEnvironment::_setup_tables ()
 319 {
 320   for (int i = 0; atmosphere_data[i][0] != -1; i++) {
 321     _temperature_degc_table->addEntry(atmosphere_data[i][0],
 322                                       atmosphere_data[i][1]);
 323     _pressure_inhg_table->addEntry(atmosphere_data[i][0],
 324                                    atmosphere_data[i][2]);
 325   }
 326 }
 327
 328 void
 329 FGEnvironment::_recalc_hdgspd ()
 330 {
 331   double angle_rad;
 332
 333   if (wind_from_east_fps == 0) {
 334     angle_rad = (wind_from_north_fps >= 0 ? SGD_PI/2 : -SGD_PI/2);
 335   } else {
 336     angle_rad = atan(wind_from_north_fps/wind_from_east_fps);
 337   }
 338   wind_from_heading_deg = angle_rad * SGD_RADIANS_TO_DEGREES;
 339   if (wind_from_east_fps >= 0)
 340     wind_from_heading_deg = 90 - wind_from_heading_deg;
 341   else
 342     wind_from_heading_deg = 270 - wind_from_heading_deg;
 343   if (angle_rad == 0)
 344     wind_speed_kt = fabs(wind_from_east_fps
 345                          * SG_METER_TO_NM * SG_FEET_TO_METER * 3600);
 346   else
 347     wind_speed_kt = (wind_from_north_fps / sin(angle_rad))
 348       * SG_METER_TO_NM * SG_FEET_TO_METER * 3600;
 349 }
 350
 351 void
 352 FGEnvironment::_recalc_ne ()
 353 {
 354   double speed_fps =
 355     wind_speed_kt * SG_NM_TO_METER * SG_METER_TO_FEET * (1.0/3600);
 356
 357   wind_from_north_fps = speed_fps *
 358     cos(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
 359   wind_from_east_fps = speed_fps *
 360     sin(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
 361 }
 362
 363 void
 364 FGEnvironment::_recalc_sl_temperature ()
 365 {
 366   // If we're in the stratosphere, leave sea-level temp alone
 367   if (elevation_ft < 38000) {
 368     temperature_sea_level_degc =
 369       (temperature_degc + 273.15)
 370       /_temperature_degc_table->interpolate(elevation_ft)
 371       - 273.15;
 372   }
 373 }
 374
 375 void
 376 FGEnvironment::_recalc_alt_temperature ()
 377 {
 378   if (elevation_ft < 38000) {
 379     temperature_degc =
 380       (temperature_sea_level_degc + 273.15) *
 381       _temperature_degc_table->interpolate(elevation_ft) - 273.15;
 382   } else {
 383     temperature_degc = -56.49;  // Stratosphere is constant
 384   }
 385 }
 386
 387 void
 388 FGEnvironment::_recalc_sl_dewpoint ()
 389 {
 390                                 // 0.2degC/1000ft
 391                                 // FIXME: this will work only for low
 392                                 // elevations
 393   dewpoint_sea_level_degc = dewpoint_degc + (elevation_ft * .0002);
 394   if (dewpoint_sea_level_degc > temperature_sea_level_degc)
 395     dewpoint_sea_level_degc = temperature_sea_level_degc;
 396 }
 397
 398 void
 399 FGEnvironment::_recalc_alt_dewpoint ()
 400 {
 401                                 // 0.2degC/1000ft
 402                                 // FIXME: this will work only for low
 403                                 // elevations
 404   dewpoint_degc = dewpoint_sea_level_degc + (elevation_ft * .0002);
 405   if (dewpoint_degc > temperature_degc)
 406     dewpoint_degc = temperature_degc;
 407 }
 408
 409 void
 410 FGEnvironment::_recalc_sl_pressure ()
 411 {
 412   pressure_sea_level_inhg =
 413     pressure_inhg / _pressure_inhg_table->interpolate(elevation_ft);
 414 }
 415
 416 void
 417 FGEnvironment::_recalc_alt_pressure ()
 418 {
 419   pressure_inhg =
 420     pressure_sea_level_inhg * _pressure_inhg_table->interpolate(elevation_ft);
 421 }
 422
 423 void
 424 FGEnvironment::_recalc_density ()
 425 {
 426   double pressure_psf = pressure_inhg * 70.7487;
 427
 428   // adjust for humidity
 429   // calculations taken from USA Today (oops!) at
 430   // http://www.usatoday.com/weather/basics/density-calculations.htm
 431   double temperature_degk = temperature_degc + 273.15;
 432   double pressure_mb = pressure_inhg * 33.86;
 433   double vapor_pressure_mb =
 434     6.11 * pow(10.0, 7.5 * dewpoint_degc / (237.7 + dewpoint_degc));
 435   double virtual_temperature_degk = temperature_degk / (1 - (vapor_pressure_mb / pressure_mb) * (1.0 - 0.622));
 436   double virtual_temperature_degr = virtual_temperature_degk * 1.8;
 437
 438   density_slugft3 = pressure_psf / (virtual_temperature_degr * 1718);
 439 }
 440
 441 // end of environment.cxx