src/Environment/ridge_lift.cxx

   1 // simulates ridge lift
   2 //
   3 // Written by Patrice Poly
   4 // Copyright (C) 2009 Patrice Poly - p.polypa@gmail.com
   5 //
   6 //
   7 // Entirely based  on the paper :
   8 // http://carrier.csi.cam.ac.uk/forsterlewis/soaring/sim/fsx/dev/sim_probe/sim_probe_paper.html
   9 // by Ian Forster-Lewis, University of Cambridge, 26th December 2007
  10 //
  11 //
  12 // This program is free software; you can redistribute it and/or
  13 // modify it under the terms of the GNU General Public License as
  14 // published by the Free Software Foundation; either version 2 of the
  15 // License, or (at your option) any later version.
  16 //
  17 // This program is distributed in the hope that it will be useful, but
  18 // WITHOUT ANY WARRANTY; without even the implied warranty of
  19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20 // General Public License for more details.
  21 //
  22 // You should have received a copy of the GNU General Public License
  23 // along with this program; if not, write to the Free Software
  24 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  25 //
  26 //
  27
  28
  29 #ifdef HAVE_CONFIG_H
  30 #  include <config.h>
  31 #endif
  32
  33 #include <Main/fg_props.hxx>
  34 #include <Main/globals.hxx>
  35 #include <Main/util.hxx>
  36 #include <Scenery/scenery.hxx>
  37 #include <string>
  38 #include <math.h>
  39 #include <simgear/sg_inlines.h>
  40
  41 using std::string;
  42
  43 #include "ridge_lift.hxx"
  44
  45 static const double BOUNDARY1_m = 40.0;
  46
  47 const double FGRidgeLift::dist_probe_m[] = { // in meters
  48      0.0,
  49    250.0,
  50    750.0,
  51   2000.0,
  52   -100.0
  53 };
  54
  55 //constructor
  56 FGRidgeLift::FGRidgeLift () :
  57   lift_factor(0.0)
  58 {
  59         strength = 0.0;
  60         timer = 0.0;
  61         for( int i = 0; i < 5; i++ )
  62                 probe_elev_m[i] = probe_lat_deg[i] = probe_lon_deg[i] = 0.0;
  63
  64         for( int i = 0; i < 4; i++ )
  65                 slope[i] = 0.0;
  66 }
  67
  68 //destructor
  69 FGRidgeLift::~FGRidgeLift()
  70 {
  71
  72 }
  73
  74 void FGRidgeLift::init(void)
  75 {
  76         _enabled_node = fgGetNode( "/environment/ridge-lift/enabled", false );
  77
  78         _ridge_lift_fps_node = fgGetNode("/environment/ridge-lift-fps", true);
  79         _surface_wind_from_deg_node =
  80                         fgGetNode("/environment/config/boundary/entry[0]/wind-from-heading-deg"
  81                         , true);
  82         _surface_wind_speed_node =
  83                         fgGetNode("/environment/config/boundary/entry[0]/wind-speed-kt"
  84                         , true);
  85         _user_longitude_node = fgGetNode("/position/longitude-deg", true);
  86         _user_latitude_node = fgGetNode("/position/latitude-deg", true);
  87         _user_altitude_agl_ft_node = fgGetNode("/position/altitude-agl-ft", true);
  88         _ground_elev_node = fgGetNode("/position/ground-elev-ft", true );
  89 }
  90
  91 void FGRidgeLift::bind() {
  92         string prop;
  93
  94         _tiedProperties.setRoot( fgGetNode("/environment/ridge-lift",true));
  95         for( int i = 0; i < 5; i++ ) {
  96                 _tiedProperties.Tie( "probe-elev-m", i, this, i, &FGRidgeLift::get_probe_elev_m );
  97                 _tiedProperties.Tie( "probe-lat-deg", i, this, i, &FGRidgeLift::get_probe_lat_deg );
  98                 _tiedProperties.Tie( "probe-lon-deg", i, this, i, &FGRidgeLift::get_probe_lon_deg );
  99         }
 100
 101         for( int i = 0; i < 4; i++ ) {
 102                 _tiedProperties.Tie( "slope", i, this, i, &FGRidgeLift::get_slope );
 103         }
 104 }
 105
 106 void FGRidgeLift::unbind() {
 107         _tiedProperties.Untie();
 108 }
 109
 110 void FGRidgeLift::update(double dt) {
 111
 112         if( dt <= SGLimitsd::min() ) // paused, do nothing but keep current lift
 113                 return;
 114
 115         if( _enabled_node && false == _enabled_node->getBoolValue() ) {
 116                 // do nothing if lift has been zeroed
 117                 if( strength != 0.0 ) {
 118                         if( strength > 0.1 ) {
 119                                 // slowly fade out strong lifts
 120                                 strength = fgGetLowPass( strength, 0, dt );
 121                         } else {
 122                                 strength = 0.0;
 123                         }
 124                         _ridge_lift_fps_node->setDoubleValue( strength );
 125                 }
 126                 return;
 127         }
 128
 129         timer -= dt;
 130         if (timer <= 0.0 ) {
 131
 132                 // probe0 is current position
 133                 probe_lat_deg[0] = _user_latitude_node->getDoubleValue();
 134                 probe_lon_deg[0] = _user_longitude_node->getDoubleValue();
 135                 probe_elev_m[0]  = _ground_elev_node->getDoubleValue() * SG_FEET_TO_METER;
 136
 137                 // position is geodetic, need geocentric for advanceRadM
 138                 SGGeod myGeodPos = SGGeod::fromDegM( probe_lon_deg[0], probe_lat_deg[0], 20000.0 );
 139                 SGGeoc myGeocPos = SGGeoc::fromGeod( myGeodPos );
 140                 double ground_wind_from_rad = _surface_wind_from_deg_node->getDoubleValue() * SG_DEGREES_TO_RADIANS;
 141
 142                 // compute the remaining probes
 143                 for (unsigned i = 1; i < sizeof(probe_elev_m)/sizeof(probe_elev_m[0]); i++) {
 144                         SGGeoc probe = myGeocPos.advanceRadM( ground_wind_from_rad, dist_probe_m[i] );
 145                         // convert to geodetic position for ground level computation
 146                         SGGeod probeGeod = SGGeod::fromGeoc( probe );
 147                         probe_lat_deg[i] = probeGeod.getLatitudeDeg();
 148                         probe_lon_deg[i] = probeGeod.getLongitudeDeg();
 149                         if (!globals->get_scenery()->get_elevation_m( probeGeod, probe_elev_m[i], NULL )) {
 150                                 // no ground found? use elevation of previous probe :-(
 151                                 probe_elev_m[i] = probe_elev_m[i-1];
 152                         }
 153                 }
 154
 155                 // slopes
 156                 double adj_slope[sizeof(slope)];
 157                 slope[0] = (probe_elev_m[0] - probe_elev_m[1]) / dist_probe_m[1];
 158                 slope[1] = (probe_elev_m[1] - probe_elev_m[2]) / dist_probe_m[2];
 159                 slope[2] = (probe_elev_m[2] - probe_elev_m[3]) / dist_probe_m[3];
 160                 slope[3] = (probe_elev_m[4] - probe_elev_m[0]) / -dist_probe_m[4];
 161
 162                 for (unsigned i = 0; i < sizeof(slope)/sizeof(slope[0]); i++)
 163                         adj_slope[i] = sin(atan(5.0 * pow ( (fabs(slope[i])),1.7) ) ) *SG_SIGN<double>(slope[i]);
 164
 165                 //adjustment
 166                 adj_slope[0] *= 0.2;
 167                 adj_slope[1] *= 0.2;
 168                 if ( adj_slope [2] < 0.0 ) {
 169                         adj_slope[2] *= 0.5;
 170                 } else {
 171                         adj_slope[2] = 0.0 ;
 172                 }
 173
 174                 if ( ( adj_slope [0] >= 0.0 ) && ( adj_slope [3] < 0.0 ) ) {
 175                         adj_slope[3] = 0.0;
 176                 } else {
 177                         adj_slope[3] *= 0.2;
 178                 }
 179                 lift_factor = adj_slope[0]+adj_slope[1]+adj_slope[2]+adj_slope[3];
 180
 181                 // restart the timer
 182                 timer = 1.0;
 183         }
 184
 185         //user altitude above ground
 186         double user_altitude_agl_m = _user_altitude_agl_ft_node->getDoubleValue() * SG_FEET_TO_METER;
 187
 188         //boundaries
 189         double boundary2_m = 130.0; // in the lift
 190         if (lift_factor < 0.0) { // in the sink
 191                 double highest_probe_temp= std::max ( probe_elev_m[1], probe_elev_m[2] );
 192                 double highest_probe_downwind_m= std::max ( highest_probe_temp, probe_elev_m[3] );
 193                 boundary2_m = highest_probe_downwind_m - probe_elev_m[0];
 194         }
 195
 196         double agl_factor;
 197         if ( user_altitude_agl_m < BOUNDARY1_m ) {
 198                 agl_factor = 0.5+0.5*user_altitude_agl_m /BOUNDARY1_m ;
 199         } else if ( user_altitude_agl_m < boundary2_m ) {
 200                 agl_factor = 1.0;
 201         } else {
 202                 agl_factor = exp(-(2 + probe_elev_m[0] / 2000) *
 203                      (user_altitude_agl_m - boundary2_m) / std::max(probe_elev_m[0],200.0));
 204         }
 205
 206         double ground_wind_speed_mps = _surface_wind_speed_node->getDoubleValue() * SG_NM_TO_METER / 3600;
 207         double lift_mps = lift_factor* ground_wind_speed_mps * agl_factor;
 208
 209         //the updraft, finally, in ft per second
 210         strength = fgGetLowPass( strength, lift_mps * SG_METER_TO_FEET, dt );
 211         _ridge_lift_fps_node->setDoubleValue( strength );
 212 }