use smallvec::SmallVec;
use sounding_base::Profile::*;
use sounding_base::{DataRow, Profile, Sounding};
use error::AnalysisError::*;
use error::*;
const FREEZING: f64 = 0.0;
#[derive(Debug, Clone, Copy)]
pub struct Layer {
pub bottom: DataRow,
pub top: DataRow,
}
pub type Layers = SmallVec<[Layer; ::VEC_SIZE]>;
impl Layer {
pub fn lapse_rate(&self) -> Result<f64> {
let top_t = self.top.temperature.ok_or(MissingValue)?;
let bottom_t = self.bottom.temperature.ok_or(MissingValue)?;
let dt = top_t - bottom_t;
let dz = self.height_thickness()?;
Ok(dt / dz * 1000.0)
}
#[cfg_attr(feature = "cargo-clippy", allow(float_cmp))]
pub fn height_thickness(&self) -> Result<f64> {
let top = self.top.height.ok_or(MissingValue)?;
let bottom = self.bottom.height.ok_or(MissingValue)?;
if top == bottom {
Err(InvalidInput)
} else {
Ok(top - bottom)
}
}
#[cfg_attr(feature = "cargo-clippy", allow(float_cmp))]
pub fn pressure_thickness(&self) -> Result<f64> {
let bottom_p = self.bottom.pressure.ok_or(MissingValue)?;
let top_p = self.top.pressure.ok_or(MissingValue)?;
if bottom_p == top_p {
Err(InvalidInput)
} else {
Ok(bottom_p - top_p)
}
}
pub fn wind_shear(&self) -> Result<(f64, f64)> {
let top_spd = self.top.speed.ok_or(MissingValue)?;
let top_dir = -(90.0 + self.top.direction.ok_or(MissingValue)?);
let bottom_spd = self.bottom.speed.ok_or(MissingValue)?;
let bottom_dir = -(90.0 + self.bottom.direction.ok_or(MissingValue)?);
let top_u = top_dir.to_radians().cos() * top_spd;
let top_v = top_dir.to_radians().sin() * top_spd;
let bottom_u = bottom_dir.to_radians().cos() * bottom_spd;
let bottom_v = bottom_dir.to_radians().sin() * bottom_spd;
let du = top_u - bottom_u;
let dv = top_v - bottom_v;
let shear_spd = du.hypot(dv);
let mut shear_dir = -(90.0 + dv.atan2(du).to_degrees());
while shear_dir < 0.0 {
shear_dir += 360.0;
}
while shear_dir > 360.0 {
shear_dir -= 360.0;
}
Ok((shear_spd, shear_dir))
}
}
#[cfg(test)]
mod layer_tests {
use super::*;
use optional::some;
use sounding_base::DataRow;
fn make_test_layer() -> Layer {
let mut bottom = DataRow::default();
bottom.pressure = some(1000.0);
bottom.temperature = some(20.0);
bottom.height = some(5.0);
bottom.speed = some(1.0);
bottom.direction = some(180.0);
let mut top = DataRow::default();
top.pressure = some(700.0);
top.temperature = some(-2.0);
top.height = some(3012.0);
top.speed = some(1.0);
top.direction = some(90.0);
Layer { bottom, top }
}
fn approx_eq(a: f64, b: f64, tol: f64) -> bool {
(a - b).abs() <= tol
}
#[test]
fn test_height_thickness() {
let lyr = make_test_layer();
println!("{:#?}", lyr);
assert!(approx_eq(
lyr.height_thickness().unwrap(),
3007.0,
::std::f64::EPSILON
));
}
#[test]
fn test_pressure_thickness() {
let lyr = make_test_layer();
println!("{:#?}", lyr);
assert!(approx_eq(
lyr.pressure_thickness().unwrap(),
300.0,
::std::f64::EPSILON
));
}
#[test]
fn test_lapse_rate() {
let lyr = make_test_layer();
println!(
"{:#?}\n\n -- \n\n {} \n\n --",
lyr,
lyr.lapse_rate().unwrap()
);
assert!(approx_eq(lyr.lapse_rate().unwrap(), -7.31626, 1.0e-5));
}
#[test]
fn test_wind_shear() {
let lyr = make_test_layer();
println!(
"{:#?}\n\n -- \n\n {:#?} \n\n --",
lyr,
lyr.wind_shear().unwrap()
);
let (speed_shear, direction_shear) = lyr.wind_shear().unwrap();
assert!(approx_eq(speed_shear, ::std::f64::consts::SQRT_2, 1.0e-5));
assert!(approx_eq(direction_shear, 45.0, 1.0e-5));
}
}
pub fn dendritic_snow_zone(snd: &Sounding) -> Result<Layers> {
temperature_layer(snd, -12.0, -18.0, 300.0)
}
pub fn hail_growth_zone(snd: &Sounding) -> Result<Layers> {
temperature_layer(snd, -10.0, -30.0, 1.0)
}
fn temperature_layer(
snd: &Sounding,
warm_side: f64,
cold_side: f64,
top_pressure: f64,
) -> Result<Layers> {
use interpolation::{linear_interp, linear_interpolate_sounding};
let mut to_return: Layers = Layers::new();
let t_profile = snd.get_profile(Temperature);
let p_profile = snd.get_profile(Pressure);
if t_profile.is_empty() || p_profile.is_empty() {
return Err(AnalysisError::MissingProfile);
}
izip!(p_profile, t_profile)
.filter_map(|pair| {
if pair.0.is_some() && pair.1.is_some() {
let (p, t) = (pair.0.unpack(), pair.1.unpack());
Some((p, t))
} else {
None
}
})
.take_while(|&(p, _)| p > top_pressure)
.fold(
Ok((None, None, None)),
|acc: Result<(Option<f64>, Option<f64>, Option<_>)>, (p, t)| {
match acc {
Ok((Some(last_p), Some(last_t), None)) => {
if last_t < cold_side && t >= cold_side && t <= warm_side {
let target_p = linear_interp(cold_side, last_t, t, last_p, p);
let bottom = linear_interpolate_sounding(snd, target_p)?;
Ok((Some(p), Some(t), Some(bottom)))
} else if last_t > warm_side && t >= cold_side && t <= warm_side {
let target_p = linear_interp(warm_side, last_t, t, last_p, p);
let bottom = linear_interpolate_sounding(snd, target_p)?;
Ok((Some(p), Some(t), Some(bottom)))
} else if (last_t < cold_side && t >= warm_side)
|| (last_t > warm_side && t <= cold_side)
{
let warm_p = linear_interp(warm_side, last_t, t, last_p, p);
let cold_p = linear_interp(cold_side, last_t, t, last_p, p);
let bottom = linear_interpolate_sounding(snd, warm_p.max(cold_p))?;
let top = linear_interpolate_sounding(snd, warm_p.min(cold_p))?;
to_return.push(Layer { bottom, top });
Ok((Some(p), Some(t), None))
} else {
Ok((Some(p), Some(t), None))
}
}
Ok((Some(last_p), Some(last_t), Some(bottom))) => {
if t < cold_side {
let target_p = linear_interp(cold_side, last_t, t, last_p, p);
let top = linear_interpolate_sounding(snd, target_p)?;
to_return.push(Layer { bottom, top });
Ok((Some(p), Some(t), None))
} else if t > warm_side {
let target_p = linear_interp(warm_side, last_t, t, last_p, p);
let top = linear_interpolate_sounding(snd, target_p)?;
to_return.push(Layer { bottom, top });
Ok((Some(p), Some(t), None))
} else {
Ok((Some(p), Some(t), Some(bottom)))
}
}
e @ Err(_) => e,
Ok((None, None, None)) => {
if t <= warm_side && t >= cold_side {
let dr = linear_interpolate_sounding(snd, p)?;
Ok((Some(p), Some(t), Some(dr)))
} else {
Ok((Some(p), Some(t), None))
}
}
_ => unreachable!(),
}
},
)
.and_then(|_| Ok(to_return))
}
pub fn warm_temperature_layer_aloft(snd: &Sounding) -> Result<Layers> {
warm_layer_aloft(snd, Temperature)
}
pub fn warm_wet_bulb_layer_aloft(snd: &Sounding) -> Result<Layers> {
warm_layer_aloft(snd, WetBulb)
}
fn warm_layer_aloft(snd: &Sounding, var: Profile) -> Result<Layers> {
use std::f64::MAX;
assert!(var == Temperature || var == WetBulb);
let mut to_return: Layers = Layers::new();
let t_profile = snd.get_profile(var);
let p_profile = snd.get_profile(Pressure);
if t_profile.is_empty() || p_profile.is_empty() {
return Err(AnalysisError::MissingProfile);
}
izip!(p_profile, t_profile)
.filter_map(|pair| {
if pair.0.is_some() && pair.1.is_some() {
let (p, t) = (pair.0.unpack(), pair.1.unpack());
Some((p, t))
} else {
None
}
})
.take_while(|&(p, _)| p > 500.0)
.fold(
Ok((MAX, MAX, None)),
|last_iter_res: Result<(_, _, _)>, (p, t)| {
let (last_p, last_t, mut bottom) = last_iter_res?;
if last_t <= FREEZING && t > FREEZING && bottom.is_none() {
let bottom_p = ::interpolation::linear_interp(FREEZING, last_t, t, last_p, p);
bottom = Some(::interpolation::linear_interpolate_sounding(snd, bottom_p)?);
}
if bottom.is_some() && last_t > FREEZING && t <= FREEZING {
let top_p = ::interpolation::linear_interp(FREEZING, last_t, t, last_p, p);
let top = ::interpolation::linear_interpolate_sounding(snd, top_p)?;
{
let bottom = bottom.unwrap();
to_return.push(Layer { bottom, top });
}
bottom = None;
}
Ok((p, t, bottom))
},
)?;
Ok(to_return)
}
pub fn cold_surface_temperature_layer(snd: &Sounding, warm_layers: &[Layer]) -> Result<Layer> {
cold_surface_layer(snd, Temperature, warm_layers)
}
fn cold_surface_layer(snd: &Sounding, var: Profile, warm_layers: &[Layer]) -> Result<Layer> {
debug_assert!(var == Temperature || var == WetBulb);
if warm_layers.is_empty() {
return Err(InvalidInput);
}
let t_profile = snd.get_profile(var);
let p_profile = snd.get_profile(Pressure);
if t_profile.is_empty() || p_profile.is_empty() {
return Err(MissingProfile);
}
izip!(0usize.., p_profile, t_profile)
.filter_map(|triplet| {
if triplet.1.is_some() && triplet.2.is_some() {
let (i, t) = (triplet.0, triplet.2.unpack());
Some((i, t))
} else {
None
}
})
.map(|(i, t)| {
if t > FREEZING {
Err(InvalidInput)
} else {
Ok(i)
}
})
.next()
.unwrap_or(Err(NoDataProfile))
.and_then(|index| snd.get_data_row(index).ok_or(InvalidInput))
.map(|bottom| Layer {
bottom,
top: warm_layers[0].bottom,
})
}
pub fn layer_agl(snd: &Sounding, meters_agl: f64) -> Result<Layer> {
use std::f64::MAX;
let tgt_elev = {
let elev = snd.get_station_info().elevation();
if elev.is_some() {
elev.unpack() + meters_agl
} else {
return Err(MissingValue);
}
};
let h_profile = snd.get_profile(GeopotentialHeight);
let p_profile = snd.get_profile(Pressure);
if h_profile.is_empty() || p_profile.is_empty() {
return Err(MissingProfile);
}
let bottom = snd.surface_as_data_row();
izip!(p_profile, h_profile)
.filter_map(|pair| {
if pair.0.is_some() && pair.1.is_some() {
let (p, h) = (pair.0.unpack(), pair.1.unpack());
Some((p, h))
} else {
None
}
})
.fold(
Ok((MAX, 0.0f64, None)),
|acc: Result<(_, _, Option<_>)>, (p, h)| {
match acc {
Ok((last_p, last_h, None)) => {
if h > tgt_elev {
let tgt_p =
::interpolation::linear_interp(tgt_elev, last_h, h, last_p, p);
Ok((MAX, MAX, Some(tgt_p)))
} else {
Ok((p, h, None))
}
}
ok @ Ok((_, _, Some(_))) => ok,
e @ Err(_) => e,
}
},
)
.and_then(|(_, _, opt)| opt.ok_or(NotEnoughData))
.and_then(|target_p| ::interpolation::linear_interpolate_sounding(snd, target_p))
.map(|top| Layer { bottom, top })
}
pub fn pressure_layer(snd: &Sounding, bottom_p: f64, top_p: f64) -> Result<Layer> {
let sfc_pressure = snd.surface_as_data_row().pressure;
if sfc_pressure.is_some() && sfc_pressure.unwrap() < bottom_p {
return Err(InvalidInput);
}
let bottom = ::interpolation::linear_interpolate_sounding(snd, bottom_p)?;
let top = ::interpolation::linear_interpolate_sounding(snd, top_p)?;
Ok(Layer { bottom, top })
}
pub fn inversions(snd: &Sounding, top_p: f64) -> Result<Layers> {
use std::f64::MAX;
let mut to_return: Layers = Layers::new();
let t_profile = snd.get_profile(Temperature);
let p_profile = snd.get_profile(Pressure);
if t_profile.is_empty() || p_profile.is_empty() {
return Err(AnalysisError::MissingProfile);
}
izip!(0usize.., p_profile, t_profile)
.filter_map(|triple| {
if triple.1.is_some() && triple.2.is_some() {
let (i, p, t) = (triple.0, triple.1.unpack(), triple.2.unpack());
Some((i, p, t))
} else {
None
}
})
.filter_map(|(i, p, t)| if p < top_p { None } else { Some((i, t)) })
.fold(
(0, MAX, None),
|(last_i, last_t, mut bottom_opt), (i, t)| {
if bottom_opt.is_none() && last_t < t {
bottom_opt = snd.get_data_row(last_i);
} else if bottom_opt.is_some() && last_t > t {
if let Some(layer) = bottom_opt.and_then(|bottom| {
snd.get_data_row(last_i)
.and_then(|top| Some(Layer { bottom, top }))
}) {
to_return.push(layer);
bottom_opt = None;
}
}
(i, t, bottom_opt)
},
);
Ok(to_return)
}
pub fn sfc_based_inversion(snd: &Sounding) -> Result<Option<Layer>> {
let t_profile = snd.get_profile(Temperature);
let p_profile = snd.get_profile(Pressure);
let h_profile = snd.get_profile(GeopotentialHeight);
if t_profile.is_empty() || p_profile.is_empty() || h_profile.is_empty() {
return Err(AnalysisError::MissingProfile);
}
izip!(0usize.., p_profile, h_profile, t_profile)
.filter_map(|tuple| {
if tuple.1.is_some() && tuple.2.is_some() && tuple.2.is_some() {
Some(tuple.0)
} else {
None
}
})
.nth(0)
.ok_or(AnalysisError::NotEnoughData)
.and_then(|index| snd.get_data_row(index).ok_or(AnalysisError::MissingValue))
.and_then(|bottom_row| {
let sfc_t = bottom_row.temperature;
if sfc_t.is_some() {
let sfc_t = sfc_t.unpack();
let val = izip!(0usize.., p_profile, t_profile, h_profile)
.filter_map(|tuple| {
if tuple.1.is_some() && tuple.2.is_some() && tuple.3.is_some() {
let (i, p, t) = (tuple.0, tuple.1.unpack(), tuple.2.unpack());
Some((i, p, t))
} else {
None
}
})
.skip(1)
.take_while(|(_, p, _)| *p > 690.0)
.filter(|(_, _, t)| *t > sfc_t)
.fold(None, |max_t_info, (i, _, t)| {
if let Some((max_t, _max_t_idx)) = max_t_info {
if t > max_t {
Some((t, i))
} else {
max_t_info
}
} else {
Some((t, i))
}
});
match val {
Some((_, idx)) => snd
.get_data_row(idx)
.ok_or(AnalysisError::MissingValue)
.and_then(|top_row| {
Ok(Some(Layer {
bottom: bottom_row,
top: top_row,
}))
}),
None => Ok(None),
}
} else {
Err(AnalysisError::MissingValue)
}
})
}