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 <Scenery/scenery.hxx>
  36 #include <string>
  37 #include <math.h>
  38
  39 using std::string;
  40
  41 #include "ridge_lift.hxx"
  42
  43 //constructor
  44 FGRidgeLift::FGRidgeLift ()
  45 {
  46         dist_probe_m[0] = 0.0; // in meters
  47         dist_probe_m[1] = 250.0;
  48         dist_probe_m[2] = 750.0;
  49         dist_probe_m[3] = 2000.0;
  50         dist_probe_m[4] = -100.0;
  51
  52         BOUNDARY1_m = 40.0; // in meters
  53         BOUNDARY2_m = 130.0;
  54
  55         PI = 4.0*atan(1.0); // pi
  56         deg2rad = (PI/180.0);
  57         rad2deg = (180.0/PI);
  58         strength = 0.0;
  59         timer = 0.0;
  60         scanned = false;
  61
  62         earth_rad_ft=20899773.07;
  63         for( int i = 0; i < sizeof(probe_elev_m)/sizeof(probe_elev_m[0]); i++ )
  64           probe_elev_m[i] = 0.0;
  65 }
  66
  67 //destructor
  68 FGRidgeLift::~FGRidgeLift()
  69 {
  70
  71 }
  72
  73 void FGRidgeLift::init(void)
  74 {
  75         _ridge_lift_fps_node =
  76                         fgGetNode("/environment/ridge-lift-fps"
  77                         , true);
  78         _surface_wind_from_deg_node =
  79                         fgGetNode("/environment/config/boundary/entry[0]/wind-from-heading-deg"
  80                         , true);
  81         _surface_wind_speed_node =
  82                         fgGetNode("/environment/config/boundary/entry[0]/wind-speed-kt"
  83                         , true);
  84         _earth_radius_node =
  85                         fgGetNode("/position/sea-level-radius-ft"
  86                         , true);
  87         _user_longitude_node =
  88                         fgGetNode("/position/longitude-deg"
  89                         , true);
  90         _user_latitude_node =
  91                         fgGetNode("/position/latitude-deg"
  92                         , true);
  93         _user_altitude_ft_node =
  94                         fgGetNode("/position/altitude-ft"
  95                         , true);
  96         _user_altitude_agl_ft_node =
  97                         fgGetNode("/position/altitude-agl-ft"
  98                         , true);
  99 }
 100
 101 void FGRidgeLift::bind() {
 102
 103         fgTie("/environment/ridge-lift/probe-elev-m[0]", this,
 104                 &FGRidgeLift::get_probe_elev_m_0); // read-only
 105         fgTie("/environment/ridge-lift/probe-elev-m[1]", this,
 106                 &FGRidgeLift::get_probe_elev_m_1); // read-only
 107         fgTie("/environment/ridge-lift/probe-elev-m[2]", this,
 108                 &FGRidgeLift::get_probe_elev_m_2); // read-only
 109         fgTie("/environment/ridge-lift/probe-elev-m[3]", this,
 110                 &FGRidgeLift::get_probe_elev_m_3); // read-only
 111         fgTie("/environment/ridge-lift/probe-elev-m[4]", this,
 112                 &FGRidgeLift::get_probe_elev_m_4); // read-only
 113
 114         fgTie("/environment/ridge-lift/probe-lat-deg[0]", this,
 115                 &FGRidgeLift::get_probe_lat_0); // read-only
 116         fgTie("/environment/ridge-lift/probe-lat-deg[1]", this,
 117                 &FGRidgeLift::get_probe_lat_1); // read-only
 118         fgTie("/environment/ridge-lift/probe-lat-deg[2]", this,
 119                 &FGRidgeLift::get_probe_lat_2); // read-only
 120         fgTie("/environment/ridge-lift/probe-lat-deg[3]", this,
 121                 &FGRidgeLift::get_probe_lat_3); // read-only
 122         fgTie("/environment/ridge-lift/probe-lat-deg[4]", this,
 123                 &FGRidgeLift::get_probe_lat_4); // read-only
 124
 125         fgTie("/environment/ridge-lift/probe-lon-deg[0]", this,
 126                 &FGRidgeLift::get_probe_lon_0); // read-only
 127         fgTie("/environment/ridge-lift/probe-lon-deg[1]", this,
 128                 &FGRidgeLift::get_probe_lon_1); // read-only
 129         fgTie("/environment/ridge-lift/probe-lon-deg[2]", this,
 130                 &FGRidgeLift::get_probe_lon_2); // read-only
 131         fgTie("/environment/ridge-lift/probe-lon-deg[3]", this,
 132                 &FGRidgeLift::get_probe_lon_3); // read-only
 133         fgTie("/environment/ridge-lift/probe-lon-deg[4]", this,
 134                 &FGRidgeLift::get_probe_lon_4); // read-only
 135
 136         fgTie("/environment/ridge-lift/slope[0]", this,
 137                 &FGRidgeLift::get_slope_0); // read-only
 138         fgTie("/environment/ridge-lift/slope[1]", this,
 139                 &FGRidgeLift::get_slope_1); // read-only
 140         fgTie("/environment/ridge-lift/slope[2]", this,
 141                 &FGRidgeLift::get_slope_2); // read-only
 142         fgTie("/environment/ridge-lift/slope[3]", this,
 143                 &FGRidgeLift::get_slope_3); // read-only
 144
 145 }
 146
 147 void FGRidgeLift::unbind() {
 148
 149         fgUntie("/environment/ridge-lift/probe-elev-m[0]");
 150         fgUntie("/environment/ridge-lift/probe-elev-m[1]");
 151         fgUntie("/environment/ridge-lift/probe-elev-m[2]");
 152         fgUntie("/environment/ridge-lift/probe-elev-m[3]");
 153         fgUntie("/environment/ridge-lift/probe-elev-m[4]");
 154
 155         fgUntie("/environment/ridge-lift/probe-lat-deg[0]");
 156         fgUntie("/environment/ridge-lift/probe-lat-deg[1]");
 157         fgUntie("/environment/ridge-lift/probe-lat-deg[2]");
 158         fgUntie("/environment/ridge-lift/probe-lat-deg[3]");
 159         fgUntie("/environment/ridge-lift/probe-lat-deg[4]");
 160
 161         fgUntie("/environment/ridge-lift/probe-lon-deg[0]");
 162         fgUntie("/environment/ridge-lift/probe-lon-deg[1]");
 163         fgUntie("/environment/ridge-lift/probe-lon-deg[2]");
 164         fgUntie("/environment/ridge-lift/probe-lon-deg[3]");
 165         fgUntie("/environment/ridge-lift/probe-lon-deg[4]");
 166
 167         fgUntie("/environment/ridge-lift/slope[0]");
 168         fgUntie("/environment/ridge-lift/slope[1]");
 169         fgUntie("/environment/ridge-lift/slope[2]");
 170         fgUntie("/environment/ridge-lift/slope[3]");
 171
 172 }
 173
 174 void FGRidgeLift::update(double dt) {
 175         Run(dt);
 176 }
 177
 178 double FGRidgeLift::sign(double x) {
 179     if (x == 0.0)
 180         return x;
 181     else
 182         return x/fabs(x);
 183 }
 184
 185 void FGRidgeLift::Run(double dt) {
 186
 187         // copy values
 188
 189         user_latitude_deg  = _user_latitude_node->getDoubleValue();
 190         user_longitude_deg = _user_longitude_node->getDoubleValue();
 191         //user_altitude_ft  = _user_altitude_ft_node->getDoubleValue();
 192
 193         if ( ( _earth_radius_node->getDoubleValue() ) > 1.0 ) {
 194         earth_rad_ft =_earth_radius_node->getDoubleValue(); }
 195         else {  earth_rad_ft=20899773.07; }
 196
 197         //earth_rad_m = earth_rad_ft * 0.3048 ;
 198         earth_rad_m = earth_rad_ft * SG_FEET_TO_METER ;
 199
 200         //get the windspeed at ground level
 201
 202         double ground_wind_from_deg = _surface_wind_from_deg_node->getDoubleValue();
 203         double ground_wind_speed_kts  = _surface_wind_speed_node->getDoubleValue();
 204         //double ground_wind_speed_mps = ground_wind_speed_kts / SG_METER_TO_FEET;
 205         double ground_wind_speed_mps = ground_wind_speed_kts / 3.2808399;
 206
 207         //double ground_wind_from_rad = (user_longitude_deg < 0.0)  ?
 208         //      PI*( ground_wind_from_deg/180.0) +PI : PI*( ground_wind_from_deg/180.0);
 209         double ground_wind_from_rad = PI*( ground_wind_from_deg/180.0);
 210
 211         // Placing the probes
 212
 213         for (int i = 0; i <= 4; i++)
 214         {
 215                         probe_lat_rad[i] = asin(sin(deg2rad*user_latitude_deg)*cos(dist_probe_m[i]/earth_rad_m)
 216                                 +cos(deg2rad*user_latitude_deg)*sin(dist_probe_m[i]/earth_rad_m)*cos(ground_wind_from_rad));
 217                 if (probe_lat_rad[i] == 0.0) {
 218                         probe_lon_rad[i] = (deg2rad*user_latitude_deg);      // probe on a pole
 219                 }
 220                 else {
 221                         probe_lon_rad[i] = fmod((deg2rad*user_longitude_deg+asin(sin(ground_wind_from_rad)
 222                                                 *sin(dist_probe_m[i]/earth_rad_m)/cos(probe_lat_rad[i]))+PI)
 223                                         ,(2.0*PI))-PI;
 224                 }
 225                 probe_lat_deg[i]= rad2deg*probe_lat_rad[i];
 226                 probe_lon_deg[i]= rad2deg*probe_lon_rad[i];
 227         }
 228
 229         // ground elevations
 230         // every second
 231
 232         timer -= dt;
 233         if (timer <= 0.0 ) {
 234                 for (int i = 0; i <= 4; i++)
 235                 {
 236                         if (globals->get_scenery()->get_elevation_m(SGGeod::fromGeodM(
 237                             SGGeod::fromRad(probe_lon_rad[i],probe_lat_rad[i]), 20000), alt, 0))
 238                         {
 239                                 if ( alt > 0.1 ) { probe_elev_m[i] =  alt; } else { probe_elev_m[i] = 0.1 ;};
 240                         }
 241                         else { probe_elev_m[i] = 0.1;};
 242                 }
 243                 timer = 1.0;
 244
 245         }
 246
 247         // slopes
 248
 249         double adj_slope[5];
 250
 251         slope[0] = (probe_elev_m[0] - probe_elev_m[1]) / dist_probe_m[1];
 252         slope[1] = (probe_elev_m[1] - probe_elev_m[2]) / dist_probe_m[2];
 253         slope[2] = (probe_elev_m[2] - probe_elev_m[3]) / dist_probe_m[3];
 254         slope[3] = (probe_elev_m[4] - probe_elev_m[0]) / -dist_probe_m[4];
 255
 256         for (int i = 0; i <= 4; i++)
 257         {
 258                 adj_slope[i] = sin(atan(5.0 * pow ( (fabs(slope[i])),1.7) ) ) *sign(slope[i]);
 259         }
 260
 261         //adjustment
 262
 263         adj_slope[0] = 0.2 * adj_slope[0];
 264         adj_slope[1] = 0.2 * adj_slope[1];
 265         if ( adj_slope [2] < 0.0 )
 266         {
 267                 adj_slope[2] = 0.5 * adj_slope[2];
 268         }
 269         else
 270         {
 271                 adj_slope[2] = 0.0 ;
 272         }
 273
 274         if ( ( adj_slope [0] >= 0.0 ) && ( adj_slope [3] < 0.0 ) )
 275         {
 276                 adj_slope[3] = 0.0;
 277         }
 278         else
 279         {
 280                 adj_slope[3] = 0.2 * adj_slope[3];
 281         }
 282
 283         double lift_factor = adj_slope[0]+adj_slope[1]+adj_slope[2]+adj_slope[3];
 284
 285         //user altitude above ground
 286
 287         user_altitude_agl_ft = _user_altitude_agl_ft_node->getDoubleValue();
 288         user_altitude_agl_m = ( user_altitude_agl_ft / SG_METER_TO_FEET );
 289
 290         //boundaries
 291         double agl_factor;
 292
 293         if (lift_factor < 0.0) // in the sink
 294         {
 295                 double highest_probe_temp= max ( probe_elev_m[1], probe_elev_m[2] );
 296                 double highest_probe_downwind_m= max ( highest_probe_temp, probe_elev_m[3] );
 297                 BOUNDARY2_m = highest_probe_downwind_m - probe_elev_m[0];
 298         }
 299         else // in the lift
 300         {
 301                 BOUNDARY2_m = 130.0;
 302         }
 303
 304         if ( user_altitude_agl_m < BOUNDARY1_m )
 305         {
 306                 agl_factor = 0.5+0.5*user_altitude_agl_m /BOUNDARY1_m ;
 307         }
 308         else if ( user_altitude_agl_m < BOUNDARY2_m )
 309         {
 310                 agl_factor = 1.0;
 311         }
 312         else
 313         {
 314                 agl_factor =  exp(-(2 + 2 * probe_elev_m[0] / 4000) *
 315                                 (user_altitude_agl_m - BOUNDARY2_m) / max(probe_elev_m[0],200.0));
 316         }
 317
 318         double lift_mps = lift_factor* ground_wind_speed_mps * agl_factor;
 319
 320         //the updraft, finally, in ft per second
 321         strength = lift_mps * SG_METER_TO_FEET ;
 322 //      if(isnan(strength)) strength=0;
 323          _ridge_lift_fps_node->setDoubleValue( strength );
 324 }
 325
 326
 327
 328
 329
 330