src/Environment/environment_ctrl.cxx

   1 // environment_ctrl.cxx -- manager for natural environment information.
   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 #ifdef HAVE_CONFIG_H
  24 #  include "config.h"
  25 #endif
  26
  27 #include <algorithm>
  28
  29 #include <simgear/debug/logstream.hxx>
  30 #include <simgear/structure/commands.hxx>
  31 #include <simgear/structure/exception.hxx>
  32
  33 #include <Airports/simple.hxx>
  34 #include <Main/fg_props.hxx>
  35 #include <Main/util.hxx>
  36
  37 #include "atmosphere.hxx"
  38 #include "fgmetar.hxx"
  39 #include "environment_ctrl.hxx"
  40
  41 using std::sort;
  42
  43 class AirportWithMetar : public FGAirport::AirportFilter {
  44 public:
  45         virtual bool passAirport(FGAirport* aApt) const {
  46                 return aApt->getMetar();
  47         }
  48
  49   // permit heliports and seaports too
  50   virtual FGPositioned::Type maxType() const
  51   { return FGPositioned::SEAPORT; }
  52 };
  53
  54 static AirportWithMetar airportWithMetarFilter;
  55 \f
  56 ////////////////////////////////////////////////////////////////////////
  57 // Implementation of FGEnvironmentCtrl abstract base class.
  58 ////////////////////////////////////////////////////////////////////////
  59
  60 FGEnvironmentCtrl::FGEnvironmentCtrl ()
  61         : _environment(0),
  62         _lon_deg(0),
  63         _lat_deg(0),
  64         _elev_ft(0)
  65 {
  66 }
  67
  68 FGEnvironmentCtrl::~FGEnvironmentCtrl ()
  69 {
  70 }
  71
  72 void
  73 FGEnvironmentCtrl::setEnvironment (FGEnvironment * environment)
  74 {
  75         _environment = environment;
  76 }
  77
  78 void
  79 FGEnvironmentCtrl::setLongitudeDeg (double lon_deg)
  80 {
  81         _lon_deg = lon_deg;
  82 }
  83
  84 void
  85 FGEnvironmentCtrl::setLatitudeDeg (double lat_deg)
  86 {
  87         _lat_deg = lat_deg;
  88 }
  89
  90 void
  91 FGEnvironmentCtrl::setElevationFt (double elev_ft)
  92 {
  93         _elev_ft = elev_ft;
  94 }
  95
  96 void
  97 FGEnvironmentCtrl::setPosition (double lon_deg, double lat_deg, double elev_ft)
  98 {
  99         _lon_deg = lon_deg;
 100         _lat_deg = lat_deg;
 101         _elev_ft = elev_ft;
 102 }
 103
 104
 105 \f
 106 ////////////////////////////////////////////////////////////////////////
 107 // Implementation of FGInterpolateEnvironmentCtrl.
 108 ////////////////////////////////////////////////////////////////////////
 109
 110
 111 FGInterpolateEnvironmentCtrl::FGInterpolateEnvironmentCtrl ()
 112 {
 113         altitude_n = fgGetNode("/position/altitude-ft", true);
 114         altitude_agl_n = fgGetNode("/position/altitude-agl-ft", true);
 115         boundary_transition_n = fgGetNode("/environment/config/boundary-transition-ft", false );
 116         boundary_n = fgGetNode("/environment/config/boundary", true );
 117         aloft_n = fgGetNode("/environment/config/aloft", true );
 118 }
 119
 120 FGInterpolateEnvironmentCtrl::~FGInterpolateEnvironmentCtrl ()
 121 {
 122         unsigned int i;
 123         for (i = 0; i < _boundary_table.size(); i++)
 124                 delete _boundary_table[i];
 125         for (i = 0; i < _aloft_table.size(); i++)
 126                 delete _aloft_table[i];
 127 }
 128
 129
 130
 131 void
 132 FGInterpolateEnvironmentCtrl::init ()
 133 {
 134         read_table( boundary_n, _boundary_table);
 135         read_table( aloft_n, _aloft_table);
 136 }
 137
 138 void
 139 FGInterpolateEnvironmentCtrl::reinit ()
 140 {
 141         init();
 142 }
 143
 144 void
 145 FGInterpolateEnvironmentCtrl::read_table (const SGPropertyNode * node, vector<bucket *> &table)
 146 {
 147         double last_altitude_ft = 0.0;
 148         double sort_required = false;
 149         size_t i;
 150
 151         for (i = 0; i < (size_t)node->nChildren(); i++) {
 152                 const SGPropertyNode * child = node->getChild(i);
 153                 if ( strcmp(child->getName(), "entry") == 0
 154                  && child->getStringValue("elevation-ft", "")[0] != '\0'
 155                  && ( child->getDoubleValue("elevation-ft") > 0.1 || i == 0 ) )
 156         {
 157                         bucket * b;
 158                         if( i < table.size() ) {
 159                                 // recycle existing bucket
 160                                 b = table[i];
 161                         } else {
 162                                 // more nodes than buckets in table, add a new one
 163                                 b = new bucket;
 164                                 table.push_back(b);
 165                         }
 166                         if (i > 0)
 167                                 b->environment.copy(table[i-1]->environment);
 168                         b->environment.read(child);
 169                         b->altitude_ft = b->environment.get_elevation_ft();
 170
 171                         // check, if altitudes are in ascending order
 172                         if( b->altitude_ft < last_altitude_ft )
 173                                 sort_required = true;
 174                         last_altitude_ft = b->altitude_ft;
 175                 }
 176         }
 177         // remove leftover buckets
 178         while( table.size() > i ) {
 179                 bucket * b = *(table.end() - 1);
 180                 delete b;
 181                 table.pop_back();
 182         }
 183
 184         if( sort_required )
 185                 sort(table.begin(), table.end(), bucket::lessThan);
 186
 187         // cleanup entries with (almost)same altitude
 188         for( vector<bucket *>::size_type n = 1; n < table.size(); n++ ) {
 189                 if( fabs(table[n]->altitude_ft - table[n-1]->altitude_ft ) < 1 ) {
 190                         SG_LOG( SG_GENERAL, SG_ALERT, "Removing duplicate altitude entry in environment config for altitude " << table[n]->altitude_ft );
 191                         table.erase( table.begin() + n );
 192                 }
 193         }
 194 }
 195
 196 void
 197 FGInterpolateEnvironmentCtrl::update (double delta_time_sec)
 198 {
 199         double altitude_ft = altitude_n->getDoubleValue();
 200         double altitude_agl_ft = altitude_agl_n->getDoubleValue();
 201         double boundary_transition =
 202                 boundary_transition_n == NULL ? 500 : boundary_transition_n->getDoubleValue();
 203
 204         int length = _boundary_table.size();
 205
 206         if (length > 0) {
 207                 // boundary table
 208                 double boundary_limit = _boundary_table[length-1]->altitude_ft;
 209                 if (boundary_limit >= altitude_agl_ft) {
 210                         do_interpolate(_boundary_table, altitude_agl_ft, _environment);
 211                         return;
 212                 } else if ((boundary_limit + boundary_transition) >= altitude_agl_ft) {
 213                         //TODO: this is 500ft above the top altitude of boundary layer
 214                         //shouldn't this be +/-250 ft off of the top altitude?
 215                         // both tables
 216                         do_interpolate(_boundary_table, altitude_agl_ft, &env1);
 217                         do_interpolate(_aloft_table, altitude_ft, &env2);
 218                         double fraction = boundary_transition > SGLimitsd::min() ?
 219                                 (altitude_agl_ft - boundary_limit) / boundary_transition : 1.0;
 220                         interpolate(&env1, &env2, fraction, _environment);
 221                         return;
 222                 }
 223         }
 224         // aloft table
 225         do_interpolate(_aloft_table, altitude_ft, _environment);
 226 }
 227
 228 void
 229 FGInterpolateEnvironmentCtrl::do_interpolate (vector<bucket *> &table, double altitude_ft, FGEnvironment * environment)
 230 {
 231         int length = table.size();
 232         if (length == 0)
 233                 return;
 234
 235         // Boundary conditions
 236         if ((length == 1) || (table[0]->altitude_ft >= altitude_ft)) {
 237                 environment->copy(table[0]->environment); // below bottom of table
 238                 return;
 239         } else if (table[length-1]->altitude_ft <= altitude_ft) {
 240                 environment->copy(table[length-1]->environment); // above top of table
 241                 return;
 242         }
 243
 244         // Search the interpolation table
 245         int layer;
 246         for ( layer = 1; // can't be below bottom layer, handled above
 247               layer < length && table[layer]->altitude_ft <= altitude_ft;
 248               layer++);
 249         FGEnvironment * env1 = &(table[layer-1]->environment);
 250         FGEnvironment * env2 = &(table[layer]->environment);
 251         // two layers of same altitude were sorted out in read_table
 252         double fraction = ((altitude_ft - table[layer-1]->altitude_ft) /
 253                           (table[layer]->altitude_ft - table[layer-1]->altitude_ft));
 254         interpolate(env1, env2, fraction, environment);
 255 }
 256
 257 bool
 258 FGInterpolateEnvironmentCtrl::bucket::operator< (const bucket &b) const
 259 {
 260         return (altitude_ft < b.altitude_ft);
 261 }
 262
 263 bool
 264 FGInterpolateEnvironmentCtrl::bucket::lessThan(bucket *a, bucket *b)
 265 {
 266         return (a->altitude_ft) < (b->altitude_ft);
 267 }
 268
 269 \f
 270 ////////////////////////////////////////////////////////////////////////
 271 // Implementation of FGMetarCtrl.
 272 ////////////////////////////////////////////////////////////////////////
 273
 274 FGMetarCtrl::FGMetarCtrl( SGSubsystem * environmentCtrl )
 275         :
 276         metar_valid(false),
 277         setup_winds_aloft(true),
 278         wind_interpolation_required(true),
 279         metar_sealevel_temperature(15.0),
 280         metar_sealevel_dewpoint(5.0),
 281         // Interpolation constant definitions.
 282         MaxWindChangeKtsSec( 0.2 ),
 283         MaxVisChangePercentSec( 0.05 ),
 284         MaxPressureChangeInHgSec( 0.0005 ), // approx 1hpa/min
 285         MaxTemperatureChangeDegcSec(10.0/60.0), // approx 10degc/min
 286         MaxCloudAltitudeChangeFtSec( 20.0 ),
 287         MaxCloudThicknessChangeFtSec( 50.0 ),
 288         MaxCloudInterpolationHeightFt( 5000.0 ),
 289         MaxCloudInterpolationDeltaFt( 4000.0 ),
 290         _environmentCtrl(environmentCtrl)
 291 {
 292         windModulator = new FGBasicWindModulator();
 293
 294         metar_base_n = fgGetNode( "/environment/metar", true );
 295         station_id_n = metar_base_n->getNode("station-id", true );
 296         station_elevation_n = metar_base_n->getNode("station-elevation-ft", true );
 297         min_visibility_n = metar_base_n->getNode("min-visibility-m", true );
 298         max_visibility_n = metar_base_n->getNode("max-visibility-m", true );
 299         base_wind_range_from_n = metar_base_n->getNode("base-wind-range-from", true );
 300         base_wind_range_to_n = metar_base_n->getNode("base-wind-range-to", true );
 301         base_wind_speed_n = metar_base_n->getNode("base-wind-speed-kt", true );
 302         base_wind_dir_n = metar_base_n->getNode("base-wind-dir-deg", true );
 303         gust_wind_speed_n = metar_base_n->getNode("gust-wind-speed-kt", true );
 304         temperature_n = metar_base_n->getNode("temperature-degc", true );
 305         dewpoint_n = metar_base_n->getNode("dewpoint-degc", true );
 306         humidity_n = metar_base_n->getNode("rel-humidity-norm", true );
 307         pressure_n = metar_base_n->getNode("pressure-inhg", true );
 308         clouds_n = metar_base_n->getNode("clouds", true );
 309         rain_n = metar_base_n->getNode("rain-norm", true );
 310         hail_n = metar_base_n->getNode("hail-norm", true );
 311         snow_n = metar_base_n->getNode("snow-norm", true );
 312         snow_cover_n = metar_base_n->getNode("snow-cover", true );
 313         magnetic_variation_n = fgGetNode( "/environment/magnetic-variation-deg", true );
 314         ground_elevation_n = fgGetNode( "/position/ground-elev-m", true );
 315         longitude_n = fgGetNode( "/position/longitude-deg", true );
 316         latitude_n = fgGetNode( "/position/latitude-deg", true );
 317         environment_clouds_n = fgGetNode("/environment/clouds");
 318
 319         boundary_wind_speed_n = fgGetNode("/environment/config/boundary/entry/wind-speed-kt", true );
 320         boundary_wind_from_heading_n = fgGetNode("/environment/config/boundary/entry/wind-from-heading-deg", true );
 321         boundary_visibility_n = fgGetNode("/environment/config/boundary/entry/visibility-m", true );
 322         boundary_sea_level_pressure_n = fgGetNode("/environment/config/boundary/entry/pressure-sea-level-inhg", true );
 323         boundary_sea_level_temperature_n = fgGetNode("/environment/config/boundary/entry/temperature-sea-level-degc", true );
 324         boundary_sea_level_dewpoint_n = fgGetNode("/environment/config/boundary/entry/dewpoint-sea-level-degc", true );
 325 }
 326
 327 FGMetarCtrl::~FGMetarCtrl ()
 328 {
 329 }
 330
 331 void FGMetarCtrl::bind ()
 332 {
 333         fgTie("/environment/metar/valid", this, &FGMetarCtrl::get_valid );
 334         fgTie("/environment/params/metar-updates-environment", this, &FGMetarCtrl::get_enabled, &FGMetarCtrl::set_enabled );
 335         fgTie("/environment/params/metar-updates-winds-aloft", this, &FGMetarCtrl::get_setup_winds_aloft, &FGMetarCtrl::set_setup_winds_aloft );
 336 }
 337
 338 void FGMetarCtrl::unbind ()
 339 {
 340         fgUntie("/environment/metar/valid");
 341         fgUntie("/environment/params/metar-updates-environment");
 342         fgUntie("/environment/params/metar-updates-winds-aloft");
 343 }
 344
 345 // use a "command" to set station temp at station elevation
 346 static void set_temp_at_altitude( double temp_degc, double altitude_ft ) {
 347         SGPropertyNode args;
 348         SGPropertyNode *node = args.getNode("temp-degc", 0, true);
 349         node->setDoubleValue( temp_degc );
 350         node = args.getNode("altitude-ft", 0, true);
 351         node->setDoubleValue( altitude_ft );
 352         globals->get_commands()->execute( altitude_ft == 0.0 ?
 353                 "set-sea-level-air-temp-degc" :
 354                 "set-outside-air-temp-degc", &args);
 355 }
 356
 357 static void set_dewpoint_at_altitude( double dewpoint_degc, double altitude_ft ) {
 358         SGPropertyNode args;
 359         SGPropertyNode *node = args.getNode("dewpoint-degc", 0, true);
 360         node->setDoubleValue( dewpoint_degc );
 361         node = args.getNode("altitude-ft", 0, true);
 362         node->setDoubleValue( altitude_ft );
 363         globals->get_commands()->execute( altitude_ft == 0.0 ?
 364                 "set-dewpoint-sea-level-air-temp-degc" :
 365                 "set-dewpoint-temp-degc", &args);
 366 }
 367
 368 /*
 369  Setup the wind nodes for a branch in the /environment/config/<branchName>/entry nodes
 370
 371  Output properties:
 372  wind-from-heading-deg
 373  wind-speed-kt
 374  turbulence/magnitude-norm
 375
 376  Input properties:
 377  wind-heading-change-deg       how many degrees does the wind direction change at this level
 378  wind-speed-change-rel         relative change of wind speed at this level
 379  turbulence/factor             factor for the calculated turbulence magnitude at this level
 380  */
 381 static void setupWindBranch( string branchName, double dir, double speed, double gust )
 382 {
 383         SGPropertyNode_ptr branch = fgGetNode("/environment/config", true)->getNode(branchName,true);
 384         vector<SGPropertyNode_ptr> entries = branch->getChildren("entry");
 385         for ( vector<SGPropertyNode_ptr>::iterator it = entries.begin(); it != entries.end(); it++) {
 386
 387                 // change wind direction as configured
 388                 double layer_dir = dir + (*it)->getDoubleValue("wind-heading-change-deg", 0.0 );
 389                 if( layer_dir >= 360.0 ) layer_dir -= 360.0;
 390                 if( layer_dir < 0.0 ) layer_dir += 360.0;
 391                 (*it)->setDoubleValue("wind-from-heading-deg", layer_dir);
 392
 393                 double layer_speed = speed*(1 + (*it)->getDoubleValue("wind-speed-change-rel", 0.0 ));
 394                 (*it)->setDoubleValue("wind-speed-kt", layer_speed );
 395
 396                 // add some turbulence
 397                 SGPropertyNode_ptr turbulence = (*it)->getNode("turbulence",true);
 398
 399                 double turbulence_norm = speed/50;
 400                 if( gust > speed ) {
 401                         turbulence_norm += (gust-speed)/25;
 402                 }
 403                 if( turbulence_norm > 1.0 ) turbulence_norm = 1.0;
 404
 405                 turbulence_norm *= turbulence->getDoubleValue("factor", 0.0 );
 406                 turbulence->setDoubleValue( "magnitude-norm", turbulence_norm );
 407         }
 408 }
 409
 410 static void setupWind( bool setup_aloft, double dir, double speed, double gust )
 411 {
 412         setupWindBranch( "boundary", dir, speed, gust );
 413         if( setup_aloft )
 414                 setupWindBranch( "aloft", dir, speed, gust );
 415 }
 416
 417 double FGMetarCtrl::interpolate_val(double currentval, double requiredval, double dval )
 418 {
 419         if (fabs(currentval - requiredval) < dval) return requiredval;
 420         if (currentval < requiredval) return (currentval + dval);
 421         if (currentval > requiredval) return (currentval - dval);
 422         return requiredval;
 423 }
 424
 425 void
 426 FGMetarCtrl::init ()
 427 {
 428         first_update = true;
 429         wind_interpolation_required = true;
 430 }
 431
 432 void
 433 FGMetarCtrl::reinit ()
 434 {
 435         init();
 436 }
 437
 438 static inline double convert_to_360( double d )
 439 {
 440         if( d < 0.0 ) return d + 360.0;
 441         if( d >= 360.0 ) return d - 360.0;
 442         return d;
 443 }
 444
 445 static inline double convert_to_180( double d )
 446 {
 447         return d > 180.0 ? d - 360.0 : d;
 448 }
 449
 450 // Return the sea level pressure for a metar observation, in inHg.
 451 // This is different from QNH because it accounts for the current
 452 // temperature at the observation point.
 453 // metarPressure in inHg
 454 // fieldHt in ft
 455 // fieldTemp in C
 456
 457 static double reducePressureSl(double metarPressure, double fieldHt,
 458                                double fieldTemp)
 459 {
 460         double elev = fieldHt * SG_FEET_TO_METER;
 461         double fieldPressure
 462                 = FGAtmo::fieldPressure(elev, metarPressure * atmodel::inHg);
 463         double slPressure = P_layer(0, elev, fieldPressure,
 464                 fieldTemp + atmodel::freezing, atmodel::ISA::lam0);
 465         return slPressure / atmodel::inHg;
 466 }
 467
 468 void
 469 FGMetarCtrl::update(double dt)
 470 {
 471         if( dt <= 0 || !metar_valid ||!enabled)
 472                 return;
 473
 474         windModulator->update(dt);
 475         // Interpolate the current configuration closer to the actual METAR
 476
 477         bool reinit_required = false;
 478         bool layer_rebuild_required = false;
 479         double station_elevation_ft = station_elevation_n->getDoubleValue();
 480
 481         if (first_update) {
 482                 double dir = base_wind_dir_n->getDoubleValue()+magnetic_variation_n->getDoubleValue();
 483                 double speed = base_wind_speed_n->getDoubleValue();
 484                 double gust = gust_wind_speed_n->getDoubleValue();
 485                 setupWind(setup_winds_aloft, dir, speed, gust);
 486
 487                 double metarvis = min_visibility_n->getDoubleValue();
 488                 fgDefaultWeatherValue("visibility-m", metarvis);
 489
 490                 set_temp_at_altitude(temperature_n->getDoubleValue(), station_elevation_ft);
 491                 set_dewpoint_at_altitude(dewpoint_n->getDoubleValue(), station_elevation_ft);
 492
 493                 double metarpressure = pressure_n->getDoubleValue();
 494                 fgDefaultWeatherValue("pressure-sea-level-inhg",
 495                         reducePressureSl(metarpressure,
 496                         station_elevation_ft,
 497                         temperature_n->getDoubleValue()));
 498
 499                 // We haven't already loaded a METAR, so apply it immediately.
 500                 vector<SGPropertyNode_ptr> layers = clouds_n->getChildren("layer");
 501                 vector<SGPropertyNode_ptr>::const_iterator layer;
 502                 vector<SGPropertyNode_ptr>::const_iterator layers_end = layers.end();
 503
 504                 int i;
 505                 for (i = 0, layer = layers.begin(); layer != layers_end; ++layer, i++) {
 506                         SGPropertyNode *target = environment_clouds_n->getChild("layer", i, true);
 507
 508                         target->setStringValue("coverage",
 509                                         (*layer)->getStringValue("coverage", "clear"));
 510                         target->setDoubleValue("elevation-ft",
 511                                         (*layer)->getDoubleValue("elevation-ft"));
 512                         target->setDoubleValue("thickness-ft",
 513                                         (*layer)->getDoubleValue("thickness-ft"));
 514                         target->setDoubleValue("span-m", 40000.0);
 515                 }
 516
 517                 first_update = false;
 518                 reinit_required = true;
 519                 layer_rebuild_required = true;
 520
 521         } else {
 522                 if( wind_interpolation_required ) {
 523                         // Generate interpolated values between the METAR and the current
 524                         // configuration.
 525
 526                         // Pick up the METAR wind values and convert them into a vector.
 527                         double metar[2];
 528                         double metar_speed = base_wind_speed_n->getDoubleValue();
 529                         double metar_heading = base_wind_dir_n->getDoubleValue()+magnetic_variation_n->getDoubleValue();
 530
 531                         metar[0] = metar_speed * sin(metar_heading * SG_DEGREES_TO_RADIANS );
 532                         metar[1] = metar_speed * cos(metar_heading * SG_DEGREES_TO_RADIANS);
 533
 534                         // Convert the current wind values and convert them into a vector
 535                         double current[2];
 536                         double speed = boundary_wind_speed_n->getDoubleValue();
 537                         double dir_from = boundary_wind_from_heading_n->getDoubleValue();;
 538
 539                         current[0] = speed * sin(dir_from * SG_DEGREES_TO_RADIANS );
 540                         current[1] = speed * cos(dir_from * SG_DEGREES_TO_RADIANS );
 541
 542                         // Determine the maximum component-wise value that the wind can change.
 543                         // First we determine the fraction in the X and Y component, then
 544                         // factor by the maximum wind change.
 545                         double x = fabs(current[0] - metar[0]);
 546                         double y = fabs(current[1] - metar[1]);
 547
 548                         // only interpolate if we have a difference
 549                         if (x + y > 0.01 ) {
 550                                 double dx = x / (x + y);
 551                                 double dy = 1 - dx;
 552
 553                                 double maxdx = dx * MaxWindChangeKtsSec;
 554                                 double maxdy = dy * MaxWindChangeKtsSec;
 555
 556                                 // Interpolate each component separately.
 557                                 current[0] = interpolate_val(current[0], metar[0], maxdx*dt);
 558                                 current[1] = interpolate_val(current[1], metar[1], maxdy*dt);
 559
 560                                 // Now convert back to polar coordinates.
 561                                 if ((fabs(current[0]) > 0.1) || (fabs(current[1]) > 0.1)) {
 562                                         // Some real wind to convert back from. Work out the speed
 563                                         // and direction value in degrees.
 564                                         speed = sqrt((current[0] * current[0]) + (current[1] * current[1]));
 565                                         dir_from = (atan2(current[0], current[1]) * SG_RADIANS_TO_DEGREES );
 566
 567                                         // Normalize the direction.
 568                                         if (dir_from < 0.0)
 569                                                 dir_from += 360.0;
 570
 571                                         SG_LOG( SG_GENERAL, SG_DEBUG, "Wind : " << dir_from << "@" << speed);
 572                                 } else {
 573                                         // Special case where there is no wind (otherwise atan2 barfs)
 574                                         speed = 0.0;
 575                                 }
 576                                 double gust = gust_wind_speed_n->getDoubleValue();
 577                                 setupWind(setup_winds_aloft, dir_from, speed, gust);
 578                                 reinit_required = true;
 579                         } else {
 580                                 wind_interpolation_required = false;
 581                         }
 582                 } else { // if(wind_interpolation_required)
 583                         // interpolation of wind vector is finished, apply wind
 584                         // variations and gusts for the boundary layer only
 585
 586
 587                         bool wind_modulated = false;
 588
 589                         // start with the main wind direction
 590                         double wind_dir = base_wind_dir_n->getDoubleValue()+magnetic_variation_n->getDoubleValue();
 591                         double min = convert_to_180(base_wind_range_from_n->getDoubleValue()+magnetic_variation_n->getDoubleValue());
 592                         double max = convert_to_180(base_wind_range_to_n->getDoubleValue()+magnetic_variation_n->getDoubleValue());
 593                         if( max > min ) {
 594                                 // if variable winds configured, modulate the wind direction
 595                                 double f = windModulator->get_direction_offset_norm();
 596                                 wind_dir = min+(max-min)*f;
 597                                 double old = convert_to_180(boundary_wind_from_heading_n->getDoubleValue());
 598                                 wind_dir = convert_to_360(fgGetLowPass(old, wind_dir, dt ));
 599                                 wind_modulated = true;
 600                         }
 601
 602                         // start with main wind speed
 603                         double wind_speed = base_wind_speed_n->getDoubleValue();
 604                         max = gust_wind_speed_n->getDoubleValue();
 605                         if( max > wind_speed ) {
 606                                 // if gusts are configured, modulate wind magnitude
 607                                 double f = windModulator->get_magnitude_factor_norm();
 608                                 wind_speed = wind_speed+(max-wind_speed)*f;
 609                                 wind_speed = fgGetLowPass(boundary_wind_speed_n->getDoubleValue(), wind_speed, dt );
 610                                 wind_modulated = true;
 611                         }
 612                         if( wind_modulated ) {
 613                                 setupWind(false, wind_dir, wind_speed, max);
 614                                 reinit_required = true;
 615                         }
 616                 }
 617
 618                 // Now handle the visibility. We convert both visibility values
 619                 // to X-values, then interpolate from there, then back to real values.
 620                 // The length_scale is fixed to 1000m, so the visibility changes by
 621                 // by MaxVisChangePercentSec or 1000m X MaxVisChangePercentSec,
 622                 // whichever is more.
 623                 double vis = boundary_visibility_n->getDoubleValue();;
 624                 double metarvis = min_visibility_n->getDoubleValue();
 625                 if( vis != metarvis ) {
 626                         double currentxval = log(1000.0 + vis);
 627                         double metarxval = log(1000.0 + metarvis);
 628
 629                         currentxval = interpolate_val(currentxval, metarxval, MaxVisChangePercentSec*dt);
 630
 631                         // Now convert back from an X-value to a straightforward visibility.
 632                         vis = exp(currentxval) - 1000.0;
 633                         fgDefaultWeatherValue("visibility-m", vis);
 634                         reinit_required = true;
 635                 }
 636
 637                 double pressure = boundary_sea_level_pressure_n->getDoubleValue();
 638                 double metarpressure = pressure_n->getDoubleValue();
 639                 double newpressure = reducePressureSl(metarpressure,
 640                         station_elevation_ft,
 641                         temperature_n->getDoubleValue());
 642                 if( pressure != newpressure ) {
 643                         pressure = interpolate_val( pressure, newpressure, MaxPressureChangeInHgSec*dt );
 644                         fgDefaultWeatherValue("pressure-sea-level-inhg", pressure);
 645                         reinit_required = true;
 646                 }
 647
 648                 {
 649                         double temperature = boundary_sea_level_temperature_n->getDoubleValue();
 650                         double dewpoint = boundary_sea_level_dewpoint_n->getDoubleValue();
 651                         if( metar_sealevel_temperature != temperature ) {
 652                                 temperature = interpolate_val( temperature, metar_sealevel_temperature, MaxTemperatureChangeDegcSec*dt );
 653                                 set_temp_at_altitude( temperature, 0.0 );
 654                         }
 655                         if( metar_sealevel_dewpoint != dewpoint ) {
 656                                 dewpoint = interpolate_val( dewpoint, metar_sealevel_dewpoint, MaxTemperatureChangeDegcSec*dt );
 657                                 set_dewpoint_at_altitude( dewpoint, 0.0 );
 658                         }
 659                 }
 660
 661                 // Set the cloud layers by interpolating over the METAR versions.
 662                 vector<SGPropertyNode_ptr> layers = clouds_n->getChildren("layer");
 663                 vector<SGPropertyNode_ptr>::const_iterator layer;
 664                 vector<SGPropertyNode_ptr>::const_iterator layers_end = layers.end();
 665
 666                 double aircraft_alt = fgGetDouble("/position/altitude-ft");
 667                 int i;
 668
 669                 for (i = 0, layer = layers.begin(); layer != layers_end; ++layer, i++) {
 670                         SGPropertyNode *target = environment_clouds_n->getChild("layer", i, true);
 671
 672                         // In the case of clouds, we want to avoid writing if nothing has
 673                         // changed, as these properties are tied to the renderer and will
 674                         // cause the clouds to be updated, reseting the texture locations.
 675
 676                         // We don't interpolate the coverage values as no-matter how we
 677                         // do it, it will be quite a sudden change of texture. Better to
 678                         // have a single change than four or five.
 679                         const char *coverage = (*layer)->getStringValue("coverage", "clear");
 680                         SGPropertyNode *cov = target->getNode("coverage", true);
 681                         if (strcmp(cov->getStringValue(), coverage) != 0) {
 682                                 cov->setStringValue(coverage);
 683                                 layer_rebuild_required = true;
 684                         }
 685
 686                         double required_alt = (*layer)->getDoubleValue("elevation-ft");
 687                         double current_alt = target->getDoubleValue("elevation-ft");
 688                         double required_thickness = (*layer)->getDoubleValue("thickness-ft");
 689                         SGPropertyNode *thickness = target->getNode("thickness-ft", true);
 690
 691                         if (current_alt < -9000 || required_alt < -9000 ||
 692                                 fabs(aircraft_alt - required_alt) > MaxCloudInterpolationHeightFt ||
 693                                 fabs(current_alt - required_alt) > MaxCloudInterpolationDeltaFt) {
 694                                 // We don't interpolate any layers that are
 695                                 //  - too far above us to be visible
 696                                 //  - too far below us to be visible
 697                                 //  - with too large a difference to make interpolation sensible
 698                                 //  - to or from -9999 (used as a placeholder)
 699                                 //  - any values that are too high above us,
 700                                 if (current_alt != required_alt)
 701                                         target->setDoubleValue("elevation-ft", required_alt);
 702
 703                                 if (thickness->getDoubleValue() != required_thickness)
 704                                         thickness->setDoubleValue(required_thickness);
 705
 706                         } else {
 707                                 // Interpolate the other values in the usual way
 708                                 if (current_alt != required_alt) {
 709                                         current_alt = interpolate_val(current_alt, required_alt, MaxCloudAltitudeChangeFtSec*dt);
 710                                         target->setDoubleValue("elevation-ft", current_alt);
 711                                 }
 712
 713                                 double current_thickness = thickness->getDoubleValue();
 714
 715                                 if (current_thickness != required_thickness) {
 716                                         current_thickness = interpolate_val(current_thickness,
 717                                                                                                  required_thickness,
 718                                                                                                  MaxCloudThicknessChangeFtSec*dt);
 719                                         thickness->setDoubleValue(current_thickness);
 720                                 }
 721                         }
 722                 }
 723         }
 724
 725         // Force an update of the 3D clouds
 726         if( layer_rebuild_required )
 727                 fgSetInt("/environment/rebuild-layers", 1 );
 728
 729         // Reinitializing of the environment controller required
 730         if( reinit_required )
 731                 _environmentCtrl->reinit();
 732 }
 733
 734 const char * FGMetarCtrl::get_metar(void) const
 735 {
 736         return metar.c_str();
 737 }
 738
 739 static const char *coverage_string[] = { "clear", "few", "scattered", "broken", "overcast" };
 740 static const double thickness_value[] = { 0, 65, 600, 750, 1000 };
 741
 742 void FGMetarCtrl::set_metar( const char * metar_string )
 743 {
 744         int i;
 745
 746         metar = metar_string;
 747
 748         SGSharedPtr<FGMetar> m;
 749         try {
 750                 m = new FGMetar( metar_string );
 751         }
 752         catch( sg_io_exception ) {
 753                 SG_LOG( SG_GENERAL, SG_WARN, "Can't get metar: " << metar_string );
 754                 metar_valid = false;
 755                 return;
 756         }
 757
 758         wind_interpolation_required = true;
 759
 760         min_visibility_n->setDoubleValue( m->getMinVisibility().getVisibility_m() );
 761         max_visibility_n->setDoubleValue( m->getMaxVisibility().getVisibility_m() );
 762
 763         const SGMetarVisibility *dirvis = m->getDirVisibility();
 764         for (i = 0; i < 8; i++, dirvis++) {
 765                 SGPropertyNode *vis = metar_base_n->getChild("visibility", i, true);
 766                 double v = dirvis->getVisibility_m();
 767
 768                 vis->setDoubleValue("min-m", v);
 769                 vis->setDoubleValue("max-m", v);
 770         }
 771
 772         base_wind_dir_n->setIntValue( m->getWindDir() );
 773         base_wind_range_from_n->setIntValue( m->getWindRangeFrom() );
 774         base_wind_range_to_n->setIntValue( m->getWindRangeTo() );
 775         base_wind_speed_n->setDoubleValue( m->getWindSpeed_kt() );
 776         gust_wind_speed_n->setDoubleValue( m->getGustSpeed_kt() );
 777         temperature_n->setDoubleValue( m->getTemperature_C() );
 778         dewpoint_n->setDoubleValue( m->getDewpoint_C() );
 779         humidity_n->setDoubleValue( m->getRelHumidity() );
 780         pressure_n->setDoubleValue( m->getPressure_inHg() );
 781
 782
 783         // get station elevation to compute cloud base
 784         double station_elevation_ft = 0;
 785         {
 786                 // 1. check the id given in the metar
 787                 FGAirport* a = FGAirport::findByIdent(m->getId());
 788
 789                 // 2. if unknown, find closest airport with metar to current position
 790                 if( a == NULL ) {
 791                         SGGeod pos = SGGeod::fromDeg(longitude_n->getDoubleValue(), latitude_n->getDoubleValue());
 792                         a = FGAirport::findClosest(pos, 10000.0, &airportWithMetarFilter);
 793                 }
 794
 795                 // 3. otherwise use ground elevation
 796                 if( a != NULL ) {
 797                         station_elevation_ft = a->getElevation();
 798                         station_id_n->setStringValue( a->ident());
 799                 } else {
 800                         station_elevation_ft = ground_elevation_n->getDoubleValue() * SG_METER_TO_FEET;
 801                         station_id_n->setStringValue( m->getId());
 802                 }
 803         }
 804
 805         station_elevation_n->setDoubleValue( station_elevation_ft );
 806
 807         {       // calculate sea level temperature and dewpoint
 808                 FGEnvironment dummy; // instantiate a dummy so we can leech a method
 809                 dummy.set_elevation_ft( station_elevation_ft );
 810                 dummy.set_temperature_degc( temperature_n->getDoubleValue() );
 811                 dummy.set_dewpoint_degc( dewpoint_n->getDoubleValue() );
 812                 metar_sealevel_temperature = dummy.get_temperature_sea_level_degc();
 813                 metar_sealevel_dewpoint = dummy.get_dewpoint_sea_level_degc();
 814         }
 815
 816         vector<SGMetarCloud> cv = m->getClouds();
 817         vector<SGMetarCloud>::const_iterator cloud, cloud_end = cv.end();
 818
 819         int layer_cnt = environment_clouds_n->getChildren("layer").size();
 820         for (i = 0, cloud = cv.begin(); i < layer_cnt; i++) {
 821
 822
 823                 const char *coverage = "clear";
 824                 double elevation = -9999.0;
 825                 double thickness = 0.0;
 826                 const double span = 40000.0;
 827
 828                 if (cloud != cloud_end) {
 829                         int c = cloud->getCoverage();
 830                         coverage = coverage_string[c];
 831                         elevation = cloud->getAltitude_ft() + station_elevation_ft;
 832                         thickness = thickness_value[c];
 833                         ++cloud;
 834                 }
 835
 836                 SGPropertyNode *layer = clouds_n->getChild("layer", i, true );
 837
 838                 // if the coverage has changed, a rebuild of the layer is needed
 839                 if( strcmp(layer->getStringValue("coverage"), coverage ) ) {
 840                         layer->setStringValue("coverage", coverage);
 841                 }
 842                 layer->setDoubleValue("elevation-ft", elevation);
 843                 layer->setDoubleValue("thickness-ft", thickness);
 844                 layer->setDoubleValue("span-m", span);
 845         }
 846
 847         rain_n->setDoubleValue(m->getRain());
 848         hail_n->setDoubleValue(m->getHail());
 849         snow_n->setDoubleValue(m->getSnow());
 850         snow_cover_n->setBoolValue(m->getSnowCover());
 851         metar_valid = true;
 852 }
 853
 854 #if defined(ENABLE_THREADS)
 855 /**
 856  * This class represents the thread of execution responsible for
 857  * fetching the metar data.
 858  */
 859 class MetarThread : public OpenThreads::Thread {
 860 public:
 861         MetarThread( FGMetarFetcher * f ) : metar_fetcher(f) {}
 862         ~MetarThread() {}
 863
 864         /**
 865          * Fetche the metar data from the NOAA.
 866          */
 867         void run();
 868
 869 private:
 870         FGMetarFetcher * metar_fetcher;
 871 };
 872
 873 void MetarThread::run()
 874 {
 875         for( ;; ) {
 876                 string airport_id = metar_fetcher->request_queue.pop();
 877
 878                 if( airport_id.size() == 0 )
 879                         break;
 880
 881                 if( metar_fetcher->_error_count > 3 ) {
 882                         SG_LOG( SG_GENERAL, SG_WARN, "Too many erros fetching METAR, thread stopped permanently.");
 883                         break;
 884                 }
 885
 886                 metar_fetcher->fetch( airport_id );
 887         }
 888 }
 889 #endif
 890
 891 FGMetarFetcher::FGMetarFetcher() :
 892 #if defined(ENABLE_THREADS)
 893         metar_thread(NULL),
 894 #endif
 895         fetch_timer(0.0),
 896         search_timer(0.0),
 897         error_timer(0.0),
 898         _stale_count(0),
 899         _error_count(0),
 900         enabled(false)
 901 {
 902         longitude_n = fgGetNode( "/position/longitude-deg", true );
 903         latitude_n  = fgGetNode( "/position/latitude-deg", true );
 904         enable_n    = fgGetNode( "/environment/params/real-world-weather-fetch", true );
 905
 906         proxy_host_n = fgGetNode("/sim/presets/proxy/host", true);
 907         proxy_port_n = fgGetNode("/sim/presets/proxy/port", true);
 908         proxy_auth_n = fgGetNode("/sim/presets/proxy/authentication", true);
 909         max_age_n    = fgGetNode("/environment/params/metar-max-age-min", true);
 910
 911         output_n         = fgGetNode("/environment/metar/data", true );
 912 #if defined(ENABLE_THREADS)
 913         metar_thread = new MetarThread(this);
 914 // FIXME: do we really need setProcessorAffinity()?
 915 //      metar_thread->setProcessorAffinity(1);
 916         metar_thread->start();
 917 #endif // ENABLE_THREADS
 918 }
 919
 920
 921 FGMetarFetcher::~FGMetarFetcher()
 922 {
 923 #if defined(ENABLE_THREADS)
 924         request_queue.push("");
 925         metar_thread->join();
 926         delete metar_thread;
 927 #endif // ENABLE_THREADS
 928 }
 929
 930 void FGMetarFetcher::init ()
 931 {
 932         fetch_timer = 0.0;
 933         search_timer = 0.0;
 934         error_timer = 0.0;
 935         _stale_count = 0;
 936         _error_count = 0;
 937         current_airport_id.clear();
 938         /* Torsten Dreyer:
 939            hack to stop startup.nas complaining if metar arrives after nasal-dir-initialized
 940            is fired. Immediately fetch and wait for the METAR before continuing. This gets the
 941            /environment/metar/xxx properties filled before nasal-dir is initialized.
 942            Maybe the runway selection should happen here to make startup.nas obsolete?
 943         */
 944         const char * startup_airport = fgGetString("/sim/startup/options/airport");
 945         if( *startup_airport ) {
 946                 FGAirport * a = FGAirport::getByIdent( startup_airport );
 947                 if( a ) {
 948                         SGGeod pos = SGGeod::fromDeg(a->getLongitude(), a->getLatitude());
 949                         a = FGAirport::findClosest(pos, 10000.0, &airportWithMetarFilter);
 950                         current_airport_id = a->getId();
 951                         fetch( current_airport_id );
 952                 }
 953         }
 954 }
 955
 956 void FGMetarFetcher::reinit ()
 957 {
 958         init();
 959 }
 960
 961 /* search for closest airport with metar every xx seconds */
 962 static const int search_interval_sec = 60;
 963
 964 /* fetch metar for airport, even if airport has not changed every xx seconds */
 965 static const int fetch_interval_sec = 900;
 966
 967 /* reset error counter after xxx seconds */
 968 static const int error_timer_sec = 3;
 969
 970 void FGMetarFetcher::update (double delta_time_sec)
 971 {
 972         fetch_timer -= delta_time_sec;
 973         search_timer -= delta_time_sec;
 974         error_timer -= delta_time_sec;
 975
 976         if( error_timer <= 0.0 ) {
 977                 error_timer = error_timer_sec;
 978                 _error_count = 0;
 979         }
 980
 981         if( enable_n->getBoolValue() == false ) {
 982                 enabled = false;
 983                 return;
 984         }
 985
 986         // we were just enabled, reset all timers to
 987         // trigger immediate metar fetch
 988         if( !enabled ) {
 989                 search_timer = 0.0;
 990                 fetch_timer = 0.0;
 991                 error_timer = error_timer_sec;
 992                 enabled = true;
 993         }
 994
 995         FGAirport * a = NULL;
 996
 997         if( search_timer <= 0.0 ) {
 998                 // search timer expired, search closest airport with metar
 999                 SGGeod pos = SGGeod::fromDeg(longitude_n->getDoubleValue(), latitude_n->getDoubleValue());
1000                 a = FGAirport::findClosest(pos, 10000.0, &airportWithMetarFilter);
1001                 search_timer = search_interval_sec;
1002         }
1003
1004         if( a == NULL )
1005                 return;
1006
1007
1008         if( a->ident() != current_airport_id || fetch_timer <= 0 ) {
1009                 // fetch timer expired or airport has changed, schedule a fetch
1010                 current_airport_id = a->ident();
1011                 fetch_timer = fetch_interval_sec;
1012 #if defined(ENABLE_THREADS)
1013                 // push this airport id into the queue for the worker thread
1014                 request_queue.push( current_airport_id );
1015 #else
1016                 // if there is no worker thread, immediately fetch the data
1017                 fetch( current_airport_id );
1018 #endif
1019         }
1020 }
1021
1022 void FGMetarFetcher::fetch( const string & id )
1023 {
1024         if( enable_n->getBoolValue() == false )
1025                 return;
1026
1027         SGSharedPtr<FGMetar> result = NULL;
1028
1029         // fetch current metar data
1030         try {
1031                 string host = proxy_host_n->getStringValue();
1032                 string auth = proxy_auth_n->getStringValue();
1033                 string port = proxy_port_n->getStringValue();
1034
1035                 result = new FGMetar( id, host, port, auth);
1036
1037                 long max_age = max_age_n->getLongValue();
1038                 long age = result->getAge_min();
1039
1040                 if (max_age && age > max_age) {
1041                         SG_LOG( SG_GENERAL, SG_WARN, "METAR data too old (" << age << " min).");
1042                         if (++_stale_count > 10) {
1043                                 _error_count = 1000;
1044                                 throw sg_io_exception("More than 10 stale METAR messages in a row." " Check your system time!");
1045                         }
1046                 } else {
1047                         _stale_count = 0;
1048                 }
1049
1050         } catch (const sg_io_exception& e) {
1051                 SG_LOG( SG_GENERAL, SG_WARN, "Error fetching live weather data: " << e.getFormattedMessage().c_str() );
1052                 result = NULL;
1053                 // remove METAR flag from the airport
1054                 FGAirport * a = FGAirport::findByIdent( id );
1055                 if( a ) a->setMetar( false );
1056                 // immediately schedule a new search
1057                 search_timer = 0.0;
1058         }
1059
1060         // write the metar to the property node, the rest is done by the methods tied to this property
1061         // don't write the metar data, if real-weather-fetch has been disabled in the meantime
1062         if( result != NULL && enable_n->getBoolValue() == true )
1063                 output_n->setStringValue( result->getData() );
1064 }
1065
1066 // end of environment_ctrl.cxx
1067