C library for Geodesics  1.50.1
geodtest.c
Go to the documentation of this file.
1 /**
2  * \file geodtest.c
3  * \brief Test suite for the geodesic routines in C
4  *
5  * Run these tests by configuring with cmake and running "make test".
6  *
7  * Copyright (c) Charles Karney (2015-2019) <charles@karney.com> and licensed
8  * under the MIT/X11 License. For more information, see
9  * https://geographiclib.sourceforge.io/
10  **********************************************************************/
11 
12 /** @cond SKIP */
13 
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 #if !defined(__cplusplus)
24 #define nullptr 0
25 #endif
26 
27 static const double wgs84_a = 6378137, wgs84_f = 1/298.257223563; /* WGS84 */
28 
29 static int checkEquals(double x, double y, double d) {
30  if (fabs(x - y) <= d)
31  return 0;
32  printf("checkEquals fails: %.7g != %.7g +/- %.7g\n", x, y, d);
33  return 1;
34 }
35 
36 static int checkNaN(double x) {
37  /* cppcheck-suppress duplicateExpression */
38  if (x != x)
39  return 0;
40  printf("checkNaN fails: %.7g\n", x);
41  return 1;
42 }
43 
44 static const int ncases = 20;
45 static const double testcases[20][12] = {
46  {35.60777, -139.44815, 111.098748429560326,
47  -11.17491, -69.95921, 129.289270889708762,
48  8935244.5604818305, 80.50729714281974, 6273170.2055303837,
49  0.16606318447386067, 0.16479116945612937, 12841384694976.432},
50  {55.52454, 106.05087, 22.020059880982801,
51  77.03196, 197.18234, 109.112041110671519,
52  4105086.1713924406, 36.892740690445894, 3828869.3344387607,
53  0.80076349608092607, 0.80101006984201008, 61674961290615.615},
54  {-21.97856, 142.59065, -32.44456876433189,
55  41.84138, 98.56635, -41.84359951440466,
56  8394328.894657671, 75.62930491011522, 6161154.5773110616,
57  0.24816339233950381, 0.24930251203627892, -6637997720646.717},
58  {-66.99028, 112.2363, 173.73491240878403,
59  -12.70631, 285.90344, 2.512956620913668,
60  11150344.2312080241, 100.278634181155759, 6289939.5670446687,
61  -0.17199490274700385, -0.17722569526345708, -121287239862139.744},
62  {-17.42761, 173.34268, -159.033557661192928,
63  -15.84784, 5.93557, -20.787484651536988,
64  16076603.1631180673, 144.640108810286253, 3732902.1583877189,
65  -0.81273638700070476, -0.81299800519154474, 97825992354058.708},
66  {32.84994, 48.28919, 150.492927788121982,
67  -56.28556, 202.29132, 48.113449399816759,
68  16727068.9438164461, 150.565799985466607, 3147838.1910180939,
69  -0.87334918086923126, -0.86505036767110637, -72445258525585.010},
70  {6.96833, 52.74123, 92.581585386317712,
71  -7.39675, 206.17291, 90.721692165923907,
72  17102477.2496958388, 154.147366239113561, 2772035.6169917581,
73  -0.89991282520302447, -0.89986892177110739, -1311796973197.995},
74  {-50.56724, -16.30485, -105.439679907590164,
75  -33.56571, -94.97412, -47.348547835650331,
76  6455670.5118668696, 58.083719495371259, 5409150.7979815838,
77  0.53053508035997263, 0.52988722644436602, 41071447902810.047},
78  {-58.93002, -8.90775, 140.965397902500679,
79  -8.91104, 133.13503, 19.255429433416599,
80  11756066.0219864627, 105.755691241406877, 6151101.2270708536,
81  -0.26548622269867183, -0.27068483874510741, -86143460552774.735},
82  {-68.82867, -74.28391, 93.774347763114881,
83  -50.63005, -8.36685, 34.65564085411343,
84  3956936.926063544, 35.572254987389284, 3708890.9544062657,
85  0.81443963736383502, 0.81420859815358342, -41845309450093.787},
86  {-10.62672, -32.0898, -86.426713286747751,
87  5.883, -134.31681, -80.473780971034875,
88  11470869.3864563009, 103.387395634504061, 6184411.6622659713,
89  -0.23138683500430237, -0.23155097622286792, 4198803992123.548},
90  {-21.76221, 166.90563, 29.319421206936428,
91  48.72884, 213.97627, 43.508671946410168,
92  9098627.3986554915, 81.963476716121964, 6299240.9166992283,
93  0.13965943368590333, 0.14152969707656796, 10024709850277.476},
94  {-19.79938, -174.47484, 71.167275780171533,
95  -11.99349, -154.35109, 65.589099775199228,
96  2319004.8601169389, 20.896611684802389, 2267960.8703918325,
97  0.93427001867125849, 0.93424887135032789, -3935477535005.785},
98  {-11.95887, -116.94513, 92.712619830452549,
99  4.57352, 7.16501, 78.64960934409585,
100  13834722.5801401374, 124.688684161089762, 5228093.177931598,
101  -0.56879356755666463, -0.56918731952397221, -9919582785894.853},
102  {-87.85331, 85.66836, -65.120313040242748,
103  66.48646, 16.09921, -4.888658719272296,
104  17286615.3147144645, 155.58592449699137, 2635887.4729110181,
105  -0.90697975771398578, -0.91095608883042767, 42667211366919.534},
106  {1.74708, 128.32011, -101.584843631173858,
107  -11.16617, 11.87109, -86.325793296437476,
108  12942901.1241347408, 116.650512484301857, 5682744.8413270572,
109  -0.44857868222697644, -0.44824490340007729, 10763055294345.653},
110  {-25.72959, -144.90758, -153.647468693117198,
111  -57.70581, -269.17879, -48.343983158876487,
112  9413446.7452453107, 84.664533838404295, 6356176.6898881281,
113  0.09492245755254703, 0.09737058264766572, 74515122850712.444},
114  {-41.22777, 122.32875, 14.285113402275739,
115  -7.57291, 130.37946, 10.805303085187369,
116  3812686.035106021, 34.34330804743883, 3588703.8812128856,
117  0.82605222593217889, 0.82572158200920196, -2456961531057.857},
118  {11.01307, 138.25278, 79.43682622782374,
119  6.62726, 247.05981, 103.708090215522657,
120  11911190.819018408, 107.341669954114577, 6070904.722786735,
121  -0.29767608923657404, -0.29785143390252321, 17121631423099.696},
122  {-29.47124, 95.14681, -163.779130441688382,
123  -27.46601, -69.15955, -15.909335945554969,
124  13487015.8381145492, 121.294026715742277, 5481428.9945736388,
125  -0.51527225545373252, -0.51556587964721788, 104679964020340.318}};
126 
127 static int testinverse() {
128  double lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12;
129  double azi1a, azi2a, s12a, a12a, m12a, M12a, M21a, S12a;
130  struct geod_geodesic g;
131  int i, result = 0;
132  geod_init(&g, wgs84_a, wgs84_f);
133  for (i = 0; i < ncases; ++i) {
134  lat1 = testcases[i][0]; lon1 = testcases[i][1]; azi1 = testcases[i][2];
135  lat2 = testcases[i][3]; lon2 = testcases[i][4]; azi2 = testcases[i][5];
136  s12 = testcases[i][6]; a12 = testcases[i][7]; m12 = testcases[i][8];
137  M12 = testcases[i][9]; M21 = testcases[i][10]; S12 = testcases[i][11];
138  a12a = geod_geninverse(&g, lat1, lon1, lat2, lon2, &s12a, &azi1a, &azi2a,
139  &m12a, &M12a, &M21a, &S12a);
140  result += checkEquals(azi1, azi1a, 1e-13);
141  result += checkEquals(azi2, azi2a, 1e-13);
142  result += checkEquals(s12, s12a, 1e-8);
143  result += checkEquals(a12, a12a, 1e-13);
144  result += checkEquals(m12, m12a, 1e-8);
145  result += checkEquals(M12, M12a, 1e-15);
146  result += checkEquals(M21, M21a, 1e-15);
147  result += checkEquals(S12, S12a, 0.1);
148  }
149  return result;
150 }
151 
152 static int testdirect() {
153  double lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12;
154  double lat2a, lon2a, azi2a, a12a, m12a, M12a, M21a, S12a;
155  struct geod_geodesic g;
156  int i, result = 0;
157  unsigned flags = GEOD_LONG_UNROLL;
158  geod_init(&g, wgs84_a, wgs84_f);
159  for (i = 0; i < ncases; ++i) {
160  lat1 = testcases[i][0]; lon1 = testcases[i][1]; azi1 = testcases[i][2];
161  lat2 = testcases[i][3]; lon2 = testcases[i][4]; azi2 = testcases[i][5];
162  s12 = testcases[i][6]; a12 = testcases[i][7]; m12 = testcases[i][8];
163  M12 = testcases[i][9]; M21 = testcases[i][10]; S12 = testcases[i][11];
164  a12a = geod_gendirect(&g, lat1, lon1, azi1, flags, s12,
165  &lat2a, &lon2a, &azi2a, nullptr, &m12a, &M12a, &M21a, &S12a);
166  result += checkEquals(lat2, lat2a, 1e-13);
167  result += checkEquals(lon2, lon2a, 1e-13);
168  result += checkEquals(azi2, azi2a, 1e-13);
169  result += checkEquals(a12, a12a, 1e-13);
170  result += checkEquals(m12, m12a, 1e-8);
171  result += checkEquals(M12, M12a, 1e-15);
172  result += checkEquals(M21, M21a, 1e-15);
173  result += checkEquals(S12, S12a, 0.1);
174  }
175  return result;
176 }
177 
178 static int testarcdirect() {
179  double lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, M12, M21, S12;
180  double lat2a, lon2a, azi2a, s12a, m12a, M12a, M21a, S12a;
181  struct geod_geodesic g;
182  int i, result = 0;
183  unsigned flags = GEOD_ARCMODE | GEOD_LONG_UNROLL;
184  geod_init(&g, wgs84_a, wgs84_f);
185  for (i = 0; i < ncases; ++i) {
186  lat1 = testcases[i][0]; lon1 = testcases[i][1]; azi1 = testcases[i][2];
187  lat2 = testcases[i][3]; lon2 = testcases[i][4]; azi2 = testcases[i][5];
188  s12 = testcases[i][6]; a12 = testcases[i][7]; m12 = testcases[i][8];
189  M12 = testcases[i][9]; M21 = testcases[i][10]; S12 = testcases[i][11];
190  geod_gendirect(&g, lat1, lon1, azi1, flags, a12,
191  &lat2a, &lon2a, &azi2a, &s12a, &m12a, &M12a, &M21a, &S12a);
192  result += checkEquals(lat2, lat2a, 1e-13);
193  result += checkEquals(lon2, lon2a, 1e-13);
194  result += checkEquals(azi2, azi2a, 1e-13);
195  result += checkEquals(s12, s12a, 1e-8);
196  result += checkEquals(m12, m12a, 1e-8);
197  result += checkEquals(M12, M12a, 1e-15);
198  result += checkEquals(M21, M21a, 1e-15);
199  result += checkEquals(S12, S12a, 0.1);
200  }
201  return result;
202 }
203 
204 static int GeodSolve0() {
205  double azi1, azi2, s12;
206  struct geod_geodesic g;
207  int result = 0;
208  geod_init(&g, wgs84_a, wgs84_f);
209  geod_inverse(&g, 40.6, -73.8, 49.01666667, 2.55, &s12, &azi1, &azi2);
210  result += checkEquals(azi1, 53.47022, 0.5e-5);
211  result += checkEquals(azi2, 111.59367, 0.5e-5);
212  result += checkEquals(s12, 5853226, 0.5);
213  return result;
214 }
215 
216 static int GeodSolve1() {
217  double lat2, lon2, azi2;
218  struct geod_geodesic g;
219  int result = 0;
220  geod_init(&g, wgs84_a, wgs84_f);
221  geod_direct(&g, 40.63972222, -73.77888889, 53.5, 5850e3,
222  &lat2, &lon2, &azi2);
223  result += checkEquals(lat2, 49.01467, 0.5e-5);
224  result += checkEquals(lon2, 2.56106, 0.5e-5);
225  result += checkEquals(azi2, 111.62947, 0.5e-5);
226  return result;
227 }
228 
229 static int GeodSolve2() {
230  /* Check fix for antipodal prolate bug found 2010-09-04 */
231  double azi1, azi2, s12;
232  struct geod_geodesic g;
233  int result = 0;
234  geod_init(&g, 6.4e6, -1/150.0);
235  geod_inverse(&g, 0.07476, 0, -0.07476, 180, &s12, &azi1, &azi2);
236  result += checkEquals(azi1, 90.00078, 0.5e-5);
237  result += checkEquals(azi2, 90.00078, 0.5e-5);
238  result += checkEquals(s12, 20106193, 0.5);
239  geod_inverse(&g, 0.1, 0, -0.1, 180, &s12, &azi1, &azi2);
240  result += checkEquals(azi1, 90.00105, 0.5e-5);
241  result += checkEquals(azi2, 90.00105, 0.5e-5);
242  result += checkEquals(s12, 20106193, 0.5);
243  return result;
244 }
245 
246 static int GeodSolve4() {
247  /* Check fix for short line bug found 2010-05-21 */
248  double s12;
249  struct geod_geodesic g;
250  int result = 0;
251  geod_init(&g, wgs84_a, wgs84_f);
252  geod_inverse(&g, 36.493349428792, 0, 36.49334942879201, .0000008,
253  &s12, nullptr, nullptr);
254  result += checkEquals(s12, 0.072, 0.5e-3);
255  return result;
256 }
257 
258 static int GeodSolve5() {
259  /* Check fix for point2=pole bug found 2010-05-03 */
260  double lat2, lon2, azi2;
261  struct geod_geodesic g;
262  int result = 0;
263  geod_init(&g, wgs84_a, wgs84_f);
264  geod_direct(&g, 0.01777745589997, 30, 0, 10e6, &lat2, &lon2, &azi2);
265  result += checkEquals(lat2, 90, 0.5e-5);
266  if (lon2 < 0) {
267  result += checkEquals(lon2, -150, 0.5e-5);
268  result += checkEquals(fabs(azi2), 180, 0.5e-5);
269  } else {
270  result += checkEquals(lon2, 30, 0.5e-5);
271  result += checkEquals(azi2, 0, 0.5e-5);
272  }
273  return result;
274 }
275 
276 static int GeodSolve6() {
277  /* Check fix for volatile sbet12a bug found 2011-06-25 (gcc 4.4.4
278  * x86 -O3). Found again on 2012-03-27 with tdm-mingw32 (g++ 4.6.1). */
279  double s12;
280  struct geod_geodesic g;
281  int result = 0;
282  geod_init(&g, wgs84_a, wgs84_f);
283  geod_inverse(&g, 88.202499451857, 0,
284  -88.202499451857, 179.981022032992859592,
285  &s12, nullptr, nullptr);
286  result += checkEquals(s12, 20003898.214, 0.5e-3);
287  geod_inverse(&g, 89.262080389218, 0,
288  -89.262080389218, 179.992207982775375662,
289  &s12, nullptr, nullptr);
290  result += checkEquals(s12, 20003925.854, 0.5e-3);
291  geod_inverse(&g, 89.333123580033, 0,
292  -89.333123580032997687, 179.99295812360148422,
293  &s12, nullptr, nullptr);
294  result += checkEquals(s12, 20003926.881, 0.5e-3);
295  return result;
296 }
297 
298 static int GeodSolve9() {
299  /* Check fix for volatile x bug found 2011-06-25 (gcc 4.4.4 x86 -O3) */
300  double s12;
301  struct geod_geodesic g;
302  int result = 0;
303  geod_init(&g, wgs84_a, wgs84_f);
304  geod_inverse(&g, 56.320923501171, 0,
305  -56.320923501171, 179.664747671772880215,
306  &s12, nullptr, nullptr);
307  result += checkEquals(s12, 19993558.287, 0.5e-3);
308  return result;
309 }
310 
311 static int GeodSolve10() {
312  /* Check fix for adjust tol1_ bug found 2011-06-25 (Visual Studio
313  * 10 rel + debug) */
314  double s12;
315  struct geod_geodesic g;
316  int result = 0;
317  geod_init(&g, wgs84_a, wgs84_f);
318  geod_inverse(&g, 52.784459512564, 0,
319  -52.784459512563990912, 179.634407464943777557,
320  &s12, nullptr, nullptr);
321  result += checkEquals(s12, 19991596.095, 0.5e-3);
322  return result;
323 }
324 
325 static int GeodSolve11() {
326  /* Check fix for bet2 = -bet1 bug found 2011-06-25 (Visual Studio
327  * 10 rel + debug) */
328  double s12;
329  struct geod_geodesic g;
330  int result = 0;
331  geod_init(&g, wgs84_a, wgs84_f);
332  geod_inverse(&g, 48.522876735459, 0,
333  -48.52287673545898293, 179.599720456223079643,
334  &s12, nullptr, nullptr);
335  result += checkEquals(s12, 19989144.774, 0.5e-3);
336  return result;
337 }
338 
339 static int GeodSolve12() {
340  /* Check fix for inverse geodesics on extreme prolate/oblate
341  * ellipsoids Reported 2012-08-29 Stefan Guenther
342  * <stefan.gunther@embl.de>; fixed 2012-10-07 */
343  double azi1, azi2, s12;
344  struct geod_geodesic g;
345  int result = 0;
346  geod_init(&g, 89.8, -1.83);
347  geod_inverse(&g, 0, 0, -10, 160, &s12, &azi1, &azi2);
348  result += checkEquals(azi1, 120.27, 1e-2);
349  result += checkEquals(azi2, 105.15, 1e-2);
350  result += checkEquals(s12, 266.7, 1e-1);
351  return result;
352 }
353 
354 static int GeodSolve14() {
355  /* Check fix for inverse ignoring lon12 = nan */
356  double azi1, azi2, s12, nan;
357  struct geod_geodesic g;
358  int result = 0;
359  {
360  double minus1 = -1;
361  /* cppcheck-suppress wrongmathcall */
362  nan = sqrt(minus1);
363  }
364  geod_init(&g, wgs84_a, wgs84_f);
365  geod_inverse(&g, 0, 0, 1, nan, &s12, &azi1, &azi2);
366  result += checkNaN(azi1);
367  result += checkNaN(azi2);
368  result += checkNaN(s12);
369  return result;
370 }
371 
372 static int GeodSolve15() {
373  /* Initial implementation of Math::eatanhe was wrong for e^2 < 0. This
374  * checks that this is fixed. */
375  double S12;
376  struct geod_geodesic g;
377  int result = 0;
378  geod_init(&g, 6.4e6, -1/150.0);
379  geod_gendirect(&g, 1, 2, 3, 0, 4, nullptr, nullptr, nullptr,
380  nullptr, nullptr, nullptr, nullptr, &S12);
381  result += checkEquals(S12, 23700, 0.5);
382  return result;
383 }
384 
385 static int GeodSolve17() {
386  /* Check fix for LONG_UNROLL bug found on 2015-05-07 */
387  double lat2, lon2, azi2;
388  struct geod_geodesic g;
389  struct geod_geodesicline l;
390  int result = 0;
391  unsigned flags = GEOD_LONG_UNROLL;
392  geod_init(&g, wgs84_a, wgs84_f);
393  geod_gendirect(&g, 40, -75, -10, flags, 2e7, &lat2, &lon2, &azi2,
394  nullptr, nullptr, nullptr, nullptr, nullptr);
395  result += checkEquals(lat2, -39, 1);
396  result += checkEquals(lon2, -254, 1);
397  result += checkEquals(azi2, -170, 1);
398  geod_lineinit(&l, &g, 40, -75, -10, 0);
399  geod_genposition(&l, flags, 2e7, &lat2, &lon2, &azi2,
400  nullptr, nullptr, nullptr, nullptr, nullptr);
401  result += checkEquals(lat2, -39, 1);
402  result += checkEquals(lon2, -254, 1);
403  result += checkEquals(azi2, -170, 1);
404  geod_direct(&g, 40, -75, -10, 2e7, &lat2, &lon2, &azi2);
405  result += checkEquals(lat2, -39, 1);
406  result += checkEquals(lon2, 105, 1);
407  result += checkEquals(azi2, -170, 1);
408  geod_position(&l, 2e7, &lat2, &lon2, &azi2);
409  result += checkEquals(lat2, -39, 1);
410  result += checkEquals(lon2, 105, 1);
411  result += checkEquals(azi2, -170, 1);
412  return result;
413 }
414 
415 static int GeodSolve26() {
416  /* Check 0/0 problem with area calculation on sphere 2015-09-08 */
417  double S12;
418  struct geod_geodesic g;
419  int result = 0;
420  geod_init(&g, 6.4e6, 0);
421  geod_geninverse(&g, 1, 2, 3, 4, nullptr, nullptr, nullptr,
422  nullptr, nullptr, nullptr, &S12);
423  result += checkEquals(S12, 49911046115.0, 0.5);
424  return result;
425 }
426 
427 static int GeodSolve28() {
428  /* Check for bad placement of assignment of r.a12 with |f| > 0.01 (bug in
429  * Java implementation fixed on 2015-05-19). */
430  double a12;
431  struct geod_geodesic g;
432  int result = 0;
433  geod_init(&g, 6.4e6, 0.1);
434  a12 = geod_gendirect(&g, 1, 2, 10, 0, 5e6, nullptr, nullptr, nullptr,
435  nullptr, nullptr, nullptr, nullptr, nullptr);
436  result += checkEquals(a12, 48.55570690, 0.5e-8);
437  return result;
438 }
439 
440 static int GeodSolve33() {
441  /* Check max(-0.0,+0.0) issues 2015-08-22 (triggered by bugs in Octave --
442  * sind(-0.0) = +0.0 -- and in some version of Visual Studio --
443  * fmod(-0.0, 360.0) = +0.0. */
444  double azi1, azi2, s12;
445  struct geod_geodesic g;
446  int result = 0;
447  geod_init(&g, wgs84_a, wgs84_f);
448  geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2);
449  result += checkEquals(azi1, 90.00000, 0.5e-5);
450  result += checkEquals(azi2, 90.00000, 0.5e-5);
451  result += checkEquals(s12, 19926189, 0.5);
452  geod_inverse(&g, 0, 0, 0, 179.5, &s12, &azi1, &azi2);
453  result += checkEquals(azi1, 55.96650, 0.5e-5);
454  result += checkEquals(azi2, 124.03350, 0.5e-5);
455  result += checkEquals(s12, 19980862, 0.5);
456  geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2);
457  result += checkEquals(azi1, 0.00000, 0.5e-5);
458  result += checkEquals(fabs(azi2), 180.00000, 0.5e-5);
459  result += checkEquals(s12, 20003931, 0.5);
460  geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2);
461  result += checkEquals(azi1, 0.00000, 0.5e-5);
462  result += checkEquals(fabs(azi2), 180.00000, 0.5e-5);
463  result += checkEquals(s12, 19893357, 0.5);
464  geod_init(&g, 6.4e6, 0);
465  geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2);
466  result += checkEquals(azi1, 90.00000, 0.5e-5);
467  result += checkEquals(azi2, 90.00000, 0.5e-5);
468  result += checkEquals(s12, 19994492, 0.5);
469  geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2);
470  result += checkEquals(azi1, 0.00000, 0.5e-5);
471  result += checkEquals(fabs(azi2), 180.00000, 0.5e-5);
472  result += checkEquals(s12, 20106193, 0.5);
473  geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2);
474  result += checkEquals(azi1, 0.00000, 0.5e-5);
475  result += checkEquals(fabs(azi2), 180.00000, 0.5e-5);
476  result += checkEquals(s12, 19994492, 0.5);
477  geod_init(&g, 6.4e6, -1/300.0);
478  geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2);
479  result += checkEquals(azi1, 90.00000, 0.5e-5);
480  result += checkEquals(azi2, 90.00000, 0.5e-5);
481  result += checkEquals(s12, 19994492, 0.5);
482  geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2);
483  result += checkEquals(azi1, 90.00000, 0.5e-5);
484  result += checkEquals(azi2, 90.00000, 0.5e-5);
485  result += checkEquals(s12, 20106193, 0.5);
486  geod_inverse(&g, 0, 0, 0.5, 180, &s12, &azi1, &azi2);
487  result += checkEquals(azi1, 33.02493, 0.5e-5);
488  result += checkEquals(azi2, 146.97364, 0.5e-5);
489  result += checkEquals(s12, 20082617, 0.5);
490  geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2);
491  result += checkEquals(azi1, 0.00000, 0.5e-5);
492  result += checkEquals(fabs(azi2), 180.00000, 0.5e-5);
493  result += checkEquals(s12, 20027270, 0.5);
494 
495  return result;
496 }
497 
498 static int GeodSolve55() {
499  /* Check fix for nan + point on equator or pole not returning all nans in
500  * Geodesic::Inverse, found 2015-09-23. */
501  double azi1, azi2, s12, nan;
502  struct geod_geodesic g;
503  int result = 0;
504  {
505  double minus1 = -1;
506  /* cppcheck-suppress wrongmathcall */
507  nan = sqrt(minus1);
508  }
509  geod_init(&g, wgs84_a, wgs84_f);
510  geod_inverse(&g, nan, 0, 0, 90, &s12, &azi1, &azi2);
511  result += checkNaN(azi1);
512  result += checkNaN(azi2);
513  result += checkNaN(s12);
514  geod_inverse(&g, nan, 0, 90, 9, &s12, &azi1, &azi2);
515  result += checkNaN(azi1);
516  result += checkNaN(azi2);
517  result += checkNaN(s12);
518  return result;
519 }
520 
521 static int GeodSolve59() {
522  /* Check for points close with longitudes close to 180 deg apart. */
523  double azi1, azi2, s12;
524  struct geod_geodesic g;
525  int result = 0;
526  geod_init(&g, wgs84_a, wgs84_f);
527  geod_inverse(&g, 5, 0.00000000000001, 10, 180, &s12, &azi1, &azi2);
528  result += checkEquals(azi1, 0.000000000000035, 1.5e-14);
529  result += checkEquals(azi2, 179.99999999999996, 1.5e-14);
530  result += checkEquals(s12, 18345191.174332713, 5e-9);
531  return result;
532 }
533 
534 static int GeodSolve61() {
535  /* Make sure small negative azimuths are west-going */
536  double lat2, lon2, azi2;
537  struct geod_geodesic g;
538  struct geod_geodesicline l;
539  int result = 0;
540  unsigned flags = GEOD_LONG_UNROLL;
541  geod_init(&g, wgs84_a, wgs84_f);
542  geod_gendirect(&g, 45, 0, -0.000000000000000003, flags, 1e7,
543  &lat2, &lon2, &azi2,
544  nullptr, nullptr, nullptr, nullptr, nullptr);
545  result += checkEquals(lat2, 45.30632, 0.5e-5);
546  result += checkEquals(lon2, -180, 0.5e-5);
547  result += checkEquals(fabs(azi2), 180, 0.5e-5);
548  geod_inverseline(&l, &g, 45, 0, 80, -0.000000000000000003, 0);
549  geod_genposition(&l, flags, 1e7, &lat2, &lon2, &azi2,
550  nullptr, nullptr, nullptr, nullptr, nullptr);
551  result += checkEquals(lat2, 45.30632, 0.5e-5);
552  result += checkEquals(lon2, -180, 0.5e-5);
553  result += checkEquals(fabs(azi2), 180, 0.5e-5);
554  return result;
555 }
556 
557 static int GeodSolve65() {
558  /* Check for bug in east-going check in GeodesicLine (needed to check for
559  * sign of 0) and sign error in area calculation due to a bogus override of
560  * the code for alp12. Found/fixed on 2015-12-19. */
561  double lat2, lon2, azi2, s12, a12, m12, M12, M21, S12;
562  struct geod_geodesic g;
563  struct geod_geodesicline l;
564  int result = 0;
565  unsigned flags = GEOD_LONG_UNROLL, caps = GEOD_ALL;
566  geod_init(&g, wgs84_a, wgs84_f);
567  geod_inverseline(&l, &g, 30, -0.000000000000000001, -31, 180, caps);
568  a12 = geod_genposition(&l, flags, 1e7,
569  &lat2, &lon2, &azi2, &s12, &m12, &M12, &M21, &S12);
570  result += checkEquals(lat2, -60.23169, 0.5e-5);
571  result += checkEquals(lon2, -0.00000, 0.5e-5);
572  result += checkEquals(fabs(azi2), 180.00000, 0.5e-5);
573  result += checkEquals(s12, 10000000, 0.5);
574  result += checkEquals(a12, 90.06544, 0.5e-5);
575  result += checkEquals(m12, 6363636, 0.5);
576  result += checkEquals(M12, -0.0012834, 0.5e-7);
577  result += checkEquals(M21, 0.0013749, 0.5e-7);
578  result += checkEquals(S12, 0, 0.5);
579  a12 = geod_genposition(&l, flags, 2e7,
580  &lat2, &lon2, &azi2, &s12, &m12, &M12, &M21, &S12);
581  result += checkEquals(lat2, -30.03547, 0.5e-5);
582  result += checkEquals(lon2, -180.00000, 0.5e-5);
583  result += checkEquals(azi2, -0.00000, 0.5e-5);
584  result += checkEquals(s12, 20000000, 0.5);
585  result += checkEquals(a12, 179.96459, 0.5e-5);
586  result += checkEquals(m12, 54342, 0.5);
587  result += checkEquals(M12, -1.0045592, 0.5e-7);
588  result += checkEquals(M21, -0.9954339, 0.5e-7);
589  result += checkEquals(S12, 127516405431022.0, 0.5);
590  return result;
591 }
592 
593 static int GeodSolve67() {
594  /* Check for InverseLine if line is slightly west of S and that s13 is
595  * correctly set. */
596  double lat2, lon2, azi2;
597  struct geod_geodesic g;
598  struct geod_geodesicline l;
599  int result = 0;
600  unsigned flags = GEOD_LONG_UNROLL;
601  geod_init(&g, wgs84_a, wgs84_f);
602  geod_inverseline(&l, &g, -5, -0.000000000000002, -10, 180, 0);
603  geod_genposition(&l, flags, 2e7, &lat2, &lon2, &azi2,
604  nullptr, nullptr, nullptr, nullptr, nullptr);
605  result += checkEquals(lat2, 4.96445, 0.5e-5);
606  result += checkEquals(lon2, -180.00000, 0.5e-5);
607  result += checkEquals(azi2, -0.00000, 0.5e-5);
608  geod_genposition(&l, flags, 0.5 * l.s13, &lat2, &lon2, &azi2,
609  nullptr, nullptr, nullptr, nullptr, nullptr);
610  result += checkEquals(lat2, -87.52461, 0.5e-5);
611  result += checkEquals(lon2, -0.00000, 0.5e-5);
612  result += checkEquals(azi2, -180.00000, 0.5e-5);
613  return result;
614 }
615 
616 static int GeodSolve71() {
617  /* Check that DirectLine sets s13. */
618  double lat2, lon2, azi2;
619  struct geod_geodesic g;
620  struct geod_geodesicline l;
621  int result = 0;
622  geod_init(&g, wgs84_a, wgs84_f);
623  geod_directline(&l, &g, 1, 2, 45, 1e7, 0);
624  geod_position(&l, 0.5 * l.s13, &lat2, &lon2, &azi2);
625  result += checkEquals(lat2, 30.92625, 0.5e-5);
626  result += checkEquals(lon2, 37.54640, 0.5e-5);
627  result += checkEquals(azi2, 55.43104, 0.5e-5);
628  return result;
629 }
630 
631 static int GeodSolve73() {
632  /* Check for backwards from the pole bug reported by Anon on 2016-02-13.
633  * This only affected the Java implementation. It was introduced in Java
634  * version 1.44 and fixed in 1.46-SNAPSHOT on 2016-01-17.
635  * Also the + sign on azi2 is a check on the normalizing of azimuths
636  * (converting -0.0 to +0.0). */
637  double lat2, lon2, azi2;
638  struct geod_geodesic g;
639  int result = 0;
640  geod_init(&g, wgs84_a, wgs84_f);
641  geod_direct(&g, 90, 10, 180, -1e6,
642  &lat2, &lon2, &azi2);
643  result += checkEquals(lat2, 81.04623, 0.5e-5);
644  result += checkEquals(lon2, -170, 0.5e-5);
645  result += azi2 == 0 ? 0 : 1;
646  result += 1/azi2 > 0 ? 0 : 1; /* Check that azi2 = +0.0 not -0.0 */
647  return result;
648 }
649 
650 static void planimeter(const struct geod_geodesic* g,
651  double points[][2], int N,
652  double* perimeter, double* area) {
653  struct geod_polygon p;
654  int i;
655  geod_polygon_init(&p, 0);
656  for (i = 0; i < N; ++i)
657  geod_polygon_addpoint(g, &p, points[i][0], points[i][1]);
658  geod_polygon_compute(g, &p, 0, 1, area, perimeter);
659 }
660 
661 static void polylength(const struct geod_geodesic* g,
662  double points[][2], int N,
663  double* perimeter) {
664  struct geod_polygon p;
665  int i;
666  geod_polygon_init(&p, 1);
667  for (i = 0; i < N; ++i)
668  geod_polygon_addpoint(g, &p, points[i][0], points[i][1]);
669  geod_polygon_compute(g, &p, 0, 1, nullptr, perimeter);
670 }
671 
672 static int GeodSolve74() {
673  /* Check fix for inaccurate areas, bug introduced in v1.46, fixed
674  * 2015-10-16. */
675  double a12, s12, azi1, azi2, m12, M12, M21, S12;
676  struct geod_geodesic g;
677  int result = 0;
678  geod_init(&g, wgs84_a, wgs84_f);
679  a12 = geod_geninverse(&g, 54.1589, 15.3872, 54.1591, 15.3877,
680  &s12, &azi1, &azi2, &m12, &M12, &M21, &S12);
681  result += checkEquals(azi1, 55.723110355, 5e-9);
682  result += checkEquals(azi2, 55.723515675, 5e-9);
683  result += checkEquals(s12, 39.527686385, 5e-9);
684  result += checkEquals(a12, 0.000355495, 5e-9);
685  result += checkEquals(m12, 39.527686385, 5e-9);
686  result += checkEquals(M12, 0.999999995, 5e-9);
687  result += checkEquals(M21, 0.999999995, 5e-9);
688  result += checkEquals(S12, 286698586.30197, 5e-4);
689  return result;
690 }
691 
692 static int GeodSolve76() {
693  /* The distance from Wellington and Salamanca (a classic failure of
694  * Vincenty) */
695  double azi1, azi2, s12;
696  struct geod_geodesic g;
697  int result = 0;
698  geod_init(&g, wgs84_a, wgs84_f);
699  geod_inverse(&g, -(41+19/60.0), 174+49/60.0, 40+58/60.0, -(5+30/60.0),
700  &s12, &azi1, &azi2);
701  result += checkEquals(azi1, 160.39137649664, 0.5e-11);
702  result += checkEquals(azi2, 19.50042925176, 0.5e-11);
703  result += checkEquals(s12, 19960543.857179, 0.5e-6);
704  return result;
705 }
706 
707 static int GeodSolve78() {
708  /* An example where the NGS calculator fails to converge */
709  double azi1, azi2, s12;
710  struct geod_geodesic g;
711  int result = 0;
712  geod_init(&g, wgs84_a, wgs84_f);
713  geod_inverse(&g, 27.2, 0.0, -27.1, 179.5, &s12, &azi1, &azi2);
714  result += checkEquals(azi1, 45.82468716758, 0.5e-11);
715  result += checkEquals(azi2, 134.22776532670, 0.5e-11);
716  result += checkEquals(s12, 19974354.765767, 0.5e-6);
717  return result;
718 }
719 
720 static int GeodSolve80() {
721  /* Some tests to add code coverage: computing scale in special cases + zero
722  * length geodesic (includes GeodSolve80 - GeodSolve83) + using an incapable
723  * line. */
724  double a12, s12, azi1, azi2, m12, M12, M21, S12;
725  struct geod_geodesic g;
726  struct geod_geodesicline l;
727  int result = 0;
728  geod_init(&g, wgs84_a, wgs84_f);
729 
730  geod_geninverse(&g, 0, 0, 0, 90, nullptr, nullptr, nullptr,
731  nullptr, &M12, &M21, nullptr);
732  result += checkEquals(M12, -0.00528427534, 0.5e-10);
733  result += checkEquals(M21, -0.00528427534, 0.5e-10);
734 
735  geod_geninverse(&g, 0, 0, 1e-6, 1e-6, nullptr, nullptr, nullptr,
736  nullptr, &M12, &M21, nullptr);
737  result += checkEquals(M12, 1, 0.5e-10);
738  result += checkEquals(M21, 1, 0.5e-10);
739 
740  a12 = geod_geninverse(&g, 20.001, 0, 20.001, 0,
741  &s12, &azi1, &azi2, &m12, &M12, &M21, &S12);
742  result += checkEquals(a12, 0, 1e-13);
743  result += checkEquals(s12, 0, 1e-8);
744  result += checkEquals(azi1, 180, 1e-13);
745  result += checkEquals(azi2, 180, 1e-13);
746  result += checkEquals(m12, 0, 1e-8);
747  result += checkEquals(M12, 1, 1e-15);
748  result += checkEquals(M21, 1, 1e-15);
749  result += checkEquals(S12, 0, 1e-10);
750 
751  a12 = geod_geninverse(&g, 90, 0, 90, 180,
752  &s12, &azi1, &azi2, &m12, &M12, &M21, &S12);
753  result += checkEquals(a12, 0, 1e-13);
754  result += checkEquals(s12, 0, 1e-8);
755  result += checkEquals(azi1, 0, 1e-13);
756  result += checkEquals(azi2, 180, 1e-13);
757  result += checkEquals(m12, 0, 1e-8);
758  result += checkEquals(M12, 1, 1e-15);
759  result += checkEquals(M21, 1, 1e-15);
760  result += checkEquals(S12, 127516405431022.0, 0.5);
761 
762  /* An incapable line which can't take distance as input */
763  geod_lineinit(&l, &g, 1, 2, 90, GEOD_LATITUDE);
764  a12 = geod_genposition(&l, 0, 1000, nullptr, nullptr, nullptr,
765  nullptr, nullptr, nullptr, nullptr, nullptr);
766  result += checkNaN(a12);
767  return result;
768 }
769 
770 static int GeodSolve84() {
771  /* Tests for python implementation to check fix for range errors with
772  * {fmod,sin,cos}(inf) (includes GeodSolve84 - GeodSolve86). */
773 
774  double lat2, lon2, azi2, inf, nan;
775  struct geod_geodesic g;
776  int result = 0;
777  geod_init(&g, wgs84_a, wgs84_f);
778  {
779  /* a round about way to set inf = 0 */
780  geod_direct(&g, 0, 0, 90, 0, &inf, nullptr, nullptr);
781  /* so that this doesn't give a compiler time error on Windows */
782  inf = 1.0/inf;
783  }
784  {
785  double minus1 = -1;
786  /* cppcheck-suppress wrongmathcall */
787  nan = sqrt(minus1);
788  }
789  geod_direct(&g, 0, 0, 90, inf, &lat2, &lon2, &azi2);
790  result += checkNaN(lat2);
791  result += checkNaN(lon2);
792  result += checkNaN(azi2);
793  geod_direct(&g, 0, 0, 90, nan, &lat2, &lon2, &azi2);
794  result += checkNaN(lat2);
795  result += checkNaN(lon2);
796  result += checkNaN(azi2);
797  geod_direct(&g, 0, 0, inf, 1000, &lat2, &lon2, &azi2);
798  result += checkNaN(lat2);
799  result += checkNaN(lon2);
800  result += checkNaN(azi2);
801  geod_direct(&g, 0, 0, nan, 1000, &lat2, &lon2, &azi2);
802  result += checkNaN(lat2);
803  result += checkNaN(lon2);
804  result += checkNaN(azi2);
805  geod_direct(&g, 0, inf, 90, 1000, &lat2, &lon2, &azi2);
806  result += lat2 == 0 ? 0 : 1;
807  result += checkNaN(lon2);
808  result += azi2 == 90 ? 0 : 1;
809  geod_direct(&g, 0, nan, 90, 1000, &lat2, &lon2, &azi2);
810  result += lat2 == 0 ? 0 : 1;
811  result += checkNaN(lon2);
812  result += azi2 == 90 ? 0 : 1;
813  geod_direct(&g, inf, 0, 90, 1000, &lat2, &lon2, &azi2);
814  result += checkNaN(lat2);
815  result += checkNaN(lon2);
816  result += checkNaN(azi2);
817  geod_direct(&g, nan, 0, 90, 1000, &lat2, &lon2, &azi2);
818  result += checkNaN(lat2);
819  result += checkNaN(lon2);
820  result += checkNaN(azi2);
821  return result;
822 }
823 
824 static int Planimeter0() {
825  /* Check fix for pole-encircling bug found 2011-03-16 */
826  double pa[4][2] = {{89, 0}, {89, 90}, {89, 180}, {89, 270}};
827  double pb[4][2] = {{-89, 0}, {-89, 90}, {-89, 180}, {-89, 270}};
828  double pc[4][2] = {{0, -1}, {-1, 0}, {0, 1}, {1, 0}};
829  double pd[3][2] = {{90, 0}, {0, 0}, {0, 90}};
830  struct geod_geodesic g;
831  double perimeter, area;
832  int result = 0;
833  geod_init(&g, wgs84_a, wgs84_f);
834 
835  planimeter(&g, pa, 4, &perimeter, &area);
836  result += checkEquals(perimeter, 631819.8745, 1e-4);
837  result += checkEquals(area, 24952305678.0, 1);
838 
839  planimeter(&g, pb, 4, &perimeter, &area);
840  result += checkEquals(perimeter, 631819.8745, 1e-4);
841  result += checkEquals(area, -24952305678.0, 1);
842 
843  planimeter(&g, pc, 4, &perimeter, &area);
844  result += checkEquals(perimeter, 627598.2731, 1e-4);
845  result += checkEquals(area, 24619419146.0, 1);
846 
847  planimeter(&g, pd, 3, &perimeter, &area);
848  result += checkEquals(perimeter, 30022685, 1);
849  result += checkEquals(area, 63758202715511.0, 1);
850 
851  polylength(&g, pd, 3, &perimeter);
852  result += checkEquals(perimeter, 20020719, 1);
853 
854  return result;
855 }
856 
857 static int Planimeter5() {
858  /* Check fix for Planimeter pole crossing bug found 2011-06-24 */
859  double points[3][2] = {{89, 0.1}, {89, 90.1}, {89, -179.9}};
860  struct geod_geodesic g;
861  double perimeter, area;
862  int result = 0;
863  geod_init(&g, wgs84_a, wgs84_f);
864  planimeter(&g, points, 3, &perimeter, &area);
865  result += checkEquals(perimeter, 539297, 1);
866  result += checkEquals(area, 12476152838.5, 1);
867  return result;
868 }
869 
870 static int Planimeter6() {
871  /* Check fix for Planimeter lon12 rounding bug found 2012-12-03 */
872  double pa[3][2] = {{9, -0.00000000000001}, {9, 180}, {9, 0}};
873  double pb[3][2] = {{9, 0.00000000000001}, {9, 0}, {9, 180}};
874  double pc[3][2] = {{9, 0.00000000000001}, {9, 180}, {9, 0}};
875  double pd[3][2] = {{9, -0.00000000000001}, {9, 0}, {9, 180}};
876  struct geod_geodesic g;
877  double perimeter, area;
878  int result = 0;
879  geod_init(&g, wgs84_a, wgs84_f);
880 
881  planimeter(&g, pa, 3, &perimeter, &area);
882  result += checkEquals(perimeter, 36026861, 1);
883  result += checkEquals(area, 0, 1);
884  planimeter(&g, pb, 3, &perimeter, &area);
885  result += checkEquals(perimeter, 36026861, 1);
886  result += checkEquals(area, 0, 1);
887  planimeter(&g, pc, 3, &perimeter, &area);
888  result += checkEquals(perimeter, 36026861, 1);
889  result += checkEquals(area, 0, 1);
890  planimeter(&g, pd, 3, &perimeter, &area);
891  result += checkEquals(perimeter, 36026861, 1);
892  result += checkEquals(area, 0, 1);
893  return result;
894 }
895 
896 static int Planimeter12() {
897  /* Area of arctic circle (not really -- adjunct to rhumb-area test) */
898  double points[2][2] = {{66.562222222, 0}, {66.562222222, 180}};
899  struct geod_geodesic g;
900  double perimeter, area;
901  int result = 0;
902  geod_init(&g, wgs84_a, wgs84_f);
903  planimeter(&g, points, 2, &perimeter, &area);
904  result += checkEquals(perimeter, 10465729, 1);
905  result += checkEquals(area, 0, 1);
906  return result;
907 }
908 
909 static int Planimeter13() {
910  /* Check encircling pole twice */
911  double points[6][2] = {{89,-360}, {89,-240}, {89,-120},
912  {89,0}, {89,120}, {89,240}};
913  struct geod_geodesic g;
914  double perimeter, area;
915  int result = 0;
916  geod_init(&g, wgs84_a, wgs84_f);
917  planimeter(&g, points, 6, &perimeter, &area);
918  result += checkEquals(perimeter, 1160741, 1);
919  result += checkEquals(area, 32415230256.0, 1);
920  return result;
921 }
922 
923 static int Planimeter15() {
924  /* Coverage tests, includes Planimeter15 - Planimeter18 (combinations of
925  * reverse and sign) + calls to testpoint, testedge, geod_polygonarea. */
926  struct geod_geodesic g;
927  struct geod_polygon p;
928  double lat[] = {2, 1, 3}, lon[] = {1, 2, 3};
929  double area, s12, azi1;
930  double r = 18454562325.45119,
931  a0 = 510065621724088.5093; /* ellipsoid area */
932  int result = 0;
933  geod_init(&g, wgs84_a, wgs84_f);
934  geod_polygon_init(&p, 0);
935  geod_polygon_addpoint(&g, &p, lat[0], lon[0]);
936  geod_polygon_addpoint(&g, &p, lat[1], lon[1]);
937  geod_polygon_testpoint(&g, &p, lat[2], lon[2], 0, 1, &area, nullptr);
938  result += checkEquals(area, r, 0.5);
939  geod_polygon_testpoint(&g, &p, lat[2], lon[2], 0, 0, &area, nullptr);
940  result += checkEquals(area, r, 0.5);
941  geod_polygon_testpoint(&g, &p, lat[2], lon[2], 1, 1, &area, nullptr);
942  result += checkEquals(area, -r, 0.5);
943  geod_polygon_testpoint(&g, &p, lat[2], lon[2], 1, 0, &area, nullptr);
944  result += checkEquals(area, a0-r, 0.5);
945  geod_inverse(&g, lat[1], lon[1], lat[2], lon[2], &s12, &azi1, nullptr);
946  geod_polygon_testedge(&g, &p, azi1, s12, 0, 1, &area, nullptr);
947  result += checkEquals(area, r, 0.5);
948  geod_polygon_testedge(&g, &p, azi1, s12, 0, 0, &area, nullptr);
949  result += checkEquals(area, r, 0.5);
950  geod_polygon_testedge(&g, &p, azi1, s12, 1, 1, &area, nullptr);
951  result += checkEquals(area, -r, 0.5);
952  geod_polygon_testedge(&g, &p, azi1, s12, 1, 0, &area, nullptr);
953  result += checkEquals(area, a0-r, 0.5);
954  geod_polygon_addpoint(&g, &p, lat[2], lon[2]);
955  geod_polygon_compute(&g, &p, 0, 1, &area, nullptr);
956  result += checkEquals(area, r, 0.5);
957  geod_polygon_compute(&g, &p, 0, 0, &area, nullptr);
958  result += checkEquals(area, r, 0.5);
959  geod_polygon_compute(&g, &p, 1, 1, &area, nullptr);
960  result += checkEquals(area, -r, 0.5);
961  geod_polygon_compute(&g, &p, 1, 0, &area, nullptr);
962  result += checkEquals(area, a0-r, 0.5);
963  geod_polygonarea(&g, lat, lon, 3, &area, nullptr);
964  result += checkEquals(area, r, 0.5);
965  return result;
966 }
967 
968 static int Planimeter19() {
969  /* Coverage tests, includes Planimeter19 - Planimeter20 (degenerate
970  * polygons) + extra cases. */
971  struct geod_geodesic g;
972  struct geod_polygon p;
973  double area, perim;
974  int result = 0;
975  geod_init(&g, wgs84_a, wgs84_f);
976  geod_polygon_init(&p, 0);
977  geod_polygon_compute(&g, &p, 0, 1, &area, &perim);
978  result += area == 0 ? 0 : 1;
979  result += perim == 0 ? 0 : 1;
980  geod_polygon_testpoint(&g, &p, 1, 1, 0, 1, &area, &perim);
981  result += area == 0 ? 0 : 1;
982  result += perim == 0 ? 0 : 1;
983  geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, &area, &perim);
984  result += checkNaN(area);
985  result += checkNaN(perim);
986  geod_polygon_addpoint(&g, &p, 1, 1);
987  geod_polygon_compute(&g, &p, 0, 1, &area, &perim);
988  result += area == 0 ? 0 : 1;
989  result += perim == 0 ? 0 : 1;
990  geod_polygon_init(&p, 1);
991  geod_polygon_compute(&g, &p, 0, 1, nullptr, &perim);
992  result += perim == 0 ? 0 : 1;
993  geod_polygon_testpoint(&g, &p, 1, 1, 0, 1, nullptr, &perim);
994  result += perim == 0 ? 0 : 1;
995  geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, nullptr, &perim);
996  result += checkNaN(perim);
997  geod_polygon_addpoint(&g, &p, 1, 1);
998  geod_polygon_compute(&g, &p, 0, 1, nullptr, &perim);
999  result += perim == 0 ? 0 : 1;
1000  geod_polygon_addpoint(&g, &p, 1, 1);
1001  geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, nullptr, &perim);
1002  result += checkEquals(perim, 1000, 1e-10);
1003  geod_polygon_testpoint(&g, &p, 2, 2, 0, 1, nullptr, &perim);
1004  result += checkEquals(perim, 156876.149, 0.5e-3);
1005  return result;
1006 }
1007 
1008 static int Planimeter21() {
1009  /* Some test to add code coverage: multiple circlings of pole (includes
1010  * Planimeter21 - Planimeter28) + invocations via testpoint and testedge. */
1011  struct geod_geodesic g;
1012  struct geod_polygon p;
1013  double area, lat = 45,
1014  a = 39.2144607176828184218, s = 8420705.40957178156285,
1015  r = 39433884866571.4277, /* Area for one circuit */
1016  a0 = 510065621724088.5093; /* Ellipsoid area */
1017  int result = 0, i;
1018  geod_init(&g, wgs84_a, wgs84_f);
1019  geod_polygon_init(&p, 0);
1020  geod_polygon_addpoint(&g, &p, lat, 60);
1021  geod_polygon_addpoint(&g, &p, lat, 180);
1022  geod_polygon_addpoint(&g, &p, lat, -60);
1023  geod_polygon_addpoint(&g, &p, lat, 60);
1024  geod_polygon_addpoint(&g, &p, lat, 180);
1025  geod_polygon_addpoint(&g, &p, lat, -60);
1026  for (i = 3; i <= 4; ++i) {
1027  geod_polygon_addpoint(&g, &p, lat, 60);
1028  geod_polygon_addpoint(&g, &p, lat, 180);
1029  geod_polygon_testpoint(&g, &p, lat, -60, 0, 1, &area, nullptr);
1030  result += checkEquals(area, i*r, 0.5);
1031  geod_polygon_testpoint(&g, &p, lat, -60, 0, 0, &area, nullptr);
1032  result += checkEquals(area, i*r, 0.5);
1033  geod_polygon_testpoint(&g, &p, lat, -60, 1, 1, &area, nullptr);
1034  result += checkEquals(area, -i*r, 0.5);
1035  geod_polygon_testpoint(&g, &p, lat, -60, 1, 0, &area, nullptr);
1036  result += checkEquals(area, -i*r + a0, 0.5);
1037  geod_polygon_testedge(&g, &p, a, s, 0, 1, &area, nullptr);
1038  result += checkEquals(area, i*r, 0.5);
1039  geod_polygon_testedge(&g, &p, a, s, 0, 0, &area, nullptr);
1040  result += checkEquals(area, i*r, 0.5);
1041  geod_polygon_testedge(&g, &p, a, s, 1, 1, &area, nullptr);
1042  result += checkEquals(area, -i*r, 0.5);
1043  geod_polygon_testedge(&g, &p, a, s, 1, 0, &area, nullptr);
1044  result += checkEquals(area, -i*r + a0, 0.5);
1045  geod_polygon_addpoint(&g, &p, lat, -60);
1046  geod_polygon_compute(&g, &p, 0, 1, &area, nullptr);
1047  result += checkEquals(area, i*r, 0.5);
1048  geod_polygon_compute(&g, &p, 0, 0, &area, nullptr);
1049  result += checkEquals(area, i*r, 0.5);
1050  geod_polygon_compute(&g, &p, 1, 1, &area, nullptr);
1051  result += checkEquals(area, -i*r, 0.5);
1052  geod_polygon_compute(&g, &p, 1, 0, &area, nullptr);
1053  result += checkEquals(area, -i*r + a0, 0.5);
1054  }
1055  return result;
1056 }
1057 
1058 static int Planimeter29() {
1059  /* Check fix to transitdirect vs transit zero handling inconsistency */
1060  struct geod_geodesic g;
1061  struct geod_polygon p;
1062  double area;
1063  int result = 0;
1064  geod_init(&g, wgs84_a, wgs84_f);
1065  geod_polygon_init(&p, 0);
1066  geod_polygon_addpoint(&g, &p, 0, 0);
1067  geod_polygon_addedge(&g, &p, 90, 1000);
1068  geod_polygon_addedge(&g, &p, 0, 1000);
1069  geod_polygon_addedge(&g, &p, -90, 1000);
1070  geod_polygon_compute(&g, &p, 0, 1, &area, nullptr);
1071  /* The area should be 1e6. Prior to the fix it was 1e6 - A/2, where
1072  * A = ellipsoid area. */
1073  result += checkEquals(area, 1000000.0, 0.01);
1074  return result;
1075 }
1076 
1077 int main() {
1078  int n = 0, i;
1079  if ((i = testinverse())) {++n; printf("testinverse fail: %d\n", i);}
1080  if ((i = testdirect())) {++n; printf("testdirect fail: %d\n", i);}
1081  if ((i = testarcdirect())) {++n; printf("testarcdirect fail: %d\n", i);}
1082  if ((i = GeodSolve0())) {++n; printf("GeodSolve0 fail: %d\n", i);}
1083  if ((i = GeodSolve1())) {++n; printf("GeodSolve1 fail: %d\n", i);}
1084  if ((i = GeodSolve2())) {++n; printf("GeodSolve2 fail: %d\n", i);}
1085  if ((i = GeodSolve4())) {++n; printf("GeodSolve4 fail: %d\n", i);}
1086  if ((i = GeodSolve5())) {++n; printf("GeodSolve5 fail: %d\n", i);}
1087  if ((i = GeodSolve6())) {++n; printf("GeodSolve6 fail: %d\n", i);}
1088  if ((i = GeodSolve9())) {++n; printf("GeodSolve9 fail: %d\n", i);}
1089  if ((i = GeodSolve10())) {++n; printf("GeodSolve10 fail: %d\n", i);}
1090  if ((i = GeodSolve11())) {++n; printf("GeodSolve11 fail: %d\n", i);}
1091  if ((i = GeodSolve12())) {++n; printf("GeodSolve12 fail: %d\n", i);}
1092  if ((i = GeodSolve14())) {++n; printf("GeodSolve14 fail: %d\n", i);}
1093  if ((i = GeodSolve15())) {++n; printf("GeodSolve15 fail: %d\n", i);}
1094  if ((i = GeodSolve17())) {++n; printf("GeodSolve17 fail: %d\n", i);}
1095  if ((i = GeodSolve26())) {++n; printf("GeodSolve26 fail: %d\n", i);}
1096  if ((i = GeodSolve28())) {++n; printf("GeodSolve28 fail: %d\n", i);}
1097  if ((i = GeodSolve33())) {++n; printf("GeodSolve33 fail: %d\n", i);}
1098  if ((i = GeodSolve55())) {++n; printf("GeodSolve55 fail: %d\n", i);}
1099  if ((i = GeodSolve59())) {++n; printf("GeodSolve59 fail: %d\n", i);}
1100  if ((i = GeodSolve61())) {++n; printf("GeodSolve61 fail: %d\n", i);}
1101  if ((i = GeodSolve65())) {++n; printf("GeodSolve65 fail: %d\n", i);}
1102  if ((i = GeodSolve67())) {++n; printf("GeodSolve67 fail: %d\n", i);}
1103  if ((i = GeodSolve71())) {++n; printf("GeodSolve71 fail: %d\n", i);}
1104  if ((i = GeodSolve73())) {++n; printf("GeodSolve73 fail: %d\n", i);}
1105  if ((i = GeodSolve74())) {++n; printf("GeodSolve74 fail: %d\n", i);}
1106  if ((i = GeodSolve76())) {++n; printf("GeodSolve76 fail: %d\n", i);}
1107  if ((i = GeodSolve78())) {++n; printf("GeodSolve78 fail: %d\n", i);}
1108  if ((i = GeodSolve80())) {++n; printf("GeodSolve80 fail: %d\n", i);}
1109  if ((i = GeodSolve84())) {++n; printf("GeodSolve84 fail: %d\n", i);}
1110  if ((i = Planimeter0())) {++n; printf("Planimeter0 fail: %d\n", i);}
1111  if ((i = Planimeter5())) {++n; printf("Planimeter5 fail: %d\n", i);}
1112  if ((i = Planimeter6())) {++n; printf("Planimeter6 fail: %d\n", i);}
1113  if ((i = Planimeter12())) {++n; printf("Planimeter12 fail: %d\n", i);}
1114  if ((i = Planimeter13())) {++n; printf("Planimeter13 fail: %d\n", i);}
1115  if ((i = Planimeter15())) {++n; printf("Planimeter15 fail: %d\n", i);}
1116  if ((i = Planimeter19())) {++n; printf("Planimeter19 fail: %d\n", i);}
1117  if ((i = Planimeter21())) {++n; printf("Planimeter21 fail: %d\n", i);}
1118  if ((i = Planimeter29())) {++n; printf("Planimeter29 fail: %d\n", i);}
1119  return n;
1120 }
1121 
1122 /** @endcond */
geod_polygon_testedge
unsigned GEOD_DLL geod_polygon_testedge(const struct geod_geodesic *g, const struct geod_polygon *p, double azi, double s, int reverse, int sign, double *pA, double *pP)
geod_polygon_addedge
void GEOD_DLL geod_polygon_addedge(const struct geod_geodesic *g, struct geod_polygon *p, double azi, double s)
geod_inverse
void GEOD_DLL geod_inverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2)
geod_polygon_compute
unsigned GEOD_DLL geod_polygon_compute(const struct geod_geodesic *g, const struct geod_polygon *p, int reverse, int sign, double *pA, double *pP)
geod_polygonarea
void GEOD_DLL geod_polygonarea(const struct geod_geodesic *g, double lats[], double lons[], int n, double *pA, double *pP)
GEOD_LONG_UNROLL
@ GEOD_LONG_UNROLL
Definition: geodesic.h:926
geod_genposition
double GEOD_DLL 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_directline
void GEOD_DLL geod_directline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, unsigned caps)
geodesic.h
API for the geodesic routines in C.
geod_lineinit
void GEOD_DLL geod_lineinit(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned caps)
geod_polygon_init
void GEOD_DLL geod_polygon_init(struct geod_polygon *p, int polylinep)
geod_inverseline
void GEOD_DLL geod_inverseline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, unsigned caps)
GEOD_ARCMODE
@ GEOD_ARCMODE
Definition: geodesic.h:925
geod_geodesic
Definition: geodesic.h:182
geod_geodesicline
Definition: geodesic.h:196
geod_gendirect
double GEOD_DLL 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)
geod_polygon::lon
double lon
Definition: geodesic.h:222
geod_polygon_testpoint
unsigned GEOD_DLL geod_polygon_testpoint(const struct geod_geodesic *g, const struct geod_polygon *p, double lat, double lon, int reverse, int sign, double *pA, double *pP)
geod_polygon::lat
double lat
Definition: geodesic.h:221
geod_geninverse
double GEOD_DLL 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_ALL
@ GEOD_ALL
Definition: geodesic.h:916
geod_direct
void GEOD_DLL geod_direct(const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, double *plat2, double *plon2, double *pazi2)
geod_position
void GEOD_DLL geod_position(const struct geod_geodesicline *l, double s12, double *plat2, double *plon2, double *pazi2)
geod_polygon_addpoint
void GEOD_DLL geod_polygon_addpoint(const struct geod_geodesic *g, struct geod_polygon *p, double lat, double lon)
geod_polygon
Definition: geodesic.h:220
geod_geodesicline::azi1
double azi1
Definition: geodesic.h:199
geod_geodesicline::caps
unsigned caps
Definition: geodesic.h:212
main
int main()
Definition: direct.c:21
GEOD_LATITUDE
@ GEOD_LATITUDE
Definition: geodesic.h:908
geod_init
void GEOD_DLL geod_init(struct geod_geodesic *g, double a, double f)