use metfor::{IntHelicityM2pS2, Knots, Meters, MetersPSec, Quantity, WindSpdDir, WindUV};
use sounding_base::Sounding;
use crate::error::*;
use crate::layers::{self, Layer};
use crate::levels::height_level;
use itertools::Itertools;
pub fn mean_wind(layer: &Layer, snd: &Sounding) -> Result<WindUV<MetersPSec>> {
use std::f64;
let height = snd.height_profile();
let wind = snd.wind_profile();
let max_hgt = if layer.top.height.is_some() {
layer.top.height.unpack()
} else {
return Err(AnalysisError::MissingProfile);
};
let min_hgt = if layer.bottom.height.is_some() {
layer.bottom.height.unpack()
} else {
return Err(AnalysisError::MissingProfile);
};
let (old_hgt, old_u, old_v, mut iu, mut iv, dz) = izip!(height, wind)
.filter_map(|(hgt, wind)| {
if hgt.is_some() && wind.is_some() {
Some((hgt.unpack(), wind.unpack()))
} else {
None
}
})
.skip_while(|&(hgt, _)| hgt < min_hgt)
.take_while(|&(hgt, _)| hgt <= max_hgt)
.map(|(hgt, wind)| {
let WindUV { u, v } = WindUV::<MetersPSec>::from(wind);
(hgt, u, v)
})
.fold(
(
Meters(f64::MAX),
MetersPSec(0.0),
MetersPSec(0.0),
MetersPSec(0.0),
MetersPSec(0.0),
Meters(0.0),
),
|acc, (hgt, u, v)| {
let (old_hgt, old_u, old_v, mut iu, mut iv, mut acc_dz) = acc;
let dz = hgt - old_hgt;
if dz < Meters(0.0) {
return (hgt, u, v, iu, iv, acc_dz);
}
iu += (u + old_u) * dz.unpack();
iv += (v + old_v) * dz.unpack();
acc_dz += dz;
(hgt, u, v, iu, iv, acc_dz)
},
);
if old_hgt == Meters(f64::MAX) {
return Err(AnalysisError::NotEnoughData);
} else if dz.unpack().abs() < std::f64::EPSILON {
iu = old_u;
iv = old_v;
} else {
iu /= 2.0 * dz.unpack();
iv /= 2.0 * dz.unpack();
}
Ok(WindUV { u: iu, v: iv })
}
pub fn sr_helicity<W>(
layer: &Layer,
storm_motion_uv_ms: W,
snd: &Sounding,
) -> Result<IntHelicityM2pS2>
where
WindUV<MetersPSec>: From<W>,
{
let height = snd.height_profile();
let wind = snd.wind_profile();
let storm_motion_uv_ms = WindUV::<MetersPSec>::from(storm_motion_uv_ms);
let bottom = layer.bottom.height.ok_or(AnalysisError::MissingValue)?;
let top = layer.top.height.ok_or(AnalysisError::MissingValue)?;
izip!(height, wind)
.filter_map(|(h, w)| {
if let (Some(h), Some(w)) = (h.into_option(), w.into_option()) {
Some((h, w))
} else {
None
}
})
.map(|(h, w)| {
let WindUV { u, v }: WindUV<MetersPSec> = From::<WindSpdDir<Knots>>::from(w);
(h, (u - storm_motion_uv_ms.u), (v - storm_motion_uv_ms.v))
})
.tuple_windows::<(_, _, _)>()
.skip_while(|(_, (h, _, _), _)| *h < bottom)
.take_while(|(_, (h, _, _), _)| *h <= top)
.map(|((h0, u0, v0), (h1, u1, v1), (h2, u2, v2))| {
let dz = Meters::from(h2 - h0).unpack();
let du = MetersPSec::from(u2 - u0).unpack() / dz;
let dv = MetersPSec::from(v2 - v0).unpack() / dz;
(h1, u1, v1, du, dv)
})
.tuple_windows::<(_, _)>()
.fold(Err(AnalysisError::NotEnoughData), |acc, (lvl0, lvl1)| {
let mut integrated_helicity: f64 = acc.unwrap_or(0.0);
let (z0, u0, v0, du0, dv0) = lvl0;
let (z1, u1, v1, du1, dv1) = lvl1;
let h0 = u0 * dv0 - v0 * du0;
let h1 = u1 * dv1 - v1 * du1;
let h = MetersPSec::from(h0 + h1).unpack() * Meters::from(z1 - z0).unpack();
integrated_helicity += h;
Ok(integrated_helicity)
})
.map(|integrated_helicity| IntHelicityM2pS2(-integrated_helicity / 2.0))
}
pub fn bunkers_storm_motion(snd: &Sounding) -> Result<(WindUV<MetersPSec>, WindUV<MetersPSec>)> {
let layer = &layers::layer_agl(snd, Meters(6000.0))?;
let WindUV {
u: mean_u,
v: mean_v,
} = mean_wind(layer, snd)?;
let WindUV {
u: shear_u,
v: shear_v,
} = bulk_shear_half_km(layer, snd)?;
const D: f64 = 7.5;
let scale = D / shear_u.unpack().hypot(shear_v.unpack());
let (delta_u, delta_v) = (shear_v * scale, -shear_u * scale);
Ok((
WindUV {
u: mean_u + delta_u,
v: mean_v + delta_v,
},
WindUV {
u: mean_u - delta_u,
v: mean_v - delta_v,
},
))
}
pub(crate) fn bulk_shear_half_km(layer: &Layer, snd: &Sounding) -> Result<WindUV<MetersPSec>> {
let bottom = layer
.bottom
.height
.into_option()
.ok_or(AnalysisError::MissingValue)?;
let top = layer
.top
.height
.into_option()
.ok_or(AnalysisError::MissingValue)?;
if top - bottom < Meters(750.0) {
return Err(AnalysisError::NotEnoughData);
}
let top_bottom_layer = height_level(bottom + Meters(500.0), snd)?;
let bottom_layer = &Layer {
top: top_bottom_layer,
bottom: layer.bottom,
};
let WindUV {
u: bottom_u,
v: bottom_v,
} = mean_wind(bottom_layer, snd)?;
let bottom_top_layer = height_level(top - Meters(500.0), snd)?;
let top_layer = &Layer {
top: layer.top,
bottom: bottom_top_layer,
};
let WindUV { u: top_u, v: top_v } = mean_wind(top_layer, snd)?;
Ok(WindUV {
u: top_u - bottom_u,
v: top_v - bottom_v,
})
}