Initial implementation of moist adiabat plotting (it sucks)

lib/xmet/skew_t.py

@@ -4,7 +4,7 @@ import cairo
 from typing import Callable
 
 from xmet.sounding import Sounding
-from xmet.thermo   import pressure_height, lapse
+from xmet.thermo   import pressure_height, lapse, moist_lapse_rate
 
 PRESSURE_MAX  = 1050 # millibar
 PRESSURE_MIN  =  100
@@ -104,26 +104,42 @@ class SkewTGraph():
                      x: float,
                      y: float,
                      start_temp: float,
-                     start_pressure: float):
+                     start_pressure: float,
+                     moist: bool=False):
         start_height = pressure_height(start_pressure)
-        sx_last = None
-        sy_last = None
+        sx_last      = None
+        sy_last      = None
+        height_last  = None
 
+        temp     = start_temp
         pressure = start_pressure
 
         while pressure >= PRESSURE_MIN:
-            height   = pressure_height(pressure)
-            temp_cur = lapse(start_temp, height - start_height)
+            height = pressure_height(pressure)
 
-            sx, sy = self.sample_to_screen(temp_cur, pressure)
+            if height_last is None:
+                height_last = height
+
+            if moist:
+                try:
+                    rate = moist_lapse_rate(temp, temp, pressure)
+                except OverflowError:
+                    break
+
+                temp = lapse(temp, height - height_last, rate)
+            else:
+                temp = lapse(start_temp, height - start_height)
+
+            sx, sy = self.sample_to_screen(temp, pressure)
 
             if sx_last is not None:
                 cr.move_to(x + sx, y + sy)
                 cr.line_to(x + sx_last, y + sy_last)
                 cr.stroke()
 
-            sx_last = sx
-            sy_last = sy
+            sx_last     = sx
+            sy_last     = sy
+            height_last = height
 
             pressure -= 10.0
 
@@ -137,6 +153,9 @@ class SkewTGraph():
         for temp in range(-140, 140, 10):
             self.draw_adiabat(cr, x, y, temp, PRESSURE_MAX)
 
+        for temp in range(-140, 140, 10):
+            self.draw_adiabat(cr, x, y, temp, PRESSURE_MAX, True)
+
         cr.restore()
 
     def draw_sounding(self,

lib/xmet/thermo.py

@@ -1,22 +1,21 @@
 LAPSE_RATE_DRY   = 9.8 / 1000 # degrees C per km
 LAPSE_RATE_MOIST = 4.0 / 1000
 
+def kelvin(c: float) -> float:
+    return 273.15 + c
+
 def vapor_pressure(dewpoint: float) -> float:
     """
     Return the pressure of vapor in a parcel of a given dewpoint.
     """
-    return 6.11 * 10 ** (
-        (7.5 * dewpoint) / (237.3 * dewpoint)
-    )
+    return 6.11 * 10 ** ((7.5 * dewpoint) / (237.3 + dewpoint))
 
 def saturated_vapor_pressure(temp: float) -> float:
     """
     Return the pressure at which the vapor and condensation of a parcel of a
     given temperature are at equilibrium.
     """
-    return 6.11 * 10 ** (
-        (7.5 * temp) / (237.3 * temp)
-    )
+    return 6.11 * 10 ** ((7.5 * temp) / (237.3 + temp))
 
 def mixing_ratio(dewpoint: float, pressure: float) -> float:
     """
@@ -56,3 +55,15 @@ def lapse(temp: float, delta: float, rate=LAPSE_RATE_DRY) -> float:
     given increase in height (in meters).
     """
     return temp - (rate * delta)
+
+def moist_lapse_rate(temp: float, dewpoint: float, pressure: float) -> float:
+    g   = 9.8076
+    Hv  = 2501000
+    r   = mixing_ratio(dewpoint, pressure)
+    Rsd = 287
+    Rsw = 461.5
+    Cpd = 1003.5
+    T   = kelvin(temp)
+
+    return g * (1 + (Hv * r) / (Rsd * T)) \
+             / (Cpd + ((Hv**2 * r) / (Rsw * T**2)))