geodesy 0.15.0

A platform for experiments with geodetic transformations and data flow
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

/// Geodesics
use crate::authoring::*;

// ----- F O R W A R D -----------------------------------------------------------------

fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
    let ellps = op.params.ellps(0);

    let n = operands.len();
    let sliced = 0..n;

    let mut successes = 0_usize;
    for i in sliced {
        let args = operands.get_coord(i);
        let origin = Coor2D::geo(args[0], args[1]);
        let azimuth = args[2].to_radians();
        let distance = args[3];

        let destination = ellps.geodesic_fwd(&origin, azimuth, distance).to_degrees();

        // No convergence?
        if destination[3] > 990.0 {
            operands.set_coord(i, &Coor4D::nan());
            continue;
        }

        let result = Coor4D([destination[1], destination[0], args[0], args[1]]);
        operands.set_coord(i, &result);
        successes += 1;
    }

    successes
}

// ----- I N V E R S E -----------------------------------------------------------------

fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
    let ellps = op.params.ellps(0);
    let reversible = op.params.boolean("reversible");

    let n = operands.len();
    let sliced = 0..n;

    let mut successes = 0_usize;
    for i in sliced {
        let mut from = Coor2D::origin();
        let mut to = Coor2D::origin();

        let coord = operands.get_coord(i);
        from[0] = coord[1].to_radians();
        from[1] = coord[0].to_radians();
        to[0] = coord[3].to_radians();
        to[1] = coord[2].to_radians();

        let mut geodesic = ellps.geodesic_inv(&from, &to).to_degrees();

        // No convergence?
        if geodesic[3] > 990.0 {
            operands.set_coord(i, &Coor4D::nan());
            continue;
        }
        geodesic[3] = (geodesic[1] + 180.0) % 360.0;

        if reversible {
            let distance = geodesic[2];
            let return_azi = geodesic[3];
            operands.set_coord(i, &Coor4D::raw(coord[2], coord[3], return_azi, distance));
            continue;
        }

        operands.set_coord(i, &geodesic);
        successes += 1;
    }

    successes
}

// ----- C O N S T R U C T O R ---------------------------------------------------------

// Example...
#[rustfmt::skip]
pub const GAMUT: [OpParameter; 3] = [
    OpParameter::Flag { key: "inv" },
    OpParameter::Flag { key: "reversible" },
    OpParameter::Text { key: "ellps", default: Some("GRS80") }
];

pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error> {
    let op = Op::basic(parameters, InnerOp(fwd), Some(InnerOp(inv)), &GAMUT)?;
    Ok(op)
}

// ----- T E S T S ---------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn geodesic() -> Result<(), Error> {
        let mut ctx = Minimal::default();

        // Approximate coordinates of Copenhagen and Paris airports
        let cph_cdg = Coor4D::raw(55., 12., 49., 2.);

        // A geodesic from Copenhagen to Paris
        let op = ctx.op("geodesic")?;
        let mut operands = [cph_cdg];
        ctx.apply(op, Inv, &mut operands)?;

        let expected = Coor4D([
            -130.1540604203936,
            -138.05257941840648,
            956066.2319619625,
            41.94742058159352,
        ]);

        assert!((operands[0][0] - expected[0]).abs() < 1e-9);
        assert!((operands[0][1] - expected[1]).abs() < 1e-9);
        assert!((operands[0][2] - expected[2]).abs() < 1e-9);
        assert!((operands[0][3] - expected[3]).abs() < 1e-9);

        // A geodesic from Copenhagen to Paris in the "reversible" format
        let op = ctx.op("geodesic reversible")?;
        let mut operands = [cph_cdg];
        ctx.apply(op, Inv, &mut operands)?;

        let expected = Coor4D([49.0, 2.0, 41.94742058159352, 956066.2319619625]);

        assert!((operands[0][0] - expected[0]).abs() < 1e-9);
        assert!((operands[0][1] - expected[1]).abs() < 1e-9);
        assert!((operands[0][2] - expected[2]).abs() < 1e-9);
        assert!((operands[0][3] - expected[3]).abs() < 1e-9);

        // And back to Copenhagen!
        ctx.apply(op, Fwd, &mut operands)?;

        assert!((operands[0][0] - cph_cdg[0]).abs() < 1e-10);
        assert!((operands[0][1] - cph_cdg[1]).abs() < 1e-10);
        assert!((operands[0][2] - cph_cdg[2]).abs() < 1e-10);
        assert!((operands[0][3] - cph_cdg[3]).abs() < 1e-10);

        Ok(())
    }
}