]> git.mxchange.org Git - flightgear.git/blob - src/Environment/environment.cxx
f99c31e091bd666f43ffe18ed041c6108f091da5
[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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, 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 #include <simgear/environment/visual_enviro.hxx>
40
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     altitude_half_to_sun_m(1000),
132     altitude_tropo_top_m(10000)
133 {
134   _setup_tables();
135   _recalc_density();
136   _recalc_relative_humidity();
137 }
138
139 FGEnvironment::FGEnvironment (const FGEnvironment &env)
140 {
141     FGEnvironment();
142     copy(env);
143 }
144
145 FGEnvironment::~FGEnvironment()
146 {
147 }
148
149 void
150 FGEnvironment::copy (const FGEnvironment &env)
151 {
152     elevation_ft = env.elevation_ft;
153     visibility_m = env.visibility_m;
154     temperature_sea_level_degc = env.temperature_sea_level_degc;
155     temperature_degc = env.temperature_degc;
156     dewpoint_sea_level_degc = env.dewpoint_sea_level_degc;
157     dewpoint_degc = env.dewpoint_degc;
158     pressure_sea_level_inhg = env.pressure_sea_level_inhg;
159     wind_from_heading_deg = env.wind_from_heading_deg;
160     wind_speed_kt = env.wind_speed_kt;
161     wind_from_north_fps = env.wind_from_north_fps;
162     wind_from_east_fps = env.wind_from_east_fps;
163     wind_from_down_fps = env.wind_from_down_fps;
164     turbulence_magnitude_norm = env.turbulence_magnitude_norm;
165     turbulence_rate_hz = env.turbulence_rate_hz;
166 }
167
168 static inline bool
169 maybe_copy_value (FGEnvironment * env, const SGPropertyNode * node,
170                   const char * name, void (FGEnvironment::*setter)(double))
171 {
172     const SGPropertyNode * child = node->getNode(name);
173                                 // fragile: depends on not being typed
174                                 // as a number
175     if (child != 0 && child->hasValue() &&
176         child->getStringValue()[0] != '\0') {
177         (env->*setter)(child->getDoubleValue());
178         return true;
179     } else {
180         return false;
181     }
182 }
183
184 void
185 FGEnvironment::read (const SGPropertyNode * node)
186 {
187     maybe_copy_value(this, node, "visibility-m",
188                      &FGEnvironment::set_visibility_m);
189
190     if (!maybe_copy_value(this, node, "temperature-sea-level-degc",
191                           &FGEnvironment::set_temperature_sea_level_degc))
192         maybe_copy_value(this, node, "temperature-degc",
193                          &FGEnvironment::set_temperature_degc);
194
195     if (!maybe_copy_value(this, node, "dewpoint-sea-level-degc",
196                           &FGEnvironment::set_dewpoint_sea_level_degc))
197         maybe_copy_value(this, node, "dewpoint-degc",
198                          &FGEnvironment::set_dewpoint_degc);
199
200     if (!maybe_copy_value(this, node, "pressure-sea-level-inhg",
201                           &FGEnvironment::set_pressure_sea_level_inhg))
202         maybe_copy_value(this, node, "pressure-inhg",
203                          &FGEnvironment::set_pressure_inhg);
204
205     maybe_copy_value(this, node, "wind-from-heading-deg",
206                      &FGEnvironment::set_wind_from_heading_deg);
207
208     maybe_copy_value(this, node, "wind-speed-kt",
209                      &FGEnvironment::set_wind_speed_kt);
210
211     maybe_copy_value(this, node, "elevation-ft",
212                      &FGEnvironment::set_elevation_ft);
213
214     maybe_copy_value(this, node, "turbulence/magnitude-norm",
215                      &FGEnvironment::set_turbulence_magnitude_norm);
216
217     maybe_copy_value(this, node, "turbulence/rate-hz",
218                      &FGEnvironment::set_turbulence_rate_hz);
219 }
220
221
222 double
223 FGEnvironment::get_visibility_m () const
224 {
225   return visibility_m;
226 }
227
228 double
229 FGEnvironment::get_temperature_sea_level_degc () const
230 {
231   return temperature_sea_level_degc;
232 }
233
234 double
235 FGEnvironment::get_temperature_degc () const
236 {
237   return temperature_degc;
238 }
239
240 double
241 FGEnvironment::get_temperature_degf () const
242 {
243   return (temperature_degc * 9.0 / 5.0) + 32.0;
244 }
245
246 double
247 FGEnvironment::get_dewpoint_sea_level_degc () const
248 {
249   return dewpoint_sea_level_degc;
250 }
251
252 double
253 FGEnvironment::get_dewpoint_degc () const
254 {
255   return dewpoint_degc;
256 }
257
258 double
259 FGEnvironment::get_pressure_sea_level_inhg () const
260 {
261   return pressure_sea_level_inhg;
262 }
263
264 double
265 FGEnvironment::get_pressure_inhg () const
266 {
267   return pressure_inhg;
268 }
269
270 double
271 FGEnvironment::get_density_slugft3 () const
272 {
273   return density_slugft3;
274 }
275
276 double
277 FGEnvironment::get_relative_humidity () const
278 {
279   return relative_humidity;
280 }
281
282 double
283 FGEnvironment::get_density_tropo_avg_kgm3 () const
284 {
285   return density_tropo_avg_kgm3;
286 }
287
288 double
289 FGEnvironment::get_altitude_half_to_sun_m () const
290 {
291   return altitude_half_to_sun_m;
292 }
293
294 double
295 FGEnvironment::get_altitude_tropo_top_m () const
296 {
297   return altitude_tropo_top_m;
298 }
299
300 double
301 FGEnvironment::get_wind_from_heading_deg () const
302 {
303   return wind_from_heading_deg;
304 }
305
306 double
307 FGEnvironment::get_wind_speed_kt () const
308 {
309   return wind_speed_kt;
310 }
311
312 double
313 FGEnvironment::get_wind_from_north_fps () const
314 {
315   return wind_from_north_fps;
316 }
317
318 double
319 FGEnvironment::get_wind_from_east_fps () const
320 {
321   return wind_from_east_fps;
322 }
323
324 double
325 FGEnvironment::get_wind_from_down_fps () const
326 {
327   return wind_from_down_fps;
328 }
329
330 double
331 FGEnvironment::get_turbulence_magnitude_norm () const
332 {
333   if( sgEnviro.get_turbulence_enable_state() )
334     if (fgGetBool("/environment/params/real-world-weather-fetch") == true)
335       return sgEnviro.get_cloud_turbulence();
336   return turbulence_magnitude_norm;
337 }
338
339 double
340 FGEnvironment::get_turbulence_rate_hz () const
341 {
342   return turbulence_rate_hz;
343 }
344
345 double
346 FGEnvironment::get_elevation_ft () const
347 {
348   return elevation_ft;
349 }
350
351 void
352 FGEnvironment::set_visibility_m (double v)
353 {
354   visibility_m = v;
355 }
356
357 void
358 FGEnvironment::set_temperature_sea_level_degc (double t)
359 {
360   temperature_sea_level_degc = t;
361   if (dewpoint_sea_level_degc > t)
362       dewpoint_sea_level_degc = t;
363   _recalc_alt_temperature();
364   _recalc_density();
365 }
366
367 void
368 FGEnvironment::set_temperature_degc (double t)
369 {
370   temperature_degc = t;
371   _recalc_sl_temperature();
372   _recalc_density();
373   _recalc_relative_humidity();
374 }
375
376 void
377 FGEnvironment::set_dewpoint_sea_level_degc (double t)
378 {
379   dewpoint_sea_level_degc = t;
380   if (temperature_sea_level_degc < t)
381       temperature_sea_level_degc = t;
382   _recalc_alt_dewpoint();
383   _recalc_density();
384 }
385
386 void
387 FGEnvironment::set_dewpoint_degc (double t)
388 {
389   dewpoint_degc = t;
390   _recalc_sl_dewpoint();
391   _recalc_density();
392   _recalc_relative_humidity();
393 }
394
395 void
396 FGEnvironment::set_pressure_sea_level_inhg (double p)
397 {
398   pressure_sea_level_inhg = p;
399   _recalc_alt_pressure();
400   _recalc_density();
401 }
402
403 void
404 FGEnvironment::set_pressure_inhg (double p)
405 {
406   pressure_inhg = p;
407   _recalc_sl_pressure();
408   _recalc_density();
409 }
410
411 void
412 FGEnvironment::set_wind_from_heading_deg (double h)
413 {
414   wind_from_heading_deg = h;
415   _recalc_ne();
416 }
417
418 void
419 FGEnvironment::set_wind_speed_kt (double s)
420 {
421   wind_speed_kt = s;
422   _recalc_ne();
423 }
424
425 void
426 FGEnvironment::set_wind_from_north_fps (double n)
427 {
428   wind_from_north_fps = n;
429   _recalc_hdgspd();
430 }
431
432 void
433 FGEnvironment::set_wind_from_east_fps (double e)
434 {
435   wind_from_east_fps = e;
436   _recalc_hdgspd();
437 }
438
439 void
440 FGEnvironment::set_wind_from_down_fps (double d)
441 {
442   wind_from_down_fps = d;
443   _recalc_hdgspd();
444 }
445
446 void
447 FGEnvironment::set_turbulence_magnitude_norm (double t)
448 {
449   turbulence_magnitude_norm = t;
450 }
451
452 void
453 FGEnvironment::set_turbulence_rate_hz (double r)
454 {
455   turbulence_rate_hz = r;
456 }
457
458 void
459 FGEnvironment::set_elevation_ft (double e)
460 {
461   elevation_ft = e;
462   _recalc_alt_temperature();
463   _recalc_alt_dewpoint();
464   _recalc_alt_pressure();
465   _recalc_density();
466   _recalc_relative_humidity();
467 }
468
469 void
470 FGEnvironment::set_altitude_half_to_sun_m (double alt)
471 {
472  altitude_half_to_sun_m = alt;
473  _recalc_density_tropo_avg_kgm3();
474 }
475
476 void
477 FGEnvironment::set_altitude_tropo_top_m (double alt)
478 {
479  altitude_tropo_top_m = alt;
480  _recalc_density_tropo_avg_kgm3();
481 }
482
483
484 void
485 FGEnvironment::_recalc_hdgspd ()
486 {
487   double angle_rad;
488
489   if (wind_from_east_fps == 0) {
490     angle_rad = (wind_from_north_fps >= 0 ? SGD_PI_2 : -SGD_PI_2);
491   } else {
492     angle_rad = atan(wind_from_north_fps/wind_from_east_fps);
493   }
494   wind_from_heading_deg = angle_rad * SGD_RADIANS_TO_DEGREES;
495   if (wind_from_east_fps >= 0)
496     wind_from_heading_deg = 90 - wind_from_heading_deg;
497   else
498     wind_from_heading_deg = 270 - wind_from_heading_deg;
499
500 #if 0
501   // FIXME: Windspeed can become negative with these formulas.
502   //        which can cause problems for animations that rely
503   //        on the windspeed property.
504   if (angle_rad == 0)
505     wind_speed_kt = fabs(wind_from_east_fps
506                          * SG_METER_TO_NM * SG_FEET_TO_METER * 3600);
507   else
508     wind_speed_kt = (wind_from_north_fps / sin(angle_rad))
509       * SG_METER_TO_NM * SG_FEET_TO_METER * 3600;
510 #else
511   wind_speed_kt = sqrt(wind_from_north_fps * wind_from_north_fps +
512                        wind_from_east_fps * wind_from_east_fps) 
513                   * SG_METER_TO_NM * SG_FEET_TO_METER * 3600;
514 #endif
515 }
516
517 void
518 FGEnvironment::_recalc_ne ()
519 {
520   double speed_fps =
521     wind_speed_kt * SG_NM_TO_METER * SG_METER_TO_FEET * (1.0/3600);
522
523   wind_from_north_fps = speed_fps *
524     cos(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
525   wind_from_east_fps = speed_fps *
526     sin(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
527 }
528
529 void
530 FGEnvironment::_recalc_sl_temperature ()
531 {
532   // If we're in the stratosphere, leave sea-level temp alone
533   if (elevation_ft < 38000) {
534     temperature_sea_level_degc = (temperature_degc + 273.15)
535         / _temperature_degc_table->interpolate(elevation_ft)
536       - 273.15;
537   }
538 }
539
540 void
541 FGEnvironment::_recalc_alt_temperature ()
542 {
543   if (elevation_ft < 38000) {
544     temperature_degc = (temperature_sea_level_degc + 273.15) *
545         _temperature_degc_table->interpolate(elevation_ft) - 273.15;
546   } else {
547     temperature_degc = -56.49;  // Stratosphere is constant
548   }
549 }
550
551 void
552 FGEnvironment::_recalc_sl_dewpoint ()
553 {
554                                 // 0.2degC/1000ft
555                                 // FIXME: this will work only for low
556                                 // elevations
557   dewpoint_sea_level_degc = dewpoint_degc + (elevation_ft * .0002);
558   if (dewpoint_sea_level_degc > temperature_sea_level_degc)
559     dewpoint_sea_level_degc = temperature_sea_level_degc;
560 }
561
562 void
563 FGEnvironment::_recalc_alt_dewpoint ()
564 {
565                                 // 0.2degC/1000ft
566                                 // FIXME: this will work only for low
567                                 // elevations
568   dewpoint_degc = dewpoint_sea_level_degc + (elevation_ft * .0002);
569   if (dewpoint_degc > temperature_degc)
570     dewpoint_degc = temperature_degc;
571 }
572
573 void
574 FGEnvironment::_recalc_sl_pressure ()
575 {
576   pressure_sea_level_inhg =
577     pressure_inhg / _pressure_inhg_table->interpolate(elevation_ft);
578 }
579
580 void
581 FGEnvironment::_recalc_alt_pressure ()
582 {
583   pressure_inhg =
584     pressure_sea_level_inhg * _pressure_inhg_table->interpolate(elevation_ft);
585 }
586
587 void
588 FGEnvironment::_recalc_density ()
589 {
590   double pressure_psf = pressure_inhg * 70.7487;
591   
592   // adjust for humidity
593   // calculations taken from USA Today (oops!) at
594   // http://www.usatoday.com/weather/basics/density-calculations.htm
595   double temperature_degk = temperature_degc + 273.15;
596   double pressure_mb = pressure_inhg * 33.86;
597   double vapor_pressure_mb =
598     6.11 * pow(10.0, 7.5 * dewpoint_degc / (237.7 + dewpoint_degc));
599   double virtual_temperature_degk = temperature_degk / (1 - (vapor_pressure_mb / pressure_mb) * (1.0 - 0.622));
600   double virtual_temperature_degr = virtual_temperature_degk * 1.8;
601
602   density_slugft3 = pressure_psf / (virtual_temperature_degr * 1718);
603   _recalc_density_tropo_avg_kgm3();
604 }
605
606 // This is used to calculate the average density on the path 
607 // of sunlight to the observer for calculating sun-color
608 void
609 FGEnvironment::_recalc_density_tropo_avg_kgm3 ()
610 {
611   double pressure_mb = pressure_inhg * 33.86;
612   double vaporpressure = 6.11 * pow(10, ((7.5 * dewpoint_degc) / (237.7 + dewpoint_degc)));
613
614   double virtual_temp = (temperature_degc + 273.15) / (1 - 0.379 * (vaporpressure/pressure_mb));
615
616   double density_half = (100 * pressure_mb * exp(-altitude_half_to_sun_m / 8000))
617       / (287.05 * virtual_temp);
618   double density_tropo = (100 * pressure_mb * exp((-1 * altitude_tropo_top_m) / 8000))
619       / ( 287.05 * virtual_temp);
620
621   density_tropo_avg_kgm3 = ((density_slugft3 * 515.379) + density_half + density_tropo) / 3;
622 }
623
624 void
625 FGEnvironment::_recalc_relative_humidity ()
626 {
627   double vaporpressure = 6.11 * pow(10, ((7.5 * dewpoint_degc) / ( 237.7 + dewpoint_degc)));
628   double sat_vaporpressure = 6.11 * pow(10, ((7.5 * temperature_degc)
629       / ( 237.7 + temperature_degc)) );
630   relative_humidity = 100 * vaporpressure / sat_vaporpressure ;
631 }
632
633
634
635 ////////////////////////////////////////////////////////////////////////
636 // Functions.
637 ////////////////////////////////////////////////////////////////////////
638
639 static inline double
640 do_interp (double a, double b, double fraction)
641 {
642     double retval = (a + ((b - a) * fraction));
643     return retval;
644 }
645
646 static inline double
647 do_interp_deg (double a, double b, double fraction)
648 {
649     a = fmod(a, 360);
650     b = fmod(b, 360);
651     if (fabs(b-a) > 180) {
652         if (a < b)
653             a += 360;
654         else
655             b += 360;
656     }
657     return fmod(do_interp(a, b, fraction), 360);
658 }
659
660 void
661 interpolate (const FGEnvironment * env1, const FGEnvironment * env2,
662              double fraction, FGEnvironment * result)
663 {
664     result->set_visibility_m
665         (do_interp(env1->get_visibility_m(),
666                    env2->get_visibility_m(),
667                    fraction));
668
669     result->set_temperature_sea_level_degc
670         (do_interp(env1->get_temperature_sea_level_degc(),
671                    env2->get_temperature_sea_level_degc(),
672                    fraction));
673
674     result->set_dewpoint_degc
675         (do_interp(env1->get_dewpoint_sea_level_degc(),
676                    env2->get_dewpoint_sea_level_degc(),
677                    fraction));
678
679     result->set_pressure_sea_level_inhg
680         (do_interp(env1->get_pressure_sea_level_inhg(),
681                    env2->get_pressure_sea_level_inhg(),
682                    fraction));
683
684     result->set_wind_from_heading_deg
685         (do_interp_deg(env1->get_wind_from_heading_deg(),
686                        env2->get_wind_from_heading_deg(),
687                        fraction));
688
689     result->set_wind_speed_kt
690         (do_interp(env1->get_wind_speed_kt(),
691                    env2->get_wind_speed_kt(),
692                    fraction));
693
694     result->set_elevation_ft
695         (do_interp(env1->get_elevation_ft(),
696                    env2->get_elevation_ft(),
697                    fraction));
698
699     result->set_turbulence_magnitude_norm
700         (do_interp(env1->get_turbulence_magnitude_norm(),
701                    env2->get_turbulence_magnitude_norm(),
702                    fraction));
703
704     result->set_turbulence_rate_hz
705         (do_interp(env1->get_turbulence_rate_hz(),
706                    env2->get_turbulence_rate_hz(),
707                    fraction));
708 }
709
710 // end of environment.cxx