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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
//! # 🧠The NavLab examples
//!
//! FFI, the Norwegian Defence Research Establishment, publishes ten worked navigation problems
//! with the [n-vector](https://www.ffi.no/en/research/n-vector) method. They're a good tour of
//! what geo3d does, so here they are, solved with geo3d.
//!
//! Every example runs as a test. The answers are checked against the ones published by FFI and by
//! the [nvector](https://github.com/pbrod/nvector) Python package's own tests.
//!
//! ## 1. Where is B, seen from A?
//!
//! A is at 1° N, 2° E, 3 m below the ellipsoid; B is at 4° N, 5° E, 6 m below. Find B's offset north,
//! east and down from A, and the direction from A to B.
//!
//! ```
//! use geo3d::{Geodetic, LocalFrame};
//!
//! let a = LocalFrame::new(Geodetic::new(1.0, 2.0, -3.0));
//! let b = Geodetic::new(4.0, 5.0, -6.0);
//!
//! let ned = a.ned_of(b); // 331 730 m north, 332 998 m east, 17 404 m down
//! let look = a.look_angles(b); // 45.1° from north, 2.1° below the horizon, 470 357 m away
//! assert!((ned.north - 331_730.234_780_89).abs() < 1e-6);
//! assert!((ned.east - 332_997.874_989_27).abs() < 1e-6);
//! assert!((ned.down - 17_404.271_361_94).abs() < 1e-6);
//! assert!((look.azimuth_deg - 45.109_263_24).abs() < 1e-8);
//! ```
//!
//! B is "down" even though both are near the surface: the Earth curves away over 470 km.
//!
//! ## 2. Where is something seen from a vehicle?
//!
//! Vehicle B is 400 m up, at the position whose n-vector points along (1, 2, 3). It's heading 10°,
//! pitched 20° nose up and rolled 30° right. Something C is 3 000 m ahead of it, 2 000 m to its
//! right and 100 m below, in its own body frame. Where is C? FFI uses the WGS72 ellipsoid here.
//!
//! ```
//! use geo3d::{Attitude, Frd, Geodetic, LocalFrame, WGS72};
//!
//! // The n-vector (1, 2, 3) is latitude atan2(3, √5) and longitude atan2(2, 1).
//! let (lat, lon) = (3f64.atan2(5f64.sqrt()).to_degrees(), 2f64.atan2(1.0).to_degrees());
//! let b = LocalFrame::on(WGS72, Geodetic::new(lat, lon, 400.0));
//!
//! let c = b.geodetic_of_body(Attitude::new(10.0, 20.0, 30.0), Frd::new(3_000.0, 2_000.0, 100.0));
//! assert!((c.lat_deg - 53.326_378_26).abs() < 1e-8 && (c.lon_deg - 63.468_123_44).abs() < 1e-8);
//! assert!((c.height_m - 406.007_196_07).abs() < 1e-6);
//! ```
//!
//! ## 3. From x, y, z to latitude and longitude
//!
//! ```
//! use geo3d::{ecef_to_geodetic, Ecef};
//!
//! let b = ecef_to_geodetic(Ecef::new(6_371e3 * 0.9, -6_371e3, 6_371e3 * 1.1));
//! // 39.3787° N, 48.0128° W, 4 702 059.834 m up
//! assert!((b.lat_deg - 39.3787).abs() < 1e-4 && (b.lon_deg + 48.0128).abs() < 1e-4);
//! assert!((b.height_m - 4_702_059.834).abs() < 1e-3);
//! ```
//!
//! ## 4. From latitude and longitude to x, y, z
//!
//! ```
//! use geo3d::{geodetic_to_ecef, Geodetic};
//!
//! let p = geodetic_to_ecef(Geodetic::new(1.0, 2.0, 3.0));
//! assert!((p.x - 6_373_290.277_218_279).abs() < 1e-8);
//! assert!((p.y - 222_560.200_674_736_52).abs() < 1e-9);
//! assert!((p.z - 110_568.827_181_785_93).abs() < 1e-9);
//! ```
//!
//! These agree with the published answer to a nanometre. In fact geo3d's are slightly closer to
//! the exact values (worked out to 50 digits) in all three coordinates.
//!
//! ## 5. The distance between two points
//!
//! A is at 88° N, 0° E and B at 89° N, 170° W. FFI uses a sphere of radius 6 371 km.
//!
//! ```
//! use geo3d::{Ellipsoid, Geodetic, WGS84};
//!
//! let (a, b) = (Geodetic::new(88.0, 0.0, 0.0), Geodetic::new(89.0, -170.0, 0.0));
//! let sphere = Ellipsoid::new(6_371e3, 0.0);
//!
//! let along = sphere.great_circle_distance_m(a, b); // 332.456 km over the surface
//! let straight = sphere.to_ecef(a).distance_to(sphere.to_ecef(b)); // 332.419 km through it
//! assert!((along - 332_456.444_11).abs() < 1e-5 && (straight - 332_418.724_86).abs() < 1e-5);
//!
//! // On the real (WGS84) shape it's 333.95 km.
//! let real = WGS84.geodesic_inverse(a, b).unwrap();
//! assert!((real.distance_m - 333_950.0).abs() < 5.0);
//! ```
//!
//! ## 6. Where was it at 16 seconds?
//!
//! An object is at 89° N, 0° E at 10 s and at 89° N, 180° E at 20 s. Where was it at 16 s?
//!
//! ```
//! use geo3d::{Geodetic, WGS84};
//!
//! let (t0, t1) = (Geodetic::new(89.0, 0.0, 0.0), Geodetic::new(89.0, 180.0, 0.0));
//! let at_16 = WGS84.great_circle_interpolate(t0, t1, (16.0 - 10.0) / (20.0 - 10.0));
//! assert!((at_16.lat_deg - 89.8).abs() < 1e-12); // 0.2° past the pole
//! assert!((at_16.lon_deg.abs() - 180.0).abs() < 1e-9);
//! ```
//!
//! FFI's answer is 89.79998°: it blends the two n-vectors in a straight line, which is quick but
//! doesn't move at a steady speed. geo3d moves along the circle at a constant speed, which gives
//! exactly 89.8°.
//!
//! ## 7. The centre of three points
//!
//! ```
//! use geo3d::{Geodetic, WGS84};
//!
//! let points = [Geodetic::new(90.0, 0.0, 0.0), Geodetic::new(60.0, 10.0, 0.0), Geodetic::new(50.0, -20.0, 0.0)];
//! let m = WGS84.mean_position(&points).unwrap(); // 67.2362° N, 6.9175° W
//! let n = WGS84.surface_normal(m);
//! assert!((n.x - 0.384_117_170_292_6).abs() < 1e-12);
//! assert!((n.y + 0.046_602_405_485_689_447).abs() < 1e-12);
//! assert!((n.z - 0.922_107_485_757_139_5).abs() < 1e-12);
//! ```
//!
//! ## 8. Travel 1 km on a bearing
//!
//! From 80° N, 90° W, head 200° for 1 000 m, on a sphere of radius 6 371 km.
//!
//! ```
//! use geo3d::{Ellipsoid, Geodetic};
//!
//! let sphere = Ellipsoid::new(6_371e3, 0.0);
//! let b = sphere.great_circle_destination(Geodetic::new(80.0, -90.0, 0.0), 200.0, 1_000.0);
//! assert!((b.lat_deg - 79.991_548_67).abs() < 1e-8 && (b.lon_deg + 90.017_698_37).abs() < 1e-8);
//! ```
//!
//! On the real shape, use [`Ellipsoid::geodesic_direct`](crate::Ellipsoid::geodesic_direct).
//!
//! ## 9. Where do two paths cross?
//!
//! Path A runs from (10° N, 20° E) to (30° N, 40° E), and path B from (50° N, 60° E) to
//! (70° N, 80° E).
//!
//! ```
//! use geo3d::{Geodetic, WGS84};
//!
//! let g = |lat, lon| Geodetic::new(lat, lon, 0.0);
//! let c = WGS84.great_circle_intersection(g(10.0, 20.0), g(30.0, 40.0), g(50.0, 60.0), g(70.0, 80.0)).unwrap();
//! assert!((c.lat_deg - 40.318_643_07).abs() < 1e-8 && (c.lon_deg - 55.901_867_88).abs() < 1e-8);
//! ```
//!
//! ## 10. How far off the path?
//!
//! A path runs due north from (0°, 0°) to (10° N, 0°). How far to the side of it is (1° N, 0.1° E)?
//!
//! ```
//! use geo3d::{Ellipsoid, Geodetic};
//!
//! let sphere = Ellipsoid::new(6_371e3, 0.0);
//! let (start, end) = (Geodetic::new(0.0, 0.0, 0.0), Geodetic::new(10.0, 0.0, 0.0));
//! let off = sphere.cross_track_distance_m(start, end, Geodetic::new(1.0, 0.1, 0.0));
//! assert!((off - 11_117.799_110_15).abs() < 1e-6); // 11.1 km to the right
//! ```
//!
//! [`Ellipsoid::along_track_distance_m`](crate::Ellipsoid::along_track_distance_m) gives how far
//! along the path that point is.