use super::{Layer, Layers};
use crate::{
error::{AnalysisError, Result},
interpolation::{linear_interp, linear_interpolate_sounding},
sounding::Sounding,
};
use itertools::{izip, Itertools};
use metfor::{
Celsius, CelsiusDiff, HectoPascal, JpKg, Kelvin, KelvinDiff, Meters, Quantity, FREEZING,
};
use optional::Optioned;
use std::iter::once;
pub fn dendritic_snow_zone(snd: &Sounding) -> Result<Layers> {
temperature_layer(snd, Celsius(-12.0), Celsius(-18.0), HectoPascal(300.0))
}
pub fn hail_growth_zone(snd: &Sounding) -> Result<Layers> {
temperature_layer(snd, Celsius(-10.0), Celsius(-30.0), HectoPascal(1.0))
}
#[inline]
fn temperature_layer(
snd: &Sounding,
warm_side: Celsius,
cold_side: Celsius,
top_pressure: HectoPascal,
) -> Result<Layers> {
let t_profile = snd.temperature_profile();
let p_profile = snd.pressure_profile();
if t_profile.is_empty() || p_profile.is_empty() {
return Err(AnalysisError::MissingProfile);
}
enum Crossing {
Enter(HectoPascal),
Exit(HectoPascal),
Over(HectoPascal, HectoPascal),
In(HectoPascal),
}
fn to_crossing_type(
pnt0: (HectoPascal, Celsius),
pnt1: (HectoPascal, Celsius),
warm_side: Celsius,
cold_side: Celsius,
) -> Option<Crossing> {
let (p0, t0) = pnt0;
let (p1, t1) = pnt1;
debug_assert!(p0 >= p1, "Pressure must decrease with height.");
if t0 < cold_side && t1 >= cold_side && t1 <= warm_side {
let cold_p = linear_interp(cold_side, t0, t1, p0, p1);
Some(Crossing::Enter(cold_p))
} else if t0 > warm_side && t1 <= warm_side && t1 >= cold_side {
let warm_p = linear_interp(warm_side, t0, t1, p0, p1);
Some(Crossing::Enter(warm_p))
} else if (t0 < cold_side && t1 > warm_side) || (t0 > warm_side && t1 < cold_side) {
let warm_p = linear_interp(warm_side, t0, t1, p0, p1);
let cold_p = linear_interp(cold_side, t0, t1, p0, p1);
Some(Crossing::Over(warm_p.max(cold_p), warm_p.min(cold_p)))
} else if t0 >= cold_side && t0 <= warm_side && t1 > warm_side {
let warm_p = linear_interp(warm_side, t0, t1, p0, p1);
Some(Crossing::Exit(warm_p))
} else if t0 >= cold_side && t0 <= warm_side && t1 < cold_side {
let cold_p = linear_interp(cold_side, t0, t1, p0, p1);
Some(Crossing::Exit(cold_p))
} else if t0 >= cold_side && t0 <= warm_side && t1 >= cold_side && t1 <= warm_side {
Some(Crossing::In(p0))
} else {
None
}
}
izip!(p_profile, t_profile)
.filter(|(p, t)| p.is_some() && t.is_some())
.map(|(p, t)| (p.unpack(), t.unpack()))
.take_while(move |&(p, _)| p > top_pressure)
.tuple_windows::<(_, _)>()
.filter_map(|(pnt0, pnt1)| to_crossing_type(pnt0, pnt1, warm_side, cold_side))
.scan(None, |bottom_p: &mut Option<_>, crossing_type: Crossing| {
match crossing_type {
Crossing::In(p) => {
if bottom_p.is_none() {
*bottom_p = Some(p);
}
Some(None)
}
Crossing::Enter(p) => {
*bottom_p = Some(p);
Some(None)
}
Crossing::Exit(top_p) => {
if let Some(bp) = bottom_p.take() {
Some(Some((bp, top_p)))
} else {
Some(None)
}
}
Crossing::Over(p0, p1) => {
debug_assert!(bottom_p.is_none());
Some(Some((p0, p1)))
}
}
})
.filter_map(|opt| opt)
.map(|(bottom_p, top_p)| {
linear_interpolate_sounding(snd, bottom_p).and_then(|bottom| {
linear_interpolate_sounding(snd, top_p).map(|top| Layer { bottom, top })
})
})
.collect()
}
pub fn warm_temperature_layer_aloft(snd: &Sounding) -> Result<Layers> {
warm_layer_aloft(snd, snd.temperature_profile())
}
pub fn warm_wet_bulb_layer_aloft(snd: &Sounding) -> Result<Layers> {
warm_layer_aloft(snd, snd.wet_bulb_profile())
}
#[inline]
fn warm_layer_aloft(snd: &Sounding, t_profile: &[Optioned<Celsius>]) -> Result<Layers> {
let p_profile = snd.pressure_profile();
if t_profile.is_empty() || p_profile.is_empty() {
return Err(AnalysisError::MissingProfile);
}
enum Crossing {
IntoWarmLayer(HectoPascal),
OutOfWarmLayer(HectoPascal),
}
fn crossing_type(
pnt0: (HectoPascal, Celsius),
pnt1: (HectoPascal, Celsius),
) -> Option<Crossing> {
let (p0, t0) = pnt0;
let (p1, t1) = pnt1;
if t0 < FREEZING && t1 >= FREEZING {
let crossing_p = linear_interp(FREEZING, t0, t1, p0, p1);
Some(Crossing::IntoWarmLayer(crossing_p))
} else if t0 >= FREEZING && t1 < FREEZING {
let crossing_p = linear_interp(FREEZING, t0, t1, p0, p1);
Some(Crossing::OutOfWarmLayer(crossing_p))
} else {
None
}
}
izip!(p_profile, t_profile)
.filter(|(p, t)| p.is_some() && t.is_some())
.map(|(p, t)| (p.unpack(), t.unpack()))
.skip_while(|&(_, t)| t > FREEZING)
.take_while(|&(p, _)| p > HectoPascal(500.0))
.tuple_windows::<(_, _)>()
.filter_map(|(pnt0, pnt1)| crossing_type(pnt0, pnt1))
.scan(
None,
|bottom: &mut Option<_>, crossing: Crossing| match crossing {
Crossing::IntoWarmLayer(bottom_p) => {
debug_assert!(bottom.is_none());
*bottom = Some(bottom_p);
Some(None)
}
Crossing::OutOfWarmLayer(top_p) => {
Some(bottom.take().map(|bottom_p| (bottom_p, top_p)))
}
},
)
.filter_map(|opt| opt)
.map(|(bottom_p, top_p)| {
linear_interpolate_sounding(snd, bottom_p).and_then(|bottom| {
linear_interpolate_sounding(snd, top_p).map(|top| Layer { bottom, top })
})
})
.collect()
}
pub fn cold_surface_temperature_layer(
snd: &Sounding,
warm_layers: &[Layer],
) -> Result<Option<Layer>> {
Ok(cold_surface_layer(
snd,
snd.temperature_profile(),
warm_layers,
))
}
fn cold_surface_layer(
snd: &Sounding,
t_profile: &[Optioned<Celsius>],
warm_layers: &[Layer],
) -> Option<Layer> {
if warm_layers.is_empty() {
return None;
}
let p_profile = snd.pressure_profile();
let layer: Option<Layer> = izip!(0usize.., p_profile, t_profile)
.filter(|(_, p, t)| p.is_some() && t.is_some())
.map(|(i, _, t)| (i, t.unpack()))
.take_while(|(_, t)| *t <= FREEZING)
.next()
.and_then(|(i, _)| snd.data_row(i))
.and_then(|bottom| warm_layers.get(0).map(|lyr| (bottom, lyr.bottom)))
.map(|(bottom, top)| Layer { bottom, top });
layer
}
pub fn warm_surface_temperature_layer(snd: &Sounding) -> Result<Option<Layer>> {
warm_surface_layer(snd, snd.temperature_profile())
}
fn warm_surface_layer(snd: &Sounding, t_profile: &[Optioned<Celsius>]) -> Result<Option<Layer>> {
let p_profile = snd.pressure_profile();
if t_profile.is_empty() || p_profile.is_empty() {
return Err(AnalysisError::NotEnoughData);
}
let iter = izip!(0usize.., p_profile, t_profile)
.filter(|(_, p, t)| p.is_some() && t.is_some())
.map(|(i, p, t)| (i, p.unpack(), t.unpack()));
let bottom = iter
.clone()
.take_while(|(_, _, t)| *t >= FREEZING)
.next()
.and_then(|(i, _, _)| snd.data_row(i));
let top = iter
.tuple_windows::<(_, _)>()
.take_while(|((_, _, t), _)| *t >= FREEZING)
.last()
.map(|((_, p0, t0), (_, p1, t1))| linear_interp(FREEZING, t0, t1, p0, p1))
.map(|target_p| linear_interpolate_sounding(snd, target_p))
.transpose()?;
Ok(bottom.and_then(|bottom| top.map(|top| Layer { bottom, top })))
}
pub fn melting_freezing_energy_area(snd: &Sounding, lyr: &Layer) -> Result<JpKg> {
let top = &lyr.top;
let bottom = &lyr.bottom;
let theta_top: Kelvin = top
.pressure
.map(|p| metfor::potential_temperature(p, FREEZING))
.ok_or(AnalysisError::MissingValue)?;
let theta_bottom: Kelvin = bottom
.pressure
.map(|p| metfor::potential_temperature(p, FREEZING))
.ok_or(AnalysisError::MissingValue)?;
let bottom_h = bottom.height.ok_or(AnalysisError::MissingValue)?;
let top_h = top.height.ok_or(AnalysisError::MissingValue)?;
let bottom_iter = once((bottom.height, bottom.temperature))
.filter(|(h, t)| h.is_some() && t.is_some())
.map(|(h, t)| (h.unpack(), t.unpack()));
let top_iter = once((top.height, top.temperature))
.filter(|(h, t)| h.is_some() && t.is_some())
.map(|(h, t)| (h.unpack(), t.unpack()));
let middle_iter = izip!(snd.height_profile(), snd.temperature_profile())
.filter(|(h, t)| h.is_some() && t.is_some())
.map(|(h, t)| (h.unpack(), t.unpack()))
.skip_while(|(h, _)| *h <= bottom_h)
.take_while(|(h, _)| *h < top_h);
let iter = bottom_iter
.chain(middle_iter)
.chain(top_iter)
.map(|(h, t)| (h, Kelvin::from(t)))
.tuple_windows::<(_, _)>();
let (sum_dz, sum_kelvin_dz) = iter.fold((Meters(0.0), 0.0), |acc, ((h0, t0), (h1, t1))| {
let (mut sum_dz, mut sum_kelvin_dz) = acc;
let dz: Meters = h1 - h0;
let dt0: KelvinDiff = t0 - Kelvin::from(FREEZING);
let dt1: KelvinDiff = t1 - Kelvin::from(FREEZING);
let kelvin_dz = (dt0 + dt1).unpack() * dz.unpack();
sum_dz += dz;
sum_kelvin_dz += kelvin_dz;
(sum_dz, sum_kelvin_dz)
});
let mean_temp = CelsiusDiff(sum_kelvin_dz / sum_dz.unpack() / 2.0);
let energy = metfor::cp * mean_temp * f64::ln(theta_top / theta_bottom);
Ok(energy)
}