use crate::gsw_internal_const::*;
use crate::gsw_internal_funcs::*;
use crate::gsw_sp_coefficients::*;
use crate::{Error, Result};
pub(crate) fn t68_from_t90(t90: f64) -> f64 {
t90 * 1.00024
}
pub fn sp_from_c(cndc: f64, t90: f64, p: f64) -> Result<f64> {
let r = if cfg!(feature = "compat") {
cndc * 0.023302418791070513
} else {
cndc / GSW_C3515
};
sp_from_r(r, t90, p)
}
#[cfg(test)]
mod test_sp_from_c {
use super::{sp_from_c, Error};
#[test]
fn zero_cndc() {
let sp = sp_from_c(0.0, 10.0, 100.0).unwrap();
assert!((sp - 0.0).abs() <= f64::EPSILON);
}
#[test]
fn nan() {
let sp = sp_from_c(f64::NAN, 1.0, 1.0);
assert!(sp.unwrap().is_nan());
let sp = sp_from_c(1.0, f64::NAN, 1.0);
assert!(sp.unwrap().is_nan());
let sp = sp_from_c(1.0, 1.0, f64::NAN);
assert!(sp.unwrap().is_nan());
}
#[test]
fn negative_cndc() {
let sp = sp_from_c(-0.1, 10.0, 100.0);
if cfg!(feature = "compat") {
assert!(sp.unwrap().is_nan());
} else {
match sp {
Err(Error::Undefined) => (),
Err(Error::NegativeSalinity) => (),
_ => panic!("can't reach other errors"),
}
}
}
}
pub fn c_from_sp(sp: f64, t90: f64, p: f64) -> Result<f64> {
Ok(GSW_C3515 * r_from_sp(sp, t90, p)?)
}
#[cfg(test)]
mod test_c_from_sp {
use super::{c_from_sp, Error};
#[test]
fn zero_sp() {
let cndc = c_from_sp(0.0, 0.0, 0.0).unwrap();
assert!((cndc - 0.000779962392516606).abs() <= f64::EPSILON);
}
#[test]
fn negative_sp() {
let cndc = c_from_sp(-0.1, 10.0, 100.0);
match cndc {
Err(Error::NegativeSalinity) => (),
_ => panic!("can't reach other errors"),
}
}
#[test]
#[cfg(not(feature = "compat"))]
fn overlimit_sp() {
let cndc = c_from_sp(42.1, 10.0, 100.0);
match cndc {
Err(Error::Undefined) => (),
_ => panic!("can't reach other errors"),
}
}
}
pub fn sp_from_r(r: f64, t90: f64, p: f64) -> Result<f64> {
let t68 = t68_from_t90(t90);
let ft68 = (t68 - 15.0) / (1.0 + K * (t68 - 15.0));
let rt_lc = C0 + (C1 + (C2 + (C3 + C4 * t68) * t68) * t68) * t68;
let rp = 1.0
+ (p * (E1 + E2 * p + E3 * p * p))
/ (1.0 + D1 * t68 + D2 * t68 * t68 + (D3 + D4 * t68) * r);
let rt = r / (rp * rt_lc);
if rt < 0.0 {
if cfg!(feature = "invalidasnan") {
return Ok(f64::NAN);
} else {
return Err(Error::Undefined);
}
}
let rtx = libm::sqrt(rt);
let mut sp = A0
+ (A1 + (A2 + (A3 + (A4 + A5 * rtx) * rtx) * rtx) * rtx) * rtx
+ ft68 * (B0 + (B1 + (B2 + (B3 + (B4 + B5 * rtx) * rtx) * rtx) * rtx) * rtx);
if sp < 2.0 {
let hill_ratio = hill_ratio_at_sp2(t90);
let x = 400.0 * rt;
let sqrty = 10.0 * rtx;
let part1 = 1.0 + x * (1.5 + x);
let part2 = 1.0 + sqrty * (1.0 + sqrty * (1.0 + sqrty));
let sp_hill_raw = sp - A0 / part1 - B0 * ft68 / part2;
sp = hill_ratio * sp_hill_raw;
}
if sp < 0.0 {
if cfg!(feature = "compat") {
Ok(0.0)
} else if cfg!(feature = "invalidasnan") {
Ok(f64::NAN)
} else {
Err(Error::NegativeSalinity)
}
} else {
Ok(sp)
}
}
#[cfg(test)]
mod test_sp_from_r {
use super::sp_from_r;
#[test]
fn nan() {
let sp = sp_from_r(f64::NAN, 1.0, 1.0);
assert!(sp.unwrap().is_nan());
let sp = sp_from_r(1.0, f64::NAN, 1.0);
assert!(sp.unwrap().is_nan());
let sp = sp_from_r(1.0, 1.0, f64::NAN);
assert!(sp.unwrap().is_nan());
}
}
#[cfg(test)]
mod test_round_trip {
use super::{r_from_sp, sp_from_r};
#[test]
fn sp_vs_r() {
let sp_values = [0.1, 10.0, 20.0, 25.0, 30.0, 35.0, 40.0];
let t_values = [0.0, 10.0, 20.0, 30.0, 40.0];
let p_values = [0.0, 10.0, 100.0, 1000.0, 2000.0, 5000.0];
let tol = 2.2e-14;
for sp in sp_values.iter() {
for t in t_values.iter() {
for p in p_values.iter() {
let sp_back = sp_from_r(r_from_sp(*sp, *t, *p).unwrap(), *t, *p).unwrap();
assert!(
(sp - sp_back).abs() <= tol,
"t:{}, p:{}, sp:{}, sp_back:{}, diff: {:+e}",
t,
p,
sp,
sp_back,
sp_back - sp
);
}
}
}
}
}
#[allow(clippy::manual_range_contains)]
pub fn r_from_sp(sp: f64, t90: f64, p: f64) -> Result<f64> {
let t68 = t68_from_t90(t90);
let ft68 = (t68 - 15.0) / (1.0 + K * (t68 - 15.0));
let x = libm::sqrt(sp);
if sp < 0.0 {
return Err(Error::NegativeSalinity);
} else if !cfg!(feature = "compat") && (sp > 42.0) {
return Err(Error::Undefined);
}
let mut rtx = if sp >= 9.0 {
P0 + x
* (P1
+ P4 * t68
+ x * (P3 + P7 * t68 + x * (P6 + P11 * t68 + x * (P10 + P16 * t68 + x * P15))))
+ t68
* (P2
+ t68
* (P5
+ x * x * (P12 + x * P17)
+ P8 * x
+ t68 * (P9 + x * (P13 + x * P18) + t68 * (P14 + P19 * x + P20 * t68))))
} else if sp >= 0.25 && sp < 9.0 {
Q0 + x
* (Q1
+ Q4 * t68
+ x * (Q3 + Q7 * t68 + x * (Q6 + Q11 * t68 + x * (Q10 + Q16 * t68 + x * Q15))))
+ t68
* (Q2
+ t68
* (Q5
+ x * x * (Q12 + x * Q17)
+ Q8 * x
+ t68 * (Q9 + x * (Q13 + x * Q18) + t68 * (Q14 + Q19 * x + Q20 * t68))))
} else if sp >= 0.003 && sp < 0.25 {
R0 + x
* (R1
+ R4 * t68
+ x * (R3 + R7 * t68 + x * (R6 + R11 * t68 + x * (R10 + R16 * t68 + x * R15))))
+ t68
* (R2
+ t68
* (R5
+ x * x * (R12 + x * R17)
+ R8 * x
+ t68 * (R9 + x * (R13 + x * R18) + t68 * (R14 + R19 * x + R20 * t68))))
} else {
U0 + x
* (U1
+ U4 * t68
+ x * (U3 + U7 * t68 + x * (U6 + U11 * t68 + x * (U10 + U16 * t68 + x * U15))))
+ t68
* (U2
+ t68
* (U5
+ x * x * (U12 + x * U17)
+ U8 * x
+ t68 * (U9 + x * (U13 + x * U18) + t68 * (U14 + U19 * x + U20 * t68))))
};
let mut dsp_drtx = A1
+ (2.0 * A2 + (3.0 * A3 + (4.0 * A4 + 5.0 * A5 * rtx) * rtx) * rtx) * rtx
+ ft68 * (B1 + (2.0 * B2 + (3.0 * B3 + (4.0 * B4 + 5.0 * B5 * rtx) * rtx) * rtx) * rtx);
if sp < 2.0 {
let x = 400.0 * (rtx * rtx);
let sqrty = 10.0 * rtx;
let part1 = 1.0 + x * (1.5 + x);
let part2 = 1.0 + sqrty * (1.0 + sqrty * (1.0 + sqrty));
let hill_ratio = hill_ratio_at_sp2(t90);
dsp_drtx = dsp_drtx
+ A0 * 800.0 * rtx * (1.5 + 2.0 * x) / (part1 * part1)
+ B0 * ft68 * (10.0 + sqrty * (20.0 + 30.0 * sqrty)) / (part2 * part2);
dsp_drtx *= hill_ratio;
}
let mut sp_est = A0
+ (A1 + (A2 + (A3 + (A4 + A5 * rtx) * rtx) * rtx) * rtx) * rtx
+ ft68 * (B0 + (B1 + (B2 + (B3 + (B4 + B5 * rtx) * rtx) * rtx) * rtx) * rtx);
if sp_est < 2.0 {
let x = 400.0 * (rtx * rtx);
let sqrty = 10.0 * rtx;
let part1 = 1.0 + x * (1.5e0 + x);
let part2 = 1.0 + sqrty * (1.0 + sqrty * (1.0 + sqrty));
let sp_hill_raw = sp_est - A0 / part1 - B0 * ft68 / part2;
let hill_ratio = hill_ratio_at_sp2(t90);
sp_est = hill_ratio * sp_hill_raw;
}
let rtx_old = rtx;
rtx = rtx_old - (sp_est - sp) / dsp_drtx;
let rtxm = 0.5 * (rtx + rtx_old);
let mut dsp_drtx = A1
+ (2.0 * A2 + (3.0 * A3 + (4.0 * A4 + 5.0 * A5 * rtxm) * rtxm) * rtxm) * rtxm
+ ft68 * (B1 + (2.0 * B2 + (3.0 * B3 + (4.0 * B4 + 5.0 * B5 * rtxm) * rtxm) * rtxm) * rtxm);
if sp_est < 2.0 {
let x = 400.0 * (rtxm * rtxm);
let sqrty = 10.0 * rtxm;
let part1 = 1.0 + x * (1.5e0 + x);
let part2 = 1.0 + sqrty * (1.0 + sqrty * (1.0 + sqrty));
dsp_drtx = dsp_drtx
+ A0 * 800.0 * rtxm * (1.5e0 + 2.0 * x) / (part1 * part1)
+ B0 * ft68 * (10.0 + sqrty * (20.0 + 30.0 * sqrty)) / (part2 * part2);
let hill_ratio = hill_ratio_at_sp2(t90);
dsp_drtx *= hill_ratio;
}
rtx = rtx_old - (sp_est - sp) / dsp_drtx;
sp_est = A0
+ (A1 + (A2 + (A3 + (A4 + A5 * rtx) * rtx) * rtx) * rtx) * rtx
+ ft68 * (B0 + (B1 + (B2 + (B3 + (B4 + B5 * rtx) * rtx) * rtx) * rtx) * rtx);
if sp_est < 2.0 {
let x = 400.0 * (rtx * rtx);
let sqrty = 10.0 * rtx;
let part1 = 1.0 + x * (1.5e0 + x);
let part2 = 1.0 + sqrty * (1.0 + sqrty * (1.0 + sqrty));
let sp_hill_raw = sp_est - A0 / part1 - B0 * ft68 / part2;
let hill_ratio = hill_ratio_at_sp2(t90);
sp_est = hill_ratio * sp_hill_raw;
}
rtx -= (sp_est - sp) / dsp_drtx;
let rt = rtx * rtx;
let aa = D3 + D4 * t68;
let bb = 1.0 + t68 * (D1 + D2 * t68);
let cc = p * (E1 + p * (E2 + E3 * p));
let rt_lc = C0 + (C1 + (C2 + (C3 + C4 * t68) * t68) * t68) * t68;
let dd = bb - aa * rt_lc * rt;
let ee = rt_lc * rt * aa * (bb + cc);
let ra = libm::sqrt(dd * dd + 4.0 * ee) - dd;
let r = 0.5 * ra / aa;
Ok(r)
}
#[cfg(test)]
mod test_r_from_sp {
use super::r_from_sp;
#[test]
fn nan() {
let r = r_from_sp(f64::NAN, 1.0, 1.0);
assert!(r.unwrap().is_nan());
let r = r_from_sp(1.0, f64::NAN, 1.0);
assert!(r.unwrap().is_nan());
let r = r_from_sp(1.0, 1.0, f64::NAN);
assert!(r.unwrap().is_nan());
}
}
pub fn sp_salinometer(rt: f64, t90: f64) -> Result<f64> {
let t68 = t68_from_t90(t90);
let ft68 = (t68 - 15.0) / (1.0 + K * (t68 - 15.0));
if rt < 0.0 {
if cfg!(feature = "invalidasnan") {
return Ok(f64::NAN);
} else {
return Err(Error::Undefined);
}
}
let rtx = libm::sqrt(rt);
let mut sp = A0
+ (A1 + (A2 + (A3 + (A4 + A5 * rtx) * rtx) * rtx) * rtx) * rtx
+ ft68 * (B0 + (B1 + (B2 + (B3 + (B4 + B5 * rtx) * rtx) * rtx) * rtx) * rtx);
if sp < 2.0 {
let hill_ratio = hill_ratio_at_sp2(t90);
let x = 400.0 * rt;
let sqrty = 10.0 * rtx;
let part1 = 1.0 + x * (1.5 + x);
let part2 = 1.0 + sqrty * (1.0 + sqrty * (1.0 + sqrty));
let sp_hill_raw = sp - A0 / part1 - B0 * ft68 / part2;
sp = hill_ratio * sp_hill_raw;
}
if sp < 0.0 {
if cfg!(feature = "compat") {
Ok(0.0)
} else if cfg!(feature = "invalidasnan") {
Ok(f64::NAN)
} else {
Err(Error::NegativeSalinity)
}
} else {
Ok(sp)
}
}
#[cfg(test)]
mod test_sp_salinometer {
use super::{sp_salinometer, Error};
#[test]
fn nan() {
let sp = sp_salinometer(f64::NAN, 1.0);
assert!(sp.unwrap().is_nan());
let sp = sp_salinometer(1.0, f64::NAN);
assert!(sp.unwrap().is_nan());
}
#[test]
fn ratio_one() {
let sp = sp_salinometer(1.0, 15.0).unwrap();
assert!((sp - 35.00000000000001).abs() <= f64::EPSILON);
}
#[test]
fn negative_rt() {
let sp = sp_salinometer(-0.1, 15.0);
if cfg!(feature = "compat") {
assert!(sp.unwrap().is_nan());
} else {
match sp {
Err(Error::Undefined) => (),
Err(Error::NegativeSalinity) => (),
_ => panic!("can't reach other errors"),
}
}
}
}