diff --git a/src/FDM/YASim/Model.cpp b/src/FDM/YASim/Model.cpp index c91fe7df6..24c2940c0 100644 --- a/src/FDM/YASim/Model.cpp +++ b/src/FDM/YASim/Model.cpp @@ -79,13 +79,18 @@ void Model::initIteration() int i; for(i=0; i<3; i++) _gyro[i] = _torque[i] = 0; + + // Need a local altitude for the wind calculation + float dummy[3]; + float alt = Math::abs(localGround(_s, dummy)); + for(i=0; i<_thrusters.size(); i++) { Thruster* t = (Thruster*)_thrusters.get(i); // Get the wind velocity at the thruster location float pos[3], v[3]; t->getPosition(pos); - localWind(pos, _s, v); + localWind(pos, _s, v, alt); t->setWind(v); t->setAir(_pressure, _temp, _rho); @@ -251,8 +256,7 @@ void Model::setGroundEffect(float* pos, float span, float mul) // (v dot _ground)-_ground[3] gives the distance AGL. void Model::setGroundPlane(double* planeNormal, double fromOrigin) { - int i; - for(i=0; i<3; i++) _ground[i] = planeNormal[i]; + for(int i=0; i<3; i++) _ground[i] = planeNormal[i]; _ground[3] = fromOrigin; } @@ -286,6 +290,14 @@ void Model::calcForces(State* s) _body.addForce(pos, thrust); } + // Get a ground plane in local coordinates. The first three + // elements are the normal vector, the final one is the distance + // from the local origin along that vector to the ground plane + // (negative for objects "above" the ground) + float ground[4]; + ground[3] = localGround(s, ground); + float alt = Math::abs(ground[3]); + // Gravity, convert to a force, then to local coordinates float grav[3]; Glue::geodUp(s->pos, grav); @@ -303,7 +315,7 @@ void Model::calcForces(State* s) // Vsurf = wind - velocity + (rot cross (cg - pos)) float vs[3], pos[3]; sf->getPosition(pos); - localWind(pos, s, vs); + localWind(pos, s, vs, alt); float force[3], torque[3]; sf->calcForce(vs, _rho, force, torque); @@ -318,7 +330,7 @@ void Model::calcForces(State* s) // Vsurf = wind - velocity + (rot cross (cg - pos)) float vs[3], pos[3]; sf->getPosition(pos); - localWind(pos, s, vs); + localWind(pos, s, vs, alt); float force[3], torque[3]; sf->calcForce(vs, _rho, force, torque); @@ -335,7 +347,7 @@ void Model::calcForces(State* s) // Vsurf = wind - velocity + (rot cross (cg - pos)) float vs[3], pos[3]; sf->getPosition(pos); - localWind(pos, s, vs); + localWind(pos, s, vs, alt); float force[3], torque[3]; sf->calcForce(vs, _rho, force, torque); @@ -347,13 +359,6 @@ void Model::calcForces(State* s) _body.addTorque(torque); } - // Get a ground plane in local coordinates. The first three - // elements are the normal vector, the final one is the distance - // from the local origin along that vector to the ground plane - // (negative for objects "above" the ground) - float ground[4]; - ground[3] = localGround(s, ground); - // Account for ground effect by multiplying the vertical force // component by an amount linear with the fraction of the wingspan // above the ground. @@ -433,19 +438,20 @@ float Model::localGround(State* s, float* out) // Calculates the airflow direction at the given point and for the // specified aircraft velocity. -void Model::localWind(float* pos, State* s, float* out) +void Model::localWind(float* pos, State* s, float* out, float alt) { float tmp[3], lwind[3], lrot[3], lv[3]; // Get a global coordinate for our local position, and calculate // turbulence. - // FIXME: modify turbulence for altitude, attenuating the vertical - // component near the ground. if(_turb) { - double gpos[3]; + double gpos[3]; float up[3]; Math::tmul33(s->orient, pos, tmp); - for(int i=0; i<3; i++) gpos[i] = s->pos[i] + tmp[i]; - _turb->getTurbulence(gpos, lwind); + for(int i=0; i<3; i++) { + gpos[i] = s->pos[i] + tmp[i]; + up[i] = _ground[i]; + } + _turb->getTurbulence(gpos, alt, up, lwind); Math::add3(_wind, lwind, lwind); } else { Math::set3(_wind, lwind); diff --git a/src/FDM/YASim/Model.hpp b/src/FDM/YASim/Model.hpp index fb5a99881..a27af5c05 100644 --- a/src/FDM/YASim/Model.hpp +++ b/src/FDM/YASim/Model.hpp @@ -74,7 +74,7 @@ private: void calcGearForce(Gear* g, float* v, float* rot, float* ground); float gearFriction(float wgt, float v, Gear* g); float localGround(State* s, float* out); - void localWind(float* pos, State* s, float* out); + void localWind(float* pos, State* s, float* out, float alt); Integrator _integrator; RigidBody _body; diff --git a/src/FDM/YASim/Turbulence.cpp b/src/FDM/YASim/Turbulence.cpp index 8862e5503..da9741f7e 100644 --- a/src/FDM/YASim/Turbulence.cpp +++ b/src/FDM/YASim/Turbulence.cpp @@ -25,7 +25,7 @@ const double MAGNITUDE_EXP = 2.0; // bandwidth to the higher frequency components. A turbulence field // will swing between maximal values over a distance of approximately // 2^(MEANINGFUL_GENS-1). -const int MEANINGFUL_GENS = 8; +const int MEANINGFUL_GENS = 7; static const float FT2M = 0.3048; @@ -130,7 +130,8 @@ void Turbulence::offset(float* offset) _off[i] += offset[i]; } -void Turbulence::getTurbulence(double* loc, float* turbOut) +void Turbulence::getTurbulence(double* loc, float alt, float* up, + float* turbOut) { // Convert to integer 2D coordinates; wrap to [0:_sz]. double a = (loc[0] + _off[0]) + (loc[2] + _off[2]); @@ -158,6 +159,20 @@ void Turbulence::getTurbulence(double* loc, float* turbOut) float avg1 = (1-a)*turb10[i] + a*turb11[i]; turbOut[i] = mag * ((1-b)*avg0 + b*avg1); } + + // Adjust for altitude effects + if(alt < 300) { + float altmul = 0.5 + (1-0.5) * (alt*(1.0/300)); + if(alt < 100) { + float vmul = alt * (1.0/100); + vmul = vmul / altmul; // pre-correct for the pending altmul + float dot = Math::dot3(turbOut, up); + float off[3]; + Math::mul3(dot * (vmul-1), up, off); + Math::add3(turbOut, off, turbOut); + } + Math::mul3(altmul, turbOut, turbOut); + } } // Associates a random number in the range [-1:1] with a given lattice diff --git a/src/FDM/YASim/Turbulence.hpp b/src/FDM/YASim/Turbulence.hpp index 606298dd4..a500245eb 100644 --- a/src/FDM/YASim/Turbulence.hpp +++ b/src/FDM/YASim/Turbulence.hpp @@ -6,7 +6,7 @@ public: ~Turbulence(); void update(double dt, double rate); void setMagnitude(double mag); - void getTurbulence(double* loc, float* turbOut); + void getTurbulence(double* loc, float alt, float* up, float* turbOut); void offset(float* dist); private: