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
  40
  41 using std::string;
  42
  43 #include "ridge_lift.hxx"
  44
  45 static string CreateIndexedPropertyName(string Property, int index)
  46 {
  47         std::stringstream str;
  48         str << index;
  49         string tmp;
  50         str >> tmp;
  51         return Property + "[" + tmp + "]";
  52 }
  53
  54 static inline double sign(double x) {
  55         return x == 0 ? 0 : x > 0 ? 1.0 : -1.0;
  56 }
  57
  58 static const    double BOUNDARY1_m = 40.0;
  59
  60 //constructor
  61 FGRidgeLift::FGRidgeLift ()
  62 {
  63         dist_probe_m[0] = 0.0; // in meters
  64         dist_probe_m[1] = 250.0;
  65         dist_probe_m[2] = 750.0;
  66         dist_probe_m[3] = 2000.0;
  67         dist_probe_m[4] = -100.0;
  68
  69         strength = 0.0;
  70         timer = 0.0;
  71 }
  72
  73 //destructor
  74 FGRidgeLift::~FGRidgeLift()
  75 {
  76
  77 }
  78
  79 void FGRidgeLift::init(void)
  80 {
  81         _ridge_lift_fps_node = fgGetNode("/environment/ridge-lift-fps", true);
  82         _surface_wind_from_deg_node =
  83                         fgGetNode("/environment/config/boundary/entry[0]/wind-from-heading-deg"
  84                         , true);
  85         _surface_wind_speed_node =
  86                         fgGetNode("/environment/config/boundary/entry[0]/wind-speed-kt"
  87                         , true);
  88         _earth_radius_node = fgGetNode("/position/sea-level-radius-ft", true);
  89         _user_longitude_node = fgGetNode("/position/longitude-deg", true);
  90         _user_latitude_node = fgGetNode("/position/latitude-deg", true);
  91         _user_altitude_ft_node = fgGetNode("/position/altitude-ft", true);
  92         _user_altitude_agl_ft_node = fgGetNode("/position/altitude-agl-ft", true);
  93 }
  94
  95 void FGRidgeLift::bind() {
  96         string prop;
  97
  98         for( int i = 0; i < 5; i++ ) {
  99                 prop = CreateIndexedPropertyName("/environment/ridge-lift/probe-elev-m", i );
 100                 fgTie( prop.c_str(), this, i, &FGRidgeLift::get_probe_elev_m); // read-only
 101
 102                 prop = CreateIndexedPropertyName("/environment/ridge-lift/probe-lat-deg", i );
 103                 fgTie( prop.c_str(), this, i, &FGRidgeLift::get_probe_lat_deg); // read-only
 104
 105                 prop = CreateIndexedPropertyName("/environment/ridge-lift/probe-lon-deg", i );
 106                 fgTie( prop.c_str(), this, i, &FGRidgeLift::get_probe_lon_deg); // read-only
 107         }
 108
 109         for( int i = 0; i < 4; i++ ) {
 110                 prop = CreateIndexedPropertyName("/environment/ridge-lift/slope", i );
 111                 fgTie( prop.c_str(), this, i, &FGRidgeLift::get_slope); // read-only
 112         }
 113 }
 114
 115 void FGRidgeLift::unbind() {
 116         string prop;
 117
 118         for( int i = 0; i < 5; i++ ) {
 119
 120                 prop = CreateIndexedPropertyName("/environment/ridge-lift/probe-elev-m", i );
 121                 fgUntie( prop.c_str() );
 122
 123                 prop = CreateIndexedPropertyName("/environment/ridge-lift/probe-lat-deg", i );
 124                 fgUntie( prop.c_str() );
 125
 126                 prop = CreateIndexedPropertyName("/environment/ridge-lift/probe-lon-deg", i );
 127                 fgUntie( prop.c_str() );
 128         }
 129
 130         for( int i = 0; i < 4; i++ ) {
 131                 prop = CreateIndexedPropertyName("/environment/ridge-lift/slope", i );
 132                 fgUntie( prop.c_str() );
 133         }
 134 }
 135
 136 void FGRidgeLift::update(double dt) {
 137
 138         //get the windspeed at ground level
 139
 140         double ground_wind_from_rad = _surface_wind_from_deg_node->getDoubleValue() * SG_DEGREES_TO_RADIANS;
 141         double ground_wind_speed_mps = _surface_wind_speed_node->getDoubleValue() * SG_NM_TO_METER / 3600;
 142
 143         timer -= dt;
 144         if (timer <= 0.0 ) {
 145                 // copy values
 146
 147                 double user_latitude_rad = _user_latitude_node->getDoubleValue() * SG_DEGREES_TO_RADIANS;
 148                 double user_longitude_rad = _user_longitude_node->getDoubleValue() * SG_DEGREES_TO_RADIANS;
 149
 150                 double earth_rad_m = _earth_radius_node->getDoubleValue() * SG_FEET_TO_METER;
 151                 if( earth_rad_m < SG_EPSILON )
 152                                 earth_rad_m = SG_EARTH_RAD * 1000;
 153
 154                 // Placing the probes
 155
 156                 for (int i = 0; i < sizeof(probe_lat_rad)/sizeof(probe_lat_rad[0]); i++) {
 157                         double probe_radius_ratio = dist_probe_m[i]/earth_rad_m;
 158
 159                         probe_lat_rad[i] = asin(sin(user_latitude_rad)*cos(probe_radius_ratio)
 160                                         +cos(user_latitude_rad)*sin(probe_radius_ratio)*cos(ground_wind_from_rad));
 161                         if (probe_lat_rad[i] < SG_EPSILON ) {
 162                                 probe_lon_rad[i] = user_latitude_rad; // probe on a pole
 163                         } else {
 164                                 probe_lon_rad[i] = fmod((user_longitude_rad+asin(sin(ground_wind_from_rad)
 165                                                         *sin(probe_radius_ratio)/cos(probe_lat_rad[i]))+SG_PI)
 166                                                 ,SGD_2PI)-SG_PI;
 167                         }
 168                         probe_lat_deg[i]= probe_lat_rad[i] * SG_RADIANS_TO_DEGREES;
 169                         probe_lon_deg[i]= probe_lon_rad[i] * SG_RADIANS_TO_DEGREES;
 170                 }
 171
 172                 for (int i = 0; i < sizeof(probe_elev_m)/sizeof(probe_elev_m[0]); i++) {
 173                         double elevation = 0;
 174                         if (globals->get_scenery()->get_elevation_m(SGGeod::fromGeodM(
 175                                 SGGeod::fromRad(probe_lon_rad[i],probe_lat_rad[i]), 20000), elevation, 0)) {
 176                                 if ( elevation > 0.1 ) {
 177                                         probe_elev_m[i] = elevation;
 178                                 } else {
 179                                         probe_elev_m[i] = 0.1 ;
 180                                 }
 181                         } else {
 182                                 probe_elev_m[i] = 0.1;
 183                         }
 184                 }
 185
 186                 // slopes
 187                 double adj_slope[4];
 188                 slope[0] = (probe_elev_m[0] - probe_elev_m[1]) / dist_probe_m[1];
 189                 slope[1] = (probe_elev_m[1] - probe_elev_m[2]) / dist_probe_m[2];
 190                 slope[2] = (probe_elev_m[2] - probe_elev_m[3]) / dist_probe_m[3];
 191                 slope[3] = (probe_elev_m[4] - probe_elev_m[0]) / -dist_probe_m[4];
 192
 193                 for (int i = 0; i < sizeof(slope)/sizeof(slope[0]); i++)
 194                         adj_slope[i] = sin(atan(5.0 * pow ( (fabs(slope[i])),1.7) ) ) *sign(slope[i]);
 195
 196                 //adjustment
 197                 adj_slope[0] = 0.2 * adj_slope[0];
 198                 adj_slope[1] = 0.2 * adj_slope[1];
 199                 if ( adj_slope [2] < 0.0 ) {
 200                         adj_slope[2] = 0.5 * adj_slope[2];
 201                 } else {
 202                         adj_slope[2] = 0.0 ;
 203                 }
 204
 205                 if ( ( adj_slope [0] >= 0.0 ) && ( adj_slope [3] < 0.0 ) ) {
 206                         adj_slope[3] = 0.0;
 207                 } else {
 208                         adj_slope[3] = 0.2 * adj_slope[3];
 209                 }
 210                 lift_factor = adj_slope[0]+adj_slope[1]+adj_slope[2]+adj_slope[3];
 211
 212                 // restart the timer
 213                 timer = 1.0;
 214         }
 215
 216         //user altitude above ground
 217         double user_altitude_agl_m = _user_altitude_agl_ft_node->getDoubleValue() * SG_FEET_TO_METER;
 218
 219         //boundaries
 220         double boundary2_m = 130.0; // in the lift
 221         if (lift_factor < 0.0) { // in the sink
 222                 double highest_probe_temp= max ( probe_elev_m[1], probe_elev_m[2] );
 223                 double highest_probe_downwind_m= max ( highest_probe_temp, probe_elev_m[3] );
 224                 boundary2_m = highest_probe_downwind_m - probe_elev_m[0];
 225         }
 226
 227         double agl_factor;
 228         if ( user_altitude_agl_m < BOUNDARY1_m ) {
 229                 agl_factor = 0.5+0.5*user_altitude_agl_m /BOUNDARY1_m ;
 230         } else if ( user_altitude_agl_m < boundary2_m ) {
 231                 agl_factor = 1.0;
 232         } else {
 233                 agl_factor = exp(-(2 + probe_elev_m[0] / 2000) *
 234                                 (user_altitude_agl_m - boundary2_m) / max(probe_elev_m[0],200.0));
 235         }
 236
 237         double lift_mps = lift_factor* ground_wind_speed_mps * agl_factor;
 238
 239         //the updraft, finally, in ft per second
 240         strength = fgGetLowPass( strength, lift_mps * SG_METER_TO_FEET, dt );
 241 //      if(isnan(strength)) strength=0;
 242          _ridge_lift_fps_node->setDoubleValue( strength );
 243 }