Reimplement draw_adiabat() with loft_parcel()

Implement loft_parcel() to iterate through pressure levels, yielding the
temperature and pressure of the parcel at each height

lib/xmet/skew_t.py

@@ -4,12 +4,9 @@ import cairo
 from typing import Callable
 
 from xmet.sounding import Sounding
-from xmet.thermo   import pressure_height, lapse, moist_lapse_rate, \
+from xmet.thermo   import pressure_height, loft_parcel, moist_lapse_rate, \
-                          LAPSE_RATE_DRY
+                          LAPSE_RATE_DRY, PRESSURE_MAX, PRESSURE_MIN, \
-
+                          PRESSURE_STEP
-PRESSURE_MAX  = 1050 # millibar
-PRESSURE_MIN  =  100
-PRESSURE_STEP =   50
 
 PRESSURE_LOG_MAX   = math.log(PRESSURE_MAX)
 PRESSURE_LOG_MIN   = math.log(PRESSURE_MIN)
@@ -107,26 +104,11 @@ class SkewTGraph():
                      start_temp: float,
                      start_pressure: float,
                      lapse_rate: Callable):
-        start_height = pressure_height(start_pressure)
+        sx_last = None
-        sx_last      = None
+        sy_last = None
-        sy_last      = None
-        height_last  = None
 
-        temp     = start_temp
+        for level in loft_parcel(start_temp, start_pressure, lapse_rate):
-        pressure = start_pressure
+            temp, pressure = level
-
-        while pressure >= PRESSURE_MIN:
-            height = pressure_height(pressure)
-
-            if height_last is None:
-                height_last = height
-
-            try:
-                rate = lapse_rate(temp, temp, pressure)
-            except OverflowError:
-                break
-
-            temp = lapse(temp, height - height_last, rate)
 
             sx, sy = self.sample_to_screen(temp, pressure)
 
@@ -135,16 +117,8 @@ class SkewTGraph():
                 cr.line_to(x + sx_last, y + sy_last)
                 cr.stroke()
 
-            sx_last     = sx
+            sx_last = sx
-            sy_last     = sy
+            sy_last = sy
-            height_last = height
-
-            if pressure == PRESSURE_MIN:
-                break
-            elif pressure - 10.0 < PRESSURE_MIN:
-                pressure = PRESSURE_MIN
-            else:
-                pressure -= 10.0
 
     def draw_adiabats(self,
                       cr: cairo.Context,

lib/xmet/thermo.py

@@ -1,6 +1,12 @@
+from typing import Callable
+
 LAPSE_RATE_DRY   = 9.8 / 1000 # degrees C per km
 LAPSE_RATE_MOIST = 4.0 / 1000
 
+PRESSURE_MAX  = 1050 # millibar
+PRESSURE_MIN  =  100
+PRESSURE_STEP =   50
+
 def kelvin(c: float) -> float:
     return 273.15 + c
 
@@ -67,3 +73,45 @@ def moist_lapse_rate(temp: float, dewpoint: float, pressure: float) -> float:
 
     return g * (1 + (Hv * r) / (Rsd * T)) \
              / (Cpd + (((Hv**2) * r) / (Rsw * (T**2))))
+
+def loft_parcel(start_temp: float,
+                start_pressure: float,
+                lapse_rate: Callable):
+    """
+    Loft a parcel of air from a given pressure, at a given temperature,
+    yielding a Tuple containing the temperature and pressure of that parcel
+    at each increment of elevation.  A Callable which, given a temperature,
+    dewpoint, pressure value, will be dispatched to obtain the lapse rate
+    for each height.
+    """
+    start_height = pressure_height(start_pressure)
+    sx_last      = None
+    sy_last      = None
+    height_last  = None
+
+    temp     = start_temp
+    pressure = start_pressure
+
+    while pressure >= PRESSURE_MIN:
+        height = pressure_height(pressure)
+
+        if height_last is None:
+            height_last = height
+
+        try:
+            rate = lapse_rate(temp, temp, pressure)
+        except OverflowError:
+            break
+
+        temp = lapse(temp, height - height_last, rate)
+
+        yield temp, pressure
+
+        height_last = height
+
+        if pressure == PRESSURE_MIN:
+            break
+        elif pressure - 10.0 < PRESSURE_MIN:
+            pressure = PRESSURE_MIN
+        else:
+            pressure -= 10.0