src/AIModel/AIThermal.cxx

   1 // FGAIThermal - FGAIBase-derived class creates an AI thermal
   2 //
   3 // Copyright (C) 2004  David P. Culp - davidculp2@comcast.net
   4 //
   5 // An attempt to refine the thermal shape and behaviour by WooT 2009
   6 //
   7 // Copyright (C) 2009 Patrice Poly ( WooT )
   8 //
   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.
  13 //
  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.
  18 //
  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.
  22
  23 #ifdef HAVE_CONFIG_H
  24 #  include <config.h>
  25 #endif
  26
  27 #include <Main/fg_props.hxx>
  28 #include <Main/globals.hxx>
  29 #include <Scenery/scenery.hxx>
  30 #include <string>
  31 #include <math.h>
  32
  33 using std::string;
  34
  35 #include "AIThermal.hxx"
  36
  37
  38 FGAIThermal::FGAIThermal() :
  39    FGAIBase(otThermal, false)
  40 {
  41    max_strength = 6.0;
  42    diameter = 0.5;
  43    strength = factor = 0.0;
  44    cycle_timer = 60*(rand()%31); // some random in the birth time
  45    ground_elev_ft = 0.0;
  46    dt_count=0.9;
  47    alt=0.0;
  48 }
  49
  50 FGAIThermal::~FGAIThermal() {
  51 }
  52
  53 void FGAIThermal::readFromScenario(SGPropertyNode* scFileNode) {
  54   if (!scFileNode)
  55     return;
  56
  57   FGAIBase::readFromScenario(scFileNode);
  58
  59   setMaxStrength(scFileNode->getDoubleValue("strength-fps", 8.0));
  60   setDiameter(scFileNode->getDoubleValue("diameter-ft", 0.0)/6076.11549);
  61   setHeight(scFileNode->getDoubleValue("height-msl", 5000.0));
  62 }
  63
  64 bool FGAIThermal::init(bool search_in_AI_path) {
  65    factor = 8.0 * max_strength / (diameter * diameter * diameter);
  66    setAltitude( height );
  67    _surface_wind_from_deg_node =
  68             fgGetNode("/environment/config/boundary/entry[0]/wind-from-heading-deg", true);
  69    _surface_wind_speed_node =
  70             fgGetNode("/environment/config/boundary/entry[0]/wind-speed-kt", true);
  71    _aloft_wind_from_deg_node =
  72             fgGetNode("/environment/config/aloft/entry[2]/wind-from-heading-deg", true);
  73    _aloft_wind_speed_node =
  74             fgGetNode("/environment/config/aloft/entry[2]/wind-speed-kt", true);
  75     do_agl_calc = 1;
  76    return FGAIBase::init(search_in_AI_path);
  77 }
  78
  79 void FGAIThermal::bind() {
  80         props->tie("position/altitude-agl-ft", // for debug and tweak
  81                 SGRawValuePointer<double>(&altitude_agl_ft));
  82         props->tie("alt-rel", // for debug and tweak
  83                 SGRawValuePointer<double>(&alt_rel));
  84         props->tie("time", // for debug and tweak
  85                 SGRawValuePointer<double>(&time));
  86         props->tie("xx", // for debug and tweak
  87                 SGRawValuePointer<double>(&xx));
  88         props->tie("is-forming", // for debug abd tweak
  89                 SGRawValuePointer<bool>(&is_forming));
  90         props->tie("is-formed", // for debug abd tweak
  91                 SGRawValuePointer<bool>(&is_formed));
  92         props->tie("is-dying", // for debug abd tweak
  93                 SGRawValuePointer<bool>(&is_dying));
  94         props->tie("is-dead", // for debug abd tweak
  95                 SGRawValuePointer<bool>(&is_dead));
  96     FGAIBase::bind();
  97 }
  98
  99 void FGAIThermal::unbind() {
 100         props->untie("position/altitude-agl-ft");
 101         props->untie("alt-rel");
 102         props->untie("time");
 103         props->untie("is-forming");
 104         props->untie("is-formed");
 105         props->untie("is-dying");
 106         props->untie("is-dead");
 107         props->untie("xx");
 108     FGAIBase::unbind();
 109 }
 110
 111
 112 void FGAIThermal::update(double dt) {
 113    FGAIBase::update(dt);
 114    Run(dt);
 115    Transform();
 116 }
 117
 118
 119
 120 //the formula to get the available portion of VUpMax depending on altitude
 121 //returns a double between 0 and 1
 122 double FGAIThermal::get_strength_fac(double alt_frac) {
 123
 124 double PI = 4.0 * atan(1.0);
 125 double fac = 0.0;
 126 if ( alt_frac <=0.0 ) { // do submarines get thermals ?
 127         fac = 0.0;
 128         }
 129 else if ( ( alt_frac>0.0 ) && (alt_frac<=0.1) ) { // ground layer
 130         fac = ( 0.1*( pow( (10.0*alt_frac),10.0) ) );
 131         }
 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) ) ;
 134         }
 135 else if ( ( alt_frac >1.0 ) && (alt_frac < 1.1 ) ) {  //above the ceiling, but not above the cloud
 136         fac = (0.5 * ( 1.0 + cos ( PI*( (-2.0*alt_frac)*5.0 ) ) ) );
 137         }
 138 else if ( alt_frac >= 1.1 ) {  //above the cloud
 139         fac = 0.0;
 140         }
 141 return fac;
 142 }
 143
 144
 145 void FGAIThermal::Run(double dt) {
 146
 147 // *****************************************
 148 // the thermal characteristics and variables
 149 // *****************************************
 150
 151 cycle_timer += dt ;
 152
 153 // time
 154
 155 // the time needed for the thermal to be completely formed
 156 double tmin1 = 5.0 ;
 157 // the time when the thermal begins to die
 158 double tmin2 = 20.0 ;
 159 // the time when the thermal is completely dead
 160 double tmin3 = 25.0;
 161 double alive_cycle_time = tmin3*60.0;
 162 //the time of the complete cycle, including a period of inactivity
 163 double tmin4 = 30.0;
 164 // some times expressed in a fraction of tmin3;
 165 double t1 = tmin1/tmin3 ;
 166 double t2 = tmin2/tmin3 ;
 167 double t3 = 1.0 ;
 168 double t4 = tmin4/tmin3;
 169 // the time elapsed since the thermal was born, in a 0-1 fraction of tmin3
 170
 171 time = cycle_timer/alive_cycle_time;
 172 //comment above and
 173 //uncomment below to freeze the time cycle
 174  time=0.5;
 175
 176 if ( time >= t4) {
 177         cycle_timer = 60*(rand()%31);
 178         }
 179
 180
 181 //the position of the thermal 'top'
 182 double thermal_foot_lat = (pos.getLatitudeDeg());
 183 double thermal_foot_lon = (pos.getLongitudeDeg());
 184
 185 //the max updraft one can expect in this thermal
 186 double MaxUpdraft=max_strength;
 187 //the max sink one can expect in this thermal, this is a negative number
 188 double MinUpdraft=-max_strength*0.25;
 189 //the fraction of MaxUpdraft one can expect at our height and time
 190 double maxstrengthavail;
 191 //max updraft at the user altitude and time
 192 double v_up_max;
 193 //min updraft at the user altitude and time, this is a negative number
 194 double v_up_min;
 195 double wind_speed;
 196
 197
 198 //max radius of the the thermal, including the sink area
 199 double Rmax = diameter/2.0;
 200 // 'shaping' of the thermal, the higher, the more conical the thermal- between 0 and 1
 201 double shaping=0.8;
 202 //the radius of the thermal at our altitude in FT, including sink
 203 double Rsink;
 204 //the relative radius of the thermal where we have updraft, between 0 an 1
 205 double r_up_frac=0.9;
 206 //radius of the thermal where we have updraft, in FT
 207 double Rup;
 208 //how far are we from the thermal center at our altitude in FEET
 209 double dist_center;
 210
 211 //the position of the center of the thermal slice at our altitude
 212 double slice_center_lon;
 213 double slice_center_lat;
 214
 215
 216
 217 // **************************************
 218 // various variables relative to the user
 219 // **************************************
 220
 221 double user_latitude  = manager->get_user_latitude();
 222 double user_longitude = manager->get_user_longitude();
 223 double user_altitude  = manager->get_user_altitude(); // MSL
 224
 225 //we need to know the thermal foot AGL altitude
 226
 227
 228 //we could do this only once, as thermal don't move
 229 //but then agl info is lost on user reset
 230 //so we only do this every 10 seconds to save cpu
 231 dt_count += dt;
 232 if (dt_count >= 10.0 ) {
 233         //double alt;
 234         if (getGroundElevationM(SGGeod::fromGeodM(pos, 20000), alt, 0)) {
 235         ground_elev_ft =  alt * SG_METER_TO_FEET;
 236         do_agl_calc = 0;
 237         altitude_agl_ft = height - ground_elev_ft ;
 238         dt_count = 0.0;
 239         }
 240 }
 241
 242 //user altitude relative to the thermal height, seen AGL from the thermal foot
 243 if ( user_altitude < 1.0 ) { user_altitude = 1.0 ;}; // an ugly way to avoid NaNs for users at alt 0
 244 double user_altitude_agl= ( user_altitude - ground_elev_ft ) ;
 245 alt_rel = user_altitude_agl / altitude_agl_ft;
 246
 247
 248
 249 //the updraft user feels !
 250 double Vup;
 251
 252 // *********************
 253 // environment variables
 254 // *********************
 255
 256 // the  wind heading at the user alt
 257 double wind_heading_rad;
 258
 259 // the "ambient" sink outside thermals
 260 double global_sink = -1.0;
 261
 262 // **************
 263 // some constants
 264 // **************
 265
 266 double PI = 4.0 * atan(1.0);
 267
 268
 269 // ******************
 270 // thermal main cycle
 271 // ******************
 272
 273 //we get the max strenght proportion we can expect at the time and altitude, formuled between 0 and 1
 274 //double xx;
 275 if (time <= t1) {
 276         xx= ( time / t1 );
 277         maxstrengthavail = xx* get_strength_fac ( alt_rel / xx );
 278
 279         is_forming=1;is_formed=0;is_dying=0;is_dead=0;
 280
 281         }
 282 else if ( (time > t1) && (time <= t2) ) {
 283         maxstrengthavail = get_strength_fac ( (alt_rel) );
 284
 285         is_forming=0;is_formed=1;is_dying=0;is_dead=0;
 286
 287         }
 288 else if ( (time > t2) && (time <= t3) ) {
 289         xx= ( ( time - t2) / (1.0 - t2) ) ;
 290         maxstrengthavail = get_strength_fac ( alt_rel - xx );
 291
 292         is_forming=0;is_formed=0;is_dying=1;is_dead=0;
 293
 294         }
 295 else {
 296         maxstrengthavail = 0.0;
 297         is_forming=0;is_formed=0;is_dying=0;is_dead=1;
 298
 299         }
 300
 301 //we get the diameter of the thermal slice at the user altitude
 302 //the thermal has a slight conic shape
 303
 304 if ( (alt_rel >= 0.0) && (alt_rel < 1.0 ) ) {
 305         Rsink = ( shaping*Rmax ) + ( (  (1.0-shaping)*Rmax*alt_rel ) / altitude_agl_ft );  // in the main thermal body
 306         }
 307 else if ( (alt_rel >=1.0) && (alt_rel < 1.1) ) {
 308         Rsink = (Rmax/2.0) * ( 1.0+ cos ( (10.0*PI*alt_rel)-(2.0*PI) ) ); // above the ceiling
 309         }
 310 else {
 311         Rsink = 0.0; // above the cloud
 312         }
 313
 314 //we get the portion of the diameter that produces lift
 315 Rup = r_up_frac * Rsink ;
 316
 317 //we now determine v_up_max and VupMin depending on our altitude
 318
 319 v_up_max = maxstrengthavail * MaxUpdraft;
 320 v_up_min = maxstrengthavail * MinUpdraft;
 321
 322 // Now we know, for current t and alt, v_up_max, VupMin, Rup, Rsink.
 323
 324 // We still need to know how far we are from the thermal center
 325
 326 // To determine the thermal inclinaison in the wind, we use a ugly approximation,
 327 // in which we say the thermal bends 20° (0.34906 rad ) for 10 kts wind.
 328 // We move the thermal foot upwind, to keep the cloud model over the "center" at ceiling level.
 329 // the displacement distance of the center of the thermal slice, at user altitude,
 330 // and relative to a hipothetical vertical thermal,  would be:
 331
 332 // get surface and 9000 ft wind
 333
 334 double ground_wind_from_deg = _surface_wind_from_deg_node->getDoubleValue();
 335 double ground_wind_speed_kts  = _surface_wind_speed_node->getDoubleValue();
 336 double aloft_wind_from_deg = _aloft_wind_from_deg_node->getDoubleValue();
 337 double aloft_wind_speed_kts  = _aloft_wind_speed_node->getDoubleValue();
 338
 339 double ground_wind_from_rad = (PI/2.0) - PI*( ground_wind_from_deg/180.0);
 340 double aloft_wind_from_rad = (PI/2.0) - PI*( aloft_wind_from_deg/180.0);
 341
 342 wind_heading_rad= PI+ 0.5*( ground_wind_from_rad + aloft_wind_from_rad );
 343
 344 wind_speed = ground_wind_speed_kts + user_altitude * ( (aloft_wind_speed_kts -ground_wind_speed_kts ) / 9000.0 );
 345
 346 double dt_center_alt = -(tan (0.034906*wind_speed)) * ( altitude_agl_ft-user_altitude_agl );
 347
 348 // now, lets find how far we are from this shifted slice
 349
 350 double dt_slice_lon_FT = ( dt_center_alt * cos ( wind_heading_rad ));
 351 double dt_slice_lat_FT = ( dt_center_alt * sin ( wind_heading_rad ));
 352
 353 double dt_slice_lon = dt_slice_lon_FT / ft_per_deg_lon;
 354 double dt_slice_lat = dt_slice_lat_FT / ft_per_deg_lat;
 355
 356 slice_center_lon = thermal_foot_lon + dt_slice_lon;
 357 slice_center_lat = thermal_foot_lat + dt_slice_lat;
 358
 359 double dist_center_lon = fabs(slice_center_lon - user_longitude)* ft_per_deg_lon;
 360 double dist_center_lat = fabs(slice_center_lat - user_latitude)* ft_per_deg_lat;
 361
 362 double dist_center_FT = sqrt ( dist_center_lon*dist_center_lon + dist_center_lat*dist_center_lat ); // feet
 363
 364 dist_center = dist_center_FT/ 6076.11549; //nautic miles
 365
 366
 367 // Now we can calculate Vup
 368
 369 if ( max_strength >=0.0 ) { // this is a thermal
 370
 371         if ( ( dist_center >= 0.0 ) && ( dist_center < Rup ) ) {  //user is in the updraft area
 372                 Vup = v_up_max * cos ( dist_center* PI/(2.0*Rup) );
 373                 }
 374         else if ( ( dist_center > Rup ) && ( dist_center <= ((Rup+Rsink)/2.0) ) ) { //user in the 1st half of the sink area
 375                 Vup = v_up_min * cos (( dist_center - ( Rup+Rsink)/2.0 ) * PI / ( 2.0* (  ( Rup+Rsink)/2.0 -Rup )));
 376                 }
 377         else if ( ( dist_center > ((Rup+Rsink)/2.0) ) && dist_center <= Rsink ) {   // user in the 2nd half of the sink area
 378                 Vup = ( global_sink + v_up_min )/2.0 + ( global_sink - v_up_min )/2.0 *cos ( (dist_center-Rsink) *PI/ ( (Rsink-Rup )/2.0) );
 379                 }
 380         else {  // outside the thermal
 381                 Vup = global_sink;
 382                 }
 383         }
 384
 385 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
 386         if ( alt_rel <=1.1 ) {
 387                 double fac =  ( 1.0 - ( 1.0 - 1.815*alt_rel)*( 1.0 - 1.815*alt_rel) );
 388                 Vup = fac * (global_sink + ( ( v_up_max - global_sink )/2.0 ) * ( 1.0+cos ( dist_center* PI / Rmax ) )) ;
 389                 }
 390         else { Vup = global_sink; }
 391 }
 392
 393 //correct for no global sink above clouds and outside thermals
 394 if ( ( (alt_rel > 1.0) && (alt_rel <1.1)) && ( dist_center > Rsink ) ) {
 395         Vup = global_sink * ( 11.0 -10.0 * alt_rel );
 396         }
 397 if ( alt_rel >= 1.1 ) {
 398         Vup = 0.0;
 399         }
 400
 401 strength = Vup;
 402 range = dist_center;
 403
 404 }
 405
 406
 407