Initial implementation of derived LCL and LFC params

lib/xmet/skew_t.py

@@ -1,12 +1,16 @@
 import math
 import cairo
 
-from typing import Callable
+from typing    import Callable
+from itertools import product
 
-from xmet.sounding import Sounding
+from xmet.list     import nearest
+from xmet.util     import cmp
+from xmet.sounding import Sounding, SoundingSample
 from xmet.thermo   import pressure_height, loft_parcel, moist_lapse_rate, \
                           LAPSE_RATE_DRY, PRESSURE_MAX, PRESSURE_MIN, \
-                          PRESSURE_STEP
+                          PRESSURE_STEP, mixing_ratio, saturated_mixing_ratio, \
+                          mixing_ratio_temp
 
 PRESSURE_LOG_MAX   = math.log(PRESSURE_MAX)
 PRESSURE_LOG_MIN   = math.log(PRESSURE_MIN)
@@ -176,6 +180,92 @@ class SkewTGraph():
         cr.stroke()
         cr.restore()
 
+    def follow_dry_adiabat(self, temp: float, pressure: float) -> dict:
+        levels = dict()
+
+        for level in loft_parcel(temp, pressure, lambda t, p: LAPSE_RATE_DRY):
+            t2, p2 = level
+
+            levels[p2] = t2
+
+        return levels
+
+    def follow_moist_adiabat(self, temp: float, pressure: float) -> dict:
+        levels = dict()
+
+        for level in loft_parcel(temp, pressure, moist_lapse_rate):
+            t2, p2 = level
+
+            levels[p2] = t2
+
+        return levels
+
+    def follow_saturated_mixing_ratio(self, temp: float, pressure: float) -> dict:
+        levels = dict()
+
+        ratio = saturated_mixing_ratio(temp, pressure)
+
+        p2 = pressure
+
+        while p2 >= PRESSURE_MIN:
+            levels[p2] = mixing_ratio_temp(ratio, p2)
+
+            p2 -= 10
+
+        return levels
+
+    def follow_sounding(self, sounding: Sounding) -> dict:
+        levels = dict()
+
+        for sample in sounding.samples:
+            levels[sample.pressure] = sample.temp
+
+        return levels
+
+    def intersect(self, series1: dict, series2: dict):
+        pairs = nearest(sorted(series1.keys(), reverse=True),
+                        sorted(series2.keys(), reverse=True))
+
+        sign_last = None
+
+        for pair in pairs:
+            v1, v2 = series1[pair[0]], series2[pair[1]]
+
+            sign = cmp(v1, v2)
+
+            if sign == 0 or (sign_last is not None and sign_last != sign):
+                return pair[0]
+
+            sign_last = sign
+
+    def find_lfc(self, temp: float, pressure: float, sounding: Sounding) -> float:
+        moist_adiabat = self.follow_moist_adiabat(temp, pressure)
+        temp_line     = self.follow_sounding(sounding)
+
+        pairs = nearest(sorted(moist_adiabat.keys(), reverse=True),
+                        sorted(temp_line.keys(),     reverse=True))
+
+        for pair in pairs:
+            v1, v2 = moist_adiabat[pair[0]], temp_line[pair[1]]
+
+            if v1 > v2:
+                return pair[0]
+
+    def derive_parameters(self, sounding: Sounding):
+        dry_adiabat = self.follow_dry_adiabat(sounding.samples[0].temp,
+                                              sounding.samples[0].pressure)
+
+        saturated_mr_line = self.follow_saturated_mixing_ratio(sounding.samples[0].dewpoint,
+                                                               sounding.samples[0].pressure)
+
+        lcl = self.intersect(dry_adiabat, saturated_mr_line)
+        lfc = self.find_lfc(dry_adiabat[lcl], lcl, sounding)
+
+        return {
+            'lcl': lcl,
+            'lfc': lfc
+        }
+
     def draw(self,
              cr: cairo.Context,
              x: float,