geodtest.c
Go to the documentation of this file.
1 
14 #include "geodesic.h"
15 #include <stdio.h>
16 #include <math.h>
17 
18 #if defined(_MSC_VER)
19 /* Squelch warnings about assignment within conditional expression */
20 # pragma warning (disable: 4706)
21 #endif
22 
23 double wgs84_a = 6378137, wgs84_f = 1/298.257223563; /* WGS84 */
24 
25 static int assertEquals(double x, double y, double d) {
26  if (fabs(x - y) <= d)
27  return 0;
28  printf("assertEquals fails: %.7g != %.7g +/- %.7g\n", x, y, d);
29  return 1;
30 }
31 
32 const int ncases = 20;
33 double testcases[20][12] = {
34  {35.60777, -139.44815, 111.098748429560326,
35  -11.17491, -69.95921, 129.289270889708762,
36  8935244.5604818305, 80.50729714281974, 6273170.2055303837,
37  0.16606318447386067, 0.16479116945612937, 12841384694976.432},
38  {55.52454, 106.05087, 22.020059880982801,
39  77.03196, 197.18234, 109.112041110671519,
40  4105086.1713924406, 36.892740690445894, 3828869.3344387607,
41  0.80076349608092607, 0.80101006984201008, 61674961290615.615},
42  {-21.97856, 142.59065, -32.44456876433189,
43  41.84138, 98.56635, -41.84359951440466,
44  8394328.894657671, 75.62930491011522, 6161154.5773110616,
45  0.24816339233950381, 0.24930251203627892, -6637997720646.717},
46  {-66.99028, 112.2363, 173.73491240878403,
47  -12.70631, 285.90344, 2.512956620913668,
48  11150344.2312080241, 100.278634181155759, 6289939.5670446687,
49  -0.17199490274700385, -0.17722569526345708, -121287239862139.744},
50  {-17.42761, 173.34268, -159.033557661192928,
51  -15.84784, 5.93557, -20.787484651536988,
52  16076603.1631180673, 144.640108810286253, 3732902.1583877189,
53  -0.81273638700070476, -0.81299800519154474, 97825992354058.708},
54  {32.84994, 48.28919, 150.492927788121982,
55  -56.28556, 202.29132, 48.113449399816759,
56  16727068.9438164461, 150.565799985466607, 3147838.1910180939,
57  -0.87334918086923126, -0.86505036767110637, -72445258525585.010},
58  {6.96833, 52.74123, 92.581585386317712,
59  -7.39675, 206.17291, 90.721692165923907,
60  17102477.2496958388, 154.147366239113561, 2772035.6169917581,
61  -0.89991282520302447, -0.89986892177110739, -1311796973197.995},
62  {-50.56724, -16.30485, -105.439679907590164,
63  -33.56571, -94.97412, -47.348547835650331,
64  6455670.5118668696, 58.083719495371259, 5409150.7979815838,
65  0.53053508035997263, 0.52988722644436602, 41071447902810.047},
66  {-58.93002, -8.90775, 140.965397902500679,
67  -8.91104, 133.13503, 19.255429433416599,
68  11756066.0219864627, 105.755691241406877, 6151101.2270708536,
69  -0.26548622269867183, -0.27068483874510741, -86143460552774.735},
70  {-68.82867, -74.28391, 93.774347763114881,
71  -50.63005, -8.36685, 34.65564085411343,
72  3956936.926063544, 35.572254987389284, 3708890.9544062657,
73  0.81443963736383502, 0.81420859815358342, -41845309450093.787},
74  {-10.62672, -32.0898, -86.426713286747751,
75  5.883, -134.31681, -80.473780971034875,
76  11470869.3864563009, 103.387395634504061, 6184411.6622659713,
77  -0.23138683500430237, -0.23155097622286792, 4198803992123.548},
78  {-21.76221, 166.90563, 29.319421206936428,
79  48.72884, 213.97627, 43.508671946410168,
80  9098627.3986554915, 81.963476716121964, 6299240.9166992283,
81  0.13965943368590333, 0.14152969707656796, 10024709850277.476},
82  {-19.79938, -174.47484, 71.167275780171533,
83  -11.99349, -154.35109, 65.589099775199228,
84  2319004.8601169389, 20.896611684802389, 2267960.8703918325,
85  0.93427001867125849, 0.93424887135032789, -3935477535005.785},
86  {-11.95887, -116.94513, 92.712619830452549,
87  4.57352, 7.16501, 78.64960934409585,
88  13834722.5801401374, 124.688684161089762, 5228093.177931598,
89  -0.56879356755666463, -0.56918731952397221, -9919582785894.853},
90  {-87.85331, 85.66836, -65.120313040242748,
91  66.48646, 16.09921, -4.888658719272296,
92  17286615.3147144645, 155.58592449699137, 2635887.4729110181,
93  -0.90697975771398578, -0.91095608883042767, 42667211366919.534},
94  {1.74708, 128.32011, -101.584843631173858,
95  -11.16617, 11.87109, -86.325793296437476,
96  12942901.1241347408, 116.650512484301857, 5682744.8413270572,
97  -0.44857868222697644, -0.44824490340007729, 10763055294345.653},
98  {-25.72959, -144.90758, -153.647468693117198,
99  -57.70581, -269.17879, -48.343983158876487,
100  9413446.7452453107, 84.664533838404295, 6356176.6898881281,
101  0.09492245755254703, 0.09737058264766572, 74515122850712.444},
102  {-41.22777, 122.32875, 14.285113402275739,
103  -7.57291, 130.37946, 10.805303085187369,
104  3812686.035106021, 34.34330804743883, 3588703.8812128856,
105  0.82605222593217889, 0.82572158200920196, -2456961531057.857},
106  {11.01307, 138.25278, 79.43682622782374,
107  6.62726, 247.05981, 103.708090215522657,
108  11911190.819018408, 107.341669954114577, 6070904.722786735,
109  -0.29767608923657404, -0.29785143390252321, 17121631423099.696},
110  {-29.47124, 95.14681, -163.779130441688382,
111  -27.46601, -69.15955, -15.909335945554969,
112  13487015.8381145492, 121.294026715742277, 5481428.9945736388,
113  -0.51527225545373252, -0.51556587964721788, 104679964020340.318}};
114 
115 static int testinverse() {
116  double lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12;
117  double azi1a, azi2a, s12a, a12a, m12a, M12a, M21a, S12a;
118  struct geod_geodesic g;
119  int i, result = 0;
120  geod_init(&g, wgs84_a, wgs84_f);
121  for (i = 0; i < ncases; ++i) {
122  lat1 = testcases[i][0]; lon1 = testcases[i][1]; azi1 = testcases[i][2];
123  lat2 = testcases[i][3]; lon2 = testcases[i][4]; azi2 = testcases[i][5];
124  s12 = testcases[i][6]; a12 = testcases[i][7]; m12 = testcases[i][8];
125  M12 = testcases[i][9]; M21 = testcases[i][10]; S12 = testcases[i][11];
126  a12a = geod_geninverse(&g, lat1, lon1, lat2, lon2, &s12a, &azi1a, &azi2a,
127  &m12a, &M12a, &M21a, &S12a);
128  result += assertEquals(azi1, azi1a, 1e-13);
129  result += assertEquals(azi2, azi2a, 1e-13);
130  result += assertEquals(s12, s12a, 1e-8);
131  result += assertEquals(a12, a12a, 1e-13);
132  result += assertEquals(m12, m12a, 1e-8);
133  result += assertEquals(M12, M12a, 1e-15);
134  result += assertEquals(M21, M21a, 1e-15);
135  result += assertEquals(S12, S12a, 0.1);
136  }
137  return result;
138 }
139 
140 static int testdirect() {
141  double lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12;
142  double lat2a, lon2a, azi2a, a12a, m12a, M12a, M21a, S12a;
143  struct geod_geodesic g;
144  int i, result = 0;
145  unsigned flags = GEOD_LONG_UNROLL;
146  geod_init(&g, wgs84_a, wgs84_f);
147  for (i = 0; i < ncases; ++i) {
148  lat1 = testcases[i][0]; lon1 = testcases[i][1]; azi1 = testcases[i][2];
149  lat2 = testcases[i][3]; lon2 = testcases[i][4]; azi2 = testcases[i][5];
150  s12 = testcases[i][6]; a12 = testcases[i][7]; m12 = testcases[i][8];
151  M12 = testcases[i][9]; M21 = testcases[i][10]; S12 = testcases[i][11];
152  a12a = geod_gendirect(&g, lat1, lon1, azi1, flags, s12,
153  &lat2a, &lon2a, &azi2a, 0,
154  &m12a, &M12a, &M21a, &S12a);
155  result += assertEquals(lat2, lat2a, 1e-13);
156  result += assertEquals(lon2, lon2a, 1e-13);
157  result += assertEquals(azi2, azi2a, 1e-13);
158  result += assertEquals(a12, a12a, 1e-13);
159  result += assertEquals(m12, m12a, 1e-8);
160  result += assertEquals(M12, M12a, 1e-15);
161  result += assertEquals(M21, M21a, 1e-15);
162  result += assertEquals(S12, S12a, 0.1);
163  }
164  return result;
165 }
166 
167 static int testarcdirect() {
168  double lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12;
169  double lat2a, lon2a, azi2a, s12a, m12a, M12a, M21a, S12a;
170  struct geod_geodesic g;
171  int i, result = 0;
172  unsigned flags = GEOD_ARCMODE | GEOD_LONG_UNROLL;
173  geod_init(&g, wgs84_a, wgs84_f);
174  for (i = 0; i < ncases; ++i) {
175  lat1 = testcases[i][0]; lon1 = testcases[i][1]; azi1 = testcases[i][2];
176  lat2 = testcases[i][3]; lon2 = testcases[i][4]; azi2 = testcases[i][5];
177  s12 = testcases[i][6]; a12 = testcases[i][7]; m12 = testcases[i][8];
178  M12 = testcases[i][9]; M21 = testcases[i][10]; S12 = testcases[i][11];
179  geod_gendirect(&g, lat1, lon1, azi1, flags, a12,
180  &lat2a, &lon2a, &azi2a, &s12a, &m12a, &M12a, &M21a, &S12a);
181  result += assertEquals(lat2, lat2a, 1e-13);
182  result += assertEquals(lon2, lon2a, 1e-13);
183  result += assertEquals(azi2, azi2a, 1e-13);
184  result += assertEquals(s12, s12a, 1e-8);
185  result += assertEquals(m12, m12a, 1e-8);
186  result += assertEquals(M12, M12a, 1e-15);
187  result += assertEquals(M21, M21a, 1e-15);
188  result += assertEquals(S12, S12a, 0.1);
189  }
190  return result;
191 }
192 
193 static int GeodSolve0() {
194  double azi1, azi2, s12;
195  struct geod_geodesic g;
196  int result = 0;
197  geod_init(&g, wgs84_a, wgs84_f);
198  geod_inverse(&g, 40.6, -73.8, 49.01666667, 2.55, &s12, &azi1, &azi2);
199  result += assertEquals(azi1, 53.47022, 0.5e-5);
200  result += assertEquals(azi2, 111.59367, 0.5e-5);
201  result += assertEquals(s12, 5853226, 0.5);
202  return result;
203 }
204 
205 static int GeodSolve1() {
206  double lat2, lon2, azi2;
207  struct geod_geodesic g;
208  int result = 0;
209  geod_init(&g, wgs84_a, wgs84_f);
210  geod_direct(&g, 40.63972222, -73.77888889, 53.5, 5850e3,
211  &lat2, &lon2, &azi2);
212  result += assertEquals(lat2, 49.01467, 0.5e-5);
213  result += assertEquals(lon2, 2.56106, 0.5e-5);
214  result += assertEquals(azi2, 111.62947, 0.5e-5);
215  return result;
216 }
217 
218 static int GeodSolve2() {
219  /* Check fix for antipodal prolate bug found 2010-09-04 */
220  double azi1, azi2, s12;
221  struct geod_geodesic g;
222  int result = 0;
223  geod_init(&g, 6.4e6, -1/150.0);
224  geod_inverse(&g, 0.07476, 0, -0.07476, 180, &s12, &azi1, &azi2);
225  result += assertEquals(azi1, 90.00078, 0.5e-5);
226  result += assertEquals(azi2, 90.00078, 0.5e-5);
227  result += assertEquals(s12, 20106193, 0.5);
228  geod_inverse(&g, 0.1, 0, -0.1, 180, &s12, &azi1, &azi2);
229  result += assertEquals(azi1, 90.00105, 0.5e-5);
230  result += assertEquals(azi2, 90.00105, 0.5e-5);
231  result += assertEquals(s12, 20106193, 0.5);
232  return result;
233 }
234 
235 static int GeodSolve4() {
236  /* Check fix for short line bug found 2010-05-21 */
237  double s12;
238  struct geod_geodesic g;
239  int result = 0;
240  geod_init(&g, wgs84_a, wgs84_f);
241  geod_inverse(&g, 36.493349428792, 0, 36.49334942879201, .0000008,
242  &s12, 0, 0);
243  result += assertEquals(s12, 0.072, 0.5e-3);
244  return result;
245 }
246 
247 static int GeodSolve5() {
248  /* Check fix for point2=pole bug found 2010-05-03 */
249  double lat2, lon2, azi2;
250  struct geod_geodesic g;
251  int result = 0;
252  geod_init(&g, wgs84_a, wgs84_f);
253  geod_direct(&g, 0.01777745589997, 30, 0, 10e6, &lat2, &lon2, &azi2);
254  result += assertEquals(lat2, 90, 0.5e-5);
255  if (lon2 < 0) {
256  result += assertEquals(lon2, -150, 0.5e-5);
257  result += assertEquals(fabs(azi2), 180, 0.5e-5);
258  } else {
259  result += assertEquals(lon2, 30, 0.5e-5);
260  result += assertEquals(azi2, 0, 0.5e-5);
261  }
262  return result;
263 }
264 
265 static int GeodSolve6() {
266  /* Check fix for volatile sbet12a bug found 2011-06-25 (gcc 4.4.4
267  * x86 -O3). Found again on 2012-03-27 with tdm-mingw32 (g++ 4.6.1). */
268  double s12;
269  struct geod_geodesic g;
270  int result = 0;
271  geod_init(&g, wgs84_a, wgs84_f);
272  geod_inverse(&g, 88.202499451857, 0,
273  -88.202499451857, 179.981022032992859592, &s12, 0, 0);
274  result += assertEquals(s12, 20003898.214, 0.5e-3);
275  geod_inverse(&g, 89.262080389218, 0,
276  -89.262080389218, 179.992207982775375662, &s12, 0, 0);
277  result += assertEquals(s12, 20003925.854, 0.5e-3);
278  geod_inverse(&g, 89.333123580033, 0,
279  -89.333123580032997687, 179.99295812360148422, &s12, 0, 0);
280  result += assertEquals(s12, 20003926.881, 0.5e-3);
281  return result;
282 }
283 
284 static int GeodSolve9() {
285  /* Check fix for volatile x bug found 2011-06-25 (gcc 4.4.4 x86 -O3) */
286  double s12;
287  struct geod_geodesic g;
288  int result = 0;
289  geod_init(&g, wgs84_a, wgs84_f);
290  geod_inverse(&g, 56.320923501171, 0,
291  -56.320923501171, 179.664747671772880215, &s12, 0, 0);
292  result += assertEquals(s12, 19993558.287, 0.5e-3);
293  return result;
294 }
295 
296 static int GeodSolve10() {
297  /* Check fix for adjust tol1_ bug found 2011-06-25 (Visual Studio
298  * 10 rel + debug) */
299  double s12;
300  struct geod_geodesic g;
301  int result = 0;
302  geod_init(&g, wgs84_a, wgs84_f);
303  geod_inverse(&g, 52.784459512564, 0,
304  -52.784459512563990912, 179.634407464943777557, &s12, 0, 0);
305  result += assertEquals(s12, 19991596.095, 0.5e-3);
306  return result;
307 }
308 
309 static int GeodSolve11() {
310  /* Check fix for bet2 = -bet1 bug found 2011-06-25 (Visual Studio
311  * 10 rel + debug) */
312  double s12;
313  struct geod_geodesic g;
314  int result = 0;
315  geod_init(&g, wgs84_a, wgs84_f);
316  geod_inverse(&g, 48.522876735459, 0,
317  -48.52287673545898293, 179.599720456223079643, &s12, 0, 0);
318  result += assertEquals(s12, 19989144.774, 0.5e-3);
319  return result;
320 }
321 
322 static int GeodSolve12() {
323  /* Check fix for inverse geodesics on extreme prolate/oblate
324  * ellipsoids Reported 2012-08-29 Stefan Guenther
325  * <stefan.gunther@embl.de>; fixed 2012-10-07 */
326  double azi1, azi2, s12;
327  struct geod_geodesic g;
328  int result = 0;
329  geod_init(&g, 89.8, -1.83);
330  geod_inverse(&g, 0, 0, -10, 160, &s12, &azi1, &azi2);
331  result += assertEquals(azi1, 120.27, 1e-2);
332  result += assertEquals(azi2, 105.15, 1e-2);
333  result += assertEquals(s12, 266.7, 1e-1);
334  return result;
335 }
336 
337 static int GeodSolve14() {
338  /* Check fix for inverse ignoring lon12 = nan */
339  double azi1, azi2, s12, nan;
340  struct geod_geodesic g;
341  int result = 0;
342  {
343  double minus1 = -1;
344  nan = sqrt(minus1);
345  }
346  geod_init(&g, wgs84_a, wgs84_f);
347  geod_inverse(&g, 0, 0, 1, nan, &s12, &azi1, &azi2);
348  result += azi1 == azi1 ? 1 : 0;
349  result += azi2 == azi2 ? 1 : 0;
350  result += s12 == s12 ? 1 : 0;
351  return result;
352 }
353 
354 static int GeodSolve15() {
355  /* Initial implementation of Math::eatanhe was wrong for e^2 < 0. This
356  * checks that this is fixed. */
357  double S12;
358  struct geod_geodesic g;
359  int result = 0;
360  geod_init(&g, 6.4e6, -1/150.0);
361  geod_gendirect(&g, 1, 2, 3, 0, 4,
362  0, 0, 0, 0, 0, 0, 0, &S12);
363  result += assertEquals(S12, 23700, 0.5);
364  return result;
365 }
366 
367 static int GeodSolve17() {
368  /* Check fix for LONG_UNROLL bug found on 2015-05-07 */
369  double lat2, lon2, azi2;
370  struct geod_geodesic g;
371  struct geod_geodesicline l;
372  int result = 0;
373  unsigned flags = GEOD_LONG_UNROLL;
374  geod_init(&g, wgs84_a, wgs84_f);
375  geod_gendirect(&g, 40, -75, -10, flags, 2e7,
376  &lat2, &lon2, &azi2, 0, 0, 0, 0, 0);
377  result += assertEquals(lat2, -39, 1);
378  result += assertEquals(lon2, -254, 1);
379  result += assertEquals(azi2, -170, 1);
380  geod_lineinit(&l, &g, 40, -75, -10, 0);
381  geod_genposition(&l, flags, 2e7, &lat2, &lon2, &azi2, 0, 0, 0, 0, 0);
382  result += assertEquals(lat2, -39, 1);
383  result += assertEquals(lon2, -254, 1);
384  result += assertEquals(azi2, -170, 1);
385  geod_direct(&g, 40, -75, -10, 2e7, &lat2, &lon2, &azi2);
386  result += assertEquals(lat2, -39, 1);
387  result += assertEquals(lon2, 105, 1);
388  result += assertEquals(azi2, -170, 1);
389  geod_position(&l, 2e7, &lat2, &lon2, &azi2);
390  result += assertEquals(lat2, -39, 1);
391  result += assertEquals(lon2, 105, 1);
392  result += assertEquals(azi2, -170, 1);
393  return result;
394 }
395 
396 static int GeodSolve26() {
397  /* Check 0/0 problem with area calculation on sphere 2015-09-08 */
398  double S12;
399  struct geod_geodesic g;
400  int result = 0;
401  geod_init(&g, 6.4e6, 0);
402  geod_geninverse(&g, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, &S12);
403  result += assertEquals(S12, 49911046115.0, 0.5);
404  return result;
405 }
406 
407 static int GeodSolve28() {
408  /* Check for bad placement of assignment of r.a12 with |f| > 0.01 (bug in
409  * Java implementation fixed on 2015-05-19). */
410  double a12;
411  struct geod_geodesic g;
412  int result = 0;
413  geod_init(&g, 6.4e6, 0.1);
414  a12 = geod_gendirect(&g, 1, 2, 10, 0, 5e6, 0, 0, 0, 0, 0, 0, 0, 0);
415  result += assertEquals(a12, 48.55570690, 0.5e-8);
416  return result;
417 }
418 
419 static int GeodSolve33() {
420  /* Check max(-0.0,+0.0) issues 2015-08-22 (triggered by bugs in Octave --
421  * sind(-0.0) = +0.0 -- and in some version of Visual Studio --
422  * fmod(-0.0, 360.0) = +0.0. */
423  double azi1, azi2, s12;
424  struct geod_geodesic g;
425  int result = 0;
426  geod_init(&g, wgs84_a, wgs84_f);
427  geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2);
428  result += assertEquals(azi1, 90.00000, 0.5e-5);
429  result += assertEquals(azi2, 90.00000, 0.5e-5);
430  result += assertEquals(s12, 19926189, 0.5);
431  geod_inverse(&g, 0, 0, 0, 179.5, &s12, &azi1, &azi2);
432  result += assertEquals(azi1, 55.96650, 0.5e-5);
433  result += assertEquals(azi2, 124.03350, 0.5e-5);
434  result += assertEquals(s12, 19980862, 0.5);
435  geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2);
436  result += assertEquals(azi1, 0.00000, 0.5e-5);
437  result += assertEquals(fabs(azi2), 180.00000, 0.5e-5);
438  result += assertEquals(s12, 20003931, 0.5);
439  geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2);
440  result += assertEquals(azi1, 0.00000, 0.5e-5);
441  result += assertEquals(fabs(azi2), 180.00000, 0.5e-5);
442  result += assertEquals(s12, 19893357, 0.5);
443  geod_init(&g, 6.4e6, 0);
444  geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2);
445  result += assertEquals(azi1, 90.00000, 0.5e-5);
446  result += assertEquals(azi2, 90.00000, 0.5e-5);
447  result += assertEquals(s12, 19994492, 0.5);
448  geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2);
449  result += assertEquals(azi1, 0.00000, 0.5e-5);
450  result += assertEquals(fabs(azi2), 180.00000, 0.5e-5);
451  result += assertEquals(s12, 20106193, 0.5);
452  geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2);
453  result += assertEquals(azi1, 0.00000, 0.5e-5);
454  result += assertEquals(fabs(azi2), 180.00000, 0.5e-5);
455  result += assertEquals(s12, 19994492, 0.5);
456  geod_init(&g, 6.4e6, -1/300.0);
457  geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2);
458  result += assertEquals(azi1, 90.00000, 0.5e-5);
459  result += assertEquals(azi2, 90.00000, 0.5e-5);
460  result += assertEquals(s12, 19994492, 0.5);
461  geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2);
462  result += assertEquals(azi1, 90.00000, 0.5e-5);
463  result += assertEquals(azi2, 90.00000, 0.5e-5);
464  result += assertEquals(s12, 20106193, 0.5);
465  geod_inverse(&g, 0, 0, 0.5, 180, &s12, &azi1, &azi2);
466  result += assertEquals(azi1, 33.02493, 0.5e-5);
467  result += assertEquals(azi2, 146.97364, 0.5e-5);
468  result += assertEquals(s12, 20082617, 0.5);
469  geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2);
470  result += assertEquals(azi1, 0.00000, 0.5e-5);
471  result += assertEquals(fabs(azi2), 180.00000, 0.5e-5);
472  result += assertEquals(s12, 20027270, 0.5);
473 
474  return result;
475 }
476 
477 static int GeodSolve55() {
478  /* Check fix for nan + point on equator or pole not returning all nans in
479  * Geodesic::Inverse, found 2015-09-23. */
480  double azi1, azi2, s12, nan;
481  struct geod_geodesic g;
482  int result = 0;
483  {
484  double minus1 = -1;
485  nan = sqrt(minus1);
486  }
487  geod_init(&g, wgs84_a, wgs84_f);
488  geod_inverse(&g, nan, 0, 0, 90, &s12, &azi1, &azi2);
489  result += azi1 == azi1 ? 1 : 0;
490  result += azi2 == azi2 ? 1 : 0;
491  result += s12 == s12 ? 1 : 0;
492  geod_inverse(&g, nan, 0, 90, 9, &s12, &azi1, &azi2);
493  result += azi1 == azi1 ? 1 : 0;
494  result += azi2 == azi2 ? 1 : 0;
495  result += s12 == s12 ? 1 : 0;
496  return result;
497 }
498 
499 static int GeodSolve59() {
500  /* Check for points close with longitudes close to 180 deg apart. */
501  double azi1, azi2, s12;
502  struct geod_geodesic g;
503  int result = 0;
504  geod_init(&g, wgs84_a, wgs84_f);
505  geod_inverse(&g, 5, 0.00000000000001, 10, 180, &s12, &azi1, &azi2);
506  result += assertEquals(azi1, 0.000000000000035, 1.5e-14);
507  result += assertEquals(azi2, 179.99999999999996, 1.5e-14);
508  result += assertEquals(s12, 18345191.174332713, 2.5e-9);
509  return result;
510 }
511 
512 static int GeodSolve61() {
513  /* Make sure small negative azimuths are west-going */
514  double lat2, lon2, azi2;
515  struct geod_geodesic g;
516  struct geod_geodesicline l;
517  int result = 0;
518  unsigned flags = GEOD_LONG_UNROLL;
519  geod_init(&g, wgs84_a, wgs84_f);
520  geod_gendirect(&g, 45, 0, -0.000000000000000003, flags, 1e7,
521  &lat2, &lon2, &azi2, 0, 0, 0, 0, 0);
522  result += assertEquals(lat2, 45.30632, 0.5e-5);
523  result += assertEquals(lon2, -180, 0.5e-5);
524  result += assertEquals(fabs(azi2), 180, 0.5e-5);
525  geod_inverseline(&l, &g, 45, 0, 80, -0.000000000000000003, 0);
526  geod_genposition(&l, flags, 1e7, &lat2, &lon2, &azi2, 0, 0, 0, 0, 0);
527  result += assertEquals(lat2, 45.30632, 0.5e-5);
528  result += assertEquals(lon2, -180, 0.5e-5);
529  result += assertEquals(fabs(azi2), 180, 0.5e-5);
530  return result;
531 }
532 
533 static int GeodSolve65() {
534  /* Check for bug in east-going check in GeodesicLine (needed to check for
535  * sign of 0) and sign error in area calculation due to a bogus override of
536  * the code for alp12. Found/fixed on 2015-12-19. */
537  double lat2, lon2, azi2, s12, a12, m12, M12, M21, S12;
538  struct geod_geodesic g;
539  struct geod_geodesicline l;
540  int result = 0;
541  unsigned flags = GEOD_LONG_UNROLL, caps = GEOD_ALL;
542  geod_init(&g, wgs84_a, wgs84_f);
543  geod_inverseline(&l, &g, 30, -0.000000000000000001, -31, 180, caps);
544  a12 = geod_genposition(&l, flags, 1e7,
545  &lat2, &lon2, &azi2, &s12, &m12, &M12, &M21, &S12);
546  result += assertEquals(lat2, -60.23169, 0.5e-5);
547  result += assertEquals(lon2, -0.00000, 0.5e-5);
548  result += assertEquals(fabs(azi2), 180.00000, 0.5e-5);
549  result += assertEquals(s12, 10000000, 0.5);
550  result += assertEquals(a12, 90.06544, 0.5e-5);
551  result += assertEquals(m12, 6363636, 0.5);
552  result += assertEquals(M12, -0.0012834, 0.5e-7);
553  result += assertEquals(M21, 0.0013749, 0.5e-7);
554  result += assertEquals(S12, 0, 0.5);
555  a12 = geod_genposition(&l, flags, 2e7,
556  &lat2, &lon2, &azi2, &s12, &m12, &M12, &M21, &S12);
557  result += assertEquals(lat2, -30.03547, 0.5e-5);
558  result += assertEquals(lon2, -180.00000, 0.5e-5);
559  result += assertEquals(azi2, -0.00000, 0.5e-5);
560  result += assertEquals(s12, 20000000, 0.5);
561  result += assertEquals(a12, 179.96459, 0.5e-5);
562  result += assertEquals(m12, 54342, 0.5);
563  result += assertEquals(M12, -1.0045592, 0.5e-7);
564  result += assertEquals(M21, -0.9954339, 0.5e-7);
565  result += assertEquals(S12, 127516405431022.0, 0.5);
566  return result;
567 }
568 
569 static int GeodSolve67() {
570  /* Check for InverseLine if line is slightly west of S and that s13 is
571  correctly set. */
572  double lat2, lon2, azi2;
573  struct geod_geodesic g;
574  struct geod_geodesicline l;
575  int result = 0;
576  unsigned flags = GEOD_LONG_UNROLL;
577  geod_init(&g, wgs84_a, wgs84_f);
578  geod_inverseline(&l, &g, -5, -0.000000000000002, -10, 180, 0);
579  geod_genposition(&l, flags, 2e7, &lat2, &lon2, &azi2, 0, 0, 0, 0, 0);
580  result += assertEquals(lat2, 4.96445, 0.5e-5);
581  result += assertEquals(lon2, -180.00000, 0.5e-5);
582  result += assertEquals(azi2, -0.00000, 0.5e-5);
583  geod_genposition(&l, flags, 0.5 * l.s13, &lat2, &lon2, &azi2, 0, 0, 0, 0, 0);
584  result += assertEquals(lat2, -87.52461, 0.5e-5);
585  result += assertEquals(lon2, -0.00000, 0.5e-5);
586  result += assertEquals(azi2, -180.00000, 0.5e-5);
587  return result;
588 }
589 
590 static int GeodSolve71() {
591  /* Check that DirectLine sets s13. */
592  double lat2, lon2, azi2;
593  struct geod_geodesic g;
594  struct geod_geodesicline l;
595  int result = 0;
596  geod_init(&g, wgs84_a, wgs84_f);
597  geod_directline(&l, &g, 1, 2, 45, 1e7, 0);
598  geod_position(&l, 0.5 * l.s13, &lat2, &lon2, &azi2);
599  result += assertEquals(lat2, 30.92625, 0.5e-5);
600  result += assertEquals(lon2, 37.54640, 0.5e-5);
601  result += assertEquals(azi2, 55.43104, 0.5e-5);
602  return result;
603 }
604 
605 static int GeodSolve73() {
606  /* Check for backwards from the pole bug reported by Anon on 2016-02-13.
607  * This only affected the Java implementation. It was introduced in Java
608  * version 1.44 and fixed in 1.46-SNAPSHOT on 2016-01-17. */
609  double lat2, lon2, azi2;
610  struct geod_geodesic g;
611  int result = 0;
612  geod_init(&g, wgs84_a, wgs84_f);
613  geod_direct(&g, 90, 10, 180, -1e6,
614  &lat2, &lon2, &azi2);
615  result += assertEquals(lat2, 81.04623, 0.5e-5);
616  result += assertEquals(lon2, -170, 0.5e-5);
617  result += assertEquals(azi2, 0, 0.5e-5);
618  return result;
619 }
620 
621 static void planimeter(const struct geod_geodesic* g,
622  double points[][2], int N,
623  double* perimeter, double* area) {
624  struct geod_polygon p;
625  int i;
626  geod_polygon_init(&p, 0);
627  for (i = 0; i < N; ++i)
628  geod_polygon_addpoint(g, &p, points[i][0], points[i][1]);
629  geod_polygon_compute(g, &p, 0, 1, area, perimeter);
630 }
631 
632 static void polylength(const struct geod_geodesic* g,
633  double points[][2], int N,
634  double* perimeter) {
635  struct geod_polygon p;
636  int i;
637  geod_polygon_init(&p, 1);
638  for (i = 0; i < N; ++i)
639  geod_polygon_addpoint(g, &p, points[i][0], points[i][1]);
640  geod_polygon_compute(g, &p, 0, 1, 0, perimeter);
641 }
642 
643 static int GeodSolve74() {
644  /* Check fix for inaccurate areas, bug introduced in v1.46, fixed
645  2015-10-16. */
646  double a12, s12, azi1, azi2, m12, M12, M21, S12;
647  struct geod_geodesic g;
648  int result = 0;
649  geod_init(&g, wgs84_a, wgs84_f);
650  a12 = geod_geninverse(&g, 54.1589, 15.3872, 54.1591, 15.3877,
651  &s12, &azi1, &azi2, &m12, &M12, &M21, &S12);
652  result += assertEquals(azi1, 55.723110355, 5e-9);
653  result += assertEquals(azi2, 55.723515675, 5e-9);
654  result += assertEquals(s12, 39.527686385, 5e-9);
655  result += assertEquals(a12, 0.000355495, 5e-9);
656  result += assertEquals(m12, 39.527686385, 5e-9);
657  result += assertEquals(M12, 0.999999995, 5e-9);
658  result += assertEquals(M21, 0.999999995, 5e-9);
659  result += assertEquals(S12, 286698586.30197, 5e-4);
660  return result;
661 }
662 
663 static int GeodSolve76() {
664  /* The distance from Wellington and Salamanca (a classic failure of
665  Vincenty) */
666  double azi1, azi2, s12;
667  struct geod_geodesic g;
668  int result = 0;
669  geod_init(&g, wgs84_a, wgs84_f);
670  geod_inverse(&g, -(41+19/60.0), 174+49/60.0, 40+58/60.0, -(5+30/60.0),
671  &s12, &azi1, &azi2);
672  result += assertEquals(azi1, 160.39137649664, 0.5e-11);
673  result += assertEquals(azi2, 19.50042925176, 0.5e-11);
674  result += assertEquals(s12, 19960543.857179, 0.5e-6);
675  return result;
676 }
677 
678 static int GeodSolve78() {
679  /* An example where the NGS calculator fails to converge */
680  double azi1, azi2, s12;
681  struct geod_geodesic g;
682  int result = 0;
683  geod_init(&g, wgs84_a, wgs84_f);
684  geod_inverse(&g, 27.2, 0.0, -27.1, 179.5, &s12, &azi1, &azi2);
685  result += assertEquals(azi1, 45.82468716758, 0.5e-11);
686  result += assertEquals(azi2, 134.22776532670, 0.5e-11);
687  result += assertEquals(s12, 19974354.765767, 0.5e-6);
688  return result;
689 }
690 
691 static int Planimeter0() {
692  /* Check fix for pole-encircling bug found 2011-03-16 */
693  double pa[4][2] = {{89, 0}, {89, 90}, {89, 180}, {89, 270}};
694  double pb[4][2] = {{-89, 0}, {-89, 90}, {-89, 180}, {-89, 270}};
695  double pc[4][2] = {{0, -1}, {-1, 0}, {0, 1}, {1, 0}};
696  double pd[3][2] = {{90, 0}, {0, 0}, {0, 90}};
697  struct geod_geodesic g;
698  double perimeter, area;
699  int result = 0;
700  geod_init(&g, wgs84_a, wgs84_f);
701 
702  planimeter(&g, pa, 4, &perimeter, &area);
703  result += assertEquals(perimeter, 631819.8745, 1e-4);
704  result += assertEquals(area, 24952305678.0, 1);
705 
706  planimeter(&g, pb, 4, &perimeter, &area);
707  result += assertEquals(perimeter, 631819.8745, 1e-4);
708  result += assertEquals(area, -24952305678.0, 1);
709 
710  planimeter(&g, pc, 4, &perimeter, &area);
711  result += assertEquals(perimeter, 627598.2731, 1e-4);
712  result += assertEquals(area, 24619419146.0, 1);
713 
714  planimeter(&g, pd, 3, &perimeter, &area);
715  result += assertEquals(perimeter, 30022685, 1);
716  result += assertEquals(area, 63758202715511.0, 1);
717 
718  polylength(&g, pd, 3, &perimeter);
719  result += assertEquals(perimeter, 20020719, 1);
720 
721  return result;
722 }
723 
724 static int Planimeter5() {
725  /* Check fix for Planimeter pole crossing bug found 2011-06-24 */
726  double points[3][2] = {{89, 0.1}, {89, 90.1}, {89, -179.9}};
727  struct geod_geodesic g;
728  double perimeter, area;
729  int result = 0;
730  geod_init(&g, wgs84_a, wgs84_f);
731  planimeter(&g, points, 3, &perimeter, &area);
732  result += assertEquals(perimeter, 539297, 1);
733  result += assertEquals(area, 12476152838.5, 1);
734  return result;
735 }
736 
737 static int Planimeter6() {
738  /* Check fix for Planimeter lon12 rounding bug found 2012-12-03 */
739  double pa[3][2] = {{9, -0.00000000000001}, {9, 180}, {9, 0}};
740  double pb[3][2] = {{9, 0.00000000000001}, {9, 0}, {9, 180}};
741  double pc[3][2] = {{9, 0.00000000000001}, {9, 180}, {9, 0}};
742  double pd[3][2] = {{9, -0.00000000000001}, {9, 0}, {9, 180}};
743  struct geod_geodesic g;
744  double perimeter, area;
745  int result = 0;
746  geod_init(&g, wgs84_a, wgs84_f);
747 
748  planimeter(&g, pa, 3, &perimeter, &area);
749  result += assertEquals(perimeter, 36026861, 1);
750  result += assertEquals(area, 0, 1);
751  planimeter(&g, pb, 3, &perimeter, &area);
752  result += assertEquals(perimeter, 36026861, 1);
753  result += assertEquals(area, 0, 1);
754  planimeter(&g, pc, 3, &perimeter, &area);
755  result += assertEquals(perimeter, 36026861, 1);
756  result += assertEquals(area, 0, 1);
757  planimeter(&g, pd, 3, &perimeter, &area);
758  result += assertEquals(perimeter, 36026861, 1);
759  result += assertEquals(area, 0, 1);
760  return result;
761 }
762 
763 static int Planimeter12() {
764  /* Area of arctic circle (not really -- adjunct to rhumb-area test) */
765  double points[2][2] = {{66.562222222, 0}, {66.562222222, 180}};
766  struct geod_geodesic g;
767  double perimeter, area;
768  int result = 0;
769  geod_init(&g, wgs84_a, wgs84_f);
770  planimeter(&g, points, 2, &perimeter, &area);
771  result += assertEquals(perimeter, 10465729, 1);
772  result += assertEquals(area, 0, 1);
773  return result;
774 }
775 
776 static int Planimeter13() {
777  /* Check encircling pole twice */
778  double points[6][2] = {{89,-360}, {89,-240}, {89,-120},
779  {89,0}, {89,120}, {89,240}};
780  struct geod_geodesic g;
781  double perimeter, area;
782  int result = 0;
783  geod_init(&g, wgs84_a, wgs84_f);
784  planimeter(&g, points, 6, &perimeter, &area);
785  result += assertEquals(perimeter, 1160741, 1);
786  result += assertEquals(area, 32415230256.0, 1);
787  return result;
788 }
789 
790 int main() {
791  int n = 0, i;
792  if ((i = testinverse())) {++n; printf("testinverse fail: %d\n", i);}
793  if ((i = testdirect())) {++n; printf("testdirect fail: %d\n", i);}
794  if ((i = testarcdirect())) {++n; printf("testarcdirect fail: %d\n", i);}
795  if ((i = GeodSolve0())) {++n; printf("GeodSolve0 fail: %d\n", i);}
796  if ((i = GeodSolve1())) {++n; printf("GeodSolve1 fail: %d\n", i);}
797  if ((i = GeodSolve2())) {++n; printf("GeodSolve2 fail: %d\n", i);}
798  if ((i = GeodSolve4())) {++n; printf("GeodSolve4 fail: %d\n", i);}
799  if ((i = GeodSolve5())) {++n; printf("GeodSolve5 fail: %d\n", i);}
800  if ((i = GeodSolve6())) {++n; printf("GeodSolve6 fail: %d\n", i);}
801  if ((i = GeodSolve9())) {++n; printf("GeodSolve9 fail: %d\n", i);}
802  if ((i = GeodSolve10())) {++n; printf("GeodSolve10 fail: %d\n", i);}
803  if ((i = GeodSolve11())) {++n; printf("GeodSolve11 fail: %d\n", i);}
804  if ((i = GeodSolve12())) {++n; printf("GeodSolve12 fail: %d\n", i);}
805  if ((i = GeodSolve14())) {++n; printf("GeodSolve14 fail: %d\n", i);}
806  if ((i = GeodSolve15())) {++n; printf("GeodSolve15 fail: %d\n", i);}
807  if ((i = GeodSolve17())) {++n; printf("GeodSolve17 fail: %d\n", i);}
808  if ((i = GeodSolve26())) {++n; printf("GeodSolve26 fail: %d\n", i);}
809  if ((i = GeodSolve28())) {++n; printf("GeodSolve28 fail: %d\n", i);}
810  if ((i = GeodSolve33())) {++n; printf("GeodSolve33 fail: %d\n", i);}
811  if ((i = GeodSolve55())) {++n; printf("GeodSolve55 fail: %d\n", i);}
812  if ((i = GeodSolve59())) {++n; printf("GeodSolve59 fail: %d\n", i);}
813  if ((i = GeodSolve61())) {++n; printf("GeodSolve61 fail: %d\n", i);}
814  if ((i = GeodSolve65())) {++n; printf("GeodSolve65 fail: %d\n", i);}
815  if ((i = GeodSolve67())) {++n; printf("GeodSolve67 fail: %d\n", i);}
816  if ((i = GeodSolve71())) {++n; printf("GeodSolve71 fail: %d\n", i);}
817  if ((i = GeodSolve73())) {++n; printf("GeodSolve73 fail: %d\n", i);}
818  if ((i = GeodSolve74())) {++n; printf("GeodSolve74 fail: %d\n", i);}
819  if ((i = GeodSolve76())) {++n; printf("GeodSolve76 fail: %d\n", i);}
820  if ((i = GeodSolve78())) {++n; printf("GeodSolve78 fail: %d\n", i);}
821  if ((i = Planimeter0())) {++n; printf("Planimeter0 fail: %d\n", i);}
822  if ((i = Planimeter5())) {++n; printf("Planimeter5 fail: %d\n", i);}
823  if ((i = Planimeter6())) {++n; printf("Planimeter6 fail: %d\n", i);}
824  if ((i = Planimeter12())) {++n; printf("Planimeter12 fail: %d\n", i);}
825  if ((i = Planimeter13())) {++n; printf("Planimeter13 fail: %d\n", i);}
826  return n;
827 }
828 
geod_inverseline
void geod_inverseline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, unsigned caps)
geod_polygon_compute
unsigned geod_polygon_compute(const struct geod_geodesic *g, const struct geod_polygon *p, int reverse, int sign, double *pA, double *pP)
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
d
static const double d[K][N]
Definition: igam.h:11
geod_direct
void geod_direct(const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, double *plat2, double *plon2, double *pazi2)
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
GEOD_LONG_UNROLL
@ GEOD_LONG_UNROLL
Definition: geodesic.h:903
geodesic.h
API for the geodesic routines in C.
GEOD_ARCMODE
@ GEOD_ARCMODE
Definition: geodesic.h:902
result
Values result
Definition: OdometryOptimize.cpp:8
main
int main(int argc, char **argv)
Definition: cmake/example_cmake_find_gtsam/main.cpp:63
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
geod_geodesic
Definition: geodesic.h:168
n
int n
Definition: BiCGSTAB_simple.cpp:1
geod_lineinit
void geod_lineinit(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned caps)
geod_init
void geod_init(struct geod_geodesic *g, double a, double f)
geod_position
void geod_position(const struct geod_geodesicline *l, double s12, double *plat2, double *plon2, double *pazi2)
geod_geodesicline
Definition: geodesic.h:182
geod_geninverse
double geod_geninverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2, double *pm12, double *pM12, double *pM21, double *pS12)
geod_polygon_addpoint
void geod_polygon_addpoint(const struct geod_geodesic *g, struct geod_polygon *p, double lat, double lon)
l
static const Line3 l(Rot3(), 1, 1)
pc
int RealScalar int RealScalar int RealScalar * pc
Definition: level1_cplx_impl.h:119
g
void g(const string &key, int i)
Definition: testBTree.cpp:41
y
Scalar * y
Definition: level1_cplx_impl.h:124
geod_directline
void geod_directline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, unsigned caps)
geod_genposition
double geod_genposition(const struct geod_geodesicline *l, unsigned flags, double s12_a12, double *plat2, double *plon2, double *pazi2, double *ps12, double *pm12, double *pM12, double *pM21, double *pS12)
geod_inverse
void geod_inverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2)
geod_polygon_init
void geod_polygon_init(struct geod_polygon *p, int polylinep)
GEOD_ALL
@ GEOD_ALL
Definition: geodesic.h:893
p
float * p
Definition: Tutorial_Map_using.cpp:9
geod_gendirect
double geod_gendirect(const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned flags, double s12_a12, double *plat2, double *plon2, double *pazi2, double *ps12, double *pm12, double *pM12, double *pM21, double *pS12)
N
#define N
Definition: igam.h:9
geod_polygon
Definition: geodesic.h:206
geod_geodesicline::caps
unsigned caps
Definition: geodesic.h:198
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9

gtsam
Author(s):
autogenerated on Sun Dec 22 2024 04:11:37