]> git.mxchange.org Git - flightgear.git/blob - src/Environment/environment.cxx
3e50bb459704e8b02f0c42e9417f6c1839783fa6
[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 #include <math.h>
29
30 #include <plib/sg.h>
31
32 #include <simgear/constants.h>
33 #include <simgear/debug/logstream.hxx>
34 #include <simgear/math/interpolater.hxx>
35 #include <simgear/props/props.hxx>
36 #include <simgear/environment/visual_enviro.hxx>
37
38 #include <Main/fg_props.hxx>
39
40 #include "environment.hxx"
41
42
43 \f
44 ////////////////////////////////////////////////////////////////////////
45 // Atmosphere model.
46 ////////////////////////////////////////////////////////////////////////
47
48 // Calculated based on the ISA standard day, as found at e.g.
49 // http://www.av8n.com/physics/altimetry.htm
50
51 // Each line of data has 3 elements:
52 //   Elevation (ft), 
53 //   temperature factor (dimensionless ratio of absolute temp), 
54 //   pressure factor (dimensionless ratio)
55 static double atmosphere_data[][3] = {
56  {  -3000.00,   1.021,  1.1133 },
57  {      0.00,   1.000,  1.0000 },
58  {   2952.76,   0.980,  0.8978 },
59  {   5905.51,   0.959,  0.8042 },
60  {   8858.27,   0.939,  0.7187 },
61  {  11811.02,   0.919,  0.6407 },
62  {  14763.78,   0.898,  0.5697 },
63  {  17716.54,   0.878,  0.5052 },
64  {  20669.29,   0.858,  0.4468 },
65  {  23622.05,   0.838,  0.3940 },
66  {  26574.80,   0.817,  0.3463 },
67  {  29527.56,   0.797,  0.3034 },
68  {  32480.31,   0.777,  0.2649 },
69  {  35433.07,   0.756,  0.2305 },
70  {  38385.83,   0.752,  0.2000 },
71  {  41338.58,   0.752,  0.1736 },
72  {  44291.34,   0.752,  0.1506 },
73  {  47244.09,   0.752,  0.1307 },
74  {  50196.85,   0.752,  0.1134 },
75  {  53149.61,   0.752,  0.0984 },
76  {  56102.36,   0.752,  0.0854 },
77  {  59055.12,   0.752,  0.0741 },
78  {  62007.87,   0.752,  0.0643 },
79  {  65000.00,   0.752,  0.0557 },
80  {  68000.00,   0.754,  0.0482 },
81  {  71000.00,   0.758,  0.0418 },
82  {  74000.00,   0.761,  0.0362 },
83  {  77000.00,   0.764,  0.0314 },
84  {  80000.00,   0.767,  0.0273 },
85  {  83000.00,   0.770,  0.0237 },
86  {  86000.00,   0.773,  0.0206 },
87  {  89000.00,   0.777,  0.0179 },
88  {  92000.00,   0.780,  0.0156 },
89  {  95000.00,   0.783,  0.0135 },
90  {  98000.00,   0.786,  0.0118 },
91  { 101000.00,   0.789,  0.0103 },
92  { -1, -1, -1 }
93 };
94
95 static SGInterpTable * _temperature_degc_table = 0;
96 static SGInterpTable * _pressure_inhg_table = 0;
97
98 static void
99 _setup_tables ()
100 {
101   if (_temperature_degc_table != 0)
102       return;
103
104   _temperature_degc_table = new SGInterpTable;
105   _pressure_inhg_table = new SGInterpTable;
106
107   for (int i = 0; atmosphere_data[i][0] != -1; i++) {
108     _temperature_degc_table->addEntry(atmosphere_data[i][0],
109                                       atmosphere_data[i][1]);
110     _pressure_inhg_table->addEntry(atmosphere_data[i][0],
111                                    atmosphere_data[i][2]);
112   }
113 }
114
115
116 \f
117 ////////////////////////////////////////////////////////////////////////
118 // Implementation of FGEnvironment.
119 ////////////////////////////////////////////////////////////////////////
120
121 FGEnvironment::FGEnvironment()
122   : elevation_ft(0),
123     visibility_m(32000),
124     temperature_sea_level_degc(15),
125     temperature_degc(15),
126     dewpoint_sea_level_degc(5), // guess
127     dewpoint_degc(5),
128     pressure_sea_level_inhg(29.92),
129     pressure_inhg(29.92),
130     turbulence_magnitude_norm(0),
131     turbulence_rate_hz(1),
132     wind_from_heading_deg(0),
133     wind_speed_kt(0),
134     wind_from_north_fps(0),
135     wind_from_east_fps(0),
136     wind_from_down_fps(0),
137     altitude_half_to_sun_m(1000),
138     altitude_tropo_top_m(10000)
139 {
140   _setup_tables();
141   _recalc_density();
142   _recalc_relative_humidity();
143 }
144
145 FGEnvironment::FGEnvironment (const FGEnvironment &env)
146 {
147     FGEnvironment();
148     copy(env);
149 }
150
151 FGEnvironment::~FGEnvironment()
152 {
153 }
154
155 void
156 FGEnvironment::copy (const FGEnvironment &env)
157 {
158     elevation_ft = env.elevation_ft;
159     visibility_m = env.visibility_m;
160     temperature_sea_level_degc = env.temperature_sea_level_degc;
161     temperature_degc = env.temperature_degc;
162     dewpoint_sea_level_degc = env.dewpoint_sea_level_degc;
163     dewpoint_degc = env.dewpoint_degc;
164     pressure_sea_level_inhg = env.pressure_sea_level_inhg;
165     wind_from_heading_deg = env.wind_from_heading_deg;
166     wind_speed_kt = env.wind_speed_kt;
167     wind_from_north_fps = env.wind_from_north_fps;
168     wind_from_east_fps = env.wind_from_east_fps;
169     wind_from_down_fps = env.wind_from_down_fps;
170     turbulence_magnitude_norm = env.turbulence_magnitude_norm;
171     turbulence_rate_hz = env.turbulence_rate_hz;
172 }
173
174 static inline bool
175 maybe_copy_value (FGEnvironment * env, const SGPropertyNode * node,
176                   const char * name, void (FGEnvironment::*setter)(double))
177 {
178     const SGPropertyNode * child = node->getNode(name);
179                                 // fragile: depends on not being typed
180                                 // as a number
181     if (child != 0 && child->hasValue() &&
182         child->getStringValue()[0] != '\0') {
183         (env->*setter)(child->getDoubleValue());
184         return true;
185     } else {
186         return false;
187     }
188 }
189
190 void
191 FGEnvironment::read (const SGPropertyNode * node)
192 {
193     maybe_copy_value(this, node, "visibility-m",
194                      &FGEnvironment::set_visibility_m);
195
196     if (!maybe_copy_value(this, node, "temperature-sea-level-degc",
197                           &FGEnvironment::set_temperature_sea_level_degc))
198         maybe_copy_value(this, node, "temperature-degc",
199                          &FGEnvironment::set_temperature_degc);
200
201     if (!maybe_copy_value(this, node, "dewpoint-sea-level-degc",
202                           &FGEnvironment::set_dewpoint_sea_level_degc))
203         maybe_copy_value(this, node, "dewpoint-degc",
204                          &FGEnvironment::set_dewpoint_degc);
205
206     if (!maybe_copy_value(this, node, "pressure-sea-level-inhg",
207                           &FGEnvironment::set_pressure_sea_level_inhg))
208         maybe_copy_value(this, node, "pressure-inhg",
209                          &FGEnvironment::set_pressure_inhg);
210
211     maybe_copy_value(this, node, "wind-from-heading-deg",
212                      &FGEnvironment::set_wind_from_heading_deg);
213
214     maybe_copy_value(this, node, "wind-speed-kt",
215                      &FGEnvironment::set_wind_speed_kt);
216
217     maybe_copy_value(this, node, "elevation-ft",
218                      &FGEnvironment::set_elevation_ft);
219
220     maybe_copy_value(this, node, "turbulence/magnitude-norm",
221                      &FGEnvironment::set_turbulence_magnitude_norm);
222
223     maybe_copy_value(this, node, "turbulence/rate-hz",
224                      &FGEnvironment::set_turbulence_rate_hz);
225 }
226
227
228 double
229 FGEnvironment::get_visibility_m () const
230 {
231   return visibility_m;
232 }
233
234 double
235 FGEnvironment::get_temperature_sea_level_degc () const
236 {
237   return temperature_sea_level_degc;
238 }
239
240 double
241 FGEnvironment::get_temperature_degc () const
242 {
243   return temperature_degc;
244 }
245
246 double
247 FGEnvironment::get_temperature_degf () const
248 {
249   return (temperature_degc * 9.0 / 5.0) + 32.0;
250 }
251
252 double
253 FGEnvironment::get_dewpoint_sea_level_degc () const
254 {
255   return dewpoint_sea_level_degc;
256 }
257
258 double
259 FGEnvironment::get_dewpoint_degc () const
260 {
261   return dewpoint_degc;
262 }
263
264 double
265 FGEnvironment::get_pressure_sea_level_inhg () const
266 {
267   return pressure_sea_level_inhg;
268 }
269
270 double
271 FGEnvironment::get_pressure_inhg () const
272 {
273   return pressure_inhg;
274 }
275
276 double
277 FGEnvironment::get_density_slugft3 () const
278 {
279   return density_slugft3;
280 }
281
282 double
283 FGEnvironment::get_relative_humidity () const
284 {
285   return relative_humidity;
286 }
287
288 double
289 FGEnvironment::get_density_tropo_avg_kgm3 () const
290 {
291   return density_tropo_avg_kgm3;
292 }
293
294 double
295 FGEnvironment::get_altitude_half_to_sun_m () const
296 {
297   return altitude_half_to_sun_m;
298 }
299
300 double
301 FGEnvironment::get_altitude_tropo_top_m () const
302 {
303   return altitude_tropo_top_m;
304 }
305
306 double
307 FGEnvironment::get_wind_from_heading_deg () const
308 {
309   return wind_from_heading_deg;
310 }
311
312 double
313 FGEnvironment::get_wind_speed_kt () const
314 {
315   return wind_speed_kt;
316 }
317
318 double
319 FGEnvironment::get_wind_from_north_fps () const
320 {
321   return wind_from_north_fps;
322 }
323
324 double
325 FGEnvironment::get_wind_from_east_fps () const
326 {
327   return wind_from_east_fps;
328 }
329
330 double
331 FGEnvironment::get_wind_from_down_fps () const
332 {
333   return wind_from_down_fps;
334 }
335
336 double
337 FGEnvironment::get_turbulence_magnitude_norm () const
338 {
339   if( sgEnviro.get_turbulence_enable_state() )
340     if (fgGetBool("/environment/params/real-world-weather-fetch") == true)
341       return sgEnviro.get_cloud_turbulence();
342   return turbulence_magnitude_norm;
343 }
344
345 double
346 FGEnvironment::get_turbulence_rate_hz () const
347 {
348   return turbulence_rate_hz;
349 }
350
351 double
352 FGEnvironment::get_elevation_ft () const
353 {
354   return elevation_ft;
355 }
356
357 void
358 FGEnvironment::set_visibility_m (double v)
359 {
360   visibility_m = v;
361 }
362
363 void
364 FGEnvironment::set_temperature_sea_level_degc (double t)
365 {
366   temperature_sea_level_degc = t;
367   if (dewpoint_sea_level_degc > t)
368       dewpoint_sea_level_degc = t;
369   _recalc_alt_temperature();
370   _recalc_density();
371 }
372
373 void
374 FGEnvironment::set_temperature_degc (double t)
375 {
376   temperature_degc = t;
377   _recalc_sl_temperature();
378   _recalc_density();
379   _recalc_relative_humidity();
380 }
381
382 void
383 FGEnvironment::set_dewpoint_sea_level_degc (double t)
384 {
385   dewpoint_sea_level_degc = t;
386   if (temperature_sea_level_degc < t)
387       temperature_sea_level_degc = t;
388   _recalc_alt_dewpoint();
389   _recalc_density();
390 }
391
392 void
393 FGEnvironment::set_dewpoint_degc (double t)
394 {
395   dewpoint_degc = t;
396   _recalc_sl_dewpoint();
397   _recalc_density();
398   _recalc_relative_humidity();
399 }
400
401 void
402 FGEnvironment::set_pressure_sea_level_inhg (double p)
403 {
404   pressure_sea_level_inhg = p;
405   _recalc_alt_pressure();
406   _recalc_density();
407 }
408
409 void
410 FGEnvironment::set_pressure_inhg (double p)
411 {
412   pressure_inhg = p;
413   _recalc_sl_pressure();
414   _recalc_density();
415 }
416
417 void
418 FGEnvironment::set_wind_from_heading_deg (double h)
419 {
420   wind_from_heading_deg = h;
421   _recalc_ne();
422 }
423
424 void
425 FGEnvironment::set_wind_speed_kt (double s)
426 {
427   wind_speed_kt = s;
428   _recalc_ne();
429 }
430
431 void
432 FGEnvironment::set_wind_from_north_fps (double n)
433 {
434   wind_from_north_fps = n;
435   _recalc_hdgspd();
436 }
437
438 void
439 FGEnvironment::set_wind_from_east_fps (double e)
440 {
441   wind_from_east_fps = e;
442   _recalc_hdgspd();
443 }
444
445 void
446 FGEnvironment::set_wind_from_down_fps (double d)
447 {
448   wind_from_down_fps = d;
449   _recalc_hdgspd();
450 }
451
452 void
453 FGEnvironment::set_turbulence_magnitude_norm (double t)
454 {
455   turbulence_magnitude_norm = t;
456 }
457
458 void
459 FGEnvironment::set_turbulence_rate_hz (double r)
460 {
461   turbulence_rate_hz = r;
462 }
463
464 void
465 FGEnvironment::set_elevation_ft (double e)
466 {
467   elevation_ft = e;
468   _recalc_alt_temperature();
469   _recalc_alt_dewpoint();
470   _recalc_alt_pressure();
471   _recalc_density();
472   _recalc_relative_humidity();
473 }
474
475 void
476 FGEnvironment::set_altitude_half_to_sun_m (double alt)
477 {
478  altitude_half_to_sun_m = alt;
479  _recalc_density_tropo_avg_kgm3();
480 }
481
482 void
483 FGEnvironment::set_altitude_tropo_top_m (double alt)
484 {
485  altitude_tropo_top_m = alt;
486  _recalc_density_tropo_avg_kgm3();
487 }
488
489
490 void
491 FGEnvironment::_recalc_hdgspd ()
492 {
493   double angle_rad;
494
495   if (wind_from_east_fps == 0) {
496     angle_rad = (wind_from_north_fps >= 0 ? SGD_PI_2 : -SGD_PI_2);
497   } else {
498     angle_rad = atan(wind_from_north_fps/wind_from_east_fps);
499   }
500   wind_from_heading_deg = angle_rad * SGD_RADIANS_TO_DEGREES;
501   if (wind_from_east_fps >= 0)
502     wind_from_heading_deg = 90 - wind_from_heading_deg;
503   else
504     wind_from_heading_deg = 270 - wind_from_heading_deg;
505
506 #if 0
507   // FIXME: Windspeed can become negative with these formulas.
508   //        which can cause problems for animations that rely
509   //        on the windspeed property.
510   if (angle_rad == 0)
511     wind_speed_kt = fabs(wind_from_east_fps
512                          * SG_METER_TO_NM * SG_FEET_TO_METER * 3600);
513   else
514     wind_speed_kt = (wind_from_north_fps / sin(angle_rad))
515       * SG_METER_TO_NM * SG_FEET_TO_METER * 3600;
516 #else
517   wind_speed_kt = sqrt(wind_from_north_fps * wind_from_north_fps +
518                        wind_from_east_fps * wind_from_east_fps) 
519                   * SG_METER_TO_NM * SG_FEET_TO_METER * 3600;
520 #endif
521 }
522
523 void
524 FGEnvironment::_recalc_ne ()
525 {
526   double speed_fps =
527     wind_speed_kt * SG_NM_TO_METER * SG_METER_TO_FEET * (1.0/3600);
528
529   wind_from_north_fps = speed_fps *
530     cos(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
531   wind_from_east_fps = speed_fps *
532     sin(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
533 }
534
535 void
536 FGEnvironment::_recalc_sl_temperature ()
537 {
538   // If we're in the stratosphere, leave sea-level temp alone
539   if (elevation_ft < 38000) {
540     temperature_sea_level_degc = (temperature_degc + 273.15)
541         / _temperature_degc_table->interpolate(elevation_ft)
542       - 273.15;
543   }
544 }
545
546 void
547 FGEnvironment::_recalc_alt_temperature ()
548 {
549   if (elevation_ft < 38000) {
550     temperature_degc = (temperature_sea_level_degc + 273.15) *
551         _temperature_degc_table->interpolate(elevation_ft) - 273.15;
552   } else {
553     temperature_degc = -56.49;  // Stratosphere is constant
554   }
555 }
556
557 void
558 FGEnvironment::_recalc_sl_dewpoint ()
559 {
560                                 // 0.2degC/1000ft
561                                 // FIXME: this will work only for low
562                                 // elevations
563   dewpoint_sea_level_degc = dewpoint_degc + (elevation_ft * .0002);
564   if (dewpoint_sea_level_degc > temperature_sea_level_degc)
565     dewpoint_sea_level_degc = temperature_sea_level_degc;
566 }
567
568 void
569 FGEnvironment::_recalc_alt_dewpoint ()
570 {
571                                 // 0.2degC/1000ft
572                                 // FIXME: this will work only for low
573                                 // elevations
574   dewpoint_degc = dewpoint_sea_level_degc + (elevation_ft * .0002);
575   if (dewpoint_degc > temperature_degc)
576     dewpoint_degc = temperature_degc;
577 }
578
579 void
580 FGEnvironment::_recalc_sl_pressure ()
581 {
582   pressure_sea_level_inhg =
583     pressure_inhg / _pressure_inhg_table->interpolate(elevation_ft);
584 }
585
586 void
587 FGEnvironment::_recalc_alt_pressure ()
588 {
589   pressure_inhg =
590     pressure_sea_level_inhg * _pressure_inhg_table->interpolate(elevation_ft);
591 }
592
593 void
594 FGEnvironment::_recalc_density ()
595 {
596   double pressure_psf = pressure_inhg * 70.7487;
597   
598   // adjust for humidity
599   // calculations taken from USA Today (oops!) at
600   // http://www.usatoday.com/weather/basics/density-calculations.htm
601   double temperature_degk = temperature_degc + 273.15;
602   double pressure_mb = pressure_inhg * 33.86;
603   double vapor_pressure_mb =
604     6.11 * pow(10.0, 7.5 * dewpoint_degc / (237.7 + dewpoint_degc));
605   double virtual_temperature_degk = temperature_degk / (1 - (vapor_pressure_mb / pressure_mb) * (1.0 - 0.622));
606   double virtual_temperature_degr = virtual_temperature_degk * 1.8;
607
608   density_slugft3 = pressure_psf / (virtual_temperature_degr * 1718);
609   _recalc_density_tropo_avg_kgm3();
610 }
611
612 // This is used to calculate the average density on the path 
613 // of sunlight to the observer for calculating sun-color
614 void
615 FGEnvironment::_recalc_density_tropo_avg_kgm3 ()
616 {
617   double pressure_mb = pressure_inhg * 33.86;
618   double vaporpressure = 6.11 * pow(10.0, ((7.5 * dewpoint_degc) / (237.7 + dewpoint_degc)));
619
620   double virtual_temp = (temperature_degc + 273.15) / (1 - 0.379 * (vaporpressure/pressure_mb));
621
622   double density_half = (100 * pressure_mb * exp(-altitude_half_to_sun_m / 8000))
623       / (287.05 * virtual_temp);
624   double density_tropo = (100 * pressure_mb * exp((-1 * altitude_tropo_top_m) / 8000))
625       / ( 287.05 * virtual_temp);
626
627   density_tropo_avg_kgm3 = ((density_slugft3 * 515.379) + density_half + density_tropo) / 3;
628 }
629
630 void
631 FGEnvironment::_recalc_relative_humidity ()
632 {
633   double vaporpressure = 6.11 * pow(10.0, ((7.5 * dewpoint_degc) / ( 237.7 + dewpoint_degc)));
634   double sat_vaporpressure = 6.11 * pow(10.0, ((7.5 * temperature_degc)
635       / ( 237.7 + temperature_degc)) );
636   relative_humidity = 100 * vaporpressure / sat_vaporpressure ;
637 }
638
639
640
641 ////////////////////////////////////////////////////////////////////////
642 // Functions.
643 ////////////////////////////////////////////////////////////////////////
644
645 static inline double
646 do_interp (double a, double b, double fraction)
647 {
648     double retval = (a + ((b - a) * fraction));
649     return retval;
650 }
651
652 static inline double
653 do_interp_deg (double a, double b, double fraction)
654 {
655     a = fmod(a, 360);
656     b = fmod(b, 360);
657     if (fabs(b-a) > 180) {
658         if (a < b)
659             a += 360;
660         else
661             b += 360;
662     }
663     return fmod(do_interp(a, b, fraction), 360);
664 }
665
666 void
667 interpolate (const FGEnvironment * env1, const FGEnvironment * env2,
668              double fraction, FGEnvironment * result)
669 {
670     result->set_visibility_m
671         (do_interp(env1->get_visibility_m(),
672                    env2->get_visibility_m(),
673                    fraction));
674
675     result->set_temperature_sea_level_degc
676         (do_interp(env1->get_temperature_sea_level_degc(),
677                    env2->get_temperature_sea_level_degc(),
678                    fraction));
679
680     result->set_dewpoint_degc
681         (do_interp(env1->get_dewpoint_sea_level_degc(),
682                    env2->get_dewpoint_sea_level_degc(),
683                    fraction));
684
685     result->set_pressure_sea_level_inhg
686         (do_interp(env1->get_pressure_sea_level_inhg(),
687                    env2->get_pressure_sea_level_inhg(),
688                    fraction));
689
690     result->set_wind_from_heading_deg
691         (do_interp_deg(env1->get_wind_from_heading_deg(),
692                        env2->get_wind_from_heading_deg(),
693                        fraction));
694
695     result->set_wind_speed_kt
696         (do_interp(env1->get_wind_speed_kt(),
697                    env2->get_wind_speed_kt(),
698                    fraction));
699
700     result->set_elevation_ft
701         (do_interp(env1->get_elevation_ft(),
702                    env2->get_elevation_ft(),
703                    fraction));
704
705     result->set_turbulence_magnitude_norm
706         (do_interp(env1->get_turbulence_magnitude_norm(),
707                    env2->get_turbulence_magnitude_norm(),
708                    fraction));
709
710     result->set_turbulence_rate_hz
711         (do_interp(env1->get_turbulence_rate_hz(),
712                    env2->get_turbulence_rate_hz(),
713                    fraction));
714 }
715
716 // end of environment.cxx