Make sure sign is preserved for turbulence direction.

diff --git a/src/FDM/JSBSim/FGAtmosphere.cpp b/src/FDM/JSBSim/FGAtmosphere.cpp
index b8098da0b588523546207cc4d8e3a06a6e4ea10c..bf02c1bc0a8f69da1f52a8c666c044e39a6910d8 100644 (file)
--- a/src/FDM/JSBSim/FGAtmosphere.cpp
+++ b/src/FDM/JSBSim/FGAtmosphere.cpp
@@ -260,6 +260,16 @@ void FGAtmosphere::Calculate(double altitude)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+// square a value, but preserve the original sign
+static inline double 
+square_signed (double value)
+{
+    if (value < 0)
+        return value * value * -1;
+    else
+        return value * value;
+}
+
 void FGAtmosphere::Turbulence(void)
 {
   switch (turbType) {
@@ -277,8 +287,11 @@ void FGAtmosphere::Turbulence(void)
     Magnitude         += MagnitudeAccel*rate*State->Getdt();
 
     vDirectiondAccelDt.Normalize();
-    vDirectiondAccelDt(eX) *= vDirectiondAccelDt(eX);
-    vDirectiondAccelDt(eY) *= vDirectiondAccelDt(eY);
+
+                                // deemphasise non-vertical forces
+    vDirectiondAccelDt(eX) = square_signed(vDirectiondAccelDt(eX));
+    vDirectiondAccelDt(eY) = square_signed(vDirectiondAccelDt(eY));
+
     vDirectionAccel += vDirectiondAccelDt*rate*TurbRate*State->Getdt();
     vDirectionAccel.Normalize();
     vDirection      += vDirectionAccel*rate*State->Getdt();