C library for Geodesics  1.52
geodesic.c
Go to the documentation of this file.
1 /**
2  * \file geodesic.c
3  * \brief Implementation of the geodesic routines in C
4  *
5  * For the full documentation see geodesic.h.
6  **********************************************************************/
7 
8 /** @cond SKIP */
9 
10 /*
11  * This is a C implementation of the geodesic algorithms described in
12  *
13  * C. F. F. Karney,
14  * Algorithms for geodesics,
15  * J. Geodesy <b>87</b>, 43--55 (2013);
16  * https://doi.org/10.1007/s00190-012-0578-z
17  * Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
18  *
19  * See the comments in geodesic.h for documentation.
20  *
21  * Copyright (c) Charles Karney (2012-2021) <charles@karney.com> and licensed
22  * under the MIT/X11 License. For more information, see
23  * https://geographiclib.sourceforge.io/
24  */
25 
26 #include "geodesic.h"
27 #include <math.h>
28 #include <limits.h>
29 #include <float.h>
30 
31 #if !defined(__cplusplus)
32 #define nullptr 0
33 #endif
34 
35 #define GEOGRAPHICLIB_GEODESIC_ORDER 6
36 #define nA1 GEOGRAPHICLIB_GEODESIC_ORDER
37 #define nC1 GEOGRAPHICLIB_GEODESIC_ORDER
38 #define nC1p GEOGRAPHICLIB_GEODESIC_ORDER
39 #define nA2 GEOGRAPHICLIB_GEODESIC_ORDER
40 #define nC2 GEOGRAPHICLIB_GEODESIC_ORDER
41 #define nA3 GEOGRAPHICLIB_GEODESIC_ORDER
42 #define nA3x nA3
43 #define nC3 GEOGRAPHICLIB_GEODESIC_ORDER
44 #define nC3x ((nC3 * (nC3 - 1)) / 2)
45 #define nC4 GEOGRAPHICLIB_GEODESIC_ORDER
46 #define nC4x ((nC4 * (nC4 + 1)) / 2)
47 #define nC (GEOGRAPHICLIB_GEODESIC_ORDER + 1)
48 
49 typedef double real;
50 typedef int boolx;
51 
52 static unsigned init = 0;
53 static const int FALSE = 0;
54 static const int TRUE = 1;
55 static unsigned digits, maxit1, maxit2;
56 static real epsilon, realmin, pi, degree, NaN,
57  tiny, tol0, tol1, tol2, tolb, xthresh;
58 
59 static void Init(void) {
60  if (!init) {
61  digits = DBL_MANT_DIG;
62  epsilon = DBL_EPSILON;
63  realmin = DBL_MIN;
64 #if defined(M_PI)
65  pi = M_PI;
66 #else
67  pi = atan2(0.0, -1.0);
68 #endif
69  maxit1 = 20;
70  maxit2 = maxit1 + digits + 10;
71  tiny = sqrt(realmin);
72  tol0 = epsilon;
73  /* Increase multiplier in defn of tol1 from 100 to 200 to fix inverse case
74  * 52.784459512564 0 -52.784459512563990912 179.634407464943777557
75  * which otherwise failed for Visual Studio 10 (Release and Debug) */
76  tol1 = 200 * tol0;
77  tol2 = sqrt(tol0);
78  /* Check on bisection interval */
79  tolb = tol0 * tol2;
80  xthresh = 1000 * tol2;
81  degree = pi/180;
82  NaN = nan("0");
83  init = 1;
84  }
85 }
86 
87 enum captype {
88  CAP_NONE = 0U,
89  CAP_C1 = 1U<<0,
90  CAP_C1p = 1U<<1,
91  CAP_C2 = 1U<<2,
92  CAP_C3 = 1U<<3,
93  CAP_C4 = 1U<<4,
94  CAP_ALL = 0x1FU,
95  OUT_ALL = 0x7F80U
96 };
97 
98 static real sq(real x) { return x * x; }
99 
100 static real sumx(real u, real v, real* t) {
101  volatile real s = u + v;
102  volatile real up = s - v;
103  volatile real vpp = s - up;
104  up -= u;
105  vpp -= v;
106  if (t) *t = -(up + vpp);
107  /* error-free sum:
108  * u + v = s + t
109  * = round(u + v) + t */
110  return s;
111 }
112 
113 static real polyval(int N, const real p[], real x) {
114  real y = N < 0 ? 0 : *p++;
115  while (--N >= 0) y = y * x + *p++;
116  return y;
117 }
118 
119 /* mimic C++ std::min and std::max */
120 static real minx(real a, real b)
121 { return (b < a) ? b : a; }
122 
123 static real maxx(real a, real b)
124 { return (a < b) ? b : a; }
125 
126 static void swapx(real* x, real* y)
127 { real t = *x; *x = *y; *y = t; }
128 
129 static void norm2(real* sinx, real* cosx) {
130 #if defined(_MSC_VER) && defined(_M_IX86)
131  /* hypot for Visual Studio (A=win32) fails monotonicity, e.g., with
132  * x = 0.6102683302836215
133  * y1 = 0.7906090004346522
134  * y2 = y1 + 1e-16
135  * the test
136  * hypot(x, y2) >= hypot(x, y1)
137  * fails. See also
138  * https://bugs.python.org/issue43088 */
139  real r = sqrt(*sinx * *sinx + *cosx * *cosx);
140 #else
141  real r = hypot(*sinx, *cosx);
142 #endif
143  *sinx /= r;
144  *cosx /= r;
145 }
146 
147 static real AngNormalize(real x) {
148  x = remainder(x, (real)(360));
149  return x != -180 ? x : 180;
150 }
151 
152 static real LatFix(real x)
153 { return fabs(x) > 90 ? NaN : x; }
154 
155 static real AngDiff(real x, real y, real* e) {
156  real t, d = AngNormalize(sumx(AngNormalize(-x), AngNormalize(y), &t));
157  /* Here y - x = d + t (mod 360), exactly, where d is in (-180,180] and
158  * abs(t) <= eps (eps = 2^-45 for doubles). The only case where the
159  * addition of t takes the result outside the range (-180,180] is d = 180
160  * and t > 0. The case, d = -180 + eps, t = -eps, can't happen, since
161  * sum would have returned the exact result in such a case (i.e., given t
162  * = 0). */
163  return sumx(d == 180 && t > 0 ? -180 : d, t, e);
164 }
165 
166 static real AngRound(real x) {
167  const real z = 1/(real)(16);
168  volatile real y;
169  if (x == 0) return 0;
170  y = fabs(x);
171  /* The compiler mustn't "simplify" z - (z - y) to y */
172  y = y < z ? z - (z - y) : y;
173  return x < 0 ? -y : y;
174 }
175 
176 static void sincosdx(real x, real* sinx, real* cosx) {
177  /* In order to minimize round-off errors, this function exactly reduces
178  * the argument to the range [-45, 45] before converting it to radians. */
179  real r, s, c; int q = 0;
180  r = remquo(x, (real)(90), &q);
181  /* now abs(r) <= 45 */
182  r *= degree;
183  /* Possibly could call the gnu extension sincos */
184  s = sin(r); c = cos(r);
185  switch ((unsigned)q & 3U) {
186  case 0U: *sinx = s; *cosx = c; break;
187  case 1U: *sinx = c; *cosx = -s; break;
188  case 2U: *sinx = -s; *cosx = -c; break;
189  default: *sinx = -c; *cosx = s; break; /* case 3U */
190  }
191  if (x != 0) { *sinx += (real)(0); *cosx += (real)(0); }
192 }
193 
194 static real atan2dx(real y, real x) {
195  /* In order to minimize round-off errors, this function rearranges the
196  * arguments so that result of atan2 is in the range [-pi/4, pi/4] before
197  * converting it to degrees and mapping the result to the correct
198  * quadrant. */
199  int q = 0; real ang;
200  if (fabs(y) > fabs(x)) { swapx(&x, &y); q = 2; }
201  if (x < 0) { x = -x; ++q; }
202  /* here x >= 0 and x >= abs(y), so angle is in [-pi/4, pi/4] */
203  ang = atan2(y, x) / degree;
204  switch (q) {
205  /* Note that atan2d(-0.0, 1.0) will return -0. However, we expect that
206  * atan2d will not be called with y = -0. If need be, include
207  *
208  * case 0: ang = 0 + ang; break;
209  */
210  case 1: ang = (y >= 0 ? 180 : -180) - ang; break;
211  case 2: ang = 90 - ang; break;
212  case 3: ang = -90 + ang; break;
213  default: break;
214  }
215  return ang;
216 }
217 
218 static void A3coeff(struct geod_geodesic* g);
219 static void C3coeff(struct geod_geodesic* g);
220 static void C4coeff(struct geod_geodesic* g);
221 static real SinCosSeries(boolx sinp,
222  real sinx, real cosx,
223  const real c[], int n);
224 static void Lengths(const struct geod_geodesic* g,
225  real eps, real sig12,
226  real ssig1, real csig1, real dn1,
227  real ssig2, real csig2, real dn2,
228  real cbet1, real cbet2,
229  real* ps12b, real* pm12b, real* pm0,
230  real* pM12, real* pM21,
231  /* Scratch area of the right size */
232  real Ca[]);
233 static real Astroid(real x, real y);
234 static real InverseStart(const struct geod_geodesic* g,
235  real sbet1, real cbet1, real dn1,
236  real sbet2, real cbet2, real dn2,
237  real lam12, real slam12, real clam12,
238  real* psalp1, real* pcalp1,
239  /* Only updated if return val >= 0 */
240  real* psalp2, real* pcalp2,
241  /* Only updated for short lines */
242  real* pdnm,
243  /* Scratch area of the right size */
244  real Ca[]);
245 static real Lambda12(const struct geod_geodesic* g,
246  real sbet1, real cbet1, real dn1,
247  real sbet2, real cbet2, real dn2,
248  real salp1, real calp1,
249  real slam120, real clam120,
250  real* psalp2, real* pcalp2,
251  real* psig12,
252  real* pssig1, real* pcsig1,
253  real* pssig2, real* pcsig2,
254  real* peps,
255  real* pdomg12,
256  boolx diffp, real* pdlam12,
257  /* Scratch area of the right size */
258  real Ca[]);
259 static real A3f(const struct geod_geodesic* g, real eps);
260 static void C3f(const struct geod_geodesic* g, real eps, real c[]);
261 static void C4f(const struct geod_geodesic* g, real eps, real c[]);
262 static real A1m1f(real eps);
263 static void C1f(real eps, real c[]);
264 static void C1pf(real eps, real c[]);
265 static real A2m1f(real eps);
266 static void C2f(real eps, real c[]);
267 static int transit(real lon1, real lon2);
268 static int transitdirect(real lon1, real lon2);
269 static void accini(real s[]);
270 static void acccopy(const real s[], real t[]);
271 static void accadd(real s[], real y);
272 static real accsum(const real s[], real y);
273 static void accneg(real s[]);
274 static void accrem(real s[], real y);
275 static real areareduceA(real area[], real area0,
276  int crossings, boolx reverse, boolx sign);
277 static real areareduceB(real area, real area0,
278  int crossings, boolx reverse, boolx sign);
279 
280 void geod_init(struct geod_geodesic* g, real a, real f) {
281  if (!init) Init();
282  g->a = a;
283  g->f = f;
284  g->f1 = 1 - g->f;
285  g->e2 = g->f * (2 - g->f);
286  g->ep2 = g->e2 / sq(g->f1); /* e2 / (1 - e2) */
287  g->n = g->f / ( 2 - g->f);
288  g->b = g->a * g->f1;
289  g->c2 = (sq(g->a) + sq(g->b) *
290  (g->e2 == 0 ? 1 :
291  (g->e2 > 0 ? atanh(sqrt(g->e2)) : atan(sqrt(-g->e2))) /
292  sqrt(fabs(g->e2))))/2; /* authalic radius squared */
293  /* The sig12 threshold for "really short". Using the auxiliary sphere
294  * solution with dnm computed at (bet1 + bet2) / 2, the relative error in the
295  * azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2. (Error
296  * measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. For a given f and
297  * sig12, the max error occurs for lines near the pole. If the old rule for
298  * computing dnm = (dn1 + dn2)/2 is used, then the error increases by a
299  * factor of 2.) Setting this equal to epsilon gives sig12 = etol2. Here
300  * 0.1 is a safety factor (error decreased by 100) and max(0.001, abs(f))
301  * stops etol2 getting too large in the nearly spherical case. */
302  g->etol2 = 0.1 * tol2 /
303  sqrt( maxx((real)(0.001), fabs(g->f)) * minx((real)(1), 1 - g->f/2) / 2 );
304 
305  A3coeff(g);
306  C3coeff(g);
307  C4coeff(g);
308 }
309 
310 static void geod_lineinit_int(struct geod_geodesicline* l,
311  const struct geod_geodesic* g,
312  real lat1, real lon1,
313  real azi1, real salp1, real calp1,
314  unsigned caps) {
315  real cbet1, sbet1, eps;
316  l->a = g->a;
317  l->f = g->f;
318  l->b = g->b;
319  l->c2 = g->c2;
320  l->f1 = g->f1;
321  /* If caps is 0 assume the standard direct calculation */
322  l->caps = (caps ? caps : GEOD_DISTANCE_IN | GEOD_LONGITUDE) |
323  /* always allow latitude and azimuth and unrolling of longitude */
325 
326  l->lat1 = LatFix(lat1);
327  l->lon1 = lon1;
328  l->azi1 = azi1;
329  l->salp1 = salp1;
330  l->calp1 = calp1;
331 
332  sincosdx(AngRound(l->lat1), &sbet1, &cbet1); sbet1 *= l->f1;
333  /* Ensure cbet1 = +epsilon at poles */
334  norm2(&sbet1, &cbet1); cbet1 = maxx(tiny, cbet1);
335  l->dn1 = sqrt(1 + g->ep2 * sq(sbet1));
336 
337  /* Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0), */
338  l->salp0 = l->salp1 * cbet1; /* alp0 in [0, pi/2 - |bet1|] */
339  /* Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following
340  * is slightly better (consider the case salp1 = 0). */
341  l->calp0 = hypot(l->calp1, l->salp1 * sbet1);
342  /* Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).
343  * sig = 0 is nearest northward crossing of equator.
344  * With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).
345  * With bet1 = pi/2, alp1 = -pi, sig1 = pi/2
346  * With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2
347  * Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).
348  * With alp0 in (0, pi/2], quadrants for sig and omg coincide.
349  * No atan2(0,0) ambiguity at poles since cbet1 = +epsilon.
350  * With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi. */
351  l->ssig1 = sbet1; l->somg1 = l->salp0 * sbet1;
352  l->csig1 = l->comg1 = sbet1 != 0 || l->calp1 != 0 ? cbet1 * l->calp1 : 1;
353  norm2(&l->ssig1, &l->csig1); /* sig1 in (-pi, pi] */
354  /* norm2(somg1, comg1); -- don't need to normalize! */
355 
356  l->k2 = sq(l->calp0) * g->ep2;
357  eps = l->k2 / (2 * (1 + sqrt(1 + l->k2)) + l->k2);
358 
359  if (l->caps & CAP_C1) {
360  real s, c;
361  l->A1m1 = A1m1f(eps);
362  C1f(eps, l->C1a);
363  l->B11 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C1a, nC1);
364  s = sin(l->B11); c = cos(l->B11);
365  /* tau1 = sig1 + B11 */
366  l->stau1 = l->ssig1 * c + l->csig1 * s;
367  l->ctau1 = l->csig1 * c - l->ssig1 * s;
368  /* Not necessary because C1pa reverts C1a
369  * B11 = -SinCosSeries(TRUE, stau1, ctau1, C1pa, nC1p); */
370  }
371 
372  if (l->caps & CAP_C1p)
373  C1pf(eps, l->C1pa);
374 
375  if (l->caps & CAP_C2) {
376  l->A2m1 = A2m1f(eps);
377  C2f(eps, l->C2a);
378  l->B21 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C2a, nC2);
379  }
380 
381  if (l->caps & CAP_C3) {
382  C3f(g, eps, l->C3a);
383  l->A3c = -l->f * l->salp0 * A3f(g, eps);
384  l->B31 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C3a, nC3-1);
385  }
386 
387  if (l->caps & CAP_C4) {
388  C4f(g, eps, l->C4a);
389  /* Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0) */
390  l->A4 = sq(l->a) * l->calp0 * l->salp0 * g->e2;
391  l->B41 = SinCosSeries(FALSE, l->ssig1, l->csig1, l->C4a, nC4);
392  }
393 
394  l->a13 = l->s13 = NaN;
395 }
396 
397 void geod_lineinit(struct geod_geodesicline* l,
398  const struct geod_geodesic* g,
399  real lat1, real lon1, real azi1, unsigned caps) {
400  real salp1, calp1;
401  azi1 = AngNormalize(azi1);
402  /* Guard against underflow in salp0 */
403  sincosdx(AngRound(azi1), &salp1, &calp1);
404  geod_lineinit_int(l, g, lat1, lon1, azi1, salp1, calp1, caps);
405 }
406 
408  const struct geod_geodesic* g,
409  real lat1, real lon1, real azi1,
410  unsigned flags, real s12_a12,
411  unsigned caps) {
412  geod_lineinit(l, g, lat1, lon1, azi1, caps);
413  geod_gensetdistance(l, flags, s12_a12);
414 }
415 
416 void geod_directline(struct geod_geodesicline* l,
417  const struct geod_geodesic* g,
418  real lat1, real lon1, real azi1,
419  real s12, unsigned caps) {
420  geod_gendirectline(l, g, lat1, lon1, azi1, GEOD_NOFLAGS, s12, caps);
421 }
422 
423 real geod_genposition(const struct geod_geodesicline* l,
424  unsigned flags, real s12_a12,
425  real* plat2, real* plon2, real* pazi2,
426  real* ps12, real* pm12,
427  real* pM12, real* pM21,
428  real* pS12) {
429  real lat2 = 0, lon2 = 0, azi2 = 0, s12 = 0,
430  m12 = 0, M12 = 0, M21 = 0, S12 = 0;
431  /* Avoid warning about uninitialized B12. */
432  real sig12, ssig12, csig12, B12 = 0, AB1 = 0;
433  real omg12, lam12, lon12;
434  real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2, dn2;
435  unsigned outmask =
436  (plat2 ? GEOD_LATITUDE : GEOD_NONE) |
437  (plon2 ? GEOD_LONGITUDE : GEOD_NONE) |
438  (pazi2 ? GEOD_AZIMUTH : GEOD_NONE) |
439  (ps12 ? GEOD_DISTANCE : GEOD_NONE) |
440  (pm12 ? GEOD_REDUCEDLENGTH : GEOD_NONE) |
441  (pM12 || pM21 ? GEOD_GEODESICSCALE : GEOD_NONE) |
442  (pS12 ? GEOD_AREA : GEOD_NONE);
443 
444  outmask &= l->caps & OUT_ALL;
445  if (!( (flags & GEOD_ARCMODE || (l->caps & (GEOD_DISTANCE_IN & OUT_ALL))) ))
446  /* Impossible distance calculation requested */
447  return NaN;
448 
449  if (flags & GEOD_ARCMODE) {
450  /* Interpret s12_a12 as spherical arc length */
451  sig12 = s12_a12 * degree;
452  sincosdx(s12_a12, &ssig12, &csig12);
453  } else {
454  /* Interpret s12_a12 as distance */
455  real
456  tau12 = s12_a12 / (l->b * (1 + l->A1m1)),
457  s = sin(tau12),
458  c = cos(tau12);
459  /* tau2 = tau1 + tau12 */
460  B12 = - SinCosSeries(TRUE,
461  l->stau1 * c + l->ctau1 * s,
462  l->ctau1 * c - l->stau1 * s,
463  l->C1pa, nC1p);
464  sig12 = tau12 - (B12 - l->B11);
465  ssig12 = sin(sig12); csig12 = cos(sig12);
466  if (fabs(l->f) > 0.01) {
467  /* Reverted distance series is inaccurate for |f| > 1/100, so correct
468  * sig12 with 1 Newton iteration. The following table shows the
469  * approximate maximum error for a = WGS_a() and various f relative to
470  * GeodesicExact.
471  * erri = the error in the inverse solution (nm)
472  * errd = the error in the direct solution (series only) (nm)
473  * errda = the error in the direct solution (series + 1 Newton) (nm)
474  *
475  * f erri errd errda
476  * -1/5 12e6 1.2e9 69e6
477  * -1/10 123e3 12e6 765e3
478  * -1/20 1110 108e3 7155
479  * -1/50 18.63 200.9 27.12
480  * -1/100 18.63 23.78 23.37
481  * -1/150 18.63 21.05 20.26
482  * 1/150 22.35 24.73 25.83
483  * 1/100 22.35 25.03 25.31
484  * 1/50 29.80 231.9 30.44
485  * 1/20 5376 146e3 10e3
486  * 1/10 829e3 22e6 1.5e6
487  * 1/5 157e6 3.8e9 280e6 */
488  real serr;
489  ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;
490  csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
491  B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
492  serr = (1 + l->A1m1) * (sig12 + (B12 - l->B11)) - s12_a12 / l->b;
493  sig12 = sig12 - serr / sqrt(1 + l->k2 * sq(ssig2));
494  ssig12 = sin(sig12); csig12 = cos(sig12);
495  /* Update B12 below */
496  }
497  }
498 
499  /* sig2 = sig1 + sig12 */
500  ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;
501  csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
502  dn2 = sqrt(1 + l->k2 * sq(ssig2));
504  if (flags & GEOD_ARCMODE || fabs(l->f) > 0.01)
505  B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
506  AB1 = (1 + l->A1m1) * (B12 - l->B11);
507  }
508  /* sin(bet2) = cos(alp0) * sin(sig2) */
509  sbet2 = l->calp0 * ssig2;
510  /* Alt: cbet2 = hypot(csig2, salp0 * ssig2); */
511  cbet2 = hypot(l->salp0, l->calp0 * csig2);
512  if (cbet2 == 0)
513  /* I.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case */
514  cbet2 = csig2 = tiny;
515  /* tan(alp0) = cos(sig2)*tan(alp2) */
516  salp2 = l->salp0; calp2 = l->calp0 * csig2; /* No need to normalize */
517 
518  if (outmask & GEOD_DISTANCE)
519  s12 = (flags & GEOD_ARCMODE) ?
520  l->b * ((1 + l->A1m1) * sig12 + AB1) :
521  s12_a12;
522 
523  if (outmask & GEOD_LONGITUDE) {
524  real E = copysign(1, l->salp0); /* east or west going? */
525  /* tan(omg2) = sin(alp0) * tan(sig2) */
526  somg2 = l->salp0 * ssig2; comg2 = csig2; /* No need to normalize */
527  /* omg12 = omg2 - omg1 */
528  omg12 = (flags & GEOD_LONG_UNROLL)
529  ? E * (sig12
530  - (atan2( ssig2, csig2) - atan2( l->ssig1, l->csig1))
531  + (atan2(E * somg2, comg2) - atan2(E * l->somg1, l->comg1)))
532  : atan2(somg2 * l->comg1 - comg2 * l->somg1,
533  comg2 * l->comg1 + somg2 * l->somg1);
534  lam12 = omg12 + l->A3c *
535  ( sig12 + (SinCosSeries(TRUE, ssig2, csig2, l->C3a, nC3-1)
536  - l->B31));
537  lon12 = lam12 / degree;
538  lon2 = (flags & GEOD_LONG_UNROLL) ? l->lon1 + lon12 :
539  AngNormalize(AngNormalize(l->lon1) + AngNormalize(lon12));
540  }
541 
542  if (outmask & GEOD_LATITUDE)
543  lat2 = atan2dx(sbet2, l->f1 * cbet2);
544 
545  if (outmask & GEOD_AZIMUTH)
546  azi2 = atan2dx(salp2, calp2);
547 
548  if (outmask & (GEOD_REDUCEDLENGTH | GEOD_GEODESICSCALE)) {
549  real
550  B22 = SinCosSeries(TRUE, ssig2, csig2, l->C2a, nC2),
551  AB2 = (1 + l->A2m1) * (B22 - l->B21),
552  J12 = (l->A1m1 - l->A2m1) * sig12 + (AB1 - AB2);
553  if (outmask & GEOD_REDUCEDLENGTH)
554  /* Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
555  * accurate cancellation in the case of coincident points. */
556  m12 = l->b * ((dn2 * (l->csig1 * ssig2) - l->dn1 * (l->ssig1 * csig2))
557  - l->csig1 * csig2 * J12);
558  if (outmask & GEOD_GEODESICSCALE) {
559  real t = l->k2 * (ssig2 - l->ssig1) * (ssig2 + l->ssig1) /
560  (l->dn1 + dn2);
561  M12 = csig12 + (t * ssig2 - csig2 * J12) * l->ssig1 / l->dn1;
562  M21 = csig12 - (t * l->ssig1 - l->csig1 * J12) * ssig2 / dn2;
563  }
564  }
565 
566  if (outmask & GEOD_AREA) {
567  real
568  B42 = SinCosSeries(FALSE, ssig2, csig2, l->C4a, nC4);
569  real salp12, calp12;
570  if (l->calp0 == 0 || l->salp0 == 0) {
571  /* alp12 = alp2 - alp1, used in atan2 so no need to normalize */
572  salp12 = salp2 * l->calp1 - calp2 * l->salp1;
573  calp12 = calp2 * l->calp1 + salp2 * l->salp1;
574  } else {
575  /* tan(alp) = tan(alp0) * sec(sig)
576  * tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)
577  * = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2)
578  * If csig12 > 0, write
579  * csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
580  * else
581  * csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
582  * No need to normalize */
583  salp12 = l->calp0 * l->salp0 *
584  (csig12 <= 0 ? l->csig1 * (1 - csig12) + ssig12 * l->ssig1 :
585  ssig12 * (l->csig1 * ssig12 / (1 + csig12) + l->ssig1));
586  calp12 = sq(l->salp0) + sq(l->calp0) * l->csig1 * csig2;
587  }
588  S12 = l->c2 * atan2(salp12, calp12) + l->A4 * (B42 - l->B41);
589  }
590 
591  /* In the pattern
592  *
593  * if ((outmask & GEOD_XX) && pYY)
594  * *pYY = YY;
595  *
596  * the second check "&& pYY" is redundant. It's there to make the CLang
597  * static analyzer happy.
598  */
599  if ((outmask & GEOD_LATITUDE) && plat2)
600  *plat2 = lat2;
601  if ((outmask & GEOD_LONGITUDE) && plon2)
602  *plon2 = lon2;
603  if ((outmask & GEOD_AZIMUTH) && pazi2)
604  *pazi2 = azi2;
605  if ((outmask & GEOD_DISTANCE) && ps12)
606  *ps12 = s12;
607  if ((outmask & GEOD_REDUCEDLENGTH) && pm12)
608  *pm12 = m12;
609  if (outmask & GEOD_GEODESICSCALE) {
610  if (pM12) *pM12 = M12;
611  if (pM21) *pM21 = M21;
612  }
613  if ((outmask & GEOD_AREA) && pS12)
614  *pS12 = S12;
615 
616  return (flags & GEOD_ARCMODE) ? s12_a12 : sig12 / degree;
617 }
618 
619 void geod_setdistance(struct geod_geodesicline* l, real s13) {
620  l->s13 = s13;
621  l->a13 = geod_genposition(l, GEOD_NOFLAGS, l->s13, nullptr, nullptr, nullptr,
622  nullptr, nullptr, nullptr, nullptr, nullptr);
623 }
624 
625 static void geod_setarc(struct geod_geodesicline* l, real a13) {
626  l->a13 = a13; l->s13 = NaN;
627  geod_genposition(l, GEOD_ARCMODE, l->a13, nullptr, nullptr, nullptr, &l->s13,
628  nullptr, nullptr, nullptr, nullptr);
629 }
630 
632  unsigned flags, real s13_a13) {
633  (flags & GEOD_ARCMODE) ?
634  geod_setarc(l, s13_a13) :
635  geod_setdistance(l, s13_a13);
636 }
637 
638 void geod_position(const struct geod_geodesicline* l, real s12,
639  real* plat2, real* plon2, real* pazi2) {
640  geod_genposition(l, FALSE, s12, plat2, plon2, pazi2,
641  nullptr, nullptr, nullptr, nullptr, nullptr);
642 }
643 
644 real geod_gendirect(const struct geod_geodesic* g,
645  real lat1, real lon1, real azi1,
646  unsigned flags, real s12_a12,
647  real* plat2, real* plon2, real* pazi2,
648  real* ps12, real* pm12, real* pM12, real* pM21,
649  real* pS12) {
650  struct geod_geodesicline l;
651  unsigned outmask =
652  (plat2 ? GEOD_LATITUDE : GEOD_NONE) |
653  (plon2 ? GEOD_LONGITUDE : GEOD_NONE) |
654  (pazi2 ? GEOD_AZIMUTH : GEOD_NONE) |
655  (ps12 ? GEOD_DISTANCE : GEOD_NONE) |
656  (pm12 ? GEOD_REDUCEDLENGTH : GEOD_NONE) |
657  (pM12 || pM21 ? GEOD_GEODESICSCALE : GEOD_NONE) |
658  (pS12 ? GEOD_AREA : GEOD_NONE);
659 
660  geod_lineinit(&l, g, lat1, lon1, azi1,
661  /* Automatically supply GEOD_DISTANCE_IN if necessary */
662  outmask |
663  ((flags & GEOD_ARCMODE) ? GEOD_NONE : GEOD_DISTANCE_IN));
664  return geod_genposition(&l, flags, s12_a12,
665  plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);
666 }
667 
668 void geod_direct(const struct geod_geodesic* g,
670  real s12,
671  real* plat2, real* plon2, real* pazi2) {
672  geod_gendirect(g, lat1, lon1, azi1, GEOD_NOFLAGS, s12, plat2, plon2, pazi2,
673  nullptr, nullptr, nullptr, nullptr, nullptr);
674 }
675 
676 static real geod_geninverse_int(const struct geod_geodesic* g,
677  real lat1, real lon1, real lat2, real lon2,
678  real* ps12,
679  real* psalp1, real* pcalp1,
680  real* psalp2, real* pcalp2,
681  real* pm12, real* pM12, real* pM21,
682  real* pS12) {
683  real s12 = 0, m12 = 0, M12 = 0, M21 = 0, S12 = 0;
684  real lon12, lon12s;
685  int latsign, lonsign, swapp;
686  real sbet1, cbet1, sbet2, cbet2, s12x = 0, m12x = 0;
687  real dn1, dn2, lam12, slam12, clam12;
688  real a12 = 0, sig12, calp1 = 0, salp1 = 0, calp2 = 0, salp2 = 0;
689  real Ca[nC];
690  boolx meridian;
691  /* somg12 > 1 marks that it needs to be calculated */
692  real omg12 = 0, somg12 = 2, comg12 = 0;
693 
694  unsigned outmask =
695  (ps12 ? GEOD_DISTANCE : GEOD_NONE) |
696  (pm12 ? GEOD_REDUCEDLENGTH : GEOD_NONE) |
697  (pM12 || pM21 ? GEOD_GEODESICSCALE : GEOD_NONE) |
698  (pS12 ? GEOD_AREA : GEOD_NONE);
699 
700  outmask &= OUT_ALL;
701  /* Compute longitude difference (AngDiff does this carefully). Result is
702  * in [-180, 180] but -180 is only for west-going geodesics. 180 is for
703  * east-going and meridional geodesics. */
704  lon12 = AngDiff(lon1, lon2, &lon12s);
705  /* Make longitude difference positive. */
706  lonsign = lon12 >= 0 ? 1 : -1;
707  /* If very close to being on the same half-meridian, then make it so. */
708  lon12 = lonsign * AngRound(lon12);
709  lon12s = AngRound((180 - lon12) - lonsign * lon12s);
710  lam12 = lon12 * degree;
711  if (lon12 > 90) {
712  sincosdx(lon12s, &slam12, &clam12);
713  clam12 = -clam12;
714  } else
715  sincosdx(lon12, &slam12, &clam12);
716 
717  /* If really close to the equator, treat as on equator. */
718  lat1 = AngRound(LatFix(lat1));
719  lat2 = AngRound(LatFix(lat2));
720  /* Swap points so that point with higher (abs) latitude is point 1
721  * If one latitude is a nan, then it becomes lat1. */
722  swapp = fabs(lat1) < fabs(lat2) ? -1 : 1;
723  if (swapp < 0) {
724  lonsign *= -1;
725  swapx(&lat1, &lat2);
726  }
727  /* Make lat1 <= 0 */
728  latsign = lat1 < 0 ? 1 : -1;
729  lat1 *= latsign;
730  lat2 *= latsign;
731  /* Now we have
732  *
733  * 0 <= lon12 <= 180
734  * -90 <= lat1 <= 0
735  * lat1 <= lat2 <= -lat1
736  *
737  * longsign, swapp, latsign register the transformation to bring the
738  * coordinates to this canonical form. In all cases, 1 means no change was
739  * made. We make these transformations so that there are few cases to
740  * check, e.g., on verifying quadrants in atan2. In addition, this
741  * enforces some symmetries in the results returned. */
742 
743  sincosdx(lat1, &sbet1, &cbet1); sbet1 *= g->f1;
744  /* Ensure cbet1 = +epsilon at poles */
745  norm2(&sbet1, &cbet1); cbet1 = maxx(tiny, cbet1);
746 
747  sincosdx(lat2, &sbet2, &cbet2); sbet2 *= g->f1;
748  /* Ensure cbet2 = +epsilon at poles */
749  norm2(&sbet2, &cbet2); cbet2 = maxx(tiny, cbet2);
750 
751  /* If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the
752  * |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is
753  * a better measure. This logic is used in assigning calp2 in Lambda12.
754  * Sometimes these quantities vanish and in that case we force bet2 = +/-
755  * bet1 exactly. An example where is is necessary is the inverse problem
756  * 48.522876735459 0 -48.52287673545898293 179.599720456223079643
757  * which failed with Visual Studio 10 (Release and Debug) */
758 
759  if (cbet1 < -sbet1) {
760  if (cbet2 == cbet1)
761  sbet2 = sbet2 < 0 ? sbet1 : -sbet1;
762  } else {
763  if (fabs(sbet2) == -sbet1)
764  cbet2 = cbet1;
765  }
766 
767  dn1 = sqrt(1 + g->ep2 * sq(sbet1));
768  dn2 = sqrt(1 + g->ep2 * sq(sbet2));
769 
770  meridian = lat1 == -90 || slam12 == 0;
771 
772  if (meridian) {
773 
774  /* Endpoints are on a single full meridian, so the geodesic might lie on
775  * a meridian. */
776 
777  real ssig1, csig1, ssig2, csig2;
778  calp1 = clam12; salp1 = slam12; /* Head to the target longitude */
779  calp2 = 1; salp2 = 0; /* At the target we're heading north */
780 
781  /* tan(bet) = tan(sig) * cos(alp) */
782  ssig1 = sbet1; csig1 = calp1 * cbet1;
783  ssig2 = sbet2; csig2 = calp2 * cbet2;
784 
785  /* sig12 = sig2 - sig1 */
786  sig12 = atan2(maxx((real)(0), csig1 * ssig2 - ssig1 * csig2),
787  csig1 * csig2 + ssig1 * ssig2);
788  Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
789  cbet1, cbet2, &s12x, &m12x, nullptr,
790  (outmask & GEOD_GEODESICSCALE) ? &M12 : nullptr,
791  (outmask & GEOD_GEODESICSCALE) ? &M21 : nullptr,
792  Ca);
793  /* Add the check for sig12 since zero length geodesics might yield m12 <
794  * 0. Test case was
795  *
796  * echo 20.001 0 20.001 0 | GeodSolve -i
797  *
798  * In fact, we will have sig12 > pi/2 for meridional geodesic which is
799  * not a shortest path. */
800  if (sig12 < 1 || m12x >= 0) {
801  /* Need at least 2, to handle 90 0 90 180 */
802  if (sig12 < 3 * tiny ||
803  // Prevent negative s12 or m12 for short lines
804  (sig12 < tol0 && (s12x < 0 || m12x < 0)))
805  sig12 = m12x = s12x = 0;
806  m12x *= g->b;
807  s12x *= g->b;
808  a12 = sig12 / degree;
809  } else
810  /* m12 < 0, i.e., prolate and too close to anti-podal */
811  meridian = FALSE;
812  }
813 
814  if (!meridian &&
815  sbet1 == 0 && /* and sbet2 == 0 */
816  /* Mimic the way Lambda12 works with calp1 = 0 */
817  (g->f <= 0 || lon12s >= g->f * 180)) {
818 
819  /* Geodesic runs along equator */
820  calp1 = calp2 = 0; salp1 = salp2 = 1;
821  s12x = g->a * lam12;
822  sig12 = omg12 = lam12 / g->f1;
823  m12x = g->b * sin(sig12);
824  if (outmask & GEOD_GEODESICSCALE)
825  M12 = M21 = cos(sig12);
826  a12 = lon12 / g->f1;
827 
828  } else if (!meridian) {
829 
830  /* Now point1 and point2 belong within a hemisphere bounded by a
831  * meridian and geodesic is neither meridional or equatorial. */
832 
833  /* Figure a starting point for Newton's method */
834  real dnm = 0;
835  sig12 = InverseStart(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
836  lam12, slam12, clam12,
837  &salp1, &calp1, &salp2, &calp2, &dnm,
838  Ca);
839 
840  if (sig12 >= 0) {
841  /* Short lines (InverseStart sets salp2, calp2, dnm) */
842  s12x = sig12 * g->b * dnm;
843  m12x = sq(dnm) * g->b * sin(sig12 / dnm);
844  if (outmask & GEOD_GEODESICSCALE)
845  M12 = M21 = cos(sig12 / dnm);
846  a12 = sig12 / degree;
847  omg12 = lam12 / (g->f1 * dnm);
848  } else {
849 
850  /* Newton's method. This is a straightforward solution of f(alp1) =
851  * lambda12(alp1) - lam12 = 0 with one wrinkle. f(alp) has exactly one
852  * root in the interval (0, pi) and its derivative is positive at the
853  * root. Thus f(alp) is positive for alp > alp1 and negative for alp <
854  * alp1. During the course of the iteration, a range (alp1a, alp1b) is
855  * maintained which brackets the root and with each evaluation of
856  * f(alp) the range is shrunk, if possible. Newton's method is
857  * restarted whenever the derivative of f is negative (because the new
858  * value of alp1 is then further from the solution) or if the new
859  * estimate of alp1 lies outside (0,pi); in this case, the new starting
860  * guess is taken to be (alp1a + alp1b) / 2. */
861  real ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0, domg12 = 0;
862  unsigned numit = 0;
863  /* Bracketing range */
864  real salp1a = tiny, calp1a = 1, salp1b = tiny, calp1b = -1;
865  boolx tripn = FALSE;
866  boolx tripb = FALSE;
867  for (; numit < maxit2; ++numit) {
868  /* the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
869  * WGS84 and random input: mean = 2.85, sd = 0.60 */
870  real dv = 0,
871  v = Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
872  slam12, clam12,
873  &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2,
874  &eps, &domg12, numit < maxit1, &dv, Ca);
875  /* Reversed test to allow escape with NaNs */
876  if (tripb || !(fabs(v) >= (tripn ? 8 : 1) * tol0)) break;
877  /* Update bracketing values */
878  if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b))
879  { salp1b = salp1; calp1b = calp1; }
880  else if (v < 0 && (numit > maxit1 || calp1/salp1 < calp1a/salp1a))
881  { salp1a = salp1; calp1a = calp1; }
882  if (numit < maxit1 && dv > 0) {
883  real
884  dalp1 = -v/dv;
885  real
886  sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
887  nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
888  if (nsalp1 > 0 && fabs(dalp1) < pi) {
889  calp1 = calp1 * cdalp1 - salp1 * sdalp1;
890  salp1 = nsalp1;
891  norm2(&salp1, &calp1);
892  /* In some regimes we don't get quadratic convergence because
893  * slope -> 0. So use convergence conditions based on epsilon
894  * instead of sqrt(epsilon). */
895  tripn = fabs(v) <= 16 * tol0;
896  continue;
897  }
898  }
899  /* Either dv was not positive or updated value was outside legal
900  * range. Use the midpoint of the bracket as the next estimate.
901  * This mechanism is not needed for the WGS84 ellipsoid, but it does
902  * catch problems with more eccentric ellipsoids. Its efficacy is
903  * such for the WGS84 test set with the starting guess set to alp1 =
904  * 90deg:
905  * the WGS84 test set: mean = 5.21, sd = 3.93, max = 24
906  * WGS84 and random input: mean = 4.74, sd = 0.99 */
907  salp1 = (salp1a + salp1b)/2;
908  calp1 = (calp1a + calp1b)/2;
909  norm2(&salp1, &calp1);
910  tripn = FALSE;
911  tripb = (fabs(salp1a - salp1) + (calp1a - calp1) < tolb ||
912  fabs(salp1 - salp1b) + (calp1 - calp1b) < tolb);
913  }
914  Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
915  cbet1, cbet2, &s12x, &m12x, nullptr,
916  (outmask & GEOD_GEODESICSCALE) ? &M12 : nullptr,
917  (outmask & GEOD_GEODESICSCALE) ? &M21 : nullptr, Ca);
918  m12x *= g->b;
919  s12x *= g->b;
920  a12 = sig12 / degree;
921  if (outmask & GEOD_AREA) {
922  /* omg12 = lam12 - domg12 */
923  real sdomg12 = sin(domg12), cdomg12 = cos(domg12);
924  somg12 = slam12 * cdomg12 - clam12 * sdomg12;
925  comg12 = clam12 * cdomg12 + slam12 * sdomg12;
926  }
927  }
928  }
929 
930  if (outmask & GEOD_DISTANCE)
931  s12 = 0 + s12x; /* Convert -0 to 0 */
932 
933  if (outmask & GEOD_REDUCEDLENGTH)
934  m12 = 0 + m12x; /* Convert -0 to 0 */
935 
936  if (outmask & GEOD_AREA) {
937  real
938  /* From Lambda12: sin(alp1) * cos(bet1) = sin(alp0) */
939  salp0 = salp1 * cbet1,
940  calp0 = hypot(calp1, salp1 * sbet1); /* calp0 > 0 */
941  real alp12;
942  if (calp0 != 0 && salp0 != 0) {
943  real
944  /* From Lambda12: tan(bet) = tan(sig) * cos(alp) */
945  ssig1 = sbet1, csig1 = calp1 * cbet1,
946  ssig2 = sbet2, csig2 = calp2 * cbet2,
947  k2 = sq(calp0) * g->ep2,
948  eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2),
949  /* Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0). */
950  A4 = sq(g->a) * calp0 * salp0 * g->e2;
951  real B41, B42;
952  norm2(&ssig1, &csig1);
953  norm2(&ssig2, &csig2);
954  C4f(g, eps, Ca);
955  B41 = SinCosSeries(FALSE, ssig1, csig1, Ca, nC4);
956  B42 = SinCosSeries(FALSE, ssig2, csig2, Ca, nC4);
957  S12 = A4 * (B42 - B41);
958  } else
959  /* Avoid problems with indeterminate sig1, sig2 on equator */
960  S12 = 0;
961 
962  if (!meridian && somg12 > 1) {
963  somg12 = sin(omg12); comg12 = cos(omg12);
964  }
965 
966  if (!meridian &&
967  /* omg12 < 3/4 * pi */
968  comg12 > -(real)(0.7071) && /* Long difference not too big */
969  sbet2 - sbet1 < (real)(1.75)) { /* Lat difference not too big */
970  /* Use tan(Gamma/2) = tan(omg12/2)
971  * * (tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2))
972  * with tan(x/2) = sin(x)/(1+cos(x)) */
973  real
974  domg12 = 1 + comg12, dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;
975  alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
976  domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );
977  } else {
978  /* alp12 = alp2 - alp1, used in atan2 so no need to normalize */
979  real
980  salp12 = salp2 * calp1 - calp2 * salp1,
981  calp12 = calp2 * calp1 + salp2 * salp1;
982  /* The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz
983  * salp12 = -0 and alp12 = -180. However this depends on the sign
984  * being attached to 0 correctly. The following ensures the correct
985  * behavior. */
986  if (salp12 == 0 && calp12 < 0) {
987  salp12 = tiny * calp1;
988  calp12 = -1;
989  }
990  alp12 = atan2(salp12, calp12);
991  }
992  S12 += g->c2 * alp12;
993  S12 *= swapp * lonsign * latsign;
994  /* Convert -0 to 0 */
995  S12 += 0;
996  }
997 
998  /* Convert calp, salp to azimuth accounting for lonsign, swapp, latsign. */
999  if (swapp < 0) {
1000  swapx(&salp1, &salp2);
1001  swapx(&calp1, &calp2);
1002  if (outmask & GEOD_GEODESICSCALE)
1003  swapx(&M12, &M21);
1004  }
1005 
1006  salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
1007  salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
1008 
1009  if (psalp1) *psalp1 = salp1;
1010  if (pcalp1) *pcalp1 = calp1;
1011  if (psalp2) *psalp2 = salp2;
1012  if (pcalp2) *pcalp2 = calp2;
1013 
1014  if (outmask & GEOD_DISTANCE)
1015  *ps12 = s12;
1016  if (outmask & GEOD_REDUCEDLENGTH)
1017  *pm12 = m12;
1018  if (outmask & GEOD_GEODESICSCALE) {
1019  if (pM12) *pM12 = M12;
1020  if (pM21) *pM21 = M21;
1021  }
1022  if (outmask & GEOD_AREA)
1023  *pS12 = S12;
1024 
1025  /* Returned value in [0, 180] */
1026  return a12;
1027 }
1028 
1029 real geod_geninverse(const struct geod_geodesic* g,
1030  real lat1, real lon1, real lat2, real lon2,
1031  real* ps12, real* pazi1, real* pazi2,
1032  real* pm12, real* pM12, real* pM21, real* pS12) {
1033  real salp1, calp1, salp2, calp2,
1034  a12 = geod_geninverse_int(g, lat1, lon1, lat2, lon2, ps12,
1035  &salp1, &calp1, &salp2, &calp2,
1036  pm12, pM12, pM21, pS12);
1037  if (pazi1) *pazi1 = atan2dx(salp1, calp1);
1038  if (pazi2) *pazi2 = atan2dx(salp2, calp2);
1039  return a12;
1040 }
1041 
1042 void geod_inverseline(struct geod_geodesicline* l,
1043  const struct geod_geodesic* g,
1044  real lat1, real lon1, real lat2, real lon2,
1045  unsigned caps) {
1046  real salp1, calp1,
1047  a12 = geod_geninverse_int(g, lat1, lon1, lat2, lon2, nullptr,
1048  &salp1, &calp1, nullptr, nullptr,
1049  nullptr, nullptr, nullptr, nullptr),
1050  azi1 = atan2dx(salp1, calp1);
1052  /* Ensure that a12 can be converted to a distance */
1053  if (caps & (OUT_ALL & GEOD_DISTANCE_IN)) caps |= GEOD_DISTANCE;
1054  geod_lineinit_int(l, g, lat1, lon1, azi1, salp1, calp1, caps);
1055  geod_setarc(l, a12);
1056 }
1057 
1058 void geod_inverse(const struct geod_geodesic* g,
1059  real lat1, real lon1, real lat2, real lon2,
1060  real* ps12, real* pazi1, real* pazi2) {
1061  geod_geninverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2,
1062  nullptr, nullptr, nullptr, nullptr);
1063 }
1064 
1065 real SinCosSeries(boolx sinp, real sinx, real cosx, const real c[], int n) {
1066  /* Evaluate
1067  * y = sinp ? sum(c[i] * sin( 2*i * x), i, 1, n) :
1068  * sum(c[i] * cos((2*i+1) * x), i, 0, n-1)
1069  * using Clenshaw summation. N.B. c[0] is unused for sin series
1070  * Approx operation count = (n + 5) mult and (2 * n + 2) add */
1071  real ar, y0, y1;
1072  c += (n + sinp); /* Point to one beyond last element */
1073  ar = 2 * (cosx - sinx) * (cosx + sinx); /* 2 * cos(2 * x) */
1074  y0 = (n & 1) ? *--c : 0; y1 = 0; /* accumulators for sum */
1075  /* Now n is even */
1076  n /= 2;
1077  while (n--) {
1078  /* Unroll loop x 2, so accumulators return to their original role */
1079  y1 = ar * y0 - y1 + *--c;
1080  y0 = ar * y1 - y0 + *--c;
1081  }
1082  return sinp
1083  ? 2 * sinx * cosx * y0 /* sin(2 * x) * y0 */
1084  : cosx * (y0 - y1); /* cos(x) * (y0 - y1) */
1085 }
1086 
1087 void Lengths(const struct geod_geodesic* g,
1088  real eps, real sig12,
1089  real ssig1, real csig1, real dn1,
1090  real ssig2, real csig2, real dn2,
1091  real cbet1, real cbet2,
1092  real* ps12b, real* pm12b, real* pm0,
1093  real* pM12, real* pM21,
1094  /* Scratch area of the right size */
1095  real Ca[]) {
1096  real m0 = 0, J12 = 0, A1 = 0, A2 = 0;
1097  real Cb[nC];
1098 
1099  /* Return m12b = (reduced length)/b; also calculate s12b = distance/b,
1100  * and m0 = coefficient of secular term in expression for reduced length. */
1101  boolx redlp = pm12b || pm0 || pM12 || pM21;
1102  if (ps12b || redlp) {
1103  A1 = A1m1f(eps);
1104  C1f(eps, Ca);
1105  if (redlp) {
1106  A2 = A2m1f(eps);
1107  C2f(eps, Cb);
1108  m0 = A1 - A2;
1109  A2 = 1 + A2;
1110  }
1111  A1 = 1 + A1;
1112  }
1113  if (ps12b) {
1114  real B1 = SinCosSeries(TRUE, ssig2, csig2, Ca, nC1) -
1115  SinCosSeries(TRUE, ssig1, csig1, Ca, nC1);
1116  /* Missing a factor of b */
1117  *ps12b = A1 * (sig12 + B1);
1118  if (redlp) {
1119  real B2 = SinCosSeries(TRUE, ssig2, csig2, Cb, nC2) -
1120  SinCosSeries(TRUE, ssig1, csig1, Cb, nC2);
1121  J12 = m0 * sig12 + (A1 * B1 - A2 * B2);
1122  }
1123  } else if (redlp) {
1124  /* Assume here that nC1 >= nC2 */
1125  int l;
1126  for (l = 1; l <= nC2; ++l)
1127  Cb[l] = A1 * Ca[l] - A2 * Cb[l];
1128  J12 = m0 * sig12 + (SinCosSeries(TRUE, ssig2, csig2, Cb, nC2) -
1129  SinCosSeries(TRUE, ssig1, csig1, Cb, nC2));
1130  }
1131  if (pm0) *pm0 = m0;
1132  if (pm12b)
1133  /* Missing a factor of b.
1134  * Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
1135  * accurate cancellation in the case of coincident points. */
1136  *pm12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1137  csig1 * csig2 * J12;
1138  if (pM12 || pM21) {
1139  real csig12 = csig1 * csig2 + ssig1 * ssig2;
1140  real t = g->ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);
1141  if (pM12)
1142  *pM12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;
1143  if (pM21)
1144  *pM21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;
1145  }
1146 }
1147 
1148 real Astroid(real x, real y) {
1149  /* Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k.
1150  * This solution is adapted from Geocentric::Reverse. */
1151  real k;
1152  real
1153  p = sq(x),
1154  q = sq(y),
1155  r = (p + q - 1) / 6;
1156  if ( !(q == 0 && r <= 0) ) {
1157  real
1158  /* Avoid possible division by zero when r = 0 by multiplying equations
1159  * for s and t by r^3 and r, resp. */
1160  S = p * q / 4, /* S = r^3 * s */
1161  r2 = sq(r),
1162  r3 = r * r2,
1163  /* The discriminant of the quadratic equation for T3. This is zero on
1164  * the evolute curve p^(1/3)+q^(1/3) = 1 */
1165  disc = S * (S + 2 * r3);
1166  real u = r;
1167  real v, uv, w;
1168  if (disc >= 0) {
1169  real T3 = S + r3, T;
1170  /* Pick the sign on the sqrt to maximize abs(T3). This minimizes loss
1171  * of precision due to cancellation. The result is unchanged because
1172  * of the way the T is used in definition of u. */
1173  T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc); /* T3 = (r * t)^3 */
1174  /* N.B. cbrt always returns the real root. cbrt(-8) = -2. */
1175  T = cbrt(T3); /* T = r * t */
1176  /* T can be zero; but then r2 / T -> 0. */
1177  u += T + (T != 0 ? r2 / T : 0);
1178  } else {
1179  /* T is complex, but the way u is defined the result is real. */
1180  real ang = atan2(sqrt(-disc), -(S + r3));
1181  /* There are three possible cube roots. We choose the root which
1182  * avoids cancellation. Note that disc < 0 implies that r < 0. */
1183  u += 2 * r * cos(ang / 3);
1184  }
1185  v = sqrt(sq(u) + q); /* guaranteed positive */
1186  /* Avoid loss of accuracy when u < 0. */
1187  uv = u < 0 ? q / (v - u) : u + v; /* u+v, guaranteed positive */
1188  w = (uv - q) / (2 * v); /* positive? */
1189  /* Rearrange expression for k to avoid loss of accuracy due to
1190  * subtraction. Division by 0 not possible because uv > 0, w >= 0. */
1191  k = uv / (sqrt(uv + sq(w)) + w); /* guaranteed positive */
1192  } else { /* q == 0 && r <= 0 */
1193  /* y = 0 with |x| <= 1. Handle this case directly.
1194  * for y small, positive root is k = abs(y)/sqrt(1-x^2) */
1195  k = 0;
1196  }
1197  return k;
1198 }
1199 
1200 real InverseStart(const struct geod_geodesic* g,
1201  real sbet1, real cbet1, real dn1,
1202  real sbet2, real cbet2, real dn2,
1203  real lam12, real slam12, real clam12,
1204  real* psalp1, real* pcalp1,
1205  /* Only updated if return val >= 0 */
1206  real* psalp2, real* pcalp2,
1207  /* Only updated for short lines */
1208  real* pdnm,
1209  /* Scratch area of the right size */
1210  real Ca[]) {
1211  real salp1 = 0, calp1 = 0, salp2 = 0, calp2 = 0, dnm = 0;
1212 
1213  /* Return a starting point for Newton's method in salp1 and calp1 (function
1214  * value is -1). If Newton's method doesn't need to be used, return also
1215  * salp2 and calp2 and function value is sig12. */
1216  real
1217  sig12 = -1, /* Return value */
1218  /* bet12 = bet2 - bet1 in [0, pi); bet12a = bet2 + bet1 in (-pi, 0] */
1219  sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
1220  cbet12 = cbet2 * cbet1 + sbet2 * sbet1;
1221  real sbet12a;
1222  boolx shortline = cbet12 >= 0 && sbet12 < (real)(0.5) &&
1223  cbet2 * lam12 < (real)(0.5);
1224  real somg12, comg12, ssig12, csig12;
1225  sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
1226  if (shortline) {
1227  real sbetm2 = sq(sbet1 + sbet2), omg12;
1228  /* sin((bet1+bet2)/2)^2
1229  * = (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2) */
1230  sbetm2 /= sbetm2 + sq(cbet1 + cbet2);
1231  dnm = sqrt(1 + g->ep2 * sbetm2);
1232  omg12 = lam12 / (g->f1 * dnm);
1233  somg12 = sin(omg12); comg12 = cos(omg12);
1234  } else {
1235  somg12 = slam12; comg12 = clam12;
1236  }
1237 
1238  salp1 = cbet2 * somg12;
1239  calp1 = comg12 >= 0 ?
1240  sbet12 + cbet2 * sbet1 * sq(somg12) / (1 + comg12) :
1241  sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
1242 
1243  ssig12 = hypot(salp1, calp1);
1244  csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
1245 
1246  if (shortline && ssig12 < g->etol2) {
1247  /* really short lines */
1248  salp2 = cbet1 * somg12;
1249  calp2 = sbet12 - cbet1 * sbet2 *
1250  (comg12 >= 0 ? sq(somg12) / (1 + comg12) : 1 - comg12);
1251  norm2(&salp2, &calp2);
1252  /* Set return value */
1253  sig12 = atan2(ssig12, csig12);
1254  } else if (fabs(g->n) > (real)(0.1) || /* No astroid calc if too eccentric */
1255  csig12 >= 0 ||
1256  ssig12 >= 6 * fabs(g->n) * pi * sq(cbet1)) {
1257  /* Nothing to do, zeroth order spherical approximation is OK */
1258  } else {
1259  /* Scale lam12 and bet2 to x, y coordinate system where antipodal point
1260  * is at origin and singular point is at y = 0, x = -1. */
1261  real y, lamscale, betscale;
1262  /* Volatile declaration needed to fix inverse case
1263  * 56.320923501171 0 -56.320923501171 179.664747671772880215
1264  * which otherwise fails with g++ 4.4.4 x86 -O3 */
1265  volatile real x;
1266  real lam12x = atan2(-slam12, -clam12); /* lam12 - pi */
1267  if (g->f >= 0) { /* In fact f == 0 does not get here */
1268  /* x = dlong, y = dlat */
1269  {
1270  real
1271  k2 = sq(sbet1) * g->ep2,
1272  eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
1273  lamscale = g->f * cbet1 * A3f(g, eps) * pi;
1274  }
1275  betscale = lamscale * cbet1;
1276 
1277  x = lam12x / lamscale;
1278  y = sbet12a / betscale;
1279  } else { /* f < 0 */
1280  /* x = dlat, y = dlong */
1281  real
1282  cbet12a = cbet2 * cbet1 - sbet2 * sbet1,
1283  bet12a = atan2(sbet12a, cbet12a);
1284  real m12b, m0;
1285  /* In the case of lon12 = 180, this repeats a calculation made in
1286  * Inverse. */
1287  Lengths(g, g->n, pi + bet12a,
1288  sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
1289  cbet1, cbet2, nullptr, &m12b, &m0, nullptr, nullptr, Ca);
1290  x = -1 + m12b / (cbet1 * cbet2 * m0 * pi);
1291  betscale = x < -(real)(0.01) ? sbet12a / x :
1292  -g->f * sq(cbet1) * pi;
1293  lamscale = betscale / cbet1;
1294  y = lam12x / lamscale;
1295  }
1296 
1297  if (y > -tol1 && x > -1 - xthresh) {
1298  /* strip near cut */
1299  if (g->f >= 0) {
1300  salp1 = minx((real)(1), -(real)(x)); calp1 = - sqrt(1 - sq(salp1));
1301  } else {
1302  calp1 = maxx((real)(x > -tol1 ? 0 : -1), (real)(x));
1303  salp1 = sqrt(1 - sq(calp1));
1304  }
1305  } else {
1306  /* Estimate alp1, by solving the astroid problem.
1307  *
1308  * Could estimate alpha1 = theta + pi/2, directly, i.e.,
1309  * calp1 = y/k; salp1 = -x/(1+k); for f >= 0
1310  * calp1 = x/(1+k); salp1 = -y/k; for f < 0 (need to check)
1311  *
1312  * However, it's better to estimate omg12 from astroid and use
1313  * spherical formula to compute alp1. This reduces the mean number of
1314  * Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12
1315  * (min 0 max 5). The changes in the number of iterations are as
1316  * follows:
1317  *
1318  * change percent
1319  * 1 5
1320  * 0 78
1321  * -1 16
1322  * -2 0.6
1323  * -3 0.04
1324  * -4 0.002
1325  *
1326  * The histogram of iterations is (m = number of iterations estimating
1327  * alp1 directly, n = number of iterations estimating via omg12, total
1328  * number of trials = 148605):
1329  *
1330  * iter m n
1331  * 0 148 186
1332  * 1 13046 13845
1333  * 2 93315 102225
1334  * 3 36189 32341
1335  * 4 5396 7
1336  * 5 455 1
1337  * 6 56 0
1338  *
1339  * Because omg12 is near pi, estimate work with omg12a = pi - omg12 */
1340  real k = Astroid(x, y);
1341  real
1342  omg12a = lamscale * ( g->f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
1343  somg12 = sin(omg12a); comg12 = -cos(omg12a);
1344  /* Update spherical estimate of alp1 using omg12 instead of lam12 */
1345  salp1 = cbet2 * somg12;
1346  calp1 = sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
1347  }
1348  }
1349  /* Sanity check on starting guess. Backwards check allows NaN through. */
1350  if (!(salp1 <= 0))
1351  norm2(&salp1, &calp1);
1352  else {
1353  salp1 = 1; calp1 = 0;
1354  }
1355 
1356  *psalp1 = salp1;
1357  *pcalp1 = calp1;
1358  if (shortline)
1359  *pdnm = dnm;
1360  if (sig12 >= 0) {
1361  *psalp2 = salp2;
1362  *pcalp2 = calp2;
1363  }
1364  return sig12;
1365 }
1366 
1367 real Lambda12(const struct geod_geodesic* g,
1368  real sbet1, real cbet1, real dn1,
1369  real sbet2, real cbet2, real dn2,
1370  real salp1, real calp1,
1371  real slam120, real clam120,
1372  real* psalp2, real* pcalp2,
1373  real* psig12,
1374  real* pssig1, real* pcsig1,
1375  real* pssig2, real* pcsig2,
1376  real* peps,
1377  real* pdomg12,
1378  boolx diffp, real* pdlam12,
1379  /* Scratch area of the right size */
1380  real Ca[]) {
1381  real salp2 = 0, calp2 = 0, sig12 = 0,
1382  ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0,
1383  domg12 = 0, dlam12 = 0;
1384  real salp0, calp0;
1385  real somg1, comg1, somg2, comg2, somg12, comg12, lam12;
1386  real B312, eta, k2;
1387 
1388  if (sbet1 == 0 && calp1 == 0)
1389  /* Break degeneracy of equatorial line. This case has already been
1390  * handled. */
1391  calp1 = -tiny;
1392 
1393  /* sin(alp1) * cos(bet1) = sin(alp0) */
1394  salp0 = salp1 * cbet1;
1395  calp0 = hypot(calp1, salp1 * sbet1); /* calp0 > 0 */
1396 
1397  /* tan(bet1) = tan(sig1) * cos(alp1)
1398  * tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1) */
1399  ssig1 = sbet1; somg1 = salp0 * sbet1;
1400  csig1 = comg1 = calp1 * cbet1;
1401  norm2(&ssig1, &csig1);
1402  /* norm2(&somg1, &comg1); -- don't need to normalize! */
1403 
1404  /* Enforce symmetries in the case abs(bet2) = -bet1. Need to be careful
1405  * about this case, since this can yield singularities in the Newton
1406  * iteration.
1407  * sin(alp2) * cos(bet2) = sin(alp0) */
1408  salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;
1409  /* calp2 = sqrt(1 - sq(salp2))
1410  * = sqrt(sq(calp0) - sq(sbet2)) / cbet2
1411  * and subst for calp0 and rearrange to give (choose positive sqrt
1412  * to give alp2 in [0, pi/2]). */
1413  calp2 = cbet2 != cbet1 || fabs(sbet2) != -sbet1 ?
1414  sqrt(sq(calp1 * cbet1) +
1415  (cbet1 < -sbet1 ?
1416  (cbet2 - cbet1) * (cbet1 + cbet2) :
1417  (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :
1418  fabs(calp1);
1419  /* tan(bet2) = tan(sig2) * cos(alp2)
1420  * tan(omg2) = sin(alp0) * tan(sig2). */
1421  ssig2 = sbet2; somg2 = salp0 * sbet2;
1422  csig2 = comg2 = calp2 * cbet2;
1423  norm2(&ssig2, &csig2);
1424  /* norm2(&somg2, &comg2); -- don't need to normalize! */
1425 
1426  /* sig12 = sig2 - sig1, limit to [0, pi] */
1427  sig12 = atan2(maxx((real)(0), csig1 * ssig2 - ssig1 * csig2),
1428  csig1 * csig2 + ssig1 * ssig2);
1429 
1430  /* omg12 = omg2 - omg1, limit to [0, pi] */
1431  somg12 = maxx((real)(0), comg1 * somg2 - somg1 * comg2);
1432  comg12 = comg1 * comg2 + somg1 * somg2;
1433  /* eta = omg12 - lam120 */
1434  eta = atan2(somg12 * clam120 - comg12 * slam120,
1435  comg12 * clam120 + somg12 * slam120);
1436  k2 = sq(calp0) * g->ep2;
1437  eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);
1438  C3f(g, eps, Ca);
1439  B312 = (SinCosSeries(TRUE, ssig2, csig2, Ca, nC3-1) -
1440  SinCosSeries(TRUE, ssig1, csig1, Ca, nC3-1));
1441  domg12 = -g->f * A3f(g, eps) * salp0 * (sig12 + B312);
1442  lam12 = eta + domg12;
1443 
1444  if (diffp) {
1445  if (calp2 == 0)
1446  dlam12 = - 2 * g->f1 * dn1 / sbet1;
1447  else {
1448  Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1449  cbet1, cbet2, nullptr, &dlam12, nullptr, nullptr, nullptr, Ca);
1450  dlam12 *= g->f1 / (calp2 * cbet2);
1451  }
1452  }
1453 
1454  *psalp2 = salp2;
1455  *pcalp2 = calp2;
1456  *psig12 = sig12;
1457  *pssig1 = ssig1;
1458  *pcsig1 = csig1;
1459  *pssig2 = ssig2;
1460  *pcsig2 = csig2;
1461  *peps = eps;
1462  *pdomg12 = domg12;
1463  if (diffp)
1464  *pdlam12 = dlam12;
1465 
1466  return lam12;
1467 }
1468 
1469 real A3f(const struct geod_geodesic* g, real eps) {
1470  /* Evaluate A3 */
1471  return polyval(nA3 - 1, g->A3x, eps);
1472 }
1473 
1474 void C3f(const struct geod_geodesic* g, real eps, real c[]) {
1475  /* Evaluate C3 coeffs
1476  * Elements c[1] through c[nC3 - 1] are set */
1477  real mult = 1;
1478  int o = 0, l;
1479  for (l = 1; l < nC3; ++l) { /* l is index of C3[l] */
1480  int m = nC3 - l - 1; /* order of polynomial in eps */
1481  mult *= eps;
1482  c[l] = mult * polyval(m, g->C3x + o, eps);
1483  o += m + 1;
1484  }
1485 }
1486 
1487 void C4f(const struct geod_geodesic* g, real eps, real c[]) {
1488  /* Evaluate C4 coeffs
1489  * Elements c[0] through c[nC4 - 1] are set */
1490  real mult = 1;
1491  int o = 0, l;
1492  for (l = 0; l < nC4; ++l) { /* l is index of C4[l] */
1493  int m = nC4 - l - 1; /* order of polynomial in eps */
1494  c[l] = mult * polyval(m, g->C4x + o, eps);
1495  o += m + 1;
1496  mult *= eps;
1497  }
1498 }
1499 
1500 /* The scale factor A1-1 = mean value of (d/dsigma)I1 - 1 */
1501 real A1m1f(real eps) {
1502  static const real coeff[] = {
1503  /* (1-eps)*A1-1, polynomial in eps2 of order 3 */
1504  1, 4, 64, 0, 256,
1505  };
1506  int m = nA1/2;
1507  real t = polyval(m, coeff, sq(eps)) / coeff[m + 1];
1508  return (t + eps) / (1 - eps);
1509 }
1510 
1511 /* The coefficients C1[l] in the Fourier expansion of B1 */
1512 void C1f(real eps, real c[]) {
1513  static const real coeff[] = {
1514  /* C1[1]/eps^1, polynomial in eps2 of order 2 */
1515  -1, 6, -16, 32,
1516  /* C1[2]/eps^2, polynomial in eps2 of order 2 */
1517  -9, 64, -128, 2048,
1518  /* C1[3]/eps^3, polynomial in eps2 of order 1 */
1519  9, -16, 768,
1520  /* C1[4]/eps^4, polynomial in eps2 of order 1 */
1521  3, -5, 512,
1522  /* C1[5]/eps^5, polynomial in eps2 of order 0 */
1523  -7, 1280,
1524  /* C1[6]/eps^6, polynomial in eps2 of order 0 */
1525  -7, 2048,
1526  };
1527  real
1528  eps2 = sq(eps),
1529  d = eps;
1530  int o = 0, l;
1531  for (l = 1; l <= nC1; ++l) { /* l is index of C1p[l] */
1532  int m = (nC1 - l) / 2; /* order of polynomial in eps^2 */
1533  c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
1534  o += m + 2;
1535  d *= eps;
1536  }
1537 }
1538 
1539 /* The coefficients C1p[l] in the Fourier expansion of B1p */
1540 void C1pf(real eps, real c[]) {
1541  static const real coeff[] = {
1542  /* C1p[1]/eps^1, polynomial in eps2 of order 2 */
1543  205, -432, 768, 1536,
1544  /* C1p[2]/eps^2, polynomial in eps2 of order 2 */
1545  4005, -4736, 3840, 12288,
1546  /* C1p[3]/eps^3, polynomial in eps2 of order 1 */
1547  -225, 116, 384,
1548  /* C1p[4]/eps^4, polynomial in eps2 of order 1 */
1549  -7173, 2695, 7680,
1550  /* C1p[5]/eps^5, polynomial in eps2 of order 0 */
1551  3467, 7680,
1552  /* C1p[6]/eps^6, polynomial in eps2 of order 0 */
1553  38081, 61440,
1554  };
1555  real
1556  eps2 = sq(eps),
1557  d = eps;
1558  int o = 0, l;
1559  for (l = 1; l <= nC1p; ++l) { /* l is index of C1p[l] */
1560  int m = (nC1p - l) / 2; /* order of polynomial in eps^2 */
1561  c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
1562  o += m + 2;
1563  d *= eps;
1564  }
1565 }
1566 
1567 /* The scale factor A2-1 = mean value of (d/dsigma)I2 - 1 */
1568 real A2m1f(real eps) {
1569  static const real coeff[] = {
1570  /* (eps+1)*A2-1, polynomial in eps2 of order 3 */
1571  -11, -28, -192, 0, 256,
1572  };
1573  int m = nA2/2;
1574  real t = polyval(m, coeff, sq(eps)) / coeff[m + 1];
1575  return (t - eps) / (1 + eps);
1576 }
1577 
1578 /* The coefficients C2[l] in the Fourier expansion of B2 */
1579 void C2f(real eps, real c[]) {
1580  static const real coeff[] = {
1581  /* C2[1]/eps^1, polynomial in eps2 of order 2 */
1582  1, 2, 16, 32,
1583  /* C2[2]/eps^2, polynomial in eps2 of order 2 */
1584  35, 64, 384, 2048,
1585  /* C2[3]/eps^3, polynomial in eps2 of order 1 */
1586  15, 80, 768,
1587  /* C2[4]/eps^4, polynomial in eps2 of order 1 */
1588  7, 35, 512,
1589  /* C2[5]/eps^5, polynomial in eps2 of order 0 */
1590  63, 1280,
1591  /* C2[6]/eps^6, polynomial in eps2 of order 0 */
1592  77, 2048,
1593  };
1594  real
1595  eps2 = sq(eps),
1596  d = eps;
1597  int o = 0, l;
1598  for (l = 1; l <= nC2; ++l) { /* l is index of C2[l] */
1599  int m = (nC2 - l) / 2; /* order of polynomial in eps^2 */
1600  c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
1601  o += m + 2;
1602  d *= eps;
1603  }
1604 }
1605 
1606 /* The scale factor A3 = mean value of (d/dsigma)I3 */
1607 void A3coeff(struct geod_geodesic* g) {
1608  static const real coeff[] = {
1609  /* A3, coeff of eps^5, polynomial in n of order 0 */
1610  -3, 128,
1611  /* A3, coeff of eps^4, polynomial in n of order 1 */
1612  -2, -3, 64,
1613  /* A3, coeff of eps^3, polynomial in n of order 2 */
1614  -1, -3, -1, 16,
1615  /* A3, coeff of eps^2, polynomial in n of order 2 */
1616  3, -1, -2, 8,
1617  /* A3, coeff of eps^1, polynomial in n of order 1 */
1618  1, -1, 2,
1619  /* A3, coeff of eps^0, polynomial in n of order 0 */
1620  1, 1,
1621  };
1622  int o = 0, k = 0, j;
1623  for (j = nA3 - 1; j >= 0; --j) { /* coeff of eps^j */
1624  int m = nA3 - j - 1 < j ? nA3 - j - 1 : j; /* order of polynomial in n */
1625  g->A3x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
1626  o += m + 2;
1627  }
1628 }
1629 
1630 /* The coefficients C3[l] in the Fourier expansion of B3 */
1631 void C3coeff(struct geod_geodesic* g) {
1632  static const real coeff[] = {
1633  /* C3[1], coeff of eps^5, polynomial in n of order 0 */
1634  3, 128,
1635  /* C3[1], coeff of eps^4, polynomial in n of order 1 */
1636  2, 5, 128,
1637  /* C3[1], coeff of eps^3, polynomial in n of order 2 */
1638  -1, 3, 3, 64,
1639  /* C3[1], coeff of eps^2, polynomial in n of order 2 */
1640  -1, 0, 1, 8,
1641  /* C3[1], coeff of eps^1, polynomial in n of order 1 */
1642  -1, 1, 4,
1643  /* C3[2], coeff of eps^5, polynomial in n of order 0 */
1644  5, 256,
1645  /* C3[2], coeff of eps^4, polynomial in n of order 1 */
1646  1, 3, 128,
1647  /* C3[2], coeff of eps^3, polynomial in n of order 2 */
1648  -3, -2, 3, 64,
1649  /* C3[2], coeff of eps^2, polynomial in n of order 2 */
1650  1, -3, 2, 32,
1651  /* C3[3], coeff of eps^5, polynomial in n of order 0 */
1652  7, 512,
1653  /* C3[3], coeff of eps^4, polynomial in n of order 1 */
1654  -10, 9, 384,
1655  /* C3[3], coeff of eps^3, polynomial in n of order 2 */
1656  5, -9, 5, 192,
1657  /* C3[4], coeff of eps^5, polynomial in n of order 0 */
1658  7, 512,
1659  /* C3[4], coeff of eps^4, polynomial in n of order 1 */
1660  -14, 7, 512,
1661  /* C3[5], coeff of eps^5, polynomial in n of order 0 */
1662  21, 2560,
1663  };
1664  int o = 0, k = 0, l, j;
1665  for (l = 1; l < nC3; ++l) { /* l is index of C3[l] */
1666  for (j = nC3 - 1; j >= l; --j) { /* coeff of eps^j */
1667  int m = nC3 - j - 1 < j ? nC3 - j - 1 : j; /* order of polynomial in n */
1668  g->C3x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
1669  o += m + 2;
1670  }
1671  }
1672 }
1673 
1674 /* The coefficients C4[l] in the Fourier expansion of I4 */
1675 void C4coeff(struct geod_geodesic* g) {
1676  static const real coeff[] = {
1677  /* C4[0], coeff of eps^5, polynomial in n of order 0 */
1678  97, 15015,
1679  /* C4[0], coeff of eps^4, polynomial in n of order 1 */
1680  1088, 156, 45045,
1681  /* C4[0], coeff of eps^3, polynomial in n of order 2 */
1682  -224, -4784, 1573, 45045,
1683  /* C4[0], coeff of eps^2, polynomial in n of order 3 */
1684  -10656, 14144, -4576, -858, 45045,
1685  /* C4[0], coeff of eps^1, polynomial in n of order 4 */
1686  64, 624, -4576, 6864, -3003, 15015,
1687  /* C4[0], coeff of eps^0, polynomial in n of order 5 */
1688  100, 208, 572, 3432, -12012, 30030, 45045,
1689  /* C4[1], coeff of eps^5, polynomial in n of order 0 */
1690  1, 9009,
1691  /* C4[1], coeff of eps^4, polynomial in n of order 1 */
1692  -2944, 468, 135135,
1693  /* C4[1], coeff of eps^3, polynomial in n of order 2 */
1694  5792, 1040, -1287, 135135,
1695  /* C4[1], coeff of eps^2, polynomial in n of order 3 */
1696  5952, -11648, 9152, -2574, 135135,
1697  /* C4[1], coeff of eps^1, polynomial in n of order 4 */
1698  -64, -624, 4576, -6864, 3003, 135135,
1699  /* C4[2], coeff of eps^5, polynomial in n of order 0 */
1700  8, 10725,
1701  /* C4[2], coeff of eps^4, polynomial in n of order 1 */
1702  1856, -936, 225225,
1703  /* C4[2], coeff of eps^3, polynomial in n of order 2 */
1704  -8448, 4992, -1144, 225225,
1705  /* C4[2], coeff of eps^2, polynomial in n of order 3 */
1706  -1440, 4160, -4576, 1716, 225225,
1707  /* C4[3], coeff of eps^5, polynomial in n of order 0 */
1708  -136, 63063,
1709  /* C4[3], coeff of eps^4, polynomial in n of order 1 */
1710  1024, -208, 105105,
1711  /* C4[3], coeff of eps^3, polynomial in n of order 2 */
1712  3584, -3328, 1144, 315315,
1713  /* C4[4], coeff of eps^5, polynomial in n of order 0 */
1714  -128, 135135,
1715  /* C4[4], coeff of eps^4, polynomial in n of order 1 */
1716  -2560, 832, 405405,
1717  /* C4[5], coeff of eps^5, polynomial in n of order 0 */
1718  128, 99099,
1719  };
1720  int o = 0, k = 0, l, j;
1721  for (l = 0; l < nC4; ++l) { /* l is index of C4[l] */
1722  for (j = nC4 - 1; j >= l; --j) { /* coeff of eps^j */
1723  int m = nC4 - j - 1; /* order of polynomial in n */
1724  g->C4x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
1725  o += m + 2;
1726  }
1727  }
1728 }
1729 
1730 int transit(real lon1, real lon2) {
1731  real lon12;
1732  /* Return 1 or -1 if crossing prime meridian in east or west direction.
1733  * Otherwise return zero. */
1734  /* Compute lon12 the same way as Geodesic::Inverse. */
1735  lon1 = AngNormalize(lon1);
1736  lon2 = AngNormalize(lon2);
1737  lon12 = AngDiff(lon1, lon2, nullptr);
1738  return lon1 <= 0 && lon2 > 0 && lon12 > 0 ? 1 :
1739  (lon2 <= 0 && lon1 > 0 && lon12 < 0 ? -1 : 0);
1740 }
1741 
1742 int transitdirect(real lon1, real lon2) {
1743  /* Compute exactly the parity of
1744  int(ceil(lon2 / 360)) - int(ceil(lon1 / 360)) */
1745  lon1 = remainder(lon1, (real)(720));
1746  lon2 = remainder(lon2, (real)(720));
1747  return ( (lon2 <= 0 && lon2 > -360 ? 1 : 0) -
1748  (lon1 <= 0 && lon1 > -360 ? 1 : 0) );
1749 }
1750 
1751 void accini(real s[]) {
1752  /* Initialize an accumulator; this is an array with two elements. */
1753  s[0] = s[1] = 0;
1754 }
1755 
1756 void acccopy(const real s[], real t[]) {
1757  /* Copy an accumulator; t = s. */
1758  t[0] = s[0]; t[1] = s[1];
1759 }
1760 
1761 void accadd(real s[], real y) {
1762  /* Add y to an accumulator. */
1763  real u, z = sumx(y, s[1], &u);
1764  s[0] = sumx(z, s[0], &s[1]);
1765  if (s[0] == 0)
1766  s[0] = u;
1767  else
1768  s[1] = s[1] + u;
1769 }
1770 
1771 real accsum(const real s[], real y) {
1772  /* Return accumulator + y (but don't add to accumulator). */
1773  real t[2];
1774  acccopy(s, t);
1775  accadd(t, y);
1776  return t[0];
1777 }
1778 
1779 void accneg(real s[]) {
1780  /* Negate an accumulator. */
1781  s[0] = -s[0]; s[1] = -s[1];
1782 }
1783 
1784 void accrem(real s[], real y) {
1785  /* Reduce to [-y/2, y/2]. */
1786  s[0] = remainder(s[0], y);
1787  accadd(s, (real)(0));
1788 }
1789 
1790 void geod_polygon_init(struct geod_polygon* p, boolx polylinep) {
1791  p->polyline = (polylinep != 0);
1792  geod_polygon_clear(p);
1793 }
1794 
1795 void geod_polygon_clear(struct geod_polygon* p) {
1796  p->lat0 = p->lon0 = p->lat = p->lon = NaN;
1797  accini(p->P);
1798  accini(p->A);
1799  p->num = p->crossings = 0;
1800 }
1801 
1802 void geod_polygon_addpoint(const struct geod_geodesic* g,
1803  struct geod_polygon* p,
1804  real lat, real lon) {
1805  lon = AngNormalize(lon);
1806  if (p->num == 0) {
1807  p->lat0 = p->lat = lat;
1808  p->lon0 = p->lon = lon;
1809  } else {
1810  real s12, S12 = 0; /* Initialize S12 to stop Visual Studio warning */
1811  geod_geninverse(g, p->lat, p->lon, lat, lon,
1812  &s12, nullptr, nullptr, nullptr, nullptr, nullptr,
1813  p->polyline ? nullptr : &S12);
1814  accadd(p->P, s12);
1815  if (!p->polyline) {
1816  accadd(p->A, S12);
1817  p->crossings += transit(p->lon, lon);
1818  }
1819  p->lat = lat; p->lon = lon;
1820  }
1821  ++p->num;
1822 }
1823 
1824 void geod_polygon_addedge(const struct geod_geodesic* g,
1825  struct geod_polygon* p,
1826  real azi, real s) {
1827  if (p->num) { /* Do nothing is num is zero */
1828  /* Initialize S12 to stop Visual Studio warning. Initialization of lat and
1829  * lon is to make CLang static analyzer happy. */
1830  real lat = 0, lon = 0, S12 = 0;
1831  geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s,
1832  &lat, &lon, nullptr,
1833  nullptr, nullptr, nullptr, nullptr,
1834  p->polyline ? nullptr : &S12);
1835  accadd(p->P, s);
1836  if (!p->polyline) {
1837  accadd(p->A, S12);
1838  p->crossings += transitdirect(p->lon, lon);
1839  }
1840  p->lat = lat; p->lon = lon;
1841  ++p->num;
1842  }
1843 }
1844 
1845 unsigned geod_polygon_compute(const struct geod_geodesic* g,
1846  const struct geod_polygon* p,
1847  boolx reverse, boolx sign,
1848  real* pA, real* pP) {
1849  real s12, S12, t[2];
1850  if (p->num < 2) {
1851  if (pP) *pP = 0;
1852  if (!p->polyline && pA) *pA = 0;
1853  return p->num;
1854  }
1855  if (p->polyline) {
1856  if (pP) *pP = p->P[0];
1857  return p->num;
1858  }
1859  geod_geninverse(g, p->lat, p->lon, p->lat0, p->lon0,
1860  &s12, nullptr, nullptr, nullptr, nullptr, nullptr, &S12);
1861  if (pP) *pP = accsum(p->P, s12);
1862  acccopy(p->A, t);
1863  accadd(t, S12);
1864  if (pA) *pA = areareduceA(t, 4 * pi * g->c2,
1865  p->crossings + transit(p->lon, p->lon0),
1866  reverse, sign);
1867  return p->num;
1868 }
1869 
1870 unsigned geod_polygon_testpoint(const struct geod_geodesic* g,
1871  const struct geod_polygon* p,
1872  real lat, real lon,
1873  boolx reverse, boolx sign,
1874  real* pA, real* pP) {
1875  real perimeter, tempsum;
1876  int crossings, i;
1877  unsigned num = p->num + 1;
1878  if (num == 1) {
1879  if (pP) *pP = 0;
1880  if (!p->polyline && pA) *pA = 0;
1881  return num;
1882  }
1883  perimeter = p->P[0];
1884  tempsum = p->polyline ? 0 : p->A[0];
1885  crossings = p->crossings;
1886  for (i = 0; i < (p->polyline ? 1 : 2); ++i) {
1887  real s12, S12 = 0; /* Initialize S12 to stop Visual Studio warning */
1888  geod_geninverse(g,
1889  i == 0 ? p->lat : lat, i == 0 ? p->lon : lon,
1890  i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon,
1891  &s12, nullptr, nullptr, nullptr, nullptr, nullptr,
1892  p->polyline ? nullptr : &S12);
1893  perimeter += s12;
1894  if (!p->polyline) {
1895  tempsum += S12;
1896  crossings += transit(i == 0 ? p->lon : lon,
1897  i != 0 ? p->lon0 : lon);
1898  }
1899  }
1900 
1901  if (pP) *pP = perimeter;
1902  if (p->polyline)
1903  return num;
1904 
1905  if (pA) *pA = areareduceB(tempsum, 4 * pi * g->c2, crossings, reverse, sign);
1906  return num;
1907 }
1908 
1909 unsigned geod_polygon_testedge(const struct geod_geodesic* g,
1910  const struct geod_polygon* p,
1911  real azi, real s,
1912  boolx reverse, boolx sign,
1913  real* pA, real* pP) {
1914  real perimeter, tempsum;
1915  int crossings;
1916  unsigned num = p->num + 1;
1917  if (num == 1) { /* we don't have a starting point! */
1918  if (pP) *pP = NaN;
1919  if (!p->polyline && pA) *pA = NaN;
1920  return 0;
1921  }
1922  perimeter = p->P[0] + s;
1923  if (p->polyline) {
1924  if (pP) *pP = perimeter;
1925  return num;
1926  }
1927 
1928  tempsum = p->A[0];
1929  crossings = p->crossings;
1930  {
1931  /* Initialization of lat, lon, and S12 is to make CLang static analyzer
1932  * happy. */
1933  real lat = 0, lon = 0, s12, S12 = 0;
1934  geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s,
1935  &lat, &lon, nullptr,
1936  nullptr, nullptr, nullptr, nullptr, &S12);
1937  tempsum += S12;
1938  crossings += transitdirect(p->lon, lon);
1939  geod_geninverse(g, lat, lon, p->lat0, p->lon0,
1940  &s12, nullptr, nullptr, nullptr, nullptr, nullptr, &S12);
1941  perimeter += s12;
1942  tempsum += S12;
1943  crossings += transit(lon, p->lon0);
1944  }
1945 
1946  if (pP) *pP = perimeter;
1947  if (pA) *pA = areareduceB(tempsum, 4 * pi * g->c2, crossings, reverse, sign);
1948  return num;
1949 }
1950 
1951 void geod_polygonarea(const struct geod_geodesic* g,
1952  real lats[], real lons[], int n,
1953  real* pA, real* pP) {
1954  int i;
1955  struct geod_polygon p;
1956  geod_polygon_init(&p, FALSE);
1957  for (i = 0; i < n; ++i)
1958  geod_polygon_addpoint(g, &p, lats[i], lons[i]);
1959  geod_polygon_compute(g, &p, FALSE, TRUE, pA, pP);
1960 }
1961 
1962 real areareduceA(real area[], real area0,
1963  int crossings, boolx reverse, boolx sign) {
1964  accrem(area, area0);
1965  if (crossings & 1)
1966  accadd(area, (area[0] < 0 ? 1 : -1) * area0/2);
1967  /* area is with the clockwise sense. If !reverse convert to
1968  * counter-clockwise convention. */
1969  if (!reverse)
1970  accneg(area);
1971  /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
1972  if (sign) {
1973  if (area[0] > area0/2)
1974  accadd(area, -area0);
1975  else if (area[0] <= -area0/2)
1976  accadd(area, +area0);
1977  } else {
1978  if (area[0] >= area0)
1979  accadd(area, -area0);
1980  else if (area[0] < 0)
1981  accadd(area, +area0);
1982  }
1983  return 0 + area[0];
1984 }
1985 
1986 real areareduceB(real area, real area0,
1987  int crossings, boolx reverse, boolx sign) {
1988  area = remainder(area, area0);
1989  if (crossings & 1)
1990  area += (area < 0 ? 1 : -1) * area0/2;
1991  /* area is with the clockwise sense. If !reverse convert to
1992  * counter-clockwise convention. */
1993  if (!reverse)
1994  area *= -1;
1995  /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
1996  if (sign) {
1997  if (area > area0/2)
1998  area -= area0;
1999  else if (area <= -area0/2)
2000  area += area0;
2001  } else {
2002  if (area >= area0)
2003  area -= area0;
2004  else if (area < 0)
2005  area += area0;
2006  }
2007  return 0 + area;
2008 }
2009 
2010 /** @endcond */
GeographicLib::Math::real real
API for the geodesic routines in C.
@ GEOD_NOFLAGS
Definition: geodesic.h:924
@ GEOD_LONG_UNROLL
Definition: geodesic.h:926
@ GEOD_ARCMODE
Definition: geodesic.h:925
@ GEOD_AZIMUTH
Definition: geodesic.h:910
@ GEOD_REDUCEDLENGTH
Definition: geodesic.h:913
@ GEOD_LONGITUDE
Definition: geodesic.h:909
@ GEOD_LATITUDE
Definition: geodesic.h:908
@ GEOD_AREA
Definition: geodesic.h:915
@ GEOD_GEODESICSCALE
Definition: geodesic.h:914
@ GEOD_DISTANCE
Definition: geodesic.h:911
@ GEOD_DISTANCE_IN
Definition: geodesic.h:912
@ GEOD_NONE
Definition: geodesic.h:907
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)
void GEOD_DLL geod_polygon_addedge(const struct geod_geodesic *g, struct geod_polygon *p, double azi, double s)
void GEOD_DLL geod_polygon_addpoint(const struct geod_geodesic *g, struct geod_polygon *p, double lat, double lon)
void GEOD_DLL geod_setdistance(struct geod_geodesicline *l, double s13)
void GEOD_DLL geod_gensetdistance(struct geod_geodesicline *l, unsigned flags, double s13_a13)
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)
void GEOD_DLL geod_polygon_clear(struct geod_polygon *p)
void GEOD_DLL geod_init(struct geod_geodesic *g, double a, double f)
void GEOD_DLL geod_directline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, unsigned caps)
void GEOD_DLL geod_inverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2)
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)
void GEOD_DLL geod_lineinit(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned caps)
unsigned GEOD_DLL geod_polygon_compute(const struct geod_geodesic *g, const struct geod_polygon *p, int reverse, int sign, double *pA, double *pP)
void GEOD_DLL geod_position(const struct geod_geodesicline *l, double s12, double *plat2, double *plon2, double *pazi2)
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)
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)
void GEOD_DLL geod_direct(const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, double *plat2, double *plon2, double *pazi2)
void GEOD_DLL geod_polygonarea(const struct geod_geodesic *g, double lats[], double lons[], int n, double *pA, double *pP)
void GEOD_DLL geod_inverseline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, unsigned caps)
void GEOD_DLL geod_gendirectline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned flags, double s12_a12, unsigned caps)
void GEOD_DLL geod_polygon_init(struct geod_polygon *p, int polylinep)
double f
Definition: geodesic.h:184
double a
Definition: geodesic.h:183
unsigned caps
Definition: geodesic.h:212
double lon
Definition: geodesic.h:222
unsigned num
Definition: geodesic.h:231
double lat
Definition: geodesic.h:221