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