diff --git a/lib/xmet/sounding.py b/lib/xmet/sounding.py index 6754a87..0884404 100644 --- a/lib/xmet/sounding.py +++ b/lib/xmet/sounding.py @@ -9,7 +9,10 @@ from xmet.list import nearest from xmet.series import Series, SeriesIntersection from xmet.thermo import follow_dry_adiabat, \ follow_moist_adiabat, \ - follow_saturated_mixing_ratio + follow_saturated_mixing_ratio, \ + pressure_height, \ + kelvin, \ + celsius LAPSE_RATE_DRY = 9.8 # degrees C per 1000m LAPSE_RATE_MOIST = 4.0 @@ -184,7 +187,41 @@ class Sounding(DatabaseTable): return series +def between(n, a, b): + return n > a and n < b + class SoundingParameters(): + @staticmethod + def cape(temp_line: Series, + moist_adiabat: Series, + lfc: tuple[float], + el: tuple[float]) -> float: + cape = 0.0 + + p_lfc, p_el = lfc[1], el[1] + + neighbors = temp_line.neighbors(moist_adiabat) + gph_last = None + + for pair in neighbors: + p_env, p_parcel = pair + + if not between(p_env, p_el, p_lfc): + continue + + t_env = kelvin(temp_line[p_env]) + t_parcel = kelvin(moist_adiabat[p_parcel]) + gph = pressure_height(p_env) + + if gph_last is not None: + gph_delta = gph - gph_last + + cape += (((t_parcel) - (t_env)) / (t_env)) * gph_delta + + gph_last = gph + + return cape + @staticmethod def from_sounding(sounding: Sounding) -> Self: temp_line = sounding.follow_temp() @@ -202,8 +239,12 @@ class SoundingParameters(): el = moist_adiabat.intersect(temp_line, SeriesIntersection.LESSER, lfc[1]) params = SoundingParameters() - params.lcl = lcl - params.lfc = lfc - params.el = el + params.lcl = lcl + params.lfc = lfc + params.el = el + params.cape = SoundingParameters.cape(temp_line, + moist_adiabat, + lfc, + el) return params