C library for Geodesics  1.44
geodesic.h
Go to the documentation of this file.
1 /**
2  * \file geodesic.h
3  * \brief Header for the geodesic routines in C
4  *
5  * This an implementation in C of the geodesic algorithms described in
6  * - C. F. F. Karney,
7  * <a href="https://dx.doi.org/10.1007/s00190-012-0578-z">
8  * Algorithms for geodesics</a>,
9  * J. Geodesy <b>87</b>, 43--55 (2013);
10  * DOI: <a href="https://dx.doi.org/10.1007/s00190-012-0578-z">
11  * 10.1007/s00190-012-0578-z</a>;
12  * addenda: <a href="http://geographiclib.sf.net/geod-addenda.html">
13  * geod-addenda.html</a>.
14  * .
15  * The principal advantages of these algorithms over previous ones (e.g.,
16  * Vincenty, 1975) are
17  * - accurate to round off for |<i>f</i>| &lt; 1/50;
18  * - the solution of the inverse problem is always found;
19  * - differential and integral properties of geodesics are computed.
20  *
21  * The shortest path between two points on the ellipsoid at (\e lat1, \e
22  * lon1) and (\e lat2, \e lon2) is called the geodesic. Its length is
23  * \e s12 and the geodesic from point 1 to point 2 has forward azimuths
24  * \e azi1 and \e azi2 at the two end points.
25  *
26  * Traditionally two geodesic problems are considered:
27  * - the direct problem -- given \e lat1, \e lon1, \e s12, and \e azi1,
28  * determine \e lat2, \e lon2, and \e azi2. This is solved by the function
29  * geod_direct().
30  * - the inverse problem -- given \e lat1, \e lon1, and \e lat2, \e lon2,
31  * determine \e s12, \e azi1, and \e azi2. This is solved by the function
32  * geod_inverse().
33  *
34  * The ellipsoid is specified by its equatorial radius \e a (typically in
35  * meters) and flattening \e f. The routines are accurate to round off with
36  * double precision arithmetic provided that |<i>f</i>| &lt; 1/50; for the
37  * WGS84 ellipsoid, the errors are less than 15 nanometers. (Reasonably
38  * accurate results are obtained for |<i>f</i>| &lt; 1/5.) For a prolate
39  * ellipsoid, specify \e f &lt; 0.
40  *
41  * The routines also calculate several other quantities of interest
42  * - \e S12 is the area between the geodesic from point 1 to point 2 and the
43  * equator; i.e., it is the area, measured counter-clockwise, of the
44  * quadrilateral with corners (\e lat1,\e lon1), (0,\e lon1), (0,\e lon2),
45  * and (\e lat2,\e lon2).
46  * - \e m12, the reduced length of the geodesic is defined such that if
47  * the initial azimuth is perturbed by \e dazi1 (radians) then the
48  * second point is displaced by \e m12 \e dazi1 in the direction
49  * perpendicular to the geodesic. On a curved surface the reduced
50  * length obeys a symmetry relation, \e m12 + \e m21 = 0. On a flat
51  * surface, we have \e m12 = \e s12.
52  * - \e M12 and \e M21 are geodesic scales. If two geodesics are
53  * parallel at point 1 and separated by a small distance \e dt, then
54  * they are separated by a distance \e M12 \e dt at point 2. \e M21
55  * is defined similarly (with the geodesics being parallel to one
56  * another at point 2). On a flat surface, we have \e M12 = \e M21
57  * = 1.
58  * - \e a12 is the arc length on the auxiliary sphere. This is a
59  * construct for converting the problem to one in spherical
60  * trigonometry. \e a12 is measured in degrees. The spherical arc
61  * length from one equator crossing to the next is always 180&deg;.
62  *
63  * If points 1, 2, and 3 lie on a single geodesic, then the following
64  * addition rules hold:
65  * - \e s13 = \e s12 + \e s23
66  * - \e a13 = \e a12 + \e a23
67  * - \e S13 = \e S12 + \e S23
68  * - \e m13 = \e m12 \e M23 + \e m23 \e M21
69  * - \e M13 = \e M12 \e M23 &minus; (1 &minus; \e M12 \e M21) \e
70  * m23 / \e m12
71  * - \e M31 = \e M32 \e M21 &minus; (1 &minus; \e M23 \e M32) \e
72  * m12 / \e m23
73  *
74  * The shortest distance returned by the solution of the inverse problem is
75  * (obviously) uniquely defined. However, in a few special cases there are
76  * multiple azimuths which yield the same shortest distance. Here is a
77  * catalog of those cases:
78  * - \e lat1 = &minus;\e lat2 (with neither point at a pole). If \e azi1 =
79  * \e azi2, the geodesic is unique. Otherwise there are two geodesics
80  * and the second one is obtained by setting [\e azi1, \e azi2] = [\e
81  * azi2, \e azi1], [\e M12, \e M21] = [\e M21, \e M12], \e S12 =
82  * &minus;\e S12. (This occurs when the longitude difference is near
83  * &plusmn;180&deg; for oblate ellipsoids.)
84  * - \e lon2 = \e lon1 &plusmn; 180&deg; (with neither point at a pole).
85  * If \e azi1 = 0&deg; or &plusmn;180&deg;, the geodesic is unique.
86  * Otherwise there are two geodesics and the second one is obtained by
87  * setting [\e azi1, \e azi2] = [&minus;\e azi1, &minus;\e azi2], \e S12
88  * = &minus;\e S12. (This occurs when \e lat2 is near &minus;\e lat1 for
89  * prolate ellipsoids.)
90  * - Points 1 and 2 at opposite poles. There are infinitely many
91  * geodesics which can be generated by setting [\e azi1, \e azi2] =
92  * [\e azi1, \e azi2] + [\e d, &minus;\e d], for arbitrary \e d. (For
93  * spheres, this prescription applies when points 1 and 2 are
94  * antipodal.)
95  * - \e s12 = 0 (coincident points). There are infinitely many geodesics
96  * which can be generated by setting [\e azi1, \e azi2] = [\e azi1, \e
97  * azi2] + [\e d, \e d], for arbitrary \e d.
98  *
99  * These routines are a simple transcription of the corresponding C++ classes
100  * in <a href="http://geographiclib.sf.net"> GeographicLib</a>. The "class
101  * data" is represented by the structs geod_geodesic, geod_geodesicline,
102  * geod_polygon and pointers to these objects are passed as initial arguments
103  * to the member functions. Most of the internal comments have been retained.
104  * However, in the process of transcription some documentation has been lost
105  * and the documentation for the C++ classes, GeographicLib::Geodesic,
106  * GeographicLib::GeodesicLine, and GeographicLib::PolygonAreaT, should be
107  * consulted. The C++ code remains the "reference implementation". Think
108  * twice about restructuring the internals of the C code since this may make
109  * porting fixes from the C++ code more difficult.
110  *
111  * Copyright (c) Charles Karney (2012-2015) <charles@karney.com> and licensed
112  * under the MIT/X11 License. For more information, see
113  * http://geographiclib.sourceforge.net/
114  *
115  * This library was distributed with
116  * <a href="../index.html">GeographicLib</a> 1.44.
117  **********************************************************************/
118 
119 #if !defined(GEODESIC_H)
120 #define GEODESIC_H 1
121 
122 /**
123  * The major version of the geodesic library. (This tracks the version of
124  * GeographicLib.)
125  **********************************************************************/
126 #define GEODESIC_VERSION_MAJOR 1
127 /**
128  * The minor version of the geodesic library. (This tracks the version of
129  * GeographicLib.)
130  **********************************************************************/
131 #define GEODESIC_VERSION_MINOR 44
132 /**
133  * The patch level of the geodesic library. (This tracks the version of
134  * GeographicLib.)
135  **********************************************************************/
136 #define GEODESIC_VERSION_PATCH 0
137 
138 /**
139  * Pack the version components into a single integer. Users should not rely on
140  * this particular packing of the components of the version number; see the
141  * documentation for GEODESIC_VERSION, below.
142  **********************************************************************/
143 #define GEODESIC_VERSION_NUM(a,b,c) ((((a) * 10000 + (b)) * 100) + (c))
144 
145 /**
146  * The version of the geodesic library as a single integer, packed as MMmmmmpp
147  * where MM is the major version, mmmm is the minor version, and pp is the
148  * patch level. Users should not rely on this particular packing of the
149  * components of the version number. Instead they should use a test such as
150  * @code{.c}
151  #if GEODESIC_VERSION >= GEODESIC_VERSION_NUM(1,40,0)
152  ...
153  #endif
154  * @endcode
155  **********************************************************************/
156 #define GEODESIC_VERSION \
157  GEODESIC_VERSION_NUM(GEODESIC_VERSION_MAJOR, \
158  GEODESIC_VERSION_MINOR, \
159  GEODESIC_VERSION_PATCH)
160 
161 #if defined(__cplusplus)
162 extern "C" {
163 #endif
164 
165  /**
166  * The struct containing information about the ellipsoid. This must be
167  * initialized by geod_init() before use.
168  **********************************************************************/
169  struct geod_geodesic {
170  double a; /**< the equatorial radius */
171  double f; /**< the flattening */
172  /**< @cond SKIP */
173  double f1, e2, ep2, n, b, c2, etol2;
174  double A3x[6], C3x[15], C4x[21];
175  /**< @endcond */
176  };
177 
178  /**
179  * The struct containing information about a single geodesic. This must be
180  * initialized by geod_lineinit() before use.
181  **********************************************************************/
183  double lat1; /**< the starting latitude */
184  double lon1; /**< the starting longitude */
185  double azi1; /**< the starting azimuth */
186  double a; /**< the equatorial radius */
187  double f; /**< the flattening */
188  /**< @cond SKIP */
189  double b, c2, f1, salp0, calp0, k2,
190  salp1, calp1, ssig1, csig1, dn1, stau1, ctau1, somg1, comg1,
191  A1m1, A2m1, A3c, B11, B21, B31, A4, B41;
192  double C1a[6+1], C1pa[6+1], C2a[6+1], C3a[6], C4a[6];
193  /**< @endcond */
194  unsigned caps; /**< the capabilities */
195  };
196 
197  /**
198  * The struct for accumulating information about a geodesic polygon. This is
199  * used for computing the perimeter and area of a polygon. This must be
200  * initialized by geod_polygon_init() before use.
201  **********************************************************************/
202  struct geod_polygon {
203  double lat; /**< the current latitude */
204  double lon; /**< the current longitude */
205  /**< @cond SKIP */
206  double lat0;
207  double lon0;
208  double A[2];
209  double P[2];
210  int polyline;
211  int crossings;
212  /**< @endcond */
213  unsigned num; /**< the number of points so far */
214  };
215 
216  /**
217  * Initialize a geod_geodesic object.
218  *
219  * @param[out] g a pointer to the object to be initialized.
220  * @param[in] a the equatorial radius (meters).
221  * @param[in] f the flattening.
222  **********************************************************************/
223  void geod_init(struct geod_geodesic* g, double a, double f);
224 
225  /**
226  * Initialize a geod_geodesicline object.
227  *
228  * @param[out] l a pointer to the object to be initialized.
229  * @param[in] g a pointer to the geod_geodesic object specifying the
230  * ellipsoid.
231  * @param[in] lat1 latitude of point 1 (degrees).
232  * @param[in] lon1 longitude of point 1 (degrees).
233  * @param[in] azi1 azimuth at point 1 (degrees).
234  * @param[in] caps bitor'ed combination of geod_mask() values specifying the
235  * capabilities the geod_geodesicline object should possess, i.e., which
236  * quantities can be returned in calls to geod_position() and
237  * geod_genposition().
238  *
239  * \e g must have been initialized with a call to geod_init(). \e lat1
240  * should be in the range [&minus;90&deg;, 90&deg;].
241  *
242  * The geod_mask values are [see geod_mask()]:
243  * - \e caps |= GEOD_LATITUDE for the latitude \e lat2; this is
244  * added automatically,
245  * - \e caps |= GEOD_LONGITUDE for the latitude \e lon2,
246  * - \e caps |= GEOD_AZIMUTH for the latitude \e azi2; this is
247  * added automatically,
248  * - \e caps |= GEOD_DISTANCE for the distance \e s12,
249  * - \e caps |= GEOD_REDUCEDLENGTH for the reduced length \e m12,
250  * - \e caps |= GEOD_GEODESICSCALE for the geodesic scales \e M12
251  * and \e M21,
252  * - \e caps |= GEOD_AREA for the area \e S12,
253  * - \e caps |= GEOD_DISTANCE_IN permits the length of the
254  * geodesic to be given in terms of \e s12; without this capability the
255  * length can only be specified in terms of arc length.
256  * .
257  * A value of \e caps = 0 is treated as GEOD_LATITUDE | GEOD_LONGITUDE |
258  * GEOD_AZIMUTH | GEOD_DISTANCE_IN (to support the solution of the "standard"
259  * direct problem).
260  **********************************************************************/
261  void geod_lineinit(struct geod_geodesicline* l,
262  const struct geod_geodesic* g,
263  double lat1, double lon1, double azi1, unsigned caps);
264 
265  /**
266  * Solve the direct geodesic problem.
267  *
268  * @param[in] g a pointer to the geod_geodesic object specifying the
269  * ellipsoid.
270  * @param[in] lat1 latitude of point 1 (degrees).
271  * @param[in] lon1 longitude of point 1 (degrees).
272  * @param[in] azi1 azimuth at point 1 (degrees).
273  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
274  * negative.
275  * @param[out] plat2 pointer to the latitude of point 2 (degrees).
276  * @param[out] plon2 pointer to the longitude of point 2 (degrees).
277  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
278  *
279  * \e g must have been initialized with a call to geod_init(). \e lat1
280  * should be in the range [&minus;90&deg;, 90&deg;]. The values of \e lon2
281  * and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;). Any of
282  * the "return" arguments \e plat2, etc., may be replaced by 0, if you do not
283  * need some quantities computed.
284  *
285  * If either point is at a pole, the azimuth is defined by keeping the
286  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;), and
287  * taking the limit &epsilon; &rarr; 0+. An arc length greater that 180&deg;
288  * signifies a geodesic which is not a shortest path. (For a prolate
289  * ellipsoid, an additional condition is necessary for a shortest path: the
290  * longitudinal extent must not exceed of 180&deg;.)
291  *
292  * Example, determine the point 10000 km NE of JFK:
293  @code{.c}
294  struct geod_geodesic g;
295  double lat, lon;
296  geod_init(&g, 6378137, 1/298.257223563);
297  geod_direct(&g, 40.64, -73.78, 45.0, 10e6, &lat, &lon, 0);
298  printf("%.5f %.5f\n", lat, lon);
299  @endcode
300  **********************************************************************/
301  void geod_direct(const struct geod_geodesic* g,
302  double lat1, double lon1, double azi1, double s12,
303  double* plat2, double* plon2, double* pazi2);
304 
305  /**
306  * Solve the inverse geodesic problem.
307  *
308  * @param[in] g a pointer to the geod_geodesic object specifying the
309  * ellipsoid.
310  * @param[in] lat1 latitude of point 1 (degrees).
311  * @param[in] lon1 longitude of point 1 (degrees).
312  * @param[in] lat2 latitude of point 2 (degrees).
313  * @param[in] lon2 longitude of point 2 (degrees).
314  * @param[out] ps12 pointer to the distance between point 1 and point 2
315  * (meters).
316  * @param[out] pazi1 pointer to the azimuth at point 1 (degrees).
317  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
318  *
319  * \e g must have been initialized with a call to geod_init(). \e lat1 and
320  * \e lat2 should be in the range [&minus;90&deg;, 90&deg;]. The values of
321  * \e azi1 and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;).
322  * Any of the "return" arguments, \e ps12, etc., may be replaced by 0, if you
323  * do not need some quantities computed.
324  *
325  * If either point is at a pole, the azimuth is defined by keeping the
326  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;), and
327  * taking the limit &epsilon; &rarr; 0+.
328  *
329  * The solution to the inverse problem is found using Newton's method. If
330  * this fails to converge (this is very unlikely in geodetic applications
331  * but does occur for very eccentric ellipsoids), then the bisection method
332  * is used to refine the solution.
333  *
334  * Example, determine the distance between JFK and Singapore Changi Airport:
335  @code{.c}
336  struct geod_geodesic g;
337  double s12;
338  geod_init(&g, 6378137, 1/298.257223563);
339  geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, 0, 0);
340  printf("%.3f\n", s12);
341  @endcode
342  **********************************************************************/
343  void geod_inverse(const struct geod_geodesic* g,
344  double lat1, double lon1, double lat2, double lon2,
345  double* ps12, double* pazi1, double* pazi2);
346 
347  /**
348  * Compute the position along a geod_geodesicline.
349  *
350  * @param[in] l a pointer to the geod_geodesicline object specifying the
351  * geodesic line.
352  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
353  * negative.
354  * @param[out] plat2 pointer to the latitude of point 2 (degrees).
355  * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires
356  * that \e l was initialized with \e caps |= GEOD_LONGITUDE.
357  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
358  *
359  * \e l must have been initialized with a call to geod_lineinit() with \e
360  * caps |= GEOD_DISTANCE_IN. The values of \e lon2 and \e azi2 returned are
361  * in the range [&minus;180&deg;, 180&deg;). Any of the "return" arguments
362  * \e plat2, etc., may be replaced by 0, if you do not need some quantities
363  * computed.
364  *
365  * Example, compute way points between JFK and Singapore Changi Airport
366  * the "obvious" way using geod_direct():
367  @code{.c}
368  struct geod_geodesic g;
369  double s12, azi1, lat[101],lon[101];
370  int i;
371  geod_init(&g, 6378137, 1/298.257223563);
372  geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0);
373  for (i = 0; i < 101; ++i) {
374  geod_direct(&g, 40.64, -73.78, azi1, i * s12 * 0.01, lat + i, lon + i, 0);
375  printf("%.5f %.5f\n", lat[i], lon[i]);
376  }
377  @endcode
378  * A faster way using geod_position():
379  @code{.c}
380  struct geod_geodesic g;
381  struct geod_geodesicline l;
382  double s12, azi1, lat[101],lon[101];
383  int i;
384  geod_init(&g, 6378137, 1/298.257223563);
385  geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0);
386  geod_lineinit(&l, &g, 40.64, -73.78, azi1, 0);
387  for (i = 0; i < 101; ++i) {
388  geod_position(&l, i * s12 * 0.01, lat + i, lon + i, 0);
389  printf("%.5f %.5f\n", lat[i], lon[i]);
390  }
391  @endcode
392  **********************************************************************/
393  void geod_position(const struct geod_geodesicline* l, double s12,
394  double* plat2, double* plon2, double* pazi2);
395 
396  /**
397  * The general direct geodesic problem.
398  *
399  * @param[in] g a pointer to the geod_geodesic object specifying the
400  * ellipsoid.
401  * @param[in] lat1 latitude of point 1 (degrees).
402  * @param[in] lon1 longitude of point 1 (degrees).
403  * @param[in] azi1 azimuth at point 1 (degrees).
404  * @param[in] flags bitor'ed combination of geod_flags(); \e flags &
405  * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags &
406  * GEOD_LONG_UNROLL "unrolls" \e lon2.
407  * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the distance
408  * between point 1 and point 2 (meters); otherwise it is the arc length
409  * between point 1 and point 2 (degrees); it can be negative.
410  * @param[out] plat2 pointer to the latitude of point 2 (degrees).
411  * @param[out] plon2 pointer to the longitude of point 2 (degrees).
412  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
413  * @param[out] ps12 pointer to the distance between point 1 and point 2
414  * (meters).
415  * @param[out] pm12 pointer to the reduced length of geodesic (meters).
416  * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
417  * point 1 (dimensionless).
418  * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
419  * point 2 (dimensionless).
420  * @param[out] pS12 pointer to the area under the geodesic
421  * (meters<sup>2</sup>).
422  * @return \e a12 arc length of between point 1 and point 2 (degrees).
423  *
424  * \e g must have been initialized with a call to geod_init(). \e lat1
425  * should be in the range [&minus;90&deg;, 90&deg;]. The function value \e
426  * a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE. Any of the "return"
427  * arguments, \e plat2, etc., may be replaced by 0, if you do not need some
428  * quantities computed.
429  *
430  * With \e flags & GEOD_LONG_UNROLL bit set, the longitude is "unrolled" so
431  * that the quantity \e lon2 &minus; \e lon1 indicates how many times and in
432  * what sense the geodesic encircles the ellipsoid.
433  **********************************************************************/
434  double geod_gendirect(const struct geod_geodesic* g,
435  double lat1, double lon1, double azi1,
436  unsigned flags, double s12_a12,
437  double* plat2, double* plon2, double* pazi2,
438  double* ps12, double* pm12, double* pM12, double* pM21,
439  double* pS12);
440 
441  /**
442  * The general inverse geodesic calculation.
443  *
444  * @param[in] g a pointer to the geod_geodesic object specifying the
445  * ellipsoid.
446  * @param[in] lat1 latitude of point 1 (degrees).
447  * @param[in] lon1 longitude of point 1 (degrees).
448  * @param[in] lat2 latitude of point 2 (degrees).
449  * @param[in] lon2 longitude of point 2 (degrees).
450  * @param[out] ps12 pointer to the distance between point 1 and point 2
451  * (meters).
452  * @param[out] pazi1 pointer to the azimuth at point 1 (degrees).
453  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
454  * @param[out] pm12 pointer to the reduced length of geodesic (meters).
455  * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
456  * point 1 (dimensionless).
457  * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
458  * point 2 (dimensionless).
459  * @param[out] pS12 pointer to the area under the geodesic
460  * (meters<sup>2</sup>).
461  * @return \e a12 arc length of between point 1 and point 2 (degrees).
462  *
463  * \e g must have been initialized with a call to geod_init(). \e lat1 and
464  * \e lat2 should be in the range [&minus;90&deg;, 90&deg;]. Any of the
465  * "return" arguments \e ps12, etc., may be replaced by 0, if you do not need
466  * some quantities computed.
467  **********************************************************************/
468  double geod_geninverse(const struct geod_geodesic* g,
469  double lat1, double lon1, double lat2, double lon2,
470  double* ps12, double* pazi1, double* pazi2,
471  double* pm12, double* pM12, double* pM21,
472  double* pS12);
473 
474  /**
475  * The general position function.
476  *
477  * @param[in] l a pointer to the geod_geodesicline object specifying the
478  * geodesic line.
479  * @param[in] flags bitor'ed combination of geod_flags(); \e flags &
480  * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags &
481  * GEOD_LONG_UNROLL "unrolls" \e lon2; if \e flags & GEOD_ARCMODE is 0,
482  * then \e l must have been initialized with \e caps |= GEOD_DISTANCE_IN.
483  * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the
484  * distance between point 1 and point 2 (meters); otherwise it is the
485  * arc length between point 1 and point 2 (degrees); it can be
486  * negative.
487  * @param[out] plat2 pointer to the latitude of point 2 (degrees).
488  * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires
489  * that \e l was initialized with \e caps |= GEOD_LONGITUDE.
490  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
491  * @param[out] ps12 pointer to the distance between point 1 and point 2
492  * (meters); requires that \e l was initialized with \e caps |=
493  * GEOD_DISTANCE.
494  * @param[out] pm12 pointer to the reduced length of geodesic (meters);
495  * requires that \e l was initialized with \e caps |= GEOD_REDUCEDLENGTH.
496  * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
497  * point 1 (dimensionless); requires that \e l was initialized with \e caps
498  * |= GEOD_GEODESICSCALE.
499  * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
500  * point 2 (dimensionless); requires that \e l was initialized with \e caps
501  * |= GEOD_GEODESICSCALE.
502  * @param[out] pS12 pointer to the area under the geodesic
503  * (meters<sup>2</sup>); requires that \e l was initialized with \e caps |=
504  * GEOD_AREA.
505  * @return \e a12 arc length of between point 1 and point 2 (degrees).
506  *
507  * \e l must have been initialized with a call to geod_lineinit() with \e
508  * caps |= GEOD_DISTANCE_IN. The value \e azi2 returned is in the range
509  * [&minus;180&deg;, 180&deg;). Any of the "return" arguments \e plat2,
510  * etc., may be replaced by 0, if you do not need some quantities
511  * computed. Requesting a value which \e l is not capable of computing
512  * is not an error; the corresponding argument will not be altered.
513  *
514  * With \e flags & GEOD_LONG_UNROLL bit set, the longitude is "unrolled" so
515  * that the quantity \e lon2 &minus; \e lon1 indicates how many times and in
516  * what sense the geodesic encircles the ellipsoid.
517  *
518  * Example, compute way points between JFK and Singapore Changi Airport
519  * using geod_genposition(). In this example, the points are evenly space in
520  * arc length (and so only approximately equally space in distance). This is
521  * faster than using geod_position() would be appropriate if drawing the path
522  * on a map.
523  @code{.c}
524  struct geod_geodesic g;
525  struct geod_geodesicline l;
526  double a12, azi1, lat[101], lon[101];
527  int i;
528  geod_init(&g, 6378137, 1/298.257223563);
529  a12 = geod_geninverse(&g, 40.64, -73.78, 1.36, 103.99,
530  0, &azi1, 0, 0, 0, 0, 0);
531  geod_lineinit(&l, &g, 40.64, -73.78, azi1, GEOD_LATITUDE | GEOD_LONGITUDE);
532  for (i = 0; i < 101; ++i) {
533  geod_genposition(&l, 1, i * a12 * 0.01,
534  lat + i, lon + i, 0, 0, 0, 0, 0, 0);
535  printf("%.5f %.5f\n", lat[i], lon[i]);
536  }
537  @endcode
538  **********************************************************************/
539  double geod_genposition(const struct geod_geodesicline* l,
540  unsigned flags, double s12_a12,
541  double* plat2, double* plon2, double* pazi2,
542  double* ps12, double* pm12,
543  double* pM12, double* pM21,
544  double* pS12);
545 
546  /**
547  * Initialize a geod_polygon object.
548  *
549  * @param[out] p a pointer to the object to be initialized.
550  * @param[in] polylinep non-zero if a polyline instead of a polygon.
551  *
552  * If \e polylinep is zero, then the sequence of vertices and edges added by
553  * geod_polygon_addpoint() and geod_polygon_addedge() define a polygon and
554  * the perimeter and area are returned by geod_polygon_compute(). If \e
555  * polylinep is non-zero, then the vertices and edges define a polyline and
556  * only the perimeter is returned by geod_polygon_compute().
557  *
558  * The area and perimeter are accumulated at two times the standard floating
559  * point precision to guard against the loss of accuracy with many-sided
560  * polygons. At any point you can ask for the perimeter and area so far.
561  *
562  * An example of the use of this function is given in the documentation for
563  * geod_polygon_compute().
564  **********************************************************************/
565  void geod_polygon_init(struct geod_polygon* p, int polylinep);
566 
567  /**
568  * Add a point to the polygon or polyline.
569  *
570  * @param[in] g a pointer to the geod_geodesic object specifying the
571  * ellipsoid.
572  * @param[in,out] p a pointer to the geod_polygon object specifying the
573  * polygon.
574  * @param[in] lat the latitude of the point (degrees).
575  * @param[in] lon the longitude of the point (degrees).
576  *
577  * \e g and \e p must have been initialized with calls to geod_init() and
578  * geod_polygon_init(), respectively. The same \e g must be used for all the
579  * points and edges in a polygon. \e lat should be in the range
580  * [&minus;90&deg;, 90&deg;].
581  *
582  * An example of the use of this function is given in the documentation for
583  * geod_polygon_compute().
584  **********************************************************************/
585  void geod_polygon_addpoint(const struct geod_geodesic* g,
586  struct geod_polygon* p,
587  double lat, double lon);
588 
589  /**
590  * Add an edge to the polygon or polyline.
591  *
592  * @param[in] g a pointer to the geod_geodesic object specifying the
593  * ellipsoid.
594  * @param[in,out] p a pointer to the geod_polygon object specifying the
595  * polygon.
596  * @param[in] azi azimuth at current point (degrees).
597  * @param[in] s distance from current point to next point (meters).
598  *
599  * \e g and \e p must have been initialized with calls to geod_init() and
600  * geod_polygon_init(), respectively. The same \e g must be used for all the
601  * points and edges in a polygon. This does nothing if no points have been
602  * added yet. The \e lat and \e lon fields of \e p give the location of the
603  * new vertex.
604  **********************************************************************/
605  void geod_polygon_addedge(const struct geod_geodesic* g,
606  struct geod_polygon* p,
607  double azi, double s);
608 
609  /**
610  * Return the results for a polygon.
611  *
612  * @param[in] g a pointer to the geod_geodesic object specifying the
613  * ellipsoid.
614  * @param[in] p a pointer to the geod_polygon object specifying the polygon.
615  * @param[in] reverse if non-zero then clockwise (instead of
616  * counter-clockwise) traversal counts as a positive area.
617  * @param[in] sign if non-zero then return a signed result for the area if
618  * the polygon is traversed in the "wrong" direction instead of returning
619  * the area for the rest of the earth.
620  * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
621  * only set if \e polyline is non-zero in the call to geod_polygon_init().
622  * @param[out] pP pointer to the perimeter of the polygon or length of the
623  * polyline (meters).
624  * @return the number of points.
625  *
626  * The area and perimeter are accumulated at two times the standard floating
627  * point precision to guard against the loss of accuracy with many-sided
628  * polygons. Only simple polygons (which are not self-intersecting) are
629  * allowed. There's no need to "close" the polygon by repeating the first
630  * vertex. Set \e pA or \e pP to zero, if you do not want the corresponding
631  * quantity returned.
632  *
633  * Example, compute the perimeter and area of the geodesic triangle with
634  * vertices (0&deg;N,0&deg;E), (0&deg;N,90&deg;E), (90&deg;N,0&deg;E).
635  @code{.c}
636  double A, P;
637  int n;
638  struct geod_geodesic g;
639  struct geod_polygon p;
640  geod_init(&g, 6378137, 1/298.257223563);
641  geod_polygon_init(&p, 0);
642 
643  geod_polygon_addpoint(&g, &p, 0, 0);
644  geod_polygon_addpoint(&g, &p, 0, 90);
645  geod_polygon_addpoint(&g, &p, 90, 0);
646  n = geod_polygon_compute(&g, &p, 0, 1, &A, &P);
647  printf("%d %.8f %.3f\n", n, P, A);
648  @endcode
649  **********************************************************************/
650  unsigned geod_polygon_compute(const struct geod_geodesic* g,
651  const struct geod_polygon* p,
652  int reverse, int sign,
653  double* pA, double* pP);
654 
655  /**
656  * Return the results assuming a tentative final test point is added;
657  * however, the data for the test point is not saved. This lets you report a
658  * running result for the perimeter and area as the user moves the mouse
659  * cursor. Ordinary floating point arithmetic is used to accumulate the data
660  * for the test point; thus the area and perimeter returned are less accurate
661  * than if geod_polygon_addpoint() and geod_polygon_compute() are used.
662  *
663  * @param[in] g a pointer to the geod_geodesic object specifying the
664  * ellipsoid.
665  * @param[in] p a pointer to the geod_polygon object specifying the polygon.
666  * @param[in] lat the latitude of the test point (degrees).
667  * @param[in] lon the longitude of the test point (degrees).
668  * @param[in] reverse if non-zero then clockwise (instead of
669  * counter-clockwise) traversal counts as a positive area.
670  * @param[in] sign if non-zero then return a signed result for the area if
671  * the polygon is traversed in the "wrong" direction instead of returning
672  * the area for the rest of the earth.
673  * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
674  * only set if \e polyline is non-zero in the call to geod_polygon_init().
675  * @param[out] pP pointer to the perimeter of the polygon or length of the
676  * polyline (meters).
677  * @return the number of points.
678  *
679  * \e lat should be in the range [&minus;90&deg;, 90&deg;].
680  **********************************************************************/
681  unsigned geod_polygon_testpoint(const struct geod_geodesic* g,
682  const struct geod_polygon* p,
683  double lat, double lon,
684  int reverse, int sign,
685  double* pA, double* pP);
686 
687  /**
688  * Return the results assuming a tentative final test point is added via an
689  * azimuth and distance; however, the data for the test point is not saved.
690  * This lets you report a running result for the perimeter and area as the
691  * user moves the mouse cursor. Ordinary floating point arithmetic is used
692  * to accumulate the data for the test point; thus the area and perimeter
693  * returned are less accurate than if geod_polygon_addedge() and
694  * geod_polygon_compute() are used.
695  *
696  * @param[in] g a pointer to the geod_geodesic object specifying the
697  * ellipsoid.
698  * @param[in] p a pointer to the geod_polygon object specifying the polygon.
699  * @param[in] azi azimuth at current point (degrees).
700  * @param[in] s distance from current point to final test point (meters).
701  * @param[in] reverse if non-zero then clockwise (instead of
702  * counter-clockwise) traversal counts as a positive area.
703  * @param[in] sign if non-zero then return a signed result for the area if
704  * the polygon is traversed in the "wrong" direction instead of returning
705  * the area for the rest of the earth.
706  * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
707  * only set if \e polyline is non-zero in the call to geod_polygon_init().
708  * @param[out] pP pointer to the perimeter of the polygon or length of the
709  * polyline (meters).
710  * @return the number of points.
711  **********************************************************************/
712  unsigned geod_polygon_testedge(const struct geod_geodesic* g,
713  const struct geod_polygon* p,
714  double azi, double s,
715  int reverse, int sign,
716  double* pA, double* pP);
717 
718  /**
719  * A simple interface for computing the area of a geodesic polygon.
720  *
721  * @param[in] g a pointer to the geod_geodesic object specifying the
722  * ellipsoid.
723  * @param[in] lats an array of latitudes of the polygon vertices (degrees).
724  * @param[in] lons an array of longitudes of the polygon vertices (degrees).
725  * @param[in] n the number of vertices.
726  * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>).
727  * @param[out] pP pointer to the perimeter of the polygon (meters).
728  *
729  * \e lats should be in the range [&minus;90&deg;, 90&deg;].
730  *
731  * Only simple polygons (which are not self-intersecting) are allowed.
732  * There's no need to "close" the polygon by repeating the first vertex. The
733  * area returned is signed with counter-clockwise traversal being treated as
734  * positive.
735  *
736  * Example, compute the area of Antarctica:
737  @code{.c}
738  double
739  lats[] = {-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
740  -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7},
741  lons[] = {-74, -102, -102, -131, -163, 163, 172, 140, 113,
742  88, 59, 25, -4, -14, -33, -46, -61};
743  struct geod_geodesic g;
744  double A, P;
745  geod_init(&g, 6378137, 1/298.257223563);
746  geod_polygonarea(&g, lats, lons, (sizeof lats) / (sizeof lats[0]), &A, &P);
747  printf("%.0f %.2f\n", A, P);
748  @endcode
749  **********************************************************************/
750  void geod_polygonarea(const struct geod_geodesic* g,
751  double lats[], double lons[], int n,
752  double* pA, double* pP);
753 
754  /**
755  * mask values for the \e caps argument to geod_lineinit().
756  **********************************************************************/
757  enum geod_mask {
758  GEOD_NONE = 0U, /**< Calculate nothing */
759  GEOD_LATITUDE = 1U<<7 | 0U, /**< Calculate latitude */
760  GEOD_LONGITUDE = 1U<<8 | 1U<<3, /**< Calculate longitude */
761  GEOD_AZIMUTH = 1U<<9 | 0U, /**< Calculate azimuth */
762  GEOD_DISTANCE = 1U<<10 | 1U<<0, /**< Calculate distance */
763  GEOD_DISTANCE_IN = 1U<<11 | 1U<<0 | 1U<<1, /**< Allow distance as input */
764  GEOD_REDUCEDLENGTH= 1U<<12 | 1U<<0 | 1U<<2, /**< Calculate reduced length */
765  GEOD_GEODESICSCALE= 1U<<13 | 1U<<0 | 1U<<2, /**< Calculate geodesic scale */
766  GEOD_AREA = 1U<<14 | 1U<<4, /**< Calculate reduced length */
767  GEOD_ALL = 0x7F80U| 0x1FU /**< Calculate everything */
768  };
769 
770  /**
771  * flag values for the \e flags argument to geod_gendirect() and
772  * geod_genposition()
773  **********************************************************************/
774  enum geod_flags {
775  GEOD_NOFLAGS = 0U, /**< No flags */
776  GEOD_ARCMODE = 1U<<0, /**< Position given in terms of arc distance */
777  GEOD_LONG_UNROLL = 1U<<15, /**< Unroll the longitude */
778  /**< @cond SKIP */
779  GEOD_LONG_NOWRAP = GEOD_LONG_UNROLL /* For backward compatibility only */
780  /**< @endcond */
781  };
782 
783 #if defined(__cplusplus)
784 }
785 #endif
786 
787 #endif
unsigned 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)
double geod_genposition(const struct geod_geodesicline *l, unsigned flags, double s12_a12, double *plat2, double *plon2, double *pazi2, double *ps12, double *pm12, double *pM12, double *pM21, double *pS12)
double lon
Definition: geodesic.h:204
void geod_polygon_addedge(const struct geod_geodesic *g, struct geod_polygon *p, double azi, double s)
unsigned num
Definition: geodesic.h:213
void geod_position(const struct geod_geodesicline *l, double s12, double *plat2, double *plon2, double *pazi2)
double f
Definition: geodesic.h:171
void geod_lineinit(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned caps)
unsigned caps
Definition: geodesic.h:194
double geod_geninverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2, double *pm12, double *pM12, double *pM21, double *pS12)
void geod_polygon_addpoint(const struct geod_geodesic *g, struct geod_polygon *p, double lat, double lon)
void geod_polygon_init(struct geod_polygon *p, int polylinep)
void geod_direct(const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, double *plat2, double *plon2, double *pazi2)
unsigned geod_polygon_compute(const struct geod_geodesic *g, const struct geod_polygon *p, int reverse, int sign, double *pA, double *pP)
void geod_polygonarea(const struct geod_geodesic *g, double lats[], double lons[], int n, double *pA, double *pP)
double a
Definition: geodesic.h:170
double geod_gendirect(const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned flags, double s12_a12, double *plat2, double *plon2, double *pazi2, double *ps12, double *pm12, double *pM12, double *pM21, double *pS12)
geod_flags
Definition: geodesic.h:774
geod_mask
Definition: geodesic.h:757
unsigned 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_inverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2)
void geod_init(struct geod_geodesic *g, double a, double f)
double lat
Definition: geodesic.h:203