#include

#include #include "AIThermal.hxx" FGAIThermal::FGAIThermal() : FGAIBase(object_type::otThermal, false) { altitude_agl_ft = 0.0; max_strength = 6.0; diameter = 0.5; strength = factor = 0.0; cycle_timer = 60*(rand()%31); // some random in the birth time ground_elev_ft = 0.0; dt_count=0.9; alt=0.0; } void FGAIThermal::readFromScenario(SGPropertyNode* scFileNode) { if (!scFileNode) return; FGAIBase::readFromScenario(scFileNode); setMaxStrength(scFileNode->getDoubleValue("strength-fps", 8.0)); setDiameter(scFileNode->getDoubleValue("diameter-ft", 0.0)/6076.11549); setHeight(scFileNode->getDoubleValue("height-msl", 5000.0)); } bool FGAIThermal::init(ModelSearchOrder searchOrder) { factor = 8.0 * max_strength / (diameter * diameter * diameter); setAltitude( height ); _surface_wind_from_deg_node = fgGetNode("/environment/config/boundary/entry[0]/wind-from-heading-deg", true); _surface_wind_speed_node = fgGetNode("/environment/config/boundary/entry[0]/wind-speed-kt", true); _aloft_wind_from_deg_node = fgGetNode("/environment/config/aloft/entry[2]/wind-from-heading-deg", true); _aloft_wind_speed_node = fgGetNode("/environment/config/aloft/entry[2]/wind-speed-kt", true); do_agl_calc = true; return FGAIBase::init(searchOrder); } void FGAIThermal::bind() { FGAIBase::bind(); tie("position/altitude-agl-ft", SGRawValuePointer(&altitude_agl_ft)); tie("alt-rel", SGRawValuePointer(&alt_rel)); tie("time", SGRawValuePointer(&time)); tie("xx", SGRawValuePointer(&xx)); tie("is-forming", SGRawValuePointer(&is_forming)); tie("is-formed", SGRawValuePointer(&is_formed)); tie("is-dying", SGRawValuePointer(&is_dying)); tie("is-dead", SGRawValuePointer(&is_dead)); } void FGAIThermal::update(double dt) { FGAIBase::update(dt); Run(dt); Transform(); } //the formula to get the available portion of VUpMax depending on altitude //returns a double between 0 and 1 double FGAIThermal::get_strength_fac(double alt_frac) { double PI = 4.0 * atan(1.0); double fac = 0.0; if ( alt_frac <=0.0 ) { // do submarines get thermals ? fac = 0.0; } else if ((alt_frac > 0.0) && (alt_frac <= 0.1)) { // ground layer fac = ( 0.1*( pow( (10.0*alt_frac),10.0) ) ); } else if ((alt_frac > 0.1) && (alt_frac <= 1.0)) { // main body of the thermal fac = 0.4175 - 0.5825* ( cos ( PI* (1.0-sqrt(alt_frac) ) +PI) ) ; } else if ((alt_frac > 1.0) && (alt_frac < 1.1)) { //above the ceiling, but not above the cloud fac = (0.5 * ( 1.0 + cos ( PI*( (-2.0*alt_frac)*5.0 ) ) ) ); } else if (alt_frac >= 1.1) { //above the cloud fac = 0.0; } return fac; } void FGAIThermal::Run(double dt) { // ***************************************** // the thermal characteristics and variables // ***************************************** cycle_timer += dt ; // time // the time needed for the thermal to be completely formed double tmin1 = 5.0 ; // the time when the thermal begins to die double tmin2 = 20.0 ; // the time when the thermal is completely dead double tmin3 = 25.0; double alive_cycle_time = tmin3*60.0; //the time of the complete cycle, including a period of inactivity double tmin4 = 30.0; // some times expressed in a fraction of tmin3; double t1 = tmin1/tmin3 ; double t2 = tmin2/tmin3 ; double t3 = 1.0 ; double t4 = tmin4/tmin3; // the time elapsed since the thermal was born, in a 0-1 fraction of tmin3 time = cycle_timer/alive_cycle_time; // comment above and uncomment below to freeze the time cycle // time=0.5; if ( time >= t4) { cycle_timer = 60*(rand()%31); } //the position of the thermal 'top' double thermal_foot_lat = (pos.getLatitudeDeg()); double thermal_foot_lon = (pos.getLongitudeDeg()); //the max updraft one can expect in this thermal double MaxUpdraft=max_strength; //the max sink one can expect in this thermal, this is a negative number double MinUpdraft=-max_strength*0.25; //the fraction of MaxUpdraft one can expect at our height and time double maxstrengthavail; //max updraft at the user altitude and time double v_up_max; //min updraft at the user altitude and time, this is a negative number double v_up_min; double wind_speed; //max radius of the the thermal, including the sink area double Rmax = diameter/2.0; // 'shaping' of the thermal, the higher, the more conical the thermal- between 0 and 1 double shaping=0.8; //the radius of the thermal at our altitude in FT, including sink double Rsink; //the relative radius of the thermal where we have updraft, between 0 an 1 double r_up_frac=0.9; //radius of the thermal where we have updraft, in FT double Rup; //how far are we from the thermal center at our altitude in FEET double dist_center; //the position of the center of the thermal slice at our altitude double slice_center_lon; double slice_center_lat; //we need to know the thermal foot AGL altitude //we could do this only once, as thermal don't move //but then agl info is lost on user reset //so we only do this every 10 seconds to save cpu dt_count += dt; if (dt_count >= 10.0 ) { if (getGroundElevationM(SGGeod::fromGeodM(pos, 20000), alt, 0)) { ground_elev_ft = alt * SG_METER_TO_FEET; do_agl_calc = false; altitude_agl_ft = height - ground_elev_ft ; dt_count = 0.0; } } //user altitude relative to the thermal height, seen AGL from the thermal foot double user_altitude = globals->get_aircraft_position().getElevationFt(); if ( user_altitude < 1.0 ) { user_altitude = 1.0 ;}; // an ugly way to avoid NaNs for users at alt 0 double user_altitude_agl= ( user_altitude - ground_elev_ft ) ; alt_rel = user_altitude_agl / altitude_agl_ft; //the updraft user feels ! double Vup; // ********************* // environment variables // ********************* // the wind heading at the user alt double wind_heading_rad; // the "ambient" sink outside thermals double global_sink = -1.0; // ************** // some constants // ************** double PI = 4.0 * atan(1.0); // ****************** // thermal main cycle // ****************** //we get the max strenght proportion we can expect at the time and altitude, formuled between 0 and 1 //double xx; if (time <= t1) { xx = ( time / t1 ); maxstrengthavail = xx* get_strength_fac ( alt_rel / xx ); is_forming = true; is_formed = false; is_dying = false; is_dead = false; } else if ((time > t1) && (time <= t2)) { maxstrengthavail = get_strength_fac ( (alt_rel) ); is_forming = false; is_formed = true; is_dying = false; is_dead = false; } else if ((time > t2) && (time <= t3)) { xx= ( ( time - t2) / (1.0 - t2) ) ; maxstrengthavail = get_strength_fac ( alt_rel - xx ); is_forming = false; is_formed = false; is_dying = true; is_dead = false; } else { maxstrengthavail = 0.0; is_forming = false; is_formed = false; is_dying = false; is_dead = true; } //we get the diameter of the thermal slice at the user altitude //the thermal has a slight conic shape if ((alt_rel >= 0.0) && (alt_rel < 1.0)) { Rsink = (shaping * Rmax) + (((1.0 - shaping) * Rmax * alt_rel) / altitude_agl_ft); // in the main thermal body } else if ((alt_rel >= 1.0) && (alt_rel < 1.1)) { Rsink = (Rmax / 2.0) * (1.0 + cos((10.0 * PI * alt_rel) - (2.0 * PI))); // above the ceiling } else { Rsink = 0.0; // above the cloud } //we get the portion of the diameter that produces lift Rup = r_up_frac * Rsink ; //we now determine v_up_max and VupMin depending on our altitude v_up_max = maxstrengthavail * MaxUpdraft; v_up_min = maxstrengthavail * MinUpdraft; // Now we know, for current t and alt, v_up_max, VupMin, Rup, Rsink. // We still need to know how far we are from the thermal center // To determine the thermal inclinaison in the wind, we use a ugly approximation, // in which we say the thermal bends 20° (0.34906 rad ) for 10 kts wind. // We move the thermal foot upwind, to keep the cloud model over the "center" at ceiling level. // the displacement distance of the center of the thermal slice, at user altitude, // and relative to a hipothetical vertical thermal, would be: // get surface and 9000 ft wind double ground_wind_from_deg = _surface_wind_from_deg_node->getDoubleValue(); double ground_wind_speed_kts = _surface_wind_speed_node->getDoubleValue(); double aloft_wind_from_deg = _aloft_wind_from_deg_node->getDoubleValue(); double aloft_wind_speed_kts = _aloft_wind_speed_node->getDoubleValue(); double ground_wind_from_rad = (PI/2.0) - PI*( ground_wind_from_deg/180.0); double aloft_wind_from_rad = (PI/2.0) - PI*( aloft_wind_from_deg/180.0); wind_heading_rad= PI+ 0.5*( ground_wind_from_rad + aloft_wind_from_rad ); wind_speed = ground_wind_speed_kts + user_altitude * ( (aloft_wind_speed_kts -ground_wind_speed_kts ) / 9000.0 ); double dt_center_alt = -(tan (0.034906*wind_speed)) * ( altitude_agl_ft-user_altitude_agl ); // now, lets find how far we are from this shifted slice double dt_slice_lon_FT = ( dt_center_alt * cos ( wind_heading_rad )); double dt_slice_lat_FT = ( dt_center_alt * sin ( wind_heading_rad )); double dt_slice_lon = dt_slice_lon_FT / ft_per_deg_lon; double dt_slice_lat = dt_slice_lat_FT / ft_per_deg_lat; slice_center_lon = thermal_foot_lon + dt_slice_lon; slice_center_lat = thermal_foot_lat + dt_slice_lat; dist_center = SGGeodesy::distanceNm(SGGeod::fromDeg(slice_center_lon, slice_center_lat), globals->get_aircraft_position()); // Now we can calculate Vup if (max_strength >= 0.0) { // this is a thermal if ((dist_center >= 0.0) && (dist_center < Rup)) { //user is in the updraft area Vup = v_up_max * cos(dist_center * PI / (2.0 * Rup)); } else if ((dist_center > Rup) && (dist_center <= ((Rup + Rsink) / 2.0))) { //user in the 1st half of the sink area Vup = v_up_min * cos((dist_center - (Rup + Rsink) / 2.0) * PI / (2.0 * ((Rup + Rsink) / 2.0 - Rup))); } else if ((dist_center > ((Rup + Rsink) / 2.0)) && dist_center <= Rsink) { // user in the 2nd half of the sink area Vup = (global_sink + v_up_min) / 2.0 + (global_sink - v_up_min) / 2.0 * cos((dist_center - Rsink) * PI / ((Rsink - Rup) / 2.0)); } else { // outside the thermal Vup = global_sink; } } 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 if (alt_rel <= 1.1) { double fac = (1.0 - (1.0 - 1.815 * alt_rel) * (1.0 - 1.815 * alt_rel)); Vup = fac * (global_sink + ((v_up_max - global_sink) / 2.0) * (1.0 + cos(dist_center * PI / Rmax))); } else { Vup = global_sink; } } //correct for no global sink above clouds and outside thermals if ( ( (alt_rel > 1.0) && (alt_rel <1.1)) && ( dist_center > Rsink ) ) { Vup = global_sink * ( 11.0 -10.0 * alt_rel ); } if ( alt_rel >= 1.1 ) { Vup = 0.0; } strength = Vup; range = dist_center; }