use optional::Optioned;
use metfor::{self, Celsius, HectoPascal, Kelvin, Quantity};
use sounding_base::{DataRow, Sounding};
use crate::error::*;
use crate::interpolation::{linear_interp, linear_interpolate_sounding};
use crate::profile::equivalent_potential_temperature;
#[derive(Debug, Clone, Copy)]
pub struct Parcel {
pub temperature: Celsius,
pub pressure: HectoPascal,
pub dew_point: Celsius,
}
impl Parcel {
pub fn theta(&self) -> Kelvin {
metfor::theta(self.pressure, self.temperature)
}
pub fn theta_e(&self) -> Result<Kelvin> {
metfor::theta_e(self.temperature, self.dew_point, self.pressure)
.ok_or(AnalysisError::MetForError)
}
pub fn mixing_ratio(&self) -> Result<f64> {
metfor::mixing_ratio(self.dew_point, self.pressure).ok_or(AnalysisError::MetForError)
}
pub fn virtual_temperature(&self) -> Result<Kelvin> {
metfor::virtual_temperature(self.temperature, self.dew_point, self.pressure)
.ok_or(AnalysisError::MetForError)
}
pub fn from_datarow(mut dr: DataRow) -> Option<Self> {
let temperature = dr.temperature.take()?;
let pressure = dr.pressure.take()?;
let dew_point = dr.dew_point.take()?;
Some(Parcel {
temperature,
pressure,
dew_point,
})
}
}
pub fn mixed_layer_parcel(snd: &Sounding) -> Result<Parcel> {
let press = snd.pressure_profile();
let t = snd.temperature_profile();
let dp = snd.dew_point_profile();
if press.is_empty() || t.is_empty() || dp.is_empty() {
return Err(AnalysisError::MissingProfile);
}
let bottom_p = press
.iter()
.filter_map(|&p| p.into_option())
.nth(0)
.ok_or(AnalysisError::NoDataProfile)?;
let (sum_p, sum_t, sum_mw) = izip!(press, t, dp)
.filter_map(|(p, t, dp)| {
if p.is_some() && t.is_some() && dp.is_some() {
Some((p.unpack(), t.unpack(), dp.unpack()))
} else {
None
}
})
.take_while(|&(p, _, _)| p >= bottom_p - HectoPascal(100.0))
.map(|(p, t, dp)| (p, dp, metfor::theta(p, t)))
.filter_map(|(p, dp, th)| metfor::mixing_ratio(dp, p).and_then(|mw| Some((p, th, mw))))
.fold((HectoPascal(0.0), 0.0f64, 0.0f64), |acc, (p, theta, mw)| {
let (sum_p, sum_t, sum_mw) = acc;
(
sum_p + p,
sum_t + theta.unpack() * p.unpack(),
sum_mw + mw * p.unpack(),
)
});
if sum_p == HectoPascal(0.0) {
return Err(AnalysisError::NotEnoughData);
}
let pressure = bottom_p;
let temperature = Celsius::from(metfor::temperature_from_theta(
Kelvin(sum_t / sum_p.unpack()),
bottom_p,
));
let dew_point = metfor::dew_point_from_p_and_mw(bottom_p, sum_mw / sum_p.unpack())
.ok_or(AnalysisError::InvalidInput)?;
Ok(Parcel {
temperature,
pressure,
dew_point,
})
}
pub fn surface_parcel(snd: &Sounding) -> Result<Parcel> {
let pressure = snd.station_pressure().ok_or(AnalysisError::MissingValue)?;
let temperature = snd.sfc_temperature().ok_or(AnalysisError::MissingValue)?;
let dew_point = snd.sfc_dew_point().ok_or(AnalysisError::MissingValue)?;
Ok(Parcel {
temperature,
pressure,
dew_point,
})
}
pub fn lowest_level_parcel(snd: &Sounding) -> Result<Parcel> {
snd.bottom_up()
.filter_map(Parcel::from_datarow)
.nth(0)
.ok_or(AnalysisError::NotEnoughData)
}
pub fn pressure_parcel(snd: &Sounding, pressure: HectoPascal) -> Result<Parcel> {
let row = linear_interpolate_sounding(snd, pressure)?;
Parcel::from_datarow(row).ok_or(AnalysisError::MissingValue)
}
pub fn most_unstable_parcel(snd: &Sounding) -> Result<Parcel> {
let temp_vec: Vec<Optioned<Kelvin>>;
let theta_e = if !snd.theta_e_profile().is_empty() {
snd.theta_e_profile()
} else {
temp_vec = equivalent_potential_temperature(snd);
&temp_vec
};
let press = snd.pressure_profile();
let top_pressure = press
.iter()
.filter_map(|p| p.into_option())
.nth(0)
.ok_or(AnalysisError::NotEnoughData)?
- HectoPascal(300.0);
let (idx, _) = izip!(0.., press, theta_e)
.take_while(|&(_, press_opt, _)| {
if press_opt.is_some() {
press_opt.unpack() >= top_pressure
} else {
true
}
})
.filter_map(|(idx, _, theta_e_opt)| {
if theta_e_opt.is_some() {
Some((idx, theta_e_opt.unpack()))
} else {
None
}
})
.fold(
(::std::usize::MAX, metfor::ABSOLUTE_ZERO),
|(max_idx, max_val), (i, theta_e)| {
if theta_e > max_val {
(i, theta_e)
} else {
(max_idx, max_val)
}
},
);
let row = snd.data_row(idx).ok_or(AnalysisError::NoDataProfile)?;
Parcel::from_datarow(row).ok_or(AnalysisError::MissingValue)
}
pub fn convective_parcel(snd: &Sounding) -> Result<Parcel> {
let Parcel {
dew_point: dp,
pressure: initial_p,
temperature: initial_t,
} = mixed_layer_parcel(snd)?;
let tgt_mw = metfor::mixing_ratio(dp, initial_p).ok_or(AnalysisError::MetForError)?;
let pressure = snd.pressure_profile();
let temperature = snd.temperature_profile();
let (tgt_p, tgt_t) = izip!(pressure, temperature)
.filter_map(|(p, t)| {
let p = p.into_option()?;
let t = t.into_option()?;
let mw = metfor::mixing_ratio(t, p)?;
Some((p, t, mw))
})
.skip_while(|(_, _, mw)| *mw <= tgt_mw)
.scan((0.0, 0.0, 0.0), |(old_p, old_t, old_mw), (p, t, mw)| {
let result = if mw <= tgt_mw && *old_mw >= tgt_mw {
let tgt_p = HectoPascal(linear_interp(tgt_mw, *old_mw, mw, *old_p, p.unpack()));
let tgt_t = Celsius(linear_interp(tgt_mw, *old_mw, mw, *old_t, t.unpack()));
Some(Some((tgt_p, tgt_t)))
} else {
Some(None)
};
*old_p = p.unpack();
*old_t = t.unpack();
*old_mw = mw;
result
})
.filter_map(|opt_opt| opt_opt)
.nth(0)
.ok_or(AnalysisError::NotEnoughData)?;
let tgt_theta = metfor::theta(tgt_p, tgt_t);
let tgt_t = Celsius::from(metfor::temperature_from_theta(tgt_theta, initial_p)).max(initial_t);
Ok(Parcel {
pressure: initial_p,
temperature: tgt_t,
dew_point: dp,
})
}