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