new temperature approximation

openap/casadi/aero_override.py

@@ -40,19 +40,11 @@ def atmos(h, dT=0):
             Air pressure (Pa), density (kg/m3), and temperature (K).
 
     """
-    dT = np.maximum(-15, np.minimum(dT, 15))
-
-    # ISA model
-    # T0_shift = T0 + dT
-    # T = np.maximum(T0_shift + beta * h, 216.65 + dT)
-    # rhotrop = rho0 * (T / T0_shift) ** 4.256848030018761
-    # dhstrat = np.maximum(0.0, h - 11000.0)
-    # rho = rhotrop * np.exp(-dhstrat / 6341.552161)
-
-    # use exponential model to avoid discontinuity at tropopause
-    T = 85.46369268 * np.exp(-0.00017235 * h) + 213.31449979 + dT
-    rho = 1.31788377 * np.exp(-0.00011107 * h) - 0.03933069
+    assert -15 < dT < 15
 
+    T0_ = T0 + dT
+    T = 0.65 * np.log(1 + np.exp(-10 * (h / 1000 - 11))) + 216.65
+    rho = rho0 * ((T0_ + beta * h) / T0_) ** 4.2559
     p = rho * R * T
     return p, rho, T