Ron Jensen: extend atmosphere tables to match environment data

diff --git a/src/FDM/YASim/Airplane.cpp b/src/FDM/YASim/Airplane.cpp
index c6b5ae9252f4a8b2eea64349db706308ef9eec01..07d434e09c7d4eff65f0fc755ecd0edf90f94562 100644 (file)
--- a/src/FDM/YASim/Airplane.cpp
+++ b/src/FDM/YASim/Airplane.cpp
@@ -28,6 +28,8 @@ const float STHRESH = 1;
 // oscillate.
 const float SOLVE_TWEAK = 0.3226;
 
+const float GRAV = 9.8f;
+
 Airplane::Airplane()
 {
     _emptyWeight = 0;
@@ -123,7 +125,7 @@ void Airplane::getPilotAccel(float* out)
 
     // Gravity
     Glue::geodUp(s->pos, out);
-    Math::mul3(-9.8f, out, out);
+    Math::mul3(GRAV, out, out);
 
     // The regular acceleration
     float tmp[3];
@@ -593,7 +595,7 @@ void Airplane::compileContactPoints()
     // Give it a spring constant such that at full compression it will
     // hold up 10 times the planes mass.  That's about right.  Yeah.
     float mass = _model.getBody()->getTotalMass();
-    float spring = (1/DIST) * 9.8f * 10.0f * mass;
+    float spring = (1/DIST) * GRAV * 10.0f * mass;
     float damp = 2 * Math::sqrt(spring * mass);
 
     int i;
@@ -785,6 +787,7 @@ void Airplane::setupWeights(bool isApproach)
 
 void Airplane::runCruise()
 {
+    __builtin_printf("runCruise()\n");
     setupState(_cruiseAoA, _cruiseSpeed,_approachGlideAngle, &_cruiseState);
     _model.setState(&_cruiseState);
     _model.setAir(_cruiseP, _cruiseT,
@@ -828,6 +831,7 @@ void Airplane::runCruise()
 
 void Airplane::runApproach()
 {
+    __builtin_printf("runApproach()\n");
     setupState(_approachAoA, _approachSpeed,_approachGlideAngle, &_approachState);
     _model.setState(&_approachState);
     _model.setAir(_approachP, _approachT,
@@ -936,7 +940,7 @@ void Airplane::solve()
        runCruise();
 
        _model.getThrust(tmp);
-        float thrust = tmp[0] + _cruiseWeight * Math::sin(_cruiseGlideAngle) * 9.81;
+        float thrust = tmp[0] + GRAV * _cruiseWeight * Math::sin(_cruiseGlideAngle);
 
        _model.getBody()->getAccel(tmp);
         Math::tmul33(_cruiseState.orient, tmp, tmp);
@@ -977,7 +981,7 @@ void Airplane::solve()
        float pitch1 = tmp[1];
 
        // Now calculate:
-       float awgt = 9.8f * _approachWeight;
+       float awgt = GRAV * _approachWeight;
 
        float dragFactor = thrust / (thrust-xforce);
        float liftFactor = awgt / (awgt+alift);
@@ -985,8 +989,10 @@ void Airplane::solve()
        float tailDelta = -pitch0 * (ARCMIN/(pitch1-pitch0));
 
         // Sanity:
-        if(dragFactor <= 0 || liftFactor <= 0)
+        if(dragFactor <= 0 || liftFactor <= 0) {
+            __builtin_printf("NEGATIVE drag %f lift %f\n", dragFactor, liftFactor);
             break;
+        }
 
         // And the elevator control in the approach.  This works just
         // like the tail incidence computation (it's solving for the
@@ -1005,7 +1011,7 @@ void Airplane::solve()
         // Now apply the values we just computed.  Note that the
         // "minor" variables are deferred until we get the lift/drag
         // numbers in the right ballpark.
-
+        __builtin_printf("Apply drag %f lift %f\n", dragFactor, liftFactor);
        applyDragFactor(dragFactor);
        applyLiftRatio(liftFactor);
 
@@ -1016,6 +1022,8 @@ void Airplane::solve()
            continue;
        }
 
+        __builtin_printf("Apply aoa %f tail %f\n", SOLVE_TWEAK*aoaDelta, SOLVE_TWEAK*tailDelta);
+
        // OK, now we can adjust the minor variables:
        _cruiseAoA += SOLVE_TWEAK*aoaDelta;
        _tailIncidence += SOLVE_TWEAK*tailDelta;
@@ -1032,6 +1040,8 @@ void Airplane::solve()
             if(abs(elevDelta) < STHRESH*0.0001)
                 break;
 
+            __builtin_printf("Apply elev %f\n", SOLVE_TWEAK*elevDelta);
+
             // Otherwise, adjust and do the next iteration
             _approachElevator.val += SOLVE_TWEAK * elevDelta;
             if(abs(_approachElevator.val) > 1) {

diff --git a/src/FDM/YASim/Atmosphere.cpp b/src/FDM/YASim/Atmosphere.cpp
index 92d2fbe56a6e3495762d2c2ffff5c96c42e134e8..5fb7347f94f42e8a5ed3472bb6fc1396c675f8b1 100644 (file)
--- a/src/FDM/YASim/Atmosphere.cpp
+++ b/src/FDM/YASim/Atmosphere.cpp
@@ -7,8 +7,10 @@ namespace yasim {
 // McCormick lists 299.16/101325/1.22500, but those don't agree with
 // R=287.  I chose to correct the temperature to 288.20, since 79F is
 // pretty hot for a "standard" atmosphere.
+// Numbers above 19000 meters calculated from src/Environment/environment.cxx
 //                             meters   kelvin      Pa   kg/m^3
-float Atmosphere::data[][4] = {{ 0.0f,     288.20f, 101325.0f, 1.22500f },
+float Atmosphere::data[][4] = {{ -900.0f, 293.91f, 111679.0f, 1.32353f },
+                               {    0.0f, 288.11f, 101325.0f, 1.22500f },
                               {   900.0f, 282.31f,  90971.0f, 1.12260f },
                               {  1800.0f, 276.46f,  81494.0f, 1.02690f },
                               {  2700.0f, 270.62f,  72835.0f, 0.93765f },
@@ -29,7 +31,20 @@ float Atmosphere::data[][4] = {{ 0.0f,     288.20f, 101325.0f, 1.22500f },
                               { 16200.0f, 216.66f,  10033.0f, 0.16133f },
                               { 17100.0f, 216.66f,   8712.0f, 0.14009f },
                               { 18000.0f, 216.66f,   7565.0f, 0.12165f },
-                              { 18900.0f, 216.66f,   6570.0f, 0.10564f }};
+                               {18900.0f, 216.66f,   6570.0f, 0.10564f },
+                               {19812.0f, 216.66f,   5644.0f, 0.09073f },
+                               {20726.0f, 217.23f,   4884.0f, 0.07831f },
+                               {21641.0f, 218.39f,   4235.0f, 0.06755f },
+                               {22555.0f, 219.25f,   3668.0f, 0.05827f },
+                               {23470.0f, 220.12f,   3182.0f, 0.05035f },
+                               {24384.0f, 220.98f,   2766.0f, 0.04360f },
+                               {25298.0f, 221.84f,   2401.0f, 0.03770f },
+                               {26213.0f, 222.71f,   2087.0f, 0.03265f },
+                               {27127.0f, 223.86f,   1814.0f, 0.02822f },
+                               {28042.0f, 224.73f,   1581.0f, 0.02450f },
+                               {28956.0f, 225.59f,   1368.0f, 0.02112f },
+                               {29870.0f, 226.45f,   1196.0f, 0.01839f },
+                               {30785.0f, 227.32f,   1044.0f, 0.01599f }};
 
 // Universal gas constant for air, in SI units.  P = R * rho * T.
 // P in pascals (N/m^2), rho is kg/m^3, T in kelvin.