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