1 // FGAIThermal - FGAIBase-derived class creates an AI thermal
3 // Original by Written by David Culp
5 // An attempt to refine the thermal shape and behaviour by WooT 2009
7 // Copyright (C) 2009 Patrice Poly ( WooT )
9 // This program is free software; you can redistribute it and/or
10 // modify it under the terms of the GNU General Public License as
11 // published by the Free Software Foundation; either version 2 of the
12 // License, or (at your option) any later version.
14 // This program is distributed in the hope that it will be useful, but
15 // WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 // General Public License for more details.
19 // You should have received a copy of the GNU General Public License
20 // along with this program; if not, write to the Free Software
21 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
27 #include <Main/fg_props.hxx>
28 #include <Main/globals.hxx>
29 #include <Scenery/scenery.hxx>
35 #include "AIThermal.hxx"
38 FGAIThermal::FGAIThermal() : FGAIBase(otThermal) {
41 strength = factor = 0.0;
42 cycle_timer = 60*(rand()%31); // some random in the birth time
48 FGAIThermal::~FGAIThermal() {
51 void FGAIThermal::readFromScenario(SGPropertyNode* scFileNode) {
55 FGAIBase::readFromScenario(scFileNode);
57 setMaxStrength(scFileNode->getDoubleValue("strength-fps", 8.0));
58 setDiameter(scFileNode->getDoubleValue("diameter-ft", 0.0)/6076.11549);
59 setHeight(scFileNode->getDoubleValue("height-msl", 5000.0));
62 bool FGAIThermal::init(bool search_in_AI_path) {
63 factor = 8.0 * max_strength / (diameter * diameter * diameter);
64 setAltitude( height );
65 _surface_wind_from_deg_node =
66 fgGetNode("/environment/config/boundary/entry[0]/wind-from-heading-deg", true);
67 _surface_wind_speed_node =
68 fgGetNode("/environment/config/boundary/entry[0]/wind-speed-kt", true);
69 _aloft_wind_from_deg_node =
70 fgGetNode("/environment/config/aloft/entry[2]/wind-from-heading-deg", true);
71 _aloft_wind_speed_node =
72 fgGetNode("/environment/config/aloft/entry[2]/wind-speed-kt", true);
74 return FGAIBase::init(search_in_AI_path);
77 void FGAIThermal::bind() {
78 props->tie("position/altitude-agl-ft", // for debug and tweak
79 SGRawValuePointer<double>(&altitude_agl_ft));
80 props->tie("alt-rel", // for debug and tweak
81 SGRawValuePointer<double>(&alt_rel));
82 props->tie("time", // for debug and tweak
83 SGRawValuePointer<double>(&time));
84 props->tie("xx", // for debug and tweak
85 SGRawValuePointer<double>(&xx));
86 props->tie("is-forming", // for debug abd tweak
87 SGRawValuePointer<bool>(&is_forming));
88 props->tie("is-formed", // for debug abd tweak
89 SGRawValuePointer<bool>(&is_formed));
90 props->tie("is-dying", // for debug abd tweak
91 SGRawValuePointer<bool>(&is_dying));
92 props->tie("is-dead", // for debug abd tweak
93 SGRawValuePointer<bool>(&is_dead));
97 void FGAIThermal::unbind() {
98 props->untie("position/altitude-agl-ft");
99 props->untie("alt-rel");
100 props->untie("time");
101 props->untie("is-forming");
102 props->untie("is-formed");
103 props->untie("is-dying");
104 props->untie("is-dead");
110 void FGAIThermal::update(double dt) {
111 FGAIBase::update(dt);
118 //the formula to get the available portion of VUpMax depending on altitude
119 //returns a double between 0 and 1
120 double FGAIThermal::get_strength_fac(double alt_frac) {
122 double PI = 4.0 * atan(1.0);
124 if ( alt_frac <=0.0 ) { // do submarines get thermals ?
128 else if ( ( alt_frac>0.0 ) && (alt_frac<=0.1) ) { // ground layer
129 fac = ( 0.1*( pow( (10.0*alt_frac),10.0) ) );
132 else if ( ( alt_frac>0.1 ) && (alt_frac<=1.0) ) { // main body of the thermal
133 fac = 0.4175 - 0.5825* ( cos ( PI* (1.0-sqrt(alt_frac) ) +PI) ) ;
136 else if ( ( alt_frac >1.0 ) && (alt_frac < 1.1 ) ) { //above the ceiling, but not above the cloud
137 fac = (0.5 * ( 1.0 + cos ( PI*( (-2.0*alt_frac)*5.0 ) ) ) );
140 else if ( alt_frac >= 1.1 ) { //above the cloud
147 void FGAIThermal::Run(double dt) {
149 // *****************************************
150 // the thermal characteristics and variables
151 // *****************************************
157 // the time needed for the thermal to be completely formed
159 // the time when the thermal begins to die
160 double tmin2 = 20.0 ;
161 // the time when the thermal is completely dead
163 double alive_cycle_time = tmin3*60.0;
164 //the time of the complete cycle, including a period of inactivity
166 // some times expressed in a fraction of tmin3;
167 double t1 = tmin1/tmin3 ;
168 double t2 = tmin2/tmin3 ;
170 double t4 = tmin4/tmin3;
171 // the time elapsed since the thermal was born, in a 0-1 fraction of tmin3
173 time = cycle_timer/alive_cycle_time;
175 //uncomment below to freeze the time cycle
179 cycle_timer = 60*(rand()%31);
183 //the position of the thermal 'top'
184 double thermal_foot_lat = (pos.getLatitudeDeg());
185 double thermal_foot_lon = (pos.getLongitudeDeg());
187 //the max updraft one can expect in this thermal
188 double MaxUpdraft=max_strength;
189 //the max sink one can expect in this thermal, this is a negative number
190 double MinUpdraft=-max_strength*0.25;
191 //the fraction of MaxUpdraft one can expect at our height and time
192 double maxstrengthavail;
193 //max updraft at the user altitude and time
195 //min updraft at the user altitude and time, this is a negative number
200 //max radius of the the thermal, including the sink area
201 double Rmax = diameter/2.0;
202 // 'shaping' of the thermal, the higher, the more conical the thermal- between 0 and 1
204 //the radius of the thermal at our altitude in FT, including sink
206 //the relative radius of the thermal where we have updraft, between 0 an 1
207 double r_up_frac=0.9;
208 //radius of the thermal where we have updraft, in FT
210 //how far are we from the thermal center at our altitude in FEET
213 //the position of the center of the thermal slice at our altitude
214 double slice_center_lon;
215 double slice_center_lat;
219 // **************************************
220 // various variables relative to the user
221 // **************************************
223 double user_latitude = manager->get_user_latitude();
224 double user_longitude = manager->get_user_longitude();
225 double user_altitude = manager->get_user_altitude(); // MSL
227 //we need to know the thermal foot AGL altitude
230 //we could do this only once, as thermal don't move
231 //but then agl info is lost on user reset
232 //so we only do this every 10 seconds to save cpu
234 if (dt_count >= 10.0 ) {
236 if (globals->get_scenery()->get_elevation_m(SGGeod::fromGeodM(pos, 20000), alt, 0)){
237 ground_elev_ft = alt * SG_METER_TO_FEET;
239 altitude_agl_ft = height - ground_elev_ft ;
244 //user altitude relative to the thermal height, seen AGL from the thermal foot
245 if ( user_altitude < 1.0 ) { user_altitude = 1.0 ;}; // an ugly way to avoid NaNs for users at alt 0
246 double user_altitude_agl= ( user_altitude - ground_elev_ft ) ;
247 alt_rel = user_altitude_agl / altitude_agl_ft;
251 //the updraft user feels !
254 // *********************
255 // environment variables
256 // *********************
258 // the windspeed at the user alt in KT
261 // the wind heading at the user alt
263 double wind_heading_deg;
264 double wind_heading_rad;
266 // the "ambient" sink outside thermals
267 double global_sink = -1.0;
273 double PI = 4.0 * atan(1.0);
276 // ******************
277 // thermal main cycle
278 // ******************
280 //we get the max strenght proportion we can expect at the time and altitude, formuled between 0 and 1
284 maxstrengthavail = xx* get_strength_fac ( alt_rel / xx );
286 is_forming=1;is_formed=0;is_dying=0;is_dead=0;
289 else if ( (time > t1) && (time <= t2) ) {
290 maxstrengthavail = get_strength_fac ( (alt_rel) );
292 is_forming=0;is_formed=1;is_dying=0;is_dead=0;
295 else if ( (time > t2) && (time <= t3) ) {
296 xx= ( ( time - t2) / (1.0 - t2) ) ;
297 maxstrengthavail = get_strength_fac ( alt_rel - xx );
299 is_forming=0;is_formed=0;is_dying=1;is_dead=0;
303 maxstrengthavail = 0.0;
304 is_forming=0;is_formed=0;is_dying=0;is_dead=1;
308 //we get the diameter of the thermal slice at the user altitude
309 //the thermal has a slight conic shape
311 if ( (alt_rel >= 0.0) && (alt_rel < 1.0 ) ) {
312 Rsink = ( shaping*Rmax ) + ( ( (1.0-shaping)*Rmax*alt_rel ) / altitude_agl_ft ); // in the main thermal body
314 else if ( (alt_rel >=1.0) && (alt_rel < 1.1) ) {
315 Rsink = (Rmax/2.0) * ( 1.0+ cos ( (10.0*PI*alt_rel)-(2.0*PI) ) ); // above the ceiling
318 Rsink = 0.0; // above the cloud
321 //we get the portion of the diameter that produces lift
322 Rup = r_up_frac * Rsink ;
324 //we now determine v_up_max and VupMin depending on our altitude
326 v_up_max = maxstrengthavail * MaxUpdraft;
327 v_up_min = maxstrengthavail * MinUpdraft;
329 // Now we know, for current t and alt, v_up_max, VupMin, Rup, Rsink.
331 // We still need to know how far we are from the thermal center
333 // To determine the thermal inclinaison in the wind, we use a ugly approximation,
334 // in which we say the thermal bends 20° (0.34906 rad ) for 10 kts wind.
335 // We move the thermal foot upwind, to keep the cloud model over the "center" at ceiling level.
336 // the displacement distance of the center of the thermal slice, at user altitude,
337 // and relative to a hipothetical vertical thermal, would be:
339 // get surface and 9000 ft wind
341 double ground_wind_from_deg = _surface_wind_from_deg_node->getDoubleValue();
342 double ground_wind_speed_kts = _surface_wind_speed_node->getDoubleValue();
343 double aloft_wind_from_deg = _aloft_wind_from_deg_node->getDoubleValue();
344 double aloft_wind_speed_kts = _aloft_wind_speed_node->getDoubleValue();
346 double ground_wind_from_rad = (PI/2.0) - PI*( ground_wind_from_deg/180.0);
347 double aloft_wind_from_rad = (PI/2.0) - PI*( aloft_wind_from_deg/180.0);
349 wind_heading_rad= PI+ 0.5*( ground_wind_from_rad + aloft_wind_from_rad );
351 wind_speed = ground_wind_speed_kts + user_altitude * ( (aloft_wind_speed_kts -ground_wind_speed_kts ) / 9000.0 );
353 double dt_center_alt = -(tan (0.034906*wind_speed)) * ( altitude_agl_ft-user_altitude_agl );
355 // now, lets find how far we are from this shifted slice
357 double dt_slice_lon_FT = ( dt_center_alt * cos ( wind_heading_rad ));
358 double dt_slice_lat_FT = ( dt_center_alt * sin ( wind_heading_rad ));
360 double dt_slice_lon = dt_slice_lon_FT / ft_per_deg_lon;
361 double dt_slice_lat = dt_slice_lat_FT / ft_per_deg_lat;
363 slice_center_lon = thermal_foot_lon + dt_slice_lon;
364 slice_center_lat = thermal_foot_lat + dt_slice_lat;
366 double dist_center_lon = fabs(slice_center_lon - user_longitude)* ft_per_deg_lon;
367 double dist_center_lat = fabs(slice_center_lat - user_latitude)* ft_per_deg_lat;
369 double dist_center_FT = sqrt ( dist_center_lon*dist_center_lon + dist_center_lat*dist_center_lat ); // feet
371 dist_center = dist_center_FT/ 6076.11549; //nautic miles
374 // Now we can calculate Vup
376 if ( max_strength >=0.0 ) { // this is a thermal
378 if ( ( dist_center >= 0.0 ) && ( dist_center < Rup ) ) { //user is in the updraft area
379 Vup = v_up_max * cos ( dist_center* PI/(2.0*Rup) );
381 else if ( ( dist_center > Rup ) && ( dist_center <= ((Rup+Rsink)/2.0) ) ) { //user in the 1st half of the sink area
382 Vup = v_up_min * cos (( dist_center - ( Rup+Rsink)/2.0 ) * PI / ( 2.0* ( ( Rup+Rsink)/2.0 -Rup )));
384 else if ( ( dist_center > ((Rup+Rsink)/2.0) ) && dist_center <= Rsink ) { // user in the 2nd half of the sink area
385 Vup = ( global_sink + v_up_min )/2.0 + ( global_sink - v_up_min )/2.0 *cos ( (dist_center-Rsink) *PI/ ( (Rsink-Rup )/2.0) );
387 else { // outside the thermal
392 else { // this is a sink, we don't want updraft on the sides, nor do we want to feel sink near or above ceiling and ground
393 if ( alt_rel <=1.1 ) {
394 double fac = ( 1.0 - ( 1.0 - 1.815*alt_rel)*( 1.0 - 1.815*alt_rel) );
395 Vup = fac * (global_sink + ( ( v_up_max - global_sink )/2.0 ) * ( 1.0+cos ( dist_center* PI / Rmax ) )) ;
397 else { Vup = global_sink; }
400 //correct for no global sink above clouds and outside thermals
401 if ( ( (alt_rel > 1.0) && (alt_rel <1.1)) && ( dist_center > Rsink ) ) {
402 Vup = global_sink * ( 11.0 -10.0 * alt_rel );
404 if ( alt_rel >= 1.1 ) {