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