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