From 3b456604c0d56aec88434bd11f51478e875c922a Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Mon, 14 Aug 2023 12:41:13 +0200 Subject: [PATCH] Refactor math and preludes, mostly to improve the doc structure --- src/bin/kp.rs | 2 +- src/ellipsoid/latitudes.rs | 4 +- src/inner_op/laea.rs | 8 +- src/inner_op/lcc.rs | 14 +-- src/lib.rs | 75 ++++++++----- src/math/ancillary.rs | 126 ++++++++++++++++++++++ src/math/angular.rs | 58 +++++++++- src/math/mod.rs | 208 ++---------------------------------- src/op/parsed_parameters.rs | 7 +- 9 files changed, 257 insertions(+), 245 deletions(-) create mode 100644 src/math/ancillary.rs diff --git a/src/bin/kp.rs b/src/bin/kp.rs index 8c69809..f517ba4 100644 --- a/src/bin/kp.rs +++ b/src/bin/kp.rs @@ -123,7 +123,7 @@ fn main() -> Result<(), anyhow::Error> { args.extend(&(["0", "0", "0", "NaN", "0"][args.len()..])); let mut b: Vec = vec![]; for e in args { - b.push(parse_sexagesimal(e)); + b.push(geodesy::math::angular::parse_sexagesimal(e)); } b[2] = options.height.unwrap_or(b[2]); b[3] = options.time.unwrap_or(b[3]); diff --git a/src/ellipsoid/latitudes.rs b/src/ellipsoid/latitudes.rs index 1a8e83f..f90207d 100644 --- a/src/ellipsoid/latitudes.rs +++ b/src/ellipsoid/latitudes.rs @@ -1,5 +1,5 @@ use super::*; -use crate::math::*; +use crate::authoring::*; // ----- Latitudes ------------------------------------------------------------- impl Ellipsoid { @@ -48,7 +48,7 @@ impl Ellipsoid { #[must_use] pub fn latitude_isometric_to_geographic(&self, isometric: f64) -> f64 { let e = self.eccentricity(); - sinhpsi_to_tanphi(isometric.sinh(), e).atan() + ancillary::sinhpsi_to_tanphi(isometric.sinh(), e).atan() } // --- Auxiliary latitudes --- diff --git a/src/inner_op/laea.rs b/src/inner_op/laea.rs index 1aebcf5..3934420 100644 --- a/src/inner_op/laea.rs +++ b/src/inner_op/laea.rs @@ -42,7 +42,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let lon = coord[0]; let (sin_lon, cos_lon) = (lon - lon_0).sin_cos(); - let q = qs(lat.sin(), e); + let q = ancillary::qs(lat.sin(), e); let rho = a * (qp + sign * q).sqrt(); coord[0] = x_0 + rho * sin_lon; @@ -60,7 +60,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let (sin_lon, cos_lon) = (lon - lon_0).sin_cos(); // Authalic latitude, 𝜉 - let xi = (qs(lat.sin(), e) / qp).asin(); + let xi = (ancillary::qs(lat.sin(), e) / qp).asin(); let (sin_xi, cos_xi) = xi.sin_cos(); let b = if oblique { @@ -224,9 +224,9 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result let (sin_phi_0, cos_phi_0) = lat_0.sin_cos(); // qs for the central parallel - let q0 = qs(sin_phi_0, e); + let q0 = ancillary::qs(sin_phi_0, e); // qs for the North Pole - let qp = qs(1.0, e); + let qp = ancillary::qs(1.0, e); // Authalic latitude of the central parallel - 𝛽₀ in the IOGP text let xi_0 = (q0 / qp).asin(); // Rq in the IOGP text diff --git a/src/inner_op/lcc.rs b/src/inner_op/lcc.rs index 4255fb6..33b94f7 100644 --- a/src/inner_op/lcc.rs +++ b/src/inner_op/lcc.rs @@ -36,7 +36,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { continue; } } else { - rho = c * crate::math::ts(phi.sin_cos(), e).powf(n); + rho = c * crate::math::ancillary::ts(phi.sin_cos(), e).powf(n); } let sc = (lam * n).sin_cos(); coord[0] = a * k_0 * rho * sc.0 + x_0; @@ -86,7 +86,7 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { } let ts0 = (rho / c).powf(1. / n); - let phi = crate::math::pj_phi2(ts0, e); + let phi = crate::math::ancillary::pj_phi2(ts0, e); if phi.is_infinite() || phi.is_nan() { coord = Coor4D::nan(); operands.set_coord(i, &coord); @@ -170,19 +170,19 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result } // Snyder (1982) eq. 12-15 - let m1 = crate::math::pj_msfn(sc, es); + let m1 = crate::math::ancillary::pj_msfn(sc, es); // Snyder (1982) eq. 7-10: exp(-𝜓) - let ml1 = crate::math::ts(sc, e); + let ml1 = crate::math::ancillary::ts(sc, e); // Secant case? if (phi1 - phi2).abs() >= EPS10 { let sc = phi2.sin_cos(); - n = (m1 / crate::math::pj_msfn(sc, es)).ln(); + n = (m1 / crate::math::ancillary::pj_msfn(sc, es)).ln(); if n == 0. { return Err(Error::General("Lcc: Invalid value for eccentricity")); } - let ml2 = crate::math::ts(sc, e); + let ml2 = crate::math::ancillary::ts(sc, e); let denom = (ml1 / ml2).ln(); if denom == 0. { return Err(Error::General("Lcc: Invalid value for eccentricity")); @@ -193,7 +193,7 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result let c = m1 * ml1.powf(-n) / n; let mut rho0 = 0.; if (lat_0.abs() - FRAC_PI_2).abs() > EPS10 { - rho0 = c * crate::math::ts(lat_0.sin_cos(), e).powf(n); + rho0 = c * crate::math::ancillary::ts(lat_0.sin_cos(), e).powf(n); } params.real.insert("c", c); diff --git a/src/lib.rs b/src/lib.rs index 6b2b40a..9f0c1c3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -6,54 +6,70 @@ mod coordinate; mod ellipsoid; mod grid; mod inner_op; -mod math; +pub mod math; mod op; // The bread-and-butter -pub use crate::context::minimal::Minimal; -pub use crate::context::Context; +// Context trait and types +pub use crate::context::minimal::Minimal; #[cfg(feature = "with_plain")] pub use crate::context::plain::Plain; +pub use crate::context::Context; +// Used to specify which operator to apply in `Context::apply(...)` +pub use crate::op::OpHandle; + +// Coordinate traits and types pub use crate::coordinate::coor2d::Coor2D; pub use crate::coordinate::coor32::Coor32; pub use crate::coordinate::coor3d::Coor3D; pub use crate::coordinate::coor4d::Coor4D; +pub use crate::coordinate::AngularUnits; pub use crate::coordinate::CoordinateMetadata; pub use crate::coordinate::CoordinateSet; + +// Additional types pub use crate::ellipsoid::Ellipsoid; +pub use crate::math::jacobian::Factors; +pub use crate::math::jacobian::Jacobian; + +// Constants - directionality in `Context::apply(...)` pub use crate::Direction::Fwd; pub use crate::Direction::Inv; -pub use crate::math::jacobian::Factors; -pub use crate::math::jacobian::Jacobian; +// PROJ interoperability +pub use crate::op::parse_proj; /// The bread-and-butter, shrink-wrapped for external use pub mod prelude { - pub use crate::context::minimal::Minimal; - pub use crate::context::Context; - + pub use crate::Context; + pub use crate::Minimal; #[cfg(feature = "with_plain")] - pub use crate::context::plain::Plain; - - pub use crate::coordinate::coor2d::Coor2D; - pub use crate::coordinate::coor32::Coor32; - pub use crate::coordinate::coor3d::Coor3D; - pub use crate::coordinate::coor4d::Coor4D; - pub use crate::coordinate::AngularUnits; - pub use crate::coordinate::CoordinateMetadata; - pub use crate::coordinate::CoordinateSet; - pub use crate::ellipsoid::Ellipsoid; - pub use crate::math::jacobian::Factors; - pub use crate::math::jacobian::Jacobian; - pub use crate::math::parse_sexagesimal; - pub use crate::op::parse_proj; - pub use crate::op::OpHandle; + pub use crate::Plain; + + pub use crate::AngularUnits; + pub use crate::CoordinateMetadata; + pub use crate::CoordinateSet; + + pub use crate::Coor2D; + pub use crate::Coor32; + pub use crate::Coor3D; + pub use crate::Coor4D; + + pub use crate::Ellipsoid; + pub use crate::Factors; + pub use crate::Jacobian; + + pub use crate::OpHandle; + pub use crate::Direction; pub use crate::Direction::Fwd; pub use crate::Direction::Inv; pub use crate::Error; + + pub use crate::parse_proj; + #[cfg(test)] pub fn some_basic_coor4dinates() -> [Coor4D; 2] { let copenhagen = Coor4D::raw(55., 12., 0., 0.); @@ -74,13 +90,9 @@ pub mod prelude { } } -/// Preamble for authoring Contexts and InnerOp modules (built-in or user defined) +/// Prelude for authoring Contexts and InnerOp modules (built-in or user defined) pub mod authoring { pub use crate::prelude::*; - pub use log::error; - pub use log::info; - pub use log::trace; - pub use log::warn; pub use crate::grid::Grid; pub use crate::inner_op::InnerOp; @@ -93,7 +105,14 @@ pub mod authoring { pub use crate::op::ParsedParameters; pub use crate::op::RawParameters; + // All new contexts are supposed to support these pub use crate::context::BUILTIN_ADAPTORS; + + // External material + pub use log::error; + pub use log::info; + pub use log::trace; + pub use log::warn; pub use std::collections::BTreeMap; } diff --git a/src/math/ancillary.rs b/src/math/ancillary.rs new file mode 100644 index 0000000..e557c60 --- /dev/null +++ b/src/math/ancillary.rs @@ -0,0 +1,126 @@ +/// The Gudermannian function (often written as gd), is the work horse for computations involving +/// the isometric latitude (i.e. the vertical coordinate of the Mercator projection) +pub mod gudermannian { + pub fn fwd(arg: f64) -> f64 { + arg.sinh().atan() + } + + pub fn inv(arg: f64) -> f64 { + arg.tan().asinh() + } +} + +/// ts is the equivalent of Charles Karney's PROJ function `pj_tsfn`. +/// It determines the function ts(phi) as defined in Snyder (1987), +/// Eq. (7-10) +/// +/// ts is the exponential of the negated isometric latitude, i.e. +/// exp(-𝜓), but evaluated in a numerically more stable way than +/// the naive ellps.isometric_latitude(...).exp() +/// +/// This version is essentially identical to Charles Karney's PROJ +/// version, including the majority of the comments. +/// +/// Inputs: +/// (sin 𝜙, cos 𝜙): trigs of geographic latitude +/// e: eccentricity of the ellipsoid +/// Output: +/// ts: exp(-𝜓) = 1 / (tan 𝜒 + sec 𝜒) +/// where 𝜓 is the isometric latitude (dimensionless) +/// and 𝜒 is the conformal latitude (radians) +/// +/// Here the isometric latitude is defined by +/// 𝜓 = log( +/// tan(𝜋/4 + 𝜙/2) * +/// ( (1 - e × sin 𝜙) / (1 + e × sin 𝜙) ) ^ (e/2) +/// ) +/// = asinh(tan 𝜙) - e × atanh(e × sin 𝜙) +/// = asinh(tan 𝜒) +/// +/// where 𝜒 is the conformal latitude +/// +pub fn ts(sincos: (f64, f64), e: f64) -> f64 { + // exp(-asinh(tan 𝜙)) + // = 1 / (tan 𝜙 + sec 𝜙) + // = cos 𝜙 / (1 + sin 𝜙) good for 𝜙 > 0 + // = (1 - sin 𝜙) / cos 𝜙 good for 𝜙 < 0 + let factor = if sincos.0 > 0. { + sincos.1 / (1. + sincos.0) + } else { + (1. - sincos.0) / sincos.1 + }; + (e * (e * sincos.0).atanh()).exp() * factor +} + +/// Snyder (1982) eq. 12-15, PROJ's pj_msfn() +pub fn pj_msfn(sincos: (f64, f64), es: f64) -> f64 { + sincos.1 / (1. - sincos.0 * sincos.0 * es).sqrt() +} + +/// Equivalent to the PROJ pj_phi2 function +pub fn pj_phi2(ts0: f64, e: f64) -> f64 { + sinhpsi_to_tanphi((1. / ts0 - ts0) / 2., e).atan() +} + +/// Snyder (1982) eq. ??, PROJ's pj_qsfn() +pub fn qs(sinphi: f64, e: f64) -> f64 { + let es = e * e; + let one_es = 1.0 - es; + + if e < 1e-7 { + return 2.0 * sinphi; + } + + let con = e * sinphi; + let div1 = 1.0 - con * con; + let div2 = 1.0 + con; + + one_es * (sinphi / div1 - (0.5 / e) * ((1. - con) / div2).ln()) +} + +/// Ancillary function for computing the inverse isometric latitude. Follows +/// [Karney, 2011](crate::Bibliography::Kar11), and the PROJ implementation +/// in proj/src/phi2.cpp. +pub fn sinhpsi_to_tanphi(taup: f64, e: f64) -> f64 { + // min iterations = 1, max iterations = 2; mean = 1.954 + const MAX_ITER: usize = 5; + + // rooteps, tol and tmax are compile time constants, but currently + // Rust cannot const-evaluate powers and roots, so we must either + // evaluate these "constants" as lazy_statics, or just swallow the + // penalty of an extra sqrt and two divisions on each call. + // If this shows unbearable, we can just also assume IEEE-64 bit + // arithmetic, and set rooteps = 0.000000014901161193847656 + let rooteps: f64 = f64::EPSILON.sqrt(); + let tol: f64 = rooteps / 10.; // the criterion for Newton's method + let tmax: f64 = 2. / rooteps; // threshold for large arg limit exact + + let e2m = 1. - e * e; + let stol = tol * taup.abs().max(1.0); + + // The initial guess. 70 corresponds to chi = 89.18 deg + let mut tau = if taup.abs() > 70. { + taup * (e * e.atanh()).exp() + } else { + taup / e2m + }; + + // Handle +/-inf, nan, and e = 1 + if (tau.abs() >= tmax) || tau.is_nan() { + return tau; + } + + for _ in 0..MAX_ITER { + 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) || tau.is_nan() { + return tau; + } + } + f64::NAN +} diff --git a/src/math/angular.rs b/src/math/angular.rs index 2be58f3..da9e564 100644 --- a/src/math/angular.rs +++ b/src/math/angular.rs @@ -1,3 +1,5 @@ +use log::warn; + /// Simplistic transformation from degrees, minutes and seconds-with-decimals /// to degrees-with-decimals. No sanity check: Sign taken from degree-component, /// minutes forced to unsigned by i16 type, but passing a negative value for @@ -63,14 +65,14 @@ pub fn dd_to_iso_dms(dd: f64) -> f64 { sign * (d * 10000. + m * 100. + s) } -/// normalize arbitrary angles to [-π, π): +/// normalize arbitrary angles to [-π, π) pub fn normalize_symmetric(angle: f64) -> f64 { use std::f64::consts::PI; let angle = (angle + PI) % (2.0 * PI); angle - PI * angle.signum() } -/// normalize arbitrary angles to [0, 2π): +/// normalize arbitrary angles to [0, 2π) pub fn normalize_positive(angle: f64) -> f64 { use std::f64::consts::PI; let angle = angle % (2.0 * PI); @@ -80,6 +82,47 @@ pub fn normalize_positive(angle: f64) -> f64 { angle } +/// Parse sexagesimal degrees, i.e. degrees, minutes and seconds in the +/// format 45:30:36, 45:30:36N,-45:30:36 etc. +pub fn parse_sexagesimal(angle: &str) -> f64 { + // Degrees, minutes, and seconds + let mut dms = [0.0, 0.0, 0.0]; + let mut angle = angle.trim(); + + // Empty? + let n = angle.len(); + if n == 0 || angle == "NaN" { + return f64::NAN; + } + + // Handle NSEW indicators + let mut postfix_sign = 1.0; + if "wWsSeEnN".contains(&angle[n - 1..]) { + if "wWsS".contains(&angle[n - 1..]) { + postfix_sign = -1.0; + } + angle = &angle[..n - 1]; + } + + // Split into as many elements as given: D, D:M, D:M:S + for (i, element) in angle.split(':').enumerate() { + if i < 3 { + if let Ok(v) = element.parse::() { + dms[i] = v; + continue; + } + } + // More than 3 elements? + warn!("Cannot parse {angle} as a real number or sexagesimal angle"); + return f64::NAN; + } + + // Sexagesimal conversion if we have more than one element. Otherwise + // decay gracefully to plain real/f64 conversion + let sign = dms[0].signum() * postfix_sign; + sign * (dms[0].abs() + (dms[1] + dms[2] / 60.0) / 60.0) +} + // ----- Tests --------------------------------------------------------------------- #[cfg(test)] @@ -107,4 +150,15 @@ mod tests { assert_eq!(iso_dm_to_dd(5530.60), -iso_dm_to_dd(-5530.60)); assert_eq!(iso_dms_to_dd(553036.), -iso_dms_to_dd(-553036.00)); } + + #[test] + fn test_parse_sexagesimal() { + assert_eq!(1.51, parse_sexagesimal("1:30:36")); + assert_eq!(-1.51, parse_sexagesimal("-1:30:36")); + assert_eq!(1.51, parse_sexagesimal("1:30:36N")); + assert_eq!(-1.51, parse_sexagesimal("1:30:36S")); + assert_eq!(1.51, parse_sexagesimal("1:30:36e")); + assert_eq!(-1.51, parse_sexagesimal("1:30:36w")); + assert!(parse_sexagesimal("q1:30:36w").is_nan()); + } } diff --git a/src/math/mod.rs b/src/math/mod.rs index 610a499..2fe9c97 100644 --- a/src/math/mod.rs +++ b/src/math/mod.rs @@ -1,208 +1,20 @@ +/// Free functions used in more than one module of the crate. +pub mod ancillary; +pub use ancillary::gudermannian; + +/// Free functions for handling and converting between +/// different representations of angles. pub mod angular; +/// Computations involving the Jacobian matrix for investigation +/// of the geometrical properties of map projections. pub mod jacobian; -pub use jacobian::Factors; -pub use jacobian::Jacobian; +/// Fourier- and Taylor series pub mod series; pub use series::fourier; - pub use series::taylor; -pub use series::taylor::fourier_coefficients; +pub use series::taylor::fourier_coefficients; pub use series::FourierCoefficients; pub use series::PolynomialCoefficients; - -use log::warn; - -/// The Gudermannian function (often written as gd), is the work horse for computations involving -/// the isometric latitude (i.e. the vertical coordinate of the Mercator projection) -pub mod gudermannian { - pub fn fwd(arg: f64) -> f64 { - arg.sinh().atan() - } - - pub fn inv(arg: f64) -> f64 { - arg.tan().asinh() - } -} - -// ----- A N C I L L A R Y F U N C T I O N S ----------------------------------------- - -// ts is the equivalent of Charles Karney's PROJ function `pj_tsfn`. -// It determines the function ts(phi) as defined in Snyder (1987), -// Eq. (7-10) -// -// ts is the exponential of the negated isometric latitude, i.e. -// exp(-𝜓), but evaluated in a numerically more stable way than -// the naive ellps.isometric_latitude(...).exp() -// -// This version is essentially identical to Charles Karney's PROJ -// version, including the majority of the comments. -// -// Inputs: -// (sin 𝜙, cos 𝜙): trigs of geographic latitude -// e: eccentricity of the ellipsoid -// Output: -// ts: exp(-𝜓) = 1 / (tan 𝜒 + sec 𝜒) -// where 𝜓 is the isometric latitude (dimensionless) -// and 𝜒 is the conformal latitude (radians) -// -// Here the isometric latitude is defined by -// 𝜓 = log( -// tan(𝜋/4 + 𝜙/2) * -// ( (1 - e × sin 𝜙) / (1 + e × sin 𝜙) ) ^ (e/2) -// ) -// = asinh(tan 𝜙) - e × atanh(e × sin 𝜙) -// = asinh(tan 𝜒) -// -// where 𝜒 is the conformal latitude -// -pub(crate) fn ts(sincos: (f64, f64), e: f64) -> f64 { - // exp(-asinh(tan 𝜙)) - // = 1 / (tan 𝜙 + sec 𝜙) - // = cos 𝜙 / (1 + sin 𝜙) good for 𝜙 > 0 - // = (1 - sin 𝜙) / cos 𝜙 good for 𝜙 < 0 - let factor = if sincos.0 > 0. { - sincos.1 / (1. + sincos.0) - } else { - (1. - sincos.0) / sincos.1 - }; - (e * (e * sincos.0).atanh()).exp() * factor -} - -// Snyder (1982) eq. 12-15, PROJ's pj_msfn() -pub(crate) fn pj_msfn(sincos: (f64, f64), es: f64) -> f64 { - sincos.1 / (1. - sincos.0 * sincos.0 * es).sqrt() -} - -// Equivalent to the PROJ pj_phi2 function -pub(crate) fn pj_phi2(ts0: f64, e: f64) -> f64 { - sinhpsi_to_tanphi((1. / ts0 - ts0) / 2., e).atan() -} - -// Snyder (1982) eq. ??, PROJ's pj_qsfn() -pub(crate) fn qs(sinphi: f64, e: f64) -> f64 { - let es = e * e; - let one_es = 1.0 - es; - - if e < 1e-7 { - return 2.0 * sinphi; - } - - let con = e * sinphi; - let div1 = 1.0 - con * con; - let div2 = 1.0 + con; - - one_es * (sinphi / div1 - (0.5 / e) * ((1. - con) / div2).ln()) -} - -// Ancillary function for computing the inverse isometric latitude. Follows -// [Karney, 2011](crate::Bibliography::Kar11), and the PROJ implementation -// in proj/src/phi2.cpp. -// Needs crate-visibility as it is also used in crate::ellipsoid::latitudes -pub(crate) fn sinhpsi_to_tanphi(taup: f64, e: f64) -> f64 { - // min iterations = 1, max iterations = 2; mean = 1.954 - const MAX_ITER: usize = 5; - - // rooteps, tol and tmax are compile time constants, but currently - // Rust cannot const-evaluate powers and roots, so we must either - // evaluate these "constants" as lazy_statics, or just swallow the - // penalty of an extra sqrt and two divisions on each call. - // If this shows unbearable, we can just also assume IEEE-64 bit - // arithmetic, and set rooteps = 0.000000014901161193847656 - let rooteps: f64 = f64::EPSILON.sqrt(); - let tol: f64 = rooteps / 10.; // the criterion for Newton's method - let tmax: f64 = 2. / rooteps; // threshold for large arg limit exact - - let e2m = 1. - e * e; - let stol = tol * taup.abs().max(1.0); - - // The initial guess. 70 corresponds to chi = 89.18 deg - let mut tau = if taup.abs() > 70. { - taup * (e * e.atanh()).exp() - } else { - taup / e2m - }; - - // Handle +/-inf, nan, and e = 1 - if (tau.abs() >= tmax) || tau.is_nan() { - return tau; - } - - for _ in 0..MAX_ITER { - 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) || tau.is_nan() { - return tau; - } - } - f64::NAN -} - -/// Parse sexagesimal degrees, i.e. degrees, minutes and seconds in the -/// format 45:30:36, 45:30:36N,-45:30:36 etc. -pub fn parse_sexagesimal(angle: &str) -> f64 { - // Degrees, minutes, and seconds - let mut dms = [0.0, 0.0, 0.0]; - let mut angle = angle.trim(); - - // Empty? - let n = angle.len(); - if n == 0 || angle == "NaN" { - return f64::NAN; - } - - // Handle NSEW indicators - let mut postfix_sign = 1.0; - if "wWsSeEnN".contains(&angle[n - 1..]) { - if "wWsS".contains(&angle[n - 1..]) { - postfix_sign = -1.0; - } - angle = &angle[..n - 1]; - } - - // Split into as many elements as given: D, D:M, D:M:S - for (i, element) in angle.split(':').enumerate() { - if i < 3 { - if let Ok(v) = element.parse::() { - dms[i] = v; - continue; - } - } - // More than 3 elements? - warn!("Cannot parse {angle} as a real number or sexagesimal angle"); - return f64::NAN; - } - - // Sexagesimal conversion if we have more than one element. Otherwise - // decay gracefully to plain real/f64 conversion - let sign = dms[0].signum() * postfix_sign; - sign * (dms[0].abs() + (dms[1] + dms[2] / 60.0) / 60.0) -} - -// ----- Tests --------------------------------------------------------------------- - -#[cfg(test)] -mod tests { - use super::*; - use crate::Error; - - #[test] - fn test_parse_sexagesimal() -> Result<(), Error> { - assert_eq!(1.51, parse_sexagesimal("1:30:36")); - assert_eq!(-1.51, parse_sexagesimal("-1:30:36")); - assert_eq!(1.51, parse_sexagesimal("1:30:36N")); - assert_eq!(-1.51, parse_sexagesimal("1:30:36S")); - assert_eq!(1.51, parse_sexagesimal("1:30:36e")); - assert_eq!(-1.51, parse_sexagesimal("1:30:36w")); - assert!(parse_sexagesimal("q1:30:36w").is_nan()); - - Ok(()) - } -} diff --git a/src/op/parsed_parameters.rs b/src/op/parsed_parameters.rs index b64daf4..076f933 100644 --- a/src/op/parsed_parameters.rs +++ b/src/op/parsed_parameters.rs @@ -6,6 +6,7 @@ // store them in the `ParsedParameters` struct, ready for use at run time. #![allow(clippy::needless_range_loop)] +use crate::math::angular; use crate::math::FourierCoefficients; use std::collections::BTreeSet; @@ -206,7 +207,7 @@ impl ParsedParameters { OpParameter::Real { key, default } => { if let Some(value) = chase(globals, &locals, key)? { - let v = crate::math::parse_sexagesimal(&value); + let v = angular::parse_sexagesimal(&value); if v.is_nan() { return Err(Error::BadParam(key.to_string(), value)); } @@ -231,7 +232,7 @@ impl ParsedParameters { let mut elements = Vec::::new(); if let Some(value) = chase(globals, &locals, key)? { for element in value.split(',') { - let v = crate::math::parse_sexagesimal(element); + let v = angular::parse_sexagesimal(element); if v.is_nan() { warn!("Cannot parse {key}:{value} as a series"); return Err(Error::BadParam(key.to_string(), value.to_string())); @@ -252,7 +253,7 @@ impl ParsedParameters { continue; } for element in value.split(',') { - let v = crate::math::parse_sexagesimal(element); + let v = angular::parse_sexagesimal(element); if v.is_nan() { warn!("Cannot parse {key}:{value} as a series"); return Err(Error::BadParam(key.to_string(), value.to_string()));