Implement rendering of dry adiabats

lib/xmet/skew_t.py

@@ -4,6 +4,8 @@ import cairo
 from typing import Callable
 
 from xmet.sounding import Sounding
+from xmet.thermo   import pressure_height, lapse, \
+                          LAPSE_RATE_MOIST, LAPSE_RATE_DRY
 
 PRESSURE_MAX  = 1050 # millibar
 PRESSURE_MIN  =  100
@@ -101,6 +103,47 @@ class SkewTGraph():
 
         cr.restore()
 
+    def draw_adiabat(self,
+                     cr: cairo.Context,
+                     x: float,
+                     y: float,
+                     start_temp: float,
+                     start_pressure: float,
+                     lapse_rate: float):
+        start_height = pressure_height(start_pressure)
+        sx_last = None
+        sy_last = None
+
+        pressure = start_pressure
+
+        while pressure >= PRESSURE_MIN:
+            height   = pressure_height(pressure)
+            temp_cur = lapse(start_temp, lapse_rate, height - start_height)
+
+            sx, sy = self.sample_to_screen(temp_cur, 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
+
+            pressure -= 10.0
+
+    def draw_adiabats(self,
+                      cr: cairo.Context,
+                      x: float,
+                      y: float):
+        cr.save()
+        cr.set_source_rgba(1.0, 0.2, 0, 0.2)
+
+        for temp in range(-140, 140, 10):
+            self.draw_adiabat(cr, x, y, temp, PRESSURE_MAX, LAPSE_RATE_DRY)
+
+        cr.restore()
+
     def draw_sounding(self,
                      cr: cairo.Context,
                      x: float,
@@ -149,6 +192,7 @@ class SkewTGraph():
 
         self.draw_isotherms(cr, x, y)
         self.draw_isobars(cr, x, y)
+        self.draw_adiabats(cr, x, y)
 
         cr.set_source_rgb(1, 0, 0)
         self.draw_sounding(cr, x, y, sounding, lambda s: s.temp)

lib/xmet/thermo.py

@@ -43,9 +43,16 @@ def lcl(temp: float, dewpoint: float) -> float:
     """
     return (temp - dewpoint) / 0.008
 
+def pressure_height(pressure: float) -> float:
+    """
+    Return the approximate altitude, in meters, for a given pressure in
+    millibar.
+    """
+    return (1 - (pressure / 1013.25) ** 0.190284) * 145366.45 * 0.3048
+
 def lapse(temp: float, rate: float, delta: float) -> float:
     """
     Return the temperature of a parcel cooled at either the dry or moist lapse
     rate for a given increase in height (in meters).
     """
-    return temp - (rate * (delta / 1000))
+    return temp - (rate * delta)