]> git.mxchange.org Git - flightgear.git/blob - src/Environment/environment.cxx
Make FlightGear work with 16bpp:
[flightgear.git] / 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 \f
45 ////////////////////////////////////////////////////////////////////////
46 // Atmosphere model.
47 ////////////////////////////////////////////////////////////////////////
48
49 // Copied from YASim Atmosphere.cxx, with m converted to ft, degK
50 // converted to degC, Pa converted to inHG, and kg/m^3 converted to
51 // slug/ft^3; they were then converted to deltas from the sea-level
52 // defaults (approx. 15degC, 29.92inHG, and 0.00237slugs/ft^3).
53
54 // Original comment from YASim:
55
56 // Copied from McCormick, who got it from "The ARDC Model Atmosphere"
57 // Note that there's an error in the text in the first entry,
58 // McCormick lists 299.16/101325/1.22500, but those don't agree with
59 // R=287.  I chose to correct the temperature to 288.20, since 79F is
60 // pretty hot for a "standard" atmosphere.
61
62 // Elevation (ft), temperature factor (degK), pressure factor (inHG)
63 static double atmosphere_data[][3] = {
64  { 0.00, 1.00, 1.000 },
65  { 2952.76, 0.98, 0.898 },
66  { 5905.51, 0.96, 0.804 },
67  { 8858.27, 0.94, 0.719 },
68  { 11811.02, 0.92, 0.641 },
69  { 14763.78, 0.90, 0.570 },
70  { 17716.54, 0.88, 0.506 },
71  { 20669.29, 0.86, 0.447 },
72  { 23622.05, 0.84, 0.394 },
73  { 26574.80, 0.82, 0.347 },
74  { 29527.56, 0.80, 0.304 },
75  { 32480.31, 0.78, 0.266 },
76  { 35433.07, 0.76, 0.231 },
77  { 38385.83, 0.75, 0.201 },
78  { 41338.58, 0.75, 0.174 },
79  { 44291.34, 0.75, 0.151 },
80  { 47244.09, 0.75, 0.131 },
81  { 50196.85, 0.75, 0.114 },
82  { 53149.61, 0.75, 0.099 },
83  { 56102.36, 0.75, 0.086 },
84  { 59055.12, 0.75, 0.075 },
85  { 62007.87, 0.75, 0.065 },
86  { -1, -1, -1 }
87 };
88
89 static SGInterpTable * _temperature_degc_table = 0;
90 static SGInterpTable * _pressure_inhg_table = 0;
91
92 static void
93 _setup_tables ()
94 {
95   if (_temperature_degc_table != 0)
96       return;
97
98   _temperature_degc_table = new SGInterpTable;
99   _pressure_inhg_table = new SGInterpTable;
100
101   for (int i = 0; atmosphere_data[i][0] != -1; i++) {
102     _temperature_degc_table->addEntry(atmosphere_data[i][0],
103                                       atmosphere_data[i][1]);
104     _pressure_inhg_table->addEntry(atmosphere_data[i][0],
105                                    atmosphere_data[i][2]);
106   }
107 }
108
109
110 \f
111 ////////////////////////////////////////////////////////////////////////
112 // Implementation of FGEnvironment.
113 ////////////////////////////////////////////////////////////////////////
114
115 FGEnvironment::FGEnvironment()
116   : elevation_ft(0),
117     visibility_m(32000),
118     temperature_sea_level_degc(15),
119     temperature_degc(15),
120     dewpoint_sea_level_degc(5), // guess
121     dewpoint_degc(5),
122     pressure_sea_level_inhg(29.92),
123     pressure_inhg(29.92),
124     turbulence_magnitude_norm(0),
125     turbulence_rate_hz(1),
126     wind_from_heading_deg(0),
127     wind_speed_kt(0),
128     wind_from_north_fps(0),
129     wind_from_east_fps(0),
130     wind_from_down_fps(0)
131 {
132   _setup_tables();
133   _recalc_density();
134 }
135
136 FGEnvironment::FGEnvironment (const FGEnvironment &env)
137 {
138     FGEnvironment();
139     copy(env);
140 }
141
142 FGEnvironment::~FGEnvironment()
143 {
144 }
145
146 void
147 FGEnvironment::copy (const FGEnvironment &env)
148 {
149     elevation_ft = env.elevation_ft;
150     visibility_m = env.visibility_m;
151     temperature_sea_level_degc = env.temperature_sea_level_degc;
152     temperature_degc = env.temperature_degc;
153     dewpoint_sea_level_degc = env.dewpoint_sea_level_degc;
154     dewpoint_degc = env.dewpoint_degc;
155     pressure_sea_level_inhg = env.pressure_sea_level_inhg;
156     wind_from_heading_deg = env.wind_from_heading_deg;
157     wind_speed_kt = env.wind_speed_kt;
158     wind_from_north_fps = env.wind_from_north_fps;
159     wind_from_east_fps = env.wind_from_east_fps;
160     wind_from_down_fps = env.wind_from_down_fps;
161     turbulence_magnitude_norm = env.turbulence_magnitude_norm;
162     turbulence_rate_hz = env.turbulence_rate_hz;
163 }
164
165 static inline bool
166 maybe_copy_value (FGEnvironment * env, const SGPropertyNode * node,
167                   const char * name, void (FGEnvironment::*setter)(double))
168 {
169     const SGPropertyNode * child = node->getNode(name);
170                                 // fragile: depends on not being typed
171                                 // as a number
172     if (child != 0 && child->hasValue() &&
173         child->getStringValue()[0] != '\0') {
174         (env->*setter)(child->getDoubleValue());
175         return true;
176     } else {
177         return false;
178     }
179 }
180
181 void
182 FGEnvironment::read (const SGPropertyNode * node)
183 {
184     maybe_copy_value(this, node, "visibility-m",
185                      &FGEnvironment::set_visibility_m);
186
187     if (!maybe_copy_value(this, node, "temperature-sea-level-degc",
188                           &FGEnvironment::set_temperature_sea_level_degc))
189         maybe_copy_value(this, node, "temperature-degc",
190                          &FGEnvironment::set_temperature_degc);
191
192     if (!maybe_copy_value(this, node, "dewpoint-sea-level-degc",
193                           &FGEnvironment::set_dewpoint_sea_level_degc))
194         maybe_copy_value(this, node, "dewpoint-degc",
195                          &FGEnvironment::set_dewpoint_degc);
196
197     if (!maybe_copy_value(this, node, "pressure-sea-level-inhg",
198                           &FGEnvironment::set_pressure_sea_level_inhg))
199         maybe_copy_value(this, node, "pressure-inhg",
200                          &FGEnvironment::set_pressure_inhg);
201
202     maybe_copy_value(this, node, "wind-from-heading-deg",
203                      &FGEnvironment::set_wind_from_heading_deg);
204
205     maybe_copy_value(this, node, "wind-speed-kt",
206                      &FGEnvironment::set_wind_speed_kt);
207
208     maybe_copy_value(this, node, "elevation-ft",
209                      &FGEnvironment::set_elevation_ft);
210
211     maybe_copy_value(this, node, "turbulence/magnitude-norm",
212                      &FGEnvironment::set_turbulence_magnitude_norm);
213
214     maybe_copy_value(this, node, "turbulence/rate-hz",
215                      &FGEnvironment::set_turbulence_rate_hz);
216 }
217
218
219 double
220 FGEnvironment::get_visibility_m () const
221 {
222   return visibility_m;
223 }
224
225 double
226 FGEnvironment::get_temperature_sea_level_degc () const
227 {
228   return temperature_sea_level_degc;
229 }
230
231 double
232 FGEnvironment::get_temperature_degc () const
233 {
234   return temperature_degc;
235 }
236
237 double
238 FGEnvironment::get_dewpoint_sea_level_degc () const
239 {
240   return dewpoint_sea_level_degc;
241 }
242
243 double
244 FGEnvironment::get_dewpoint_degc () const
245 {
246   return dewpoint_degc;
247 }
248
249 double
250 FGEnvironment::get_pressure_sea_level_inhg () const
251 {
252   return pressure_sea_level_inhg;
253 }
254
255 double
256 FGEnvironment::get_pressure_inhg () const
257 {
258   return pressure_inhg;
259 }
260
261 double
262 FGEnvironment::get_density_slugft3 () const
263 {
264   return density_slugft3;
265 }
266
267 double
268 FGEnvironment::get_wind_from_heading_deg () const
269 {
270   return wind_from_heading_deg;
271 }
272
273 double
274 FGEnvironment::get_wind_speed_kt () const
275 {
276   return wind_speed_kt;
277 }
278
279 double
280 FGEnvironment::get_wind_from_north_fps () const
281 {
282   return wind_from_north_fps;
283 }
284
285 double
286 FGEnvironment::get_wind_from_east_fps () const
287 {
288   return wind_from_east_fps;
289 }
290
291 double
292 FGEnvironment::get_wind_from_down_fps () const
293 {
294   return wind_from_down_fps;
295 }
296
297 double
298 FGEnvironment::get_turbulence_magnitude_norm () const
299 {
300   return turbulence_magnitude_norm;
301 }
302
303 double
304 FGEnvironment::get_turbulence_rate_hz () const
305 {
306   return turbulence_rate_hz;
307 }
308
309 double
310 FGEnvironment::get_elevation_ft () const
311 {
312   return elevation_ft;
313 }
314
315 void
316 FGEnvironment::set_visibility_m (double v)
317 {
318   visibility_m = v;
319 }
320
321 void
322 FGEnvironment::set_temperature_sea_level_degc (double t)
323 {
324   temperature_sea_level_degc = t;
325   if (dewpoint_sea_level_degc > t)
326       dewpoint_sea_level_degc = t;
327   _recalc_alt_temperature();
328   _recalc_density();
329 }
330
331 void
332 FGEnvironment::set_temperature_degc (double t)
333 {
334   temperature_degc = t;
335   _recalc_sl_temperature();
336   _recalc_density();
337 }
338
339 void
340 FGEnvironment::set_dewpoint_sea_level_degc (double t)
341 {
342   dewpoint_sea_level_degc = t;
343   if (temperature_sea_level_degc < t)
344       temperature_sea_level_degc = t;
345   _recalc_alt_dewpoint();
346   _recalc_density();
347 }
348
349 void
350 FGEnvironment::set_dewpoint_degc (double t)
351 {
352   dewpoint_degc = t;
353   _recalc_sl_dewpoint();
354   _recalc_density();
355 }
356
357 void
358 FGEnvironment::set_pressure_sea_level_inhg (double p)
359 {
360   pressure_sea_level_inhg = p;
361   _recalc_alt_pressure();
362   _recalc_density();
363 }
364
365 void
366 FGEnvironment::set_pressure_inhg (double p)
367 {
368   pressure_inhg = p;
369   _recalc_sl_pressure();
370   _recalc_density();
371 }
372
373 void
374 FGEnvironment::set_wind_from_heading_deg (double h)
375 {
376   wind_from_heading_deg = h;
377   _recalc_ne();
378 }
379
380 void
381 FGEnvironment::set_wind_speed_kt (double s)
382 {
383   wind_speed_kt = s;
384   _recalc_ne();
385 }
386
387 void
388 FGEnvironment::set_wind_from_north_fps (double n)
389 {
390   wind_from_north_fps = n;
391   _recalc_hdgspd();
392 }
393
394 void
395 FGEnvironment::set_wind_from_east_fps (double e)
396 {
397   wind_from_east_fps = e;
398   _recalc_hdgspd();
399 }
400
401 void
402 FGEnvironment::set_wind_from_down_fps (double d)
403 {
404   wind_from_down_fps = d;
405   _recalc_hdgspd();
406 }
407
408 void
409 FGEnvironment::set_turbulence_magnitude_norm (double t)
410 {
411   turbulence_magnitude_norm = t;
412 }
413
414 void
415 FGEnvironment::set_turbulence_rate_hz (double r)
416 {
417   turbulence_rate_hz = r;
418 }
419
420 void
421 FGEnvironment::set_elevation_ft (double e)
422 {
423   elevation_ft = e;
424   _recalc_alt_temperature();
425   _recalc_alt_dewpoint();
426   _recalc_alt_pressure();
427   _recalc_density();
428 }
429
430 void
431 FGEnvironment::_recalc_hdgspd ()
432 {
433   double angle_rad;
434
435   if (wind_from_east_fps == 0) {
436     angle_rad = (wind_from_north_fps >= 0 ? SGD_PI/2 : -SGD_PI/2);
437   } else {
438     angle_rad = atan(wind_from_north_fps/wind_from_east_fps);
439   }
440   wind_from_heading_deg = angle_rad * SGD_RADIANS_TO_DEGREES;
441   if (wind_from_east_fps >= 0)
442     wind_from_heading_deg = 90 - wind_from_heading_deg;
443   else
444     wind_from_heading_deg = 270 - wind_from_heading_deg;
445
446 #if 0
447   // FIXME: Windspeed can become negative with these formulas.
448   //        which can cause problems for animations that rely
449   //        on the windspeed property.
450   if (angle_rad == 0)
451     wind_speed_kt = fabs(wind_from_east_fps
452                          * SG_METER_TO_NM * SG_FEET_TO_METER * 3600);
453   else
454     wind_speed_kt = (wind_from_north_fps / sin(angle_rad))
455       * SG_METER_TO_NM * SG_FEET_TO_METER * 3600;
456 #else
457   wind_speed_kt = sqrt(wind_from_north_fps * wind_from_north_fps +
458                        wind_from_east_fps * wind_from_east_fps) 
459                   * SG_METER_TO_NM * SG_FEET_TO_METER * 3600;
460 #endif
461 }
462
463 void
464 FGEnvironment::_recalc_ne ()
465 {
466   double speed_fps =
467     wind_speed_kt * SG_NM_TO_METER * SG_METER_TO_FEET * (1.0/3600);
468
469   wind_from_north_fps = speed_fps *
470     cos(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
471   wind_from_east_fps = speed_fps *
472     sin(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
473 }
474
475 void
476 FGEnvironment::_recalc_sl_temperature ()
477 {
478   // If we're in the stratosphere, leave sea-level temp alone
479   if (elevation_ft < 38000) {
480     temperature_sea_level_degc =
481       (temperature_degc + 273.15)
482       /_temperature_degc_table->interpolate(elevation_ft)
483       - 273.15;
484   }
485 }
486
487 void
488 FGEnvironment::_recalc_alt_temperature ()
489 {
490   if (elevation_ft < 38000) {
491     temperature_degc =
492       (temperature_sea_level_degc + 273.15) *
493       _temperature_degc_table->interpolate(elevation_ft) - 273.15;
494   } else {
495     temperature_degc = -56.49;  // Stratosphere is constant
496   }
497 }
498
499 void
500 FGEnvironment::_recalc_sl_dewpoint ()
501 {
502                                 // 0.2degC/1000ft
503                                 // FIXME: this will work only for low
504                                 // elevations
505   dewpoint_sea_level_degc = dewpoint_degc + (elevation_ft * .0002);
506   if (dewpoint_sea_level_degc > temperature_sea_level_degc)
507     dewpoint_sea_level_degc = temperature_sea_level_degc;
508 }
509
510 void
511 FGEnvironment::_recalc_alt_dewpoint ()
512 {
513                                 // 0.2degC/1000ft
514                                 // FIXME: this will work only for low
515                                 // elevations
516   dewpoint_degc = dewpoint_sea_level_degc + (elevation_ft * .0002);
517   if (dewpoint_degc > temperature_degc)
518     dewpoint_degc = temperature_degc;
519 }
520
521 void
522 FGEnvironment::_recalc_sl_pressure ()
523 {
524   pressure_sea_level_inhg =
525     pressure_inhg / _pressure_inhg_table->interpolate(elevation_ft);
526 }
527
528 void
529 FGEnvironment::_recalc_alt_pressure ()
530 {
531   pressure_inhg =
532     pressure_sea_level_inhg * _pressure_inhg_table->interpolate(elevation_ft);
533 }
534
535 void
536 FGEnvironment::_recalc_density ()
537 {
538   double pressure_psf = pressure_inhg * 70.7487;
539   
540   // adjust for humidity
541   // calculations taken from USA Today (oops!) at
542   // http://www.usatoday.com/weather/basics/density-calculations.htm
543   double temperature_degk = temperature_degc + 273.15;
544   double pressure_mb = pressure_inhg * 33.86;
545   double vapor_pressure_mb =
546     6.11 * pow(10.0, 7.5 * dewpoint_degc / (237.7 + dewpoint_degc));
547   double virtual_temperature_degk = temperature_degk / (1 - (vapor_pressure_mb / pressure_mb) * (1.0 - 0.622));
548   double virtual_temperature_degr = virtual_temperature_degk * 1.8;
549
550   density_slugft3 = pressure_psf / (virtual_temperature_degr * 1718);
551 }
552
553
554 \f
555 ////////////////////////////////////////////////////////////////////////
556 // Functions.
557 ////////////////////////////////////////////////////////////////////////
558
559 static inline double
560 do_interp (double a, double b, double fraction)
561 {
562     double retval = (a + ((b - a) * fraction));
563     return retval;
564 }
565
566 static inline double
567 do_interp_deg (double a, double b, double fraction)
568 {
569     a = fmod(a, 360);
570     b = fmod(b, 360);
571     if (fabs(b-a) > 180) {
572         if (a < b)
573             a += 360;
574         else
575             b += 360;
576     }
577     return fmod(do_interp(a, b, fraction), 360);
578 }
579
580 void
581 interpolate (const FGEnvironment * env1, const FGEnvironment * env2,
582              double fraction, FGEnvironment * result)
583 {
584     result->set_visibility_m
585         (do_interp(env1->get_visibility_m(),
586                    env2->get_visibility_m(),
587                    fraction));
588
589     result->set_temperature_sea_level_degc
590         (do_interp(env1->get_temperature_sea_level_degc(),
591                    env2->get_temperature_sea_level_degc(),
592                    fraction));
593
594     result->set_dewpoint_degc
595         (do_interp(env1->get_dewpoint_sea_level_degc(),
596                    env2->get_dewpoint_sea_level_degc(),
597                    fraction));
598
599     result->set_pressure_sea_level_inhg
600         (do_interp(env1->get_pressure_sea_level_inhg(),
601                    env2->get_pressure_sea_level_inhg(),
602                    fraction));
603
604     result->set_wind_from_heading_deg
605         (do_interp_deg(env1->get_wind_from_heading_deg(),
606                        env2->get_wind_from_heading_deg(),
607                        fraction));
608
609     result->set_wind_speed_kt
610         (do_interp(env1->get_wind_speed_kt(),
611                    env2->get_wind_speed_kt(),
612                    fraction));
613
614     result->set_elevation_ft
615         (do_interp(env1->get_elevation_ft(),
616                    env2->get_elevation_ft(),
617                    fraction));
618
619     result->set_turbulence_magnitude_norm
620         (do_interp(env1->get_turbulence_magnitude_norm(),
621                    env2->get_turbulence_magnitude_norm(),
622                    fraction));
623
624     result->set_turbulence_rate_hz
625         (do_interp(env1->get_turbulence_rate_hz(),
626                    env2->get_turbulence_rate_hz(),
627                    fraction));
628 }
629
630 // end of environment.cxx