1 // simulates ridge lift
3 // Written by Patrice Poly
4 // Copyright (C) 2009 Patrice Poly - p.polypa@gmail.com
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
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.
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.
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.
33 #include <Main/fg_props.hxx>
34 #include <Main/globals.hxx>
35 #include <Main/util.hxx>
36 #include <Scenery/scenery.hxx>
43 #include "ridge_lift.hxx"
45 static string CreateIndexedPropertyName(string Property, int index)
47 std::stringstream str;
51 return Property + "[" + tmp + "]";
54 static const double BOUNDARY1_m = 40.0;
57 FGRidgeLift::FGRidgeLift ()
59 dist_probe_m[0] = 0.0; // in meters
60 dist_probe_m[1] = 250.0;
61 dist_probe_m[2] = 750.0;
62 dist_probe_m[3] = 2000.0;
63 dist_probe_m[4] = -100.0;
70 FGRidgeLift::~FGRidgeLift()
75 void FGRidgeLift::init(void)
77 _ridge_lift_fps_node = fgGetNode("/environment/ridge-lift-fps", true);
78 _surface_wind_from_deg_node =
79 fgGetNode("/environment/config/boundary/entry[0]/wind-from-heading-deg"
81 _surface_wind_speed_node =
82 fgGetNode("/environment/config/boundary/entry[0]/wind-speed-kt"
84 _earth_radius_node = fgGetNode("/position/sea-level-radius-ft", true);
85 _user_longitude_node = fgGetNode("/position/longitude-deg", true);
86 _user_latitude_node = fgGetNode("/position/latitude-deg", true);
87 _user_altitude_ft_node = fgGetNode("/position/altitude-ft", true);
88 _user_altitude_agl_ft_node = fgGetNode("/position/altitude-agl-ft", true);
91 void FGRidgeLift::bind() {
94 for( int i = 0; i < 5; i++ ) {
95 prop = CreateIndexedPropertyName("/environment/ridge-lift/probe-elev-m", i );
96 fgTie( prop.c_str(), this, i, &FGRidgeLift::get_probe_elev_m); // read-only
98 prop = CreateIndexedPropertyName("/environment/ridge-lift/probe-lat-deg", i );
99 fgTie( prop.c_str(), this, i, &FGRidgeLift::get_probe_lat_deg); // read-only
101 prop = CreateIndexedPropertyName("/environment/ridge-lift/probe-lon-deg", i );
102 fgTie( prop.c_str(), this, i, &FGRidgeLift::get_probe_lon_deg); // read-only
105 for( int i = 0; i < 4; i++ ) {
106 prop = CreateIndexedPropertyName("/environment/ridge-lift/slope", i );
107 fgTie( prop.c_str(), this, i, &FGRidgeLift::get_slope); // read-only
111 void FGRidgeLift::unbind() {
114 for( int i = 0; i < 5; i++ ) {
116 prop = CreateIndexedPropertyName("/environment/ridge-lift/probe-elev-m", i );
117 fgUntie( prop.c_str() );
119 prop = CreateIndexedPropertyName("/environment/ridge-lift/probe-lat-deg", i );
120 fgUntie( prop.c_str() );
122 prop = CreateIndexedPropertyName("/environment/ridge-lift/probe-lon-deg", i );
123 fgUntie( prop.c_str() );
126 for( int i = 0; i < 4; i++ ) {
127 prop = CreateIndexedPropertyName("/environment/ridge-lift/slope", i );
128 fgUntie( prop.c_str() );
132 void FGRidgeLift::update(double dt) {
134 //get the windspeed at ground level
136 double ground_wind_from_rad = _surface_wind_from_deg_node->getDoubleValue() * SG_DEGREES_TO_RADIANS;
137 double ground_wind_speed_mps = _surface_wind_speed_node->getDoubleValue() * SG_NM_TO_METER / 3600;
143 double user_latitude_rad = _user_latitude_node->getDoubleValue() * SG_DEGREES_TO_RADIANS;
144 double user_longitude_rad = _user_longitude_node->getDoubleValue() * SG_DEGREES_TO_RADIANS;
146 double earth_rad_m = _earth_radius_node->getDoubleValue() * SG_FEET_TO_METER;
147 if( earth_rad_m < SG_EPSILON )
148 earth_rad_m = SG_EARTH_RAD * 1000;
150 // Placing the probes
152 for (int i = 0; i < sizeof(probe_lat_rad)/sizeof(probe_lat_rad[0]); i++) {
153 double probe_radius_ratio = dist_probe_m[i]/earth_rad_m;
155 probe_lat_rad[i] = asin(sin(user_latitude_rad)*cos(probe_radius_ratio)
156 +cos(user_latitude_rad)*sin(probe_radius_ratio)*cos(ground_wind_from_rad));
157 if (probe_lat_rad[i] < SG_EPSILON ) {
158 probe_lon_rad[i] = user_latitude_rad; // probe on a pole
160 probe_lon_rad[i] = fmod((user_longitude_rad+asin(sin(ground_wind_from_rad)
161 *sin(probe_radius_ratio)/cos(probe_lat_rad[i]))+SG_PI)
164 probe_lat_deg[i]= probe_lat_rad[i] * SG_RADIANS_TO_DEGREES;
165 probe_lon_deg[i]= probe_lon_rad[i] * SG_RADIANS_TO_DEGREES;
168 for (int i = 0; i < sizeof(probe_elev_m)/sizeof(probe_elev_m[0]); i++) {
169 if (globals->get_scenery()->get_elevation_m(SGGeod::fromGeodM(
170 SGGeod::fromRad(probe_lon_rad[i],probe_lat_rad[i]), 20000), alt, 0)) {
172 probe_elev_m[i] = alt;
174 probe_elev_m[i] = 0.1 ;
177 probe_elev_m[i] = 0.1;
183 slope[0] = (probe_elev_m[0] - probe_elev_m[1]) / dist_probe_m[1];
184 slope[1] = (probe_elev_m[1] - probe_elev_m[2]) / dist_probe_m[2];
185 slope[2] = (probe_elev_m[2] - probe_elev_m[3]) / dist_probe_m[3];
186 slope[3] = (probe_elev_m[4] - probe_elev_m[0]) / -dist_probe_m[4];
188 for (int i = 0; i < sizeof(slope)/sizeof(slope[0]); i++)
189 adj_slope[i] = sin(atan(5.0 * pow ( (fabs(slope[i])),1.7) ) ) *sign(slope[i]);
192 adj_slope[0] = 0.2 * adj_slope[0];
193 adj_slope[1] = 0.2 * adj_slope[1];
194 if ( adj_slope [2] < 0.0 ) {
195 adj_slope[2] = 0.5 * adj_slope[2];
200 if ( ( adj_slope [0] >= 0.0 ) && ( adj_slope [3] < 0.0 ) ) {
203 adj_slope[3] = 0.2 * adj_slope[3];
205 lift_factor = adj_slope[0]+adj_slope[1]+adj_slope[2]+adj_slope[3];
211 //user altitude above ground
212 double user_altitude_agl_m = _user_altitude_agl_ft_node->getDoubleValue() * SG_FEET_TO_METER;
215 double boundary2_m = 130.0; // in the lift
216 if (lift_factor < 0.0) { // in the sink
217 double highest_probe_temp= max ( probe_elev_m[1], probe_elev_m[2] );
218 double highest_probe_downwind_m= max ( highest_probe_temp, probe_elev_m[3] );
219 boundary2_m = highest_probe_downwind_m - probe_elev_m[0];
223 if ( user_altitude_agl_m < BOUNDARY1_m ) {
224 agl_factor = 0.5+0.5*user_altitude_agl_m /BOUNDARY1_m ;
225 } else if ( user_altitude_agl_m < boundary2_m ) {
228 agl_factor = exp(-(2 + probe_elev_m[0] / 2000) *
229 (user_altitude_agl_m - boundary2_m) / max(probe_elev_m[0],200.0));
232 double lift_mps = lift_factor* ground_wind_speed_mps * agl_factor;
234 //the updraft, finally, in ft per second
235 strength = fgGetLowPass( strength, lift_mps * SG_METER_TO_FEET, dt );
236 // if(isnan(strength)) strength=0;
237 _ridge_lift_fps_node->setDoubleValue( strength );