const TWO_PI: f64 = std::f64::consts::TAU;
pub(crate) fn m(sinφ: f64, cosφ: f64, e_sq: f64) -> f64 {
cosφ / (1. - e_sq * sinφ * sinφ).sqrt()
}
pub(crate) fn t(cosφ: f64, sinφ: f64, e: f64) -> f64 {
(e * (e * sinφ).atanh()).exp()
* if sinφ > 0. {
cosφ / (1. + sinφ)
} else {
(1. - sinφ) / cosφ
}
}
pub(crate) fn phi2(ts0: f64, e: f64) -> Option<f64> {
let phi2 = sinhpsi2tanphi((1. / ts0 - ts0) / 2., e)?.atan();
Some(phi2)
}
pub(crate) fn sinhpsi2tanphi(taup: f64, e: f64) -> Option<f64> {
const MAX_ITER: usize = 5;
let root_eps: f64 = f64::EPSILON.sqrt();
let tol: f64 = root_eps / 10.; let tmax: f64 = 2. / root_eps; let e2m: f64 = 1. - e * e;
let stol: f64 = tol * 1.0_f64.max(taup.abs());
let mut tau = if taup.abs() > 70. {
taup * (e * e.atanh()).exp()
} else {
taup / e2m
};
if tau.abs() >= tmax {
return Some(tau);
}
let mut count = MAX_ITER;
while count > 0 {
let tau1 = (1. + tau * tau).sqrt();
let sig = (e * (e * tau / tau1).atanh()).sinh();
let taupa = (1. + sig * sig).sqrt() * tau - sig * tau1;
let dtau =
(taup - taupa) * (1. + e2m * (tau * tau)) / (e2m * tau1 * (1. + taupa * taupa).sqrt());
tau += dtau;
if dtau.abs() < stol {
return Some(tau);
}
count -= 1;
}
None
}
pub(crate) fn normalize_longitude(lambda: f64) -> f64 {
(lambda + std::f64::consts::PI).rem_euclid(TWO_PI) - std::f64::consts::PI
}