From 7e3531ef5cae0c0e9ac325de8712f3b67a3e5e5f Mon Sep 17 00:00:00 2001 From: torsten Date: Mon, 20 Apr 2009 14:20:05 +0000 Subject: [PATCH] Implementation of ridge lift from Patrice Poly based on an algorithm of Ian Forster-Lewis --- src/AIModel/AIManager.cxx | 8 +- src/AIModel/AIManager.hxx | 2 +- src/AIModel/AIThermal.cxx | 375 +++++++++++++++++++++++++--- src/AIModel/AIThermal.hxx | 46 +++- src/Environment/Makefile.am | 3 +- src/Environment/environment.cxx | 36 +++ src/Environment/environment.hxx | 7 + src/Environment/environment_mgr.cxx | 14 ++ src/Environment/ridge_lift.cxx | 317 +++++++++++++++++++++++ src/Environment/ridge_lift.hxx | 125 ++++++++++ src/Main/fg_init.cxx | 8 + 11 files changed, 892 insertions(+), 49 deletions(-) create mode 100644 src/Environment/ridge_lift.cxx create mode 100644 src/Environment/ridge_lift.hxx diff --git a/src/AIModel/AIManager.cxx b/src/AIModel/AIManager.cxx index b01a1ac10..c23785041 100644 --- a/src/AIModel/AIManager.cxx +++ b/src/AIModel/AIManager.cxx @@ -65,7 +65,7 @@ FGAIManager::init() { if (!enabled) return; - wind_from_down_node = fgGetNode("/environment/wind-from-down-fps", true); + thermal_lift_node = fgGetNode("/environment/thermal-lift-fps", true); wind_from_east_node = fgGetNode("/environment/wind-from-east-fps",true); wind_from_north_node = fgGetNode("/environment/wind-from-north-fps",true); @@ -172,7 +172,7 @@ FGAIManager::update(double dt) { } } - wind_from_down_node->setDoubleValue( strength ); // for thermals + thermal_lift_node->setDoubleValue( strength ); // for thermals } void @@ -258,6 +258,8 @@ FGAIManager::processThermal( FGAIThermal* thermal ) { } + + void FGAIManager::processScenario( const string &filename ) { @@ -397,7 +399,7 @@ FGAIManager::calcCollision(double alt, double lat, double lon, double fuse_range tgt_ht[type] += fuse_range; if (fabs(tgt_alt - alt) > tgt_ht[type] || type == FGAIBase::otBallistic - || type == FGAIBase::otStorm || type == FGAIBase::otThermal) { + || type == FGAIBase::otStorm || type == FGAIBase::otThermal ) { SG_LOG(SG_GENERAL, SG_DEBUG, "AIManager: skipping " << fabs(tgt_alt - alt) << " " diff --git a/src/AIModel/AIManager.hxx b/src/AIModel/AIManager.hxx index 8b3582b39..3f6dd9c2f 100644 --- a/src/AIModel/AIManager.hxx +++ b/src/AIModel/AIManager.hxx @@ -100,7 +100,7 @@ private: double calcRange(double lat, double lon, double lat2, double lon2)const; SGPropertyNode_ptr root; - SGPropertyNode_ptr wind_from_down_node; + SGPropertyNode_ptr thermal_lift_node; SGPropertyNode_ptr user_latitude_node; SGPropertyNode_ptr user_longitude_node; SGPropertyNode_ptr user_altitude_node; diff --git a/src/AIModel/AIThermal.cxx b/src/AIModel/AIThermal.cxx index e9ad927a1..da9987fb8 100644 --- a/src/AIModel/AIThermal.cxx +++ b/src/AIModel/AIThermal.cxx @@ -1,8 +1,10 @@ // FGAIThermal - FGAIBase-derived class creates an AI thermal // -// Written by David Culp, started Feb 2004. +// Original by Written by David Culp // -// Copyright (C) 2004 David P. Culp - davidculp2@comcast.net +// An attempt to refine the thermal shape and behaviour by WooT 2009 +// +// Copyright (C) 2009 Patrice Poly ( WooT ) // // This program is free software; you can redistribute it and/or // modify it under the terms of the GNU General Public License as @@ -37,6 +39,10 @@ FGAIThermal::FGAIThermal() : FGAIBase(otThermal) { 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; } FGAIThermal::~FGAIThermal() { @@ -49,21 +55,54 @@ void FGAIThermal::readFromScenario(SGPropertyNode* scFileNode) { 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)); + setDiameter(scFileNode->getDoubleValue("diameter-ft", 0.0)/6076.11549); + setHeight(scFileNode->getDoubleValue("height-msl", 5000.0)); } bool FGAIThermal::init(bool search_in_AI_path) { 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 = 1; return FGAIBase::init(search_in_AI_path); } void FGAIThermal::bind() { + props->tie("position/altitude-agl-ft", // for debug and tweak + SGRawValuePointer(&altitude_agl_ft)); + props->tie("alt-rel", // for debug and tweak + SGRawValuePointer(&alt_rel)); + props->tie("time", // for debug and tweak + SGRawValuePointer(&time)); + props->tie("xx", // for debug and tweak + SGRawValuePointer(&xx)); + props->tie("is-forming", // for debug abd tweak + SGRawValuePointer(&is_forming)); + props->tie("is-formed", // for debug abd tweak + SGRawValuePointer(&is_formed)); + props->tie("is-dying", // for debug abd tweak + SGRawValuePointer(&is_dying)); + props->tie("is-dead", // for debug abd tweak + SGRawValuePointer(&is_dead)); FGAIBase::bind(); } void FGAIThermal::unbind() { + props->untie("position/altitude-agl-ft"); + props->untie("alt-rel"); + props->untie("time"); + props->untie("is-forming"); + props->untie("is-formed"); + props->untie("is-dying"); + props->untie("is-dead"); + props->untie("xx"); FGAIBase::unbind(); } @@ -75,41 +114,301 @@ void FGAIThermal::update(double dt) { } + +//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; +if ( alt_frac <=0.0 ) { // do submarines get thermals ? + fac = 0.0; + return fac; + } +else if ( ( alt_frac>0.0 ) && (alt_frac<=0.1) ) { // ground layer + fac = ( 0.1*( pow( (10.0*alt_frac),10.0) ) ); + return fac; + } +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) ) ; + return fac; + } +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 ) ) ) ); + return fac; + } +else if ( alt_frac >= 1.1 ) { //above the cloud + fac = 0.0; + return fac; + } +} + + void FGAIThermal::Run(double dt) { - //###########################// - // do calculations for range // - //###########################// - - // copy values from the AIManager - double user_latitude = manager->get_user_latitude(); - double user_longitude = manager->get_user_longitude(); - double user_altitude = manager->get_user_altitude(); - - // calculate range to target in feet and nautical miles - double lat_range = fabs(pos.getLatitudeDeg() - user_latitude) * ft_per_deg_lat; - double lon_range = fabs(pos.getLongitudeDeg() - user_longitude) * ft_per_deg_lon; - double range_ft = sqrt(lat_range*lat_range + lon_range*lon_range); - range = range_ft / 6076.11549; - - // Calculate speed of rising air if within range. - // Air vertical speed is maximum at center of thermal, - // and decreases to zero at the edge (as distance cubed). - if (range < (diameter * 0.5)) { - strength = max_strength - ( range * range * range * factor ); - } else { - strength = 0.0; - } - - // Stop lift at the top of the thermal (smoothly) - if (user_altitude > (height + 100.0)) { - strength = 0.0; - } - else if (user_altitude < height) { - // do nothing - } - else { - strength -= (strength * (user_altitude - height) * 0.01); - } +// ***************************************** +// 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; + + + +// ************************************** +// various variables relative to the user +// ************************************** + +double user_latitude = manager->get_user_latitude(); +double user_longitude = manager->get_user_longitude(); +double user_altitude = manager->get_user_altitude(); // MSL + +//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 ) { + //double alt; + if (globals->get_scenery()->get_elevation_m(SGGeod::fromGeodM(pos, 20000), alt, 0)){ + ground_elev_ft = alt * SG_METER_TO_FEET; + do_agl_calc = 0; + altitude_agl_ft = height - ground_elev_ft ; + dt_count = 0.0; + } } +//user altitude relative to the thermal height, seen AGL from the thermal foot +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 windspeed at the user alt in KT +double windspeed; + +// the wind heading at the user alt +double wind_heading; +double wind_heading_deg; +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=1;is_formed=0;is_dying=0;is_dead=0; + + } +else if ( (time > t1) && (time <= t2) ) { + maxstrengthavail = get_strength_fac ( (alt_rel) ); + + is_forming=0;is_formed=1;is_dying=0;is_dead=0; + + } +else if ( (time > t2) && (time <= t3) ) { + xx= ( ( time - t2) / (1.0 - t2) ) ; + maxstrengthavail = get_strength_fac ( alt_rel - xx ); + + is_forming=0;is_formed=0;is_dying=1;is_dead=0; + + } +else { + maxstrengthavail = 0.0; + is_forming=0;is_formed=0;is_dying=0;is_dead=1; + + } + +//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; + +double dist_center_lon = fabs(slice_center_lon - user_longitude)* ft_per_deg_lon; +double dist_center_lat = fabs(slice_center_lat - user_latitude)* ft_per_deg_lat; + +double dist_center_FT = sqrt ( dist_center_lon*dist_center_lon + dist_center_lat*dist_center_lat ); // feet + +dist_center = dist_center_FT/ 6076.11549; //nautic miles + + +// 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; + +} + + + diff --git a/src/AIModel/AIThermal.hxx b/src/AIModel/AIThermal.hxx index 731ac81ce..675a8aec7 100644 --- a/src/AIModel/AIThermal.hxx +++ b/src/AIModel/AIThermal.hxx @@ -1,8 +1,10 @@ -// FGAIThermal - AIBase derived class creates an AI thunderstorm +// FGAIThermal - FGAIBase-derived class creates an AI thermal // -// Written by David Culp, started Feb 2004. +// Original by Written by David Culp // -// Copyright (C) 2004 David P. Culp - davidculp2@comcast.net +// An attempt to refine the thermal shape and behaviour by WooT 2009 +// +// Copyright (C) 2009 Patrice Poly ( WooT ) // // This program is free software; you can redistribute it and/or // modify it under the terms of the GNU General Public License as @@ -42,22 +44,54 @@ public: virtual void unbind(); virtual void update(double dt); - inline void setMaxStrength( double s ) { max_strength = s; }; + inline void setMaxStrength( double s ) { max_strength = s; }; inline void setDiameter( double d ) { diameter = d; }; inline void setHeight( double h ) { height = h; }; - inline double getStrength() const { return strength; }; + inline void setMaxUpdraft( double lift ) { v_up_max = lift; }; + inline void setMinUpdraft( double sink ) { v_up_min = sink; }; + inline void setR_up_frac( double r ) { r_up_frac = r; }; + + inline double getStrength() const { return strength; }; inline double getDiameter() const { return diameter; }; inline double getHeight() const { return height; }; + inline double getV_up_max() const { return v_up_max; }; + inline double getV_up_min() const { return v_up_min; }; + inline double getR_up_frac() const { return r_up_frac; }; virtual const char* getTypeString(void) const { return "thermal"; } + void getGroundElev(double dt); + + private: void Run(double dt); - double max_strength; + double get_strength_fac(double alt_frac); + double max_strength; double strength; double diameter; double height; double factor; + double alt_rel; + double alt; + double v_up_max; + double v_up_min; + double r_up_frac; + double cycle_timer; + double dt_count; + double time; + double xx; + double ground_elev_ft; // ground level in ft + double altitude_agl_ft; // altitude above ground in feet + bool do_agl_calc; + bool is_forming; + bool is_formed; + bool is_dying; + bool is_dead; + SGPropertyNode_ptr _surface_wind_from_deg_node; + SGPropertyNode_ptr _surface_wind_speed_node; + SGPropertyNode_ptr _aloft_wind_from_deg_node; + SGPropertyNode_ptr _aloft_wind_speed_node; + }; diff --git a/src/Environment/Makefile.am b/src/Environment/Makefile.am index 0debf38ac..63b40f6ec 100644 --- a/src/Environment/Makefile.am +++ b/src/Environment/Makefile.am @@ -10,6 +10,7 @@ libEnvironment_a_SOURCES = \ environment_ctrl.cxx environment_ctrl.hxx \ fgmetar.cxx fgmetar.hxx fgclouds.cxx fgclouds.hxx \ atmosphere.cxx atmosphere.hxx \ - precipitation_mgr.cxx precipitation_mgr.hxx + precipitation_mgr.cxx precipitation_mgr.hxx \ + ridge_lift.cxx ridge_lift.hxx INCLUDES = -I$(top_srcdir) -I$(top_srcdir)/src diff --git a/src/Environment/environment.cxx b/src/Environment/environment.cxx index dfb8824ac..7a516b139 100644 --- a/src/Environment/environment.cxx +++ b/src/Environment/environment.cxx @@ -133,6 +133,8 @@ void FGEnvironment::_init() wind_from_north_fps = 0; wind_from_east_fps = 0; wind_from_down_fps = 0; + thermal_lift_fps = 0; + ridge_lift_fps= 0; altitude_half_to_sun_m = 1000; altitude_tropo_top_m = 10000; _setup_tables(); @@ -171,6 +173,8 @@ FGEnvironment::copy (const FGEnvironment &env) wind_from_north_fps = env.wind_from_north_fps; wind_from_east_fps = env.wind_from_east_fps; wind_from_down_fps = env.wind_from_down_fps; + thermal_lift_fps = env.thermal_lift_fps; + ridge_lift_fps= env.ridge_lift_fps; turbulence_magnitude_norm = env.turbulence_magnitude_norm; turbulence_rate_hz = env.turbulence_rate_hz; } @@ -337,6 +341,18 @@ FGEnvironment::get_wind_from_down_fps () const return wind_from_down_fps; } +double +FGEnvironment::get_thermal_lift_fps () const +{ + return thermal_lift_fps; +} + +double +FGEnvironment::get_ridge_lift_fps () const +{ + return ridge_lift_fps; +} + double FGEnvironment::get_turbulence_magnitude_norm () const { @@ -453,6 +469,20 @@ FGEnvironment::set_wind_from_down_fps (double d) _recalc_hdgspd(); } +void +FGEnvironment::set_thermal_lift_fps (double th) +{ + thermal_lift_fps = th; + _recalc_updraft(); +} + +void +FGEnvironment::set_ridge_lift_fps (double ri) +{ + ridge_lift_fps = ri; + _recalc_updraft(); +} + void FGEnvironment::set_turbulence_magnitude_norm (double t) { @@ -536,6 +566,12 @@ FGEnvironment::_recalc_ne () sin(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS); } +void +FGEnvironment::_recalc_updraft () +{ + wind_from_down_fps = thermal_lift_fps + ridge_lift_fps ; +} + void FGEnvironment::_recalc_sl_temperature () { diff --git a/src/Environment/environment.hxx b/src/Environment/environment.hxx index 4d82685aa..c4dc8fb84 100644 --- a/src/Environment/environment.hxx +++ b/src/Environment/environment.hxx @@ -70,6 +70,8 @@ public: virtual double get_wind_from_north_fps () const; virtual double get_wind_from_east_fps () const; virtual double get_wind_from_down_fps () const; + virtual double get_thermal_lift_fps () const; + virtual double get_ridge_lift_fps () const; virtual double get_turbulence_magnitude_norm () const; virtual double get_turbulence_rate_hz () const; @@ -88,6 +90,8 @@ public: virtual void set_wind_from_north_fps (double n); virtual void set_wind_from_east_fps (double e); virtual void set_wind_from_down_fps (double d); + virtual void set_thermal_lift_fps (double th); + virtual void set_ridge_lift_fps (double ri); virtual void set_turbulence_magnitude_norm (double t); virtual void set_turbulence_rate_hz (double t); @@ -101,6 +105,7 @@ private: void _init(); void _recalc_hdgspd (); void _recalc_ne (); + void _recalc_updraft (); void _recalc_sl_temperature (); void _recalc_alt_temperature (); @@ -139,6 +144,8 @@ private: double wind_from_north_fps; double wind_from_east_fps; double wind_from_down_fps; + double thermal_lift_fps; + double ridge_lift_fps; }; diff --git a/src/Environment/environment_mgr.cxx b/src/Environment/environment_mgr.cxx index 0ff1abb86..7b4f91e1f 100644 --- a/src/Environment/environment_mgr.cxx +++ b/src/Environment/environment_mgr.cxx @@ -145,6 +145,16 @@ FGEnvironmentMgr::bind () &FGEnvironment::get_wind_from_down_fps, &FGEnvironment::set_wind_from_down_fps); fgSetArchivable("/environment/wind-from-down-fps"); + + fgTie("/environment/thermal-lift-fps", _environment, + &FGEnvironment::get_thermal_lift_fps, + &FGEnvironment::set_thermal_lift_fps); + fgSetArchivable("/environment/thermal-lift-fps"); + fgTie("/environment/ridge-lift-fps", _environment, + &FGEnvironment::get_ridge_lift_fps, + &FGEnvironment::set_ridge_lift_fps); + fgSetArchivable("/environment/ridge-lift-fps"); + fgTie("/environment/turbulence/magnitude-norm", _environment, &FGEnvironment::get_turbulence_magnitude_norm, &FGEnvironment::set_turbulence_magnitude_norm); @@ -227,6 +237,10 @@ FGEnvironmentMgr::unbind () fgUntie("/environment/wind-from-north-fps"); fgUntie("/environment/wind-from-east-fps"); fgUntie("/environment/wind-from-down-fps"); + + fgUntie("/environment/thermal-lift-fps"); + fgUntie("/environment/ridge-lift-fps"); + fgUntie("/environment/atmosphere/altitude-half-to-sun"); fgUntie("/environment/atmosphere/altitude-troposphere-top"); for (int i = 0; i < MAX_CLOUD_LAYERS; i++) { diff --git a/src/Environment/ridge_lift.cxx b/src/Environment/ridge_lift.cxx new file mode 100644 index 000000000..011fe006f --- /dev/null +++ b/src/Environment/ridge_lift.cxx @@ -0,0 +1,317 @@ +// simulates ridge lift +// +// Written by Patrice Poly +// Copyright (C) 2009 Patrice Poly - p.polypa@gmail.com +// +// +// Entirely based on the paper : +// http://carrier.csi.cam.ac.uk/forsterlewis/soaring/sim/fsx/dev/sim_probe/sim_probe_paper.html +// by Ian Forster-Lewis, University of Cambridge, 26th December 2007 +// +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// +// + + +#ifdef HAVE_CONFIG_H +# include +#endif + +#include
+#include
+#include +#include +#include + +using std::string; + +#include "ridge_lift.hxx" + +//constructor +FGRidgeLift::FGRidgeLift () +{ + dist_probe_m[0] = 0.0; // in meters + dist_probe_m[1] = 250.0; + dist_probe_m[2] = 750.0; + dist_probe_m[3] = 2000.0; + dist_probe_m[4] = -100.0; + + BOUNDARY1_m = 40.0; // in meters + BOUNDARY2_m = 130.0; + + PI = 4.0*atan(1.0); // pi + deg2rad = (PI/180.0); + rad2deg = (180.0/PI); + strength = 0.0; + timer = 0.0; + scanned = false; + + earth_rad_ft=20899773.07; +} + +//destructor +FGRidgeLift::~FGRidgeLift() +{ + +} + +void FGRidgeLift::init(void) +{ + _ridge_lift_fps_node = + fgGetNode("/environment/ridge-lift-fps" + , true); + _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); + _earth_radius_node = + fgGetNode("/position/sea-level-radius-ft" + , true); + _user_longitude_node = + fgGetNode("/position/longitude-deg" + , true); + _user_latitude_node = + fgGetNode("/position/latitude-deg" + , true); + _user_altitude_ft_node = + fgGetNode("/position/altitude-ft" + , true); + _user_altitude_agl_ft_node = + fgGetNode("/position/altitude-agl-ft" + , true); +} + +void FGRidgeLift::bind() { + + fgTie("/environment/ridge-lift/probe-elev-m[0]", this, + &FGRidgeLift::get_probe_elev_m_0); // read-only + fgTie("/environment/ridge-lift/probe-elev-m[1]", this, + &FGRidgeLift::get_probe_elev_m_1); // read-only + fgTie("/environment/ridge-lift/probe-elev-m[2]", this, + &FGRidgeLift::get_probe_elev_m_2); // read-only + fgTie("/environment/ridge-lift/probe-elev-m[3]", this, + &FGRidgeLift::get_probe_elev_m_3); // read-only + fgTie("/environment/ridge-lift/probe-elev-m[4]", this, + &FGRidgeLift::get_probe_elev_m_4); // read-only + + fgTie("/environment/ridge-lift/probe-lat-deg[0]", this, + &FGRidgeLift::get_probe_lat_0); // read-only + fgTie("/environment/ridge-lift/probe-lat-deg[1]", this, + &FGRidgeLift::get_probe_lat_1); // read-only + fgTie("/environment/ridge-lift/probe-lat-deg[2]", this, + &FGRidgeLift::get_probe_lat_2); // read-only + fgTie("/environment/ridge-lift/probe-lat-deg[3]", this, + &FGRidgeLift::get_probe_lat_3); // read-only + fgTie("/environment/ridge-lift/probe-lat-deg[4]", this, + &FGRidgeLift::get_probe_lat_4); // read-only + + fgTie("/environment/ridge-lift/probe-lon-deg[0]", this, + &FGRidgeLift::get_probe_lon_0); // read-only + fgTie("/environment/ridge-lift/probe-lon-deg[1]", this, + &FGRidgeLift::get_probe_lon_1); // read-only + fgTie("/environment/ridge-lift/probe-lon-deg[2]", this, + &FGRidgeLift::get_probe_lon_2); // read-only + fgTie("/environment/ridge-lift/probe-lon-deg[3]", this, + &FGRidgeLift::get_probe_lon_3); // read-only + fgTie("/environment/ridge-lift/probe-lon-deg[4]", this, + &FGRidgeLift::get_probe_lon_4); // read-only + + fgTie("/environment/ridge-lift/slope[0]", this, + &FGRidgeLift::get_slope_0); // read-only + fgTie("/environment/ridge-lift/slope[1]", this, + &FGRidgeLift::get_slope_1); // read-only + fgTie("/environment/ridge-lift/slope[2]", this, + &FGRidgeLift::get_slope_2); // read-only + fgTie("/environment/ridge-lift/slope[3]", this, + &FGRidgeLift::get_slope_3); // read-only + +} + +void FGRidgeLift::unbind() { + + fgUntie("/environment/ridge-lift/probe-elev-m[0]"); + fgUntie("/environment/ridge-lift/probe-elev-m[1]"); + fgUntie("/environment/ridge-lift/probe-elev-m[2]"); + fgUntie("/environment/ridge-lift/probe-elev-m[3]"); + fgUntie("/environment/ridge-lift/probe-elev-m[4]"); + + fgUntie("/environment/ridge-lift/probe-lat-deg[0]"); + fgUntie("/environment/ridge-lift/probe-lat-deg[1]"); + fgUntie("/environment/ridge-lift/probe-lat-deg[2]"); + fgUntie("/environment/ridge-lift/probe-lat-deg[3]"); + fgUntie("/environment/ridge-lift/probe-lat-deg[4]"); + + fgUntie("/environment/ridge-lift/probe-lon-deg[0]"); + fgUntie("/environment/ridge-lift/probe-lon-deg[1]"); + fgUntie("/environment/ridge-lift/probe-lon-deg[2]"); + fgUntie("/environment/ridge-lift/probe-lon-deg[3]"); + fgUntie("/environment/ridge-lift/probe-lon-deg[4]"); + + fgUntie("/environment/ridge-lift/slope[0]"); + fgUntie("/environment/ridge-lift/slope[1]"); + fgUntie("/environment/ridge-lift/slope[2]"); + fgUntie("/environment/ridge-lift/slope[3]"); + +} + +void FGRidgeLift::update(double dt) { + Run(dt); +} + +double FGRidgeLift::sign(double x) { + if (x == 0.0) + return x; + else + return x/fabs(x); +} + +void FGRidgeLift::Run(double dt) { + + // copy values + + user_latitude_deg = _user_latitude_node->getDoubleValue(); + user_longitude_deg = _user_longitude_node->getDoubleValue(); + //user_altitude_ft = _user_altitude_ft_node->getDoubleValue(); + + if ( ( _earth_radius_node->getDoubleValue() ) > 1.0 ) { + earth_rad_ft =_earth_radius_node->getDoubleValue(); } + else { earth_rad_ft=20899773.07; } + + //earth_rad_m = earth_rad_ft * 0.3048 ; + earth_rad_m = earth_rad_ft * SG_FEET_TO_METER ; + + //get the windspeed at ground level + + double ground_wind_from_deg = _surface_wind_from_deg_node->getDoubleValue(); + double ground_wind_speed_kts = _surface_wind_speed_node->getDoubleValue(); + //double ground_wind_speed_mps = ground_wind_speed_kts / SG_METER_TO_FEET; + double ground_wind_speed_mps = ground_wind_speed_kts / 3.2808399; + + double ground_wind_from_rad = (user_longitude_deg < 0.0) ? + PI*( ground_wind_from_deg/180.0) +PI : PI*( ground_wind_from_deg/180.0); + + // Placing the probes + + for (int i = 0; i <= 4; i++) + { + probe_lat_rad[i] = asin(sin(deg2rad*user_latitude_deg)*cos(dist_probe_m[i]/earth_rad_m) + +cos(deg2rad*user_latitude_deg)*sin(dist_probe_m[i]/earth_rad_m)*cos(ground_wind_from_rad)); + if (probe_lat_rad[i] == 0.0) { + probe_lon_rad[i] = (deg2rad*user_latitude_deg); // probe on a pole + } + else { + probe_lon_rad[i] = fmod((deg2rad*user_longitude_deg+asin(sin(ground_wind_from_rad) + *sin(dist_probe_m[i]/earth_rad_m)/cos(probe_lon_rad[i]))+PI) + ,(2.0*PI))-PI; + } + probe_lat_deg[i]= rad2deg*probe_lat_rad[i]; + probe_lon_deg[i]= rad2deg*probe_lon_rad[i]; + } + + // ground elevations + // every second + + timer += dt; + if (timer >= 1.0 ) { + scanned = true; + for (int i = 0; i <= 4; i++) + { + if (globals->get_scenery()->get_elevation_m(SGGeod::fromGeodM( + SGGeod::fromRad(probe_lon_rad[i],probe_lat_rad[i]), 20000), alt, 0)); + { + probe_elev_m[i] = alt; + } + } + timer = 0.0; + + } + + // slopes + + double adj_slope[5]; + + slope[0] = (probe_elev_m[0] - probe_elev_m[1]) / dist_probe_m[1]; + slope[1] = (probe_elev_m[1] - probe_elev_m[2]) / dist_probe_m[2]; + slope[2] = (probe_elev_m[2] - probe_elev_m[3]) / dist_probe_m[3]; + slope[3] = (probe_elev_m[4] - probe_elev_m[0]) / -dist_probe_m[4]; + + for (int i = 0; i <= 4; i++) + { + adj_slope[i] = sin(atan(5.0 * pow ( (abs(slope[i])),1.7) ) ) *sign(slope[i]); + } + + //adjustment + + adj_slope[0] = 0.2 * adj_slope[0]; + adj_slope[1] = 0.2 * adj_slope[1]; + if ( adj_slope [2] < 0.0 ) + { + adj_slope[2] = 0.5 * adj_slope[2]; + } + else + { + adj_slope[2] = 0.0 ; + } + + if ( ( adj_slope [0] >= 0.0 ) && ( adj_slope [3] < 0.0 ) ) + { + adj_slope[3] = 0.0; + } + else + { + adj_slope[3] = 0.2 * adj_slope[3]; + } + + double lift_factor = adj_slope[0]+adj_slope[1]+adj_slope[2]+adj_slope[3]; + + //user altitude above ground + + user_altitude_agl_ft = _user_altitude_agl_ft_node->getDoubleValue(); + user_altitude_agl_m = ( user_altitude_agl_ft / SG_METER_TO_FEET ); + + //boundaries + double agl_factor; + + if ( user_altitude_agl_m < BOUNDARY1_m ) + { + agl_factor = 0.5+0.5*user_altitude_agl_m /BOUNDARY1_m ; + } + else if ( user_altitude_agl_m < BOUNDARY2_m ) + { + agl_factor = 1.0; + } + else + { + agl_factor = exp(-(2 + 2 * probe_elev_m[0] / 4000) * + (user_altitude_agl_m - BOUNDARY2_m) / max(probe_elev_m[0],200.0)); + } + + double lift_mps = lift_factor* ground_wind_speed_mps * agl_factor; + + //the updraft, finally, in ft per second + strength = lift_mps * SG_METER_TO_FEET ; + + _ridge_lift_fps_node->setDoubleValue( strength ); + +} + + + + + + diff --git a/src/Environment/ridge_lift.hxx b/src/Environment/ridge_lift.hxx new file mode 100644 index 000000000..c6f49ecf2 --- /dev/null +++ b/src/Environment/ridge_lift.hxx @@ -0,0 +1,125 @@ +// simulates ridge lift +// +// Written by Patrice Poly +// Copyright (C) 2009 Patrice Poly - p.polypa@gmail.com +// +// +// Entirely based on the paper : +// http://carrier.csi.cam.ac.uk/forsterlewis/soaring/sim/fsx/dev/sim_probe/sim_probe_paper.html +// by Ian Forster-Lewis, University of Cambridge, 26th December 2007 +// +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// + + +#ifndef _FG_RidgeLift_HXX +#define _FG_RidgeLift_HXX + +#ifdef HAVE_CONFIG +# include +#endif + + +#include +using std::string; + + +class FGRidgeLift : public SGSubsystem { +public: + + FGRidgeLift(); + ~FGRidgeLift(); + + virtual void bind(); + virtual void unbind(); + virtual void update(double dt); + virtual void init(); + + inline double getStrength() const { return strength; }; + + inline double get_probe_elev_m_0() const { return probe_elev_m[0]; }; + inline double get_probe_elev_m_1() const { return probe_elev_m[1]; }; + inline double get_probe_elev_m_2() const { return probe_elev_m[2]; }; + inline double get_probe_elev_m_3() const { return probe_elev_m[3]; }; + inline double get_probe_elev_m_4() const { return probe_elev_m[4]; }; + + inline double get_probe_lat_0() const { return probe_lat_deg[0]; }; + inline double get_probe_lat_1() const { return probe_lat_deg[1]; }; + inline double get_probe_lat_2() const { return probe_lat_deg[2]; }; + inline double get_probe_lat_3() const { return probe_lat_deg[3]; }; + inline double get_probe_lat_4() const { return probe_lat_deg[4]; }; + + inline double get_probe_lon_0() const { return probe_lon_deg[0]; }; + inline double get_probe_lon_1() const { return probe_lon_deg[1]; }; + inline double get_probe_lon_2() const { return probe_lon_deg[2]; }; + inline double get_probe_lon_3() const { return probe_lon_deg[3]; }; + inline double get_probe_lon_4() const { return probe_lon_deg[4]; }; + + inline double get_slope_0() const { return slope[0]; }; + inline double get_slope_1() const { return slope[1]; }; + inline double get_slope_2() const { return slope[2]; }; + inline double get_slope_3() const { return slope[3]; }; + +private: + //void init(); + void Run(double dt); + + double dist_probe_m[5]; + double BOUNDARY1_m; + double BOUNDARY2_m; + double PI; // pi + double deg2rad; + double rad2deg; + + bool scanned; + + double strength; + double timer; + + double probe_lat_rad[5]; + double probe_lon_rad[5]; + + double probe_lat_deg[5]; + double probe_lon_deg[5]; + + double alt; + double probe_elev_m[5]; + + double slope[4]; + double earth_rad_ft; + double earth_rad_m; + double user_latitude_deg; + double user_longitude_deg; + //double user_altitude; + double user_altitude_agl_ft; + double user_altitude_agl_m; + + double sign(double x); + + SGPropertyNode_ptr _ridge_lift_fps_node; + + SGPropertyNode_ptr _surface_wind_from_deg_node; + SGPropertyNode_ptr _surface_wind_speed_node; + + SGPropertyNode_ptr _user_altitude_ft_node; + SGPropertyNode_ptr _user_altitude_agl_ft_node; + SGPropertyNode_ptr _earth_radius_node; + SGPropertyNode_ptr _user_longitude_node; + SGPropertyNode_ptr _user_latitude_node; + +}; + +#endif // _FG_RidgeLift_HXX \ No newline at end of file diff --git a/src/Main/fg_init.cxx b/src/Main/fg_init.cxx index 97c8e9821..164446cb4 100644 --- a/src/Main/fg_init.cxx +++ b/src/Main/fg_init.cxx @@ -115,6 +115,7 @@ #include #include +#include #include "fg_init.hxx" #include "fg_io.hxx" @@ -1508,6 +1509,13 @@ bool fgInitSubsystems() { // Initialize the weather modeling subsystem globals->add_subsystem("environment", new FGEnvironmentMgr); + //////////////////////////////////////////////////////////////////// + // Initialize the ridge lift simulation. + //////////////////////////////////////////////////////////////////// + + // Initialize the ridgelift subsystem + globals->add_subsystem("ridgelift", new FGRidgeLift); + //////////////////////////////////////////////////////////////////// // Initialize the aircraft systems and instrumentation (before the -- 2.39.5