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