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 
Scalar * y
void geod_directline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, unsigned caps)
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)
int main(int argc, char **argv)
int RealScalar int RealScalar int RealScalar * pc
int n
#define N
Definition: gksort.c:12
Real fabs(const Real &a)
void geod_position(const struct geod_geodesicline *l, double s12, double *plat2, double *plon2, double *pazi2)
void geod_lineinit(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned caps)
void g(const string &key, int i)
Definition: testBTree.cpp:41
unsigned caps
Definition: geodesic.h:198
static const Line3 l(Rot3(), 1, 1)
Values result
void geod_inverseline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, unsigned caps)
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)
void geod_polygon_addpoint(const struct geod_geodesic *g, struct geod_polygon *p, double lat, double lon)
Array< double, 1, 3 > e(1./3., 0.5, 2.)
void geod_polygon_init(struct geod_polygon *p, int polylinep)
void geod_direct(const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, double *plat2, double *plon2, double *pazi2)
unsigned geod_polygon_compute(const struct geod_geodesic *g, const struct geod_polygon *p, int reverse, int sign, double *pA, double *pP)
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)
float * p
void geod_inverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2)
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
void geod_init(struct geod_geodesic *g, double a, double f)
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
API for the geodesic routines in C.


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:18