C library for Geodesics  1.43
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.43.
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 43
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
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;]; \e lon1 and \e azi1
241  * should be in the range [&minus;540&deg;, 540&deg;).
242  *
243  * The geod_mask values are [see geod_mask()]:
244  * - \e caps |= GEOD_LATITUDE for the latitude \e lat2; this is
245  * added automatically,
246  * - \e caps |= GEOD_LONGITUDE for the latitude \e lon2,
247  * - \e caps |= GEOD_AZIMUTH for the latitude \e azi2; this is
248  * added automatically,
249  * - \e caps |= GEOD_DISTANCE for the distance \e s12,
250  * - \e caps |= GEOD_REDUCEDLENGTH for the reduced length \e m12,
251  * - \e caps |= GEOD_GEODESICSCALE for the geodesic scales \e M12
252  * and \e M21,
253  * - \e caps |= GEOD_AREA for the area \e S12,
254  * - \e caps |= GEOD_DISTANCE_IN permits the length of the
255  * geodesic to be given in terms of \e s12; without this capability the
256  * length can only be specified in terms of arc length.
257  * .
258  * A value of \e caps = 0 is treated as GEOD_LATITUDE | GEOD_LONGITUDE |
259  * GEOD_AZIMUTH | GEOD_DISTANCE_IN (to support the solution of the "standard"
260  * direct problem).
261  **********************************************************************/
262  void geod_lineinit(struct geod_geodesicline* l,
263  const struct geod_geodesic* g,
264  double lat1, double lon1, double azi1, unsigned caps);
265 
266  /**
267  * Solve the direct geodesic problem.
268  *
269  * @param[in] g a pointer to the geod_geodesic object specifying the
270  * ellipsoid.
271  * @param[in] lat1 latitude of point 1 (degrees).
272  * @param[in] lon1 longitude of point 1 (degrees).
273  * @param[in] azi1 azimuth at point 1 (degrees).
274  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
275  * negative.
276  * @param[out] plat2 pointer to the latitude of point 2 (degrees).
277  * @param[out] plon2 pointer to the longitude of point 2 (degrees).
278  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
279  *
280  * \e g must have been initialized with a call to geod_init(). \e lat1
281  * should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e azi1
282  * should be in the range [&minus;540&deg;, 540&deg;). The values of \e lon2
283  * and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;). Any of
284  * the "return" arguments \e plat2, etc., may be replaced by 0, if you do not
285  * need some quantities computed.
286  *
287  * If either point is at a pole, the azimuth is defined by keeping the
288  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;), and
289  * taking the limit &epsilon; &rarr; 0+. An arc length greater that 180&deg;
290  * signifies a geodesic which is not a shortest path. (For a prolate
291  * ellipsoid, an additional condition is necessary for a shortest path: the
292  * longitudinal extent must not exceed of 180&deg;.)
293  *
294  * Example, determine the point 10000 km NE of JFK:
295  @code
296  struct geod_geodesic g;
297  double lat, lon;
298  geod_init(&g, 6378137, 1/298.257223563);
299  geod_direct(&g, 40.64, -73.78, 45.0, 10e6, &lat, &lon, 0);
300  printf("%.5f %.5f\n", lat, lon);
301  @endcode
302  **********************************************************************/
303  void geod_direct(const struct geod_geodesic* g,
304  double lat1, double lon1, double azi1, double s12,
305  double* plat2, double* plon2, double* pazi2);
306 
307  /**
308  * Solve the inverse geodesic problem.
309  *
310  * @param[in] g a pointer to the geod_geodesic object specifying the
311  * ellipsoid.
312  * @param[in] lat1 latitude of point 1 (degrees).
313  * @param[in] lon1 longitude of point 1 (degrees).
314  * @param[in] lat2 latitude of point 2 (degrees).
315  * @param[in] lon2 longitude of point 2 (degrees).
316  * @param[out] ps12 pointer to the distance between point 1 and point 2
317  * (meters).
318  * @param[out] pazi1 pointer to the azimuth at point 1 (degrees).
319  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
320  *
321  * \e g must have been initialized with a call to geod_init(). \e lat1
322  * and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and
323  * \e lon2 should be in the range [&minus;540&deg;, 540&deg;). The values of
324  * \e azi1 and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;).
325  * Any of the "return" arguments \e ps12, etc., may be replaced by 0, if you
326  * do not need some quantities computed.
327  *
328  * If either point is at a pole, the azimuth is defined by keeping the
329  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;), and
330  * taking the limit &epsilon; &rarr; 0+.
331  *
332  * The solution to the inverse problem is found using Newton's method. If
333  * this fails to converge (this is very unlikely in geodetic applications
334  * but does occur for very eccentric ellipsoids), then the bisection method
335  * is used to refine the solution.
336  *
337  * Example, determine the distance between JFK and Singapore Changi Airport:
338  @code
339  struct geod_geodesic g;
340  double s12;
341  geod_init(&g, 6378137, 1/298.257223563);
342  geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, 0, 0);
343  printf("%.3f\n", s12);
344  @endcode
345  **********************************************************************/
346  void geod_inverse(const struct geod_geodesic* g,
347  double lat1, double lon1, double lat2, double lon2,
348  double* ps12, double* pazi1, double* pazi2);
349 
350  /**
351  * Compute the position along a geod_geodesicline.
352  *
353  * @param[in] l a pointer to the geod_geodesicline object specifying the
354  * geodesic line.
355  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
356  * negative.
357  * @param[out] plat2 pointer to the latitude of point 2 (degrees).
358  * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires
359  * that \e l was initialized with \e caps |= GEOD_LONGITUDE.
360  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
361  *
362  * \e l must have been initialized with a call to geod_lineinit() with \e
363  * caps |= GEOD_DISTANCE_IN. The values of \e lon2 and \e azi2 returned are
364  * in the range [&minus;180&deg;, 180&deg;). Any of the "return" arguments
365  * \e plat2, etc., may be replaced by 0, if you do not need some quantities
366  * computed.
367  *
368  * Example, compute way points between JFK and Singapore Changi Airport
369  * the "obvious" way using geod_direct():
370  @code
371  struct geod_geodesic g;
372  double s12, azi1, lat[101],lon[101];
373  int i;
374  geod_init(&g, 6378137, 1/298.257223563);
375  geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0);
376  for (i = 0; i < 101; ++i) {
377  geod_direct(&g, 40.64, -73.78, azi1, i * s12 * 0.01, lat + i, lon + i, 0);
378  printf("%.5f %.5f\n", lat[i], lon[i]);
379  }
380  @endcode
381  * A faster way using geod_position():
382  @code
383  struct geod_geodesic g;
384  struct geod_geodesicline l;
385  double s12, azi1, lat[101],lon[101];
386  int i;
387  geod_init(&g, 6378137, 1/298.257223563);
388  geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0);
389  geod_lineinit(&l, &g, 40.64, -73.78, azi1, 0);
390  for (i = 0; i < 101; ++i) {
391  geod_position(&l, i * s12 * 0.01, lat + i, lon + i, 0);
392  printf("%.5f %.5f\n", lat[i], lon[i]);
393  }
394  @endcode
395  **********************************************************************/
396  void geod_position(const struct geod_geodesicline* l, double s12,
397  double* plat2, double* plon2, double* pazi2);
398 
399  /**
400  * The general direct geodesic problem.
401  *
402  * @param[in] g a pointer to the geod_geodesic object specifying the
403  * ellipsoid.
404  * @param[in] lat1 latitude of point 1 (degrees).
405  * @param[in] lon1 longitude of point 1 (degrees).
406  * @param[in] azi1 azimuth at point 1 (degrees).
407  * @param[in] flags bitor'ed combination of geod_flags(); \e flags &
408  * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags &
409  * GEOD_LONG_UNROLL "unrolls" \e lon2.
410  * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the distance
411  * between point 1 and point 2 (meters); otherwise it is the arc length
412  * between point 1 and point 2 (degrees); it can be negative.
413  * @param[out] plat2 pointer to the latitude of point 2 (degrees).
414  * @param[out] plon2 pointer to the longitude of point 2 (degrees).
415  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
416  * @param[out] ps12 pointer to the distance between point 1 and point 2
417  * (meters).
418  * @param[out] pm12 pointer to the reduced length of geodesic (meters).
419  * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
420  * point 1 (dimensionless).
421  * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
422  * point 2 (dimensionless).
423  * @param[out] pS12 pointer to the area under the geodesic
424  * (meters<sup>2</sup>).
425  * @return \e a12 arc length of between point 1 and point 2 (degrees).
426  *
427  * \e g must have been initialized with a call to geod_init(). \e lat1
428  * should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e azi1
429  * should be in the range [&minus;540&deg;, 540&deg;). The function
430  * value \e a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE. Any of the
431  * "return" arguments \e plat2, etc., may be replaced by 0, if you do not
432  * need some quantities computed.
433  *
434  * With \e flags & GEOD_LONG_UNROLL bit set, the longitude is "unrolled" so
435  * that the quantity \e lon2 &minus; \e lon1 indicates how many times and in
436  * what sense the geodesic encircles the ellipsoid. Because \e lon2 might be
437  * outside the normal allowed range for longitudes, [&minus;540&deg;,
438  * 540&deg;), be sure to normalize it, e.g., with fmod(\e lon2, 360.0) before
439  * using it in subsequent calculations
440  **********************************************************************/
441  double geod_gendirect(const struct geod_geodesic* g,
442  double lat1, double lon1, double azi1,
443  unsigned flags, double s12_a12,
444  double* plat2, double* plon2, double* pazi2,
445  double* ps12, double* pm12, double* pM12, double* pM21,
446  double* pS12);
447 
448  /**
449  * The general inverse geodesic calculation.
450  *
451  * @param[in] g a pointer to the geod_geodesic object specifying the
452  * ellipsoid.
453  * @param[in] lat1 latitude of point 1 (degrees).
454  * @param[in] lon1 longitude of point 1 (degrees).
455  * @param[in] lat2 latitude of point 2 (degrees).
456  * @param[in] lon2 longitude of point 2 (degrees).
457  * @param[out] ps12 pointer to the distance between point 1 and point 2
458  * (meters).
459  * @param[out] pazi1 pointer to the azimuth at point 1 (degrees).
460  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
461  * @param[out] pm12 pointer to the reduced length of geodesic (meters).
462  * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
463  * point 1 (dimensionless).
464  * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
465  * point 2 (dimensionless).
466  * @param[out] pS12 pointer to the area under the geodesic
467  * (meters<sup>2</sup>).
468  * @return \e a12 arc length of between point 1 and point 2 (degrees).
469  *
470  * \e g must have been initialized with a call to geod_init(). \e lat1
471  * and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and
472  * \e lon2 should be in the range [&minus;540&deg;, 540&deg;). Any of the
473  * "return" arguments \e ps12, etc., may be replaced by 0, if you do not need
474  * some quantities computed.
475  **********************************************************************/
476  double geod_geninverse(const struct geod_geodesic* g,
477  double lat1, double lon1, double lat2, double lon2,
478  double* ps12, double* pazi1, double* pazi2,
479  double* pm12, double* pM12, double* pM21,
480  double* pS12);
481 
482  /**
483  * The general position function.
484  *
485  * @param[in] l a pointer to the geod_geodesicline object specifying the
486  * geodesic line.
487  * @param[in] flags bitor'ed combination of geod_flags(); \e flags &
488  * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags &
489  * GEOD_LONG_UNROLL "unrolls" \e lon2; if \e flags & GEOD_ARCMODE is 0,
490  * then \e l must have been initialized with \e caps |= GEOD_DISTANCE_IN.
491  * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the
492  * distance between point 1 and point 2 (meters); otherwise it is the
493  * arc length between point 1 and point 2 (degrees); it can be
494  * negative.
495  * @param[out] plat2 pointer to the latitude of point 2 (degrees).
496  * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires
497  * that \e l was initialized with \e caps |= GEOD_LONGITUDE.
498  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
499  * @param[out] ps12 pointer to the distance between point 1 and point 2
500  * (meters); requires that \e l was initialized with \e caps |=
501  * GEOD_DISTANCE.
502  * @param[out] pm12 pointer to the reduced length of geodesic (meters);
503  * requires that \e l was initialized with \e caps |= GEOD_REDUCEDLENGTH.
504  * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
505  * point 1 (dimensionless); requires that \e l was initialized with \e caps
506  * |= GEOD_GEODESICSCALE.
507  * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
508  * point 2 (dimensionless); requires that \e l was initialized with \e caps
509  * |= GEOD_GEODESICSCALE.
510  * @param[out] pS12 pointer to the area under the geodesic
511  * (meters<sup>2</sup>); requires that \e l was initialized with \e caps |=
512  * GEOD_AREA.
513  * @return \e a12 arc length of between point 1 and point 2 (degrees).
514  *
515  * \e l must have been initialized with a call to geod_lineinit() with \e
516  * caps |= GEOD_DISTANCE_IN. The value \e azi2 returned is in the range
517  * [&minus;180&deg;, 180&deg;). Any of the "return" arguments \e plat2,
518  * etc., may be replaced by 0, if you do not need some quantities
519  * computed. Requesting a value which \e l is not capable of computing
520  * is not an error; the corresponding argument will not be altered.
521  *
522  * With \e flags & GEOD_LONG_UNROLL bit set, the longitude is "unrolled" so
523  * that the quantity \e lon2 &minus; \e lon1 indicates how many times and in
524  * what sense the geodesic encircles the ellipsoid. Because \e lon2 might be
525  * outside the normal allowed range for longitudes, [&minus;540&deg;,
526  * 540&deg;), be sure to normalize it, e.g., with fmod(\e lon2, 360.0) before
527  * using it in subsequent calculations
528  *
529  * Example, compute way points between JFK and Singapore Changi Airport
530  * using geod_genposition(). In this example, the points are evenly space in
531  * arc length (and so only approximately equally space in distance). This is
532  * faster than using geod_position() would be appropriate if drawing the path
533  * on a map.
534  @code
535  struct geod_geodesic g;
536  struct geod_geodesicline l;
537  double a12, azi1, lat[101], lon[101];
538  int i;
539  geod_init(&g, 6378137, 1/298.257223563);
540  a12 = geod_geninverse(&g, 40.64, -73.78, 1.36, 103.99,
541  0, &azi1, 0, 0, 0, 0, 0);
542  geod_lineinit(&l, &g, 40.64, -73.78, azi1, GEOD_LATITUDE | GEOD_LONGITUDE);
543  for (i = 0; i < 101; ++i) {
544  geod_genposition(&l, 1, i * a12 * 0.01,
545  lat + i, lon + i, 0, 0, 0, 0, 0, 0);
546  printf("%.5f %.5f\n", lat[i], lon[i]);
547  }
548  @endcode
549  **********************************************************************/
550  double geod_genposition(const struct geod_geodesicline* l,
551  unsigned flags, double s12_a12,
552  double* plat2, double* plon2, double* pazi2,
553  double* ps12, double* pm12,
554  double* pM12, double* pM21,
555  double* pS12);
556 
557  /**
558  * Initialize a geod_polygon object.
559  *
560  * @param[out] p a pointer to the object to be initialized.
561  * @param[in] polylinep non-zero if a polyline instead of a polygon.
562  *
563  * If \e polylinep is zero, then the sequence of vertices and edges added by
564  * geod_polygon_addpoint() and geod_polygon_addedge() define a polygon and
565  * the perimeter and area are returned by geod_polygon_compute(). If \e
566  * polylinep is non-zero, then the vertices and edges define a polyline and
567  * only the perimeter is returned by geod_polygon_compute().
568  *
569  * The area and perimeter are accumulated at two times the standard floating
570  * point precision to guard against the loss of accuracy with many-sided
571  * polygons. At any point you can ask for the perimeter and area so far.
572  *
573  * An example of the use of this function is given in the documentation for
574  * geod_polygon_compute().
575  **********************************************************************/
576  void geod_polygon_init(struct geod_polygon* p, int polylinep);
577 
578  /**
579  * Add a point to the polygon or polyline.
580  *
581  * @param[in] g a pointer to the geod_geodesic object specifying the
582  * ellipsoid.
583  * @param[in,out] p a pointer to the geod_polygon object specifying the
584  * polygon.
585  * @param[in] lat the latitude of the point (degrees).
586  * @param[in] lon the longitude of the point (degrees).
587  *
588  * \e g and \e p must have been initialized with calls to geod_init() and
589  * geod_polygon_init(), respectively. The same \e g must be used for all the
590  * points and edges in a polygon. \e lat should be in the range
591  * [&minus;90&deg;, 90&deg;] and \e lon should be in the range
592  * [&minus;540&deg;, 540&deg;).
593  *
594  * An example of the use of this function is given in the documentation for
595  * geod_polygon_compute().
596  **********************************************************************/
597  void geod_polygon_addpoint(const struct geod_geodesic* g,
598  struct geod_polygon* p,
599  double lat, double lon);
600 
601  /**
602  * Add an edge to the polygon or polyline.
603  *
604  * @param[in] g a pointer to the geod_geodesic object specifying the
605  * ellipsoid.
606  * @param[in,out] p a pointer to the geod_polygon object specifying the
607  * polygon.
608  * @param[in] azi azimuth at current point (degrees).
609  * @param[in] s distance from current point to next point (meters).
610  *
611  * \e g and \e p must have been initialized with calls to geod_init() and
612  * geod_polygon_init(), respectively. The same \e g must be used for all the
613  * points and edges in a polygon. \e azi should be in the range
614  * [&minus;540&deg;, 540&deg;). This does nothing if no points have been
615  * added yet. The \e lat and \e lon fields of \e p give the location of
616  * the new vertex.
617  **********************************************************************/
618  void geod_polygon_addedge(const struct geod_geodesic* g,
619  struct geod_polygon* p,
620  double azi, double s);
621 
622  /**
623  * Return the results for a polygon.
624  *
625  * @param[in] g a pointer to the geod_geodesic object specifying the
626  * ellipsoid.
627  * @param[in] p a pointer to the geod_polygon object specifying the polygon.
628  * @param[in] reverse if non-zero then clockwise (instead of
629  * counter-clockwise) traversal counts as a positive area.
630  * @param[in] sign if non-zero then return a signed result for the area if
631  * the polygon is traversed in the "wrong" direction instead of returning
632  * the area for the rest of the earth.
633  * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
634  * only set if \e polyline is non-zero in the call to geod_polygon_init().
635  * @param[out] pP pointer to the perimeter of the polygon or length of the
636  * polyline (meters).
637  * @return the number of points.
638  *
639  * The area and perimeter are accumulated at two times the standard floating
640  * point precision to guard against the loss of accuracy with many-sided
641  * polygons. Only simple polygons (which are not self-intersecting) are
642  * allowed. There's no need to "close" the polygon by repeating the first
643  * vertex. Set \e pA or \e pP to zero, if you do not want the corresponding
644  * quantity returned.
645  *
646  * Example, compute the perimeter and area of the geodesic triangle with
647  * vertices (0&deg;N,0&deg;E), (0&deg;N,90&deg;E), (90&deg;N,0&deg;E).
648  @code
649  double A, P;
650  int n;
651  struct geod_geodesic g;
652  struct geod_polygon p;
653  geod_init(&g, 6378137, 1/298.257223563);
654  geod_polygon_init(&p, 0);
655 
656  geod_polygon_addpoint(&g, &p, 0, 0);
657  geod_polygon_addpoint(&g, &p, 0, 90);
658  geod_polygon_addpoint(&g, &p, 90, 0);
659  n = geod_polygon_compute(&g, &p, 0, 1, &A, &P);
660  printf("%d %.8f %.3f\n", n, P, A);
661  @endcode
662  **********************************************************************/
663  unsigned geod_polygon_compute(const struct geod_geodesic* g,
664  const struct geod_polygon* p,
665  int reverse, int sign,
666  double* pA, double* pP);
667 
668  /**
669  * Return the results assuming a tentative final test point is added;
670  * however, the data for the test point is not saved. This lets you report a
671  * running result for the perimeter and area as the user moves the mouse
672  * cursor. Ordinary floating point arithmetic is used to accumulate the data
673  * for the test point; thus the area and perimeter returned are less accurate
674  * than if geod_polygon_addpoint() and geod_polygon_compute() are used.
675  *
676  * @param[in] g a pointer to the geod_geodesic object specifying the
677  * ellipsoid.
678  * @param[in] p a pointer to the geod_polygon object specifying the polygon.
679  * @param[in] lat the latitude of the test point (degrees).
680  * @param[in] lon the longitude of the test point (degrees).
681  * @param[in] reverse if non-zero then clockwise (instead of
682  * counter-clockwise) traversal counts as a positive area.
683  * @param[in] sign if non-zero then return a signed result for the area if
684  * the polygon is traversed in the "wrong" direction instead of returning
685  * the area for the rest of the earth.
686  * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
687  * only set if \e polyline is non-zero in the call to geod_polygon_init().
688  * @param[out] pP pointer to the perimeter of the polygon or length of the
689  * polyline (meters).
690  * @return the number of points.
691  *
692  * \e lat should be in the range [&minus;90&deg;, 90&deg;] and \e
693  * lon should be in the range [&minus;540&deg;, 540&deg;).
694  **********************************************************************/
695  unsigned geod_polygon_testpoint(const struct geod_geodesic* g,
696  const struct geod_polygon* p,
697  double lat, double lon,
698  int reverse, int sign,
699  double* pA, double* pP);
700 
701  /**
702  * Return the results assuming a tentative final test point is added via an
703  * azimuth and distance; however, the data for the test point is not saved.
704  * This lets you report a running result for the perimeter and area as the
705  * user moves the mouse cursor. Ordinary floating point arithmetic is used
706  * to accumulate the data for the test point; thus the area and perimeter
707  * returned are less accurate than if geod_polygon_addedge() and
708  * geod_polygon_compute() are used.
709  *
710  * @param[in] g a pointer to the geod_geodesic object specifying the
711  * ellipsoid.
712  * @param[in] p a pointer to the geod_polygon object specifying the polygon.
713  * @param[in] azi azimuth at current point (degrees).
714  * @param[in] s distance from current point to final test point (meters).
715  * @param[in] reverse if non-zero then clockwise (instead of
716  * counter-clockwise) traversal counts as a positive area.
717  * @param[in] sign if non-zero then return a signed result for the area if
718  * the polygon is traversed in the "wrong" direction instead of returning
719  * the area for the rest of the earth.
720  * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
721  * only set if \e polyline is non-zero in the call to geod_polygon_init().
722  * @param[out] pP pointer to the perimeter of the polygon or length of the
723  * polyline (meters).
724  * @return the number of points.
725  *
726  * \e azi should be in the range [&minus;540&deg;, 540&deg;).
727  **********************************************************************/
728  unsigned geod_polygon_testedge(const struct geod_geodesic* g,
729  const struct geod_polygon* p,
730  double azi, double s,
731  int reverse, int sign,
732  double* pA, double* pP);
733 
734  /**
735  * A simple interface for computing the area of a geodesic polygon.
736  *
737  * @param[in] g a pointer to the geod_geodesic object specifying the
738  * ellipsoid.
739  * @param[in] lats an array of latitudes of the polygon vertices (degrees).
740  * @param[in] lons an array of longitudes of the polygon vertices (degrees).
741  * @param[in] n the number of vertices.
742  * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>).
743  * @param[out] pP pointer to the perimeter of the polygon (meters).
744  *
745  * \e lats should be in the range [&minus;90&deg;, 90&deg;]; \e lons should
746  * be in the range [&minus;540&deg;, 540&deg;).
747  *
748  * Only simple polygons (which are not self-intersecting) are allowed.
749  * There's no need to "close" the polygon by repeating the first vertex. The
750  * area returned is signed with counter-clockwise traversal being treated as
751  * positive.
752  *
753  * Example, compute the area of Antarctica:
754  @code
755  double
756  lats[] = {-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
757  -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7},
758  lons[] = {-74, -102, -102, -131, -163, 163, 172, 140, 113,
759  88, 59, 25, -4, -14, -33, -46, -61};
760  struct geod_geodesic g;
761  double A, P;
762  geod_init(&g, 6378137, 1/298.257223563);
763  geod_polygonarea(&g, lats, lons, (sizeof lats) / (sizeof lats[0]), &A, &P);
764  printf("%.0f %.2f\n", A, P);
765  @endcode
766  **********************************************************************/
767  void geod_polygonarea(const struct geod_geodesic* g,
768  double lats[], double lons[], int n,
769  double* pA, double* pP);
770 
771  /**
772  * mask values for the \e caps argument to geod_lineinit().
773  **********************************************************************/
774  enum geod_mask {
775  GEOD_NONE = 0U, /**< Calculate nothing */
776  GEOD_LATITUDE = 1U<<7 | 0U, /**< Calculate latitude */
777  GEOD_LONGITUDE = 1U<<8 | 1U<<3, /**< Calculate longitude */
778  GEOD_AZIMUTH = 1U<<9 | 0U, /**< Calculate azimuth */
779  GEOD_DISTANCE = 1U<<10 | 1U<<0, /**< Calculate distance */
780  GEOD_DISTANCE_IN = 1U<<11 | 1U<<0 | 1U<<1, /**< Allow distance as input */
781  GEOD_REDUCEDLENGTH= 1U<<12 | 1U<<0 | 1U<<2, /**< Calculate reduced length */
782  GEOD_GEODESICSCALE= 1U<<13 | 1U<<0 | 1U<<2, /**< Calculate geodesic scale */
783  GEOD_AREA = 1U<<14 | 1U<<4, /**< Calculate reduced length */
784  GEOD_ALL = 0x7F80U| 0x1FU /**< Calculate everything */
785  };
786 
787  /**
788  * flag values for the \e flags argument to geod_gendirect() and
789  * geod_genposition()
790  **********************************************************************/
791  enum geod_flags {
792  GEOD_NOFLAGS = 0U, /**< No flags */
793  GEOD_ARCMODE = 1U<<0, /**< Position given in terms of arc distance */
794  GEOD_LONG_UNROLL = 1U<<15, /**< Unroll the longitude */
795  /**< @cond SKIP */
796  GEOD_LONG_NOWRAP = GEOD_LONG_UNROLL /* For backward compatibility only */
797  /**< @endcond */
798  };
799 
800 #if defined(__cplusplus)
801 }
802 #endif
803 
804 #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:791
geod_mask
Definition: geodesic.h:774
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