src/Environment/environment_ctrl.cxx

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