geodesy/inner_op/
geodesic.rs

1/// Geodesics
2use crate::authoring::*;
3
4// ----- F O R W A R D -----------------------------------------------------------------
5
6fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
7    let ellps = op.params.ellps(0);
8
9    let n = operands.len();
10    let sliced = 0..n;
11
12    let mut successes = 0_usize;
13    for i in sliced {
14        let args = operands.get_coord(i);
15        let origin = Coor2D::geo(args[0], args[1]);
16        let azimuth = args[2].to_radians();
17        let distance = args[3];
18
19        let destination = ellps.geodesic_fwd(&origin, azimuth, distance).to_degrees();
20
21        // No convergence?
22        if destination[3] > 990.0 {
23            operands.set_coord(i, &Coor4D::nan());
24            continue;
25        }
26
27        let result = Coor4D([destination[1], destination[0], args[0], args[1]]);
28        operands.set_coord(i, &result);
29        successes += 1;
30    }
31
32    successes
33}
34
35// ----- I N V E R S E -----------------------------------------------------------------
36
37fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
38    let ellps = op.params.ellps(0);
39    let reversible = op.params.boolean("reversible");
40
41    let n = operands.len();
42    let sliced = 0..n;
43
44    let mut successes = 0_usize;
45    for i in sliced {
46        let mut from = Coor2D::origin();
47        let mut to = Coor2D::origin();
48
49        let coord = operands.get_coord(i);
50        from[0] = coord[1].to_radians();
51        from[1] = coord[0].to_radians();
52        to[0] = coord[3].to_radians();
53        to[1] = coord[2].to_radians();
54
55        let mut geodesic = ellps.geodesic_inv(&from, &to).to_degrees();
56
57        // No convergence?
58        if geodesic[3] > 990.0 {
59            operands.set_coord(i, &Coor4D::nan());
60            continue;
61        }
62        geodesic[3] = (geodesic[1] + 180.0) % 360.0;
63
64        if reversible {
65            let distance = geodesic[2];
66            let return_azi = geodesic[3];
67            operands.set_coord(i, &Coor4D::raw(coord[2], coord[3], return_azi, distance));
68            continue;
69        }
70
71        operands.set_coord(i, &geodesic);
72        successes += 1;
73    }
74
75    successes
76}
77
78// ----- C O N S T R U C T O R ---------------------------------------------------------
79
80// Example...
81#[rustfmt::skip]
82pub const GAMUT: [OpParameter; 3] = [
83    OpParameter::Flag { key: "inv" },
84    OpParameter::Flag { key: "reversible" },
85    OpParameter::Text { key: "ellps", default: Some("GRS80") }
86];
87
88pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error> {
89    let op = Op::basic(parameters, InnerOp(fwd), Some(InnerOp(inv)), &GAMUT)?;
90    Ok(op)
91}
92
93// ----- T E S T S ---------------------------------------------------------------------
94
95#[cfg(test)]
96mod tests {
97    use super::*;
98
99    #[test]
100    fn geodesic() -> Result<(), Error> {
101        let mut ctx = Minimal::default();
102
103        // Approximate coordinates of Copenhagen and Paris airports
104        let cph_cdg = Coor4D::raw(55., 12., 49., 2.);
105
106        // A geodesic from Copenhagen to Paris
107        let op = ctx.op("geodesic")?;
108        let mut operands = [cph_cdg];
109        ctx.apply(op, Inv, &mut operands)?;
110
111        let expected = Coor4D([
112            -130.1540604203936,
113            -138.05257941840648,
114            956066.2319619625,
115            41.94742058159352,
116        ]);
117
118        assert!((operands[0][0] - expected[0]).abs() < 1e-9);
119        assert!((operands[0][1] - expected[1]).abs() < 1e-9);
120        assert!((operands[0][2] - expected[2]).abs() < 1e-9);
121        assert!((operands[0][3] - expected[3]).abs() < 1e-9);
122
123        // A geodesic from Copenhagen to Paris in the "reversible" format
124        let op = ctx.op("geodesic reversible")?;
125        let mut operands = [cph_cdg];
126        ctx.apply(op, Inv, &mut operands)?;
127
128        let expected = Coor4D([49.0, 2.0, 41.94742058159352, 956066.2319619625]);
129
130        assert!((operands[0][0] - expected[0]).abs() < 1e-9);
131        assert!((operands[0][1] - expected[1]).abs() < 1e-9);
132        assert!((operands[0][2] - expected[2]).abs() < 1e-9);
133        assert!((operands[0][3] - expected[3]).abs() < 1e-9);
134
135        // And back to Copenhagen!
136        ctx.apply(op, Fwd, &mut operands)?;
137
138        assert!((operands[0][0] - cph_cdg[0]).abs() < 1e-10);
139        assert!((operands[0][1] - cph_cdg[1]).abs() < 1e-10);
140        assert!((operands[0][2] - cph_cdg[2]).abs() < 1e-10);
141        assert!((operands[0][3] - cph_cdg[3]).abs() < 1e-10);
142
143        Ok(())
144    }
145}
geodesy/inner_op/geodesic.rs

geodesy/inner_op/
geodesic.rs
