Attenuate turbulence near the ground. Tweak turbulence numbers.

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(int i=0; i<3; i++) _ground[i] = planeNormal[i];
-    for(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];
+        for(int i=0; i<3; i++) {
-        _turb->getTurbulence(gpos, lwind);
+            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);

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;

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

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: