C library for Geodesics  1.42
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-2014) <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.42.
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 42
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_NOWRAP prevents the value of \e lon2 being wrapped into
410  * the range [&minus;180&deg;, 180&deg;).
411  * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the distance
412  * between point 1 and point 2 (meters); otherwise it is the arc length
413  * between point 1 and point 2 (degrees); it can be negative.
414  * @param[out] plat2 pointer to the latitude of point 2 (degrees).
415  * @param[out] plon2 pointer to the longitude of point 2 (degrees).
416  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
417  * @param[out] ps12 pointer to the distance between point 1 and point 2
418  * (meters).
419  * @param[out] pm12 pointer to the reduced length of geodesic (meters).
420  * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
421  * point 1 (dimensionless).
422  * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
423  * point 2 (dimensionless).
424  * @param[out] pS12 pointer to the area under the geodesic
425  * (meters<sup>2</sup>).
426  * @return \e a12 arc length of between point 1 and point 2 (degrees).
427  *
428  * \e g must have been initialized with a call to geod_init(). \e lat1
429  * should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e azi1
430  * should be in the range [&minus;540&deg;, 540&deg;). The function
431  * value \e a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE. Any of the
432  * "return" arguments \e plat2, etc., may be replaced by 0, if you do not
433  * need some quantities computed.
434  *
435  * With \e flags & GEOD_LONG_NOWRAP bit set, the quantity \e lon2 &minus;
436  * \e lon1 indicates how many times the geodesic wrapped around the
437  * ellipsoid. Because \e lon2 might be outside the normal allowed range
438  * for longitudes, [&minus;540&deg;, 540&deg;), be sure to normalize it,
439  * e.g., with fmod(\e lon2, 360.0) before using it in subsequent
440  * calculations
441  **********************************************************************/
442  double geod_gendirect(const struct geod_geodesic* g,
443  double lat1, double lon1, double azi1,
444  unsigned flags, double s12_a12,
445  double* plat2, double* plon2, double* pazi2,
446  double* ps12, double* pm12, double* pM12, double* pM21,
447  double* pS12);
448 
449  /**
450  * The general inverse geodesic calculation.
451  *
452  * @param[in] g a pointer to the geod_geodesic object specifying the
453  * ellipsoid.
454  * @param[in] lat1 latitude of point 1 (degrees).
455  * @param[in] lon1 longitude of point 1 (degrees).
456  * @param[in] lat2 latitude of point 2 (degrees).
457  * @param[in] lon2 longitude of point 2 (degrees).
458  * @param[out] ps12 pointer to the distance between point 1 and point 2
459  * (meters).
460  * @param[out] pazi1 pointer to the azimuth at point 1 (degrees).
461  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
462  * @param[out] pm12 pointer to the reduced length of geodesic (meters).
463  * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
464  * point 1 (dimensionless).
465  * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
466  * point 2 (dimensionless).
467  * @param[out] pS12 pointer to the area under the geodesic
468  * (meters<sup>2</sup>).
469  * @return \e a12 arc length of between point 1 and point 2 (degrees).
470  *
471  * \e g must have been initialized with a call to geod_init(). \e lat1
472  * and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and
473  * \e lon2 should be in the range [&minus;540&deg;, 540&deg;). Any of the
474  * "return" arguments \e ps12, etc., may be replaced by 0, if you do not need
475  * some quantities computed.
476  **********************************************************************/
477  double geod_geninverse(const struct geod_geodesic* g,
478  double lat1, double lon1, double lat2, double lon2,
479  double* ps12, double* pazi1, double* pazi2,
480  double* pm12, double* pM12, double* pM21,
481  double* pS12);
482 
483  /**
484  * The general position function.
485  *
486  * @param[in] l a pointer to the geod_geodesicline object specifying the
487  * geodesic line.
488  * @param[in] flags bitor'ed combination of geod_flags(); \e flags &
489  * GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags &
490  * GEOD_LONG_NOWRAP prevents the value of \e lon2 being wrapped into
491  * the range [&minus;180&deg;, 180&deg;); if \e flags & GEOD_ARCMODE is
492  * 0, then \e l must have been initialized with \e caps |=
493  * GEOD_DISTANCE_IN.
494  * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the
495  * distance between point 1 and point 2 (meters); otherwise it is the
496  * arc length between point 1 and point 2 (degrees); it can be
497  * negative.
498  * @param[out] plat2 pointer to the latitude of point 2 (degrees).
499  * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires
500  * that \e l was initialized with \e caps |= GEOD_LONGITUDE.
501  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
502  * @param[out] ps12 pointer to the distance between point 1 and point 2
503  * (meters); requires that \e l was initialized with \e caps |=
504  * GEOD_DISTANCE.
505  * @param[out] pm12 pointer to the reduced length of geodesic (meters);
506  * requires that \e l was initialized with \e caps |= GEOD_REDUCEDLENGTH.
507  * @param[out] pM12 pointer to the geodesic scale of point 2 relative to
508  * point 1 (dimensionless); requires that \e l was initialized with \e caps
509  * |= GEOD_GEODESICSCALE.
510  * @param[out] pM21 pointer to the geodesic scale of point 1 relative to
511  * point 2 (dimensionless); requires that \e l was initialized with \e caps
512  * |= GEOD_GEODESICSCALE.
513  * @param[out] pS12 pointer to the area under the geodesic
514  * (meters<sup>2</sup>); requires that \e l was initialized with \e caps |=
515  * GEOD_AREA.
516  * @return \e a12 arc length of between point 1 and point 2 (degrees).
517  *
518  * \e l must have been initialized with a call to geod_lineinit() with \e
519  * caps |= GEOD_DISTANCE_IN. The value \e azi2 returned is in the range
520  * [&minus;180&deg;, 180&deg;). Any of the "return" arguments \e plat2,
521  * etc., may be replaced by 0, if you do not need some quantities
522  * computed. Requesting a value which \e l is not capable of computing
523  * is not an error; the corresponding argument will not be altered.
524  *
525  * With \e flags & GEOD_LONG_NOWRAP bit set, the quantity \e lon2 &minus;
526  * \e lon1 indicates how many times the geodesic wrapped around the
527  * ellipsoid. Because \e lon2 might be outside the normal allowed range
528  * for longitudes, [&minus;540&deg;, 540&deg;), be sure to normalize it,
529  * e.g., with fmod(\e lon2, 360.0) before using it in subsequent
530  * calculations
531  *
532  * Example, compute way points between JFK and Singapore Changi Airport
533  * using geod_genposition(). In this example, the points are evenly space in
534  * arc length (and so only approximately equally space in distance). This is
535  * faster than using geod_position() would be appropriate if drawing the path
536  * on a map.
537  @code
538  struct geod_geodesic g;
539  struct geod_geodesicline l;
540  double a12, azi1, lat[101], lon[101];
541  int i;
542  geod_init(&g, 6378137, 1/298.257223563);
543  a12 = geod_geninverse(&g, 40.64, -73.78, 1.36, 103.99,
544  0, &azi1, 0, 0, 0, 0, 0);
545  geod_lineinit(&l, &g, 40.64, -73.78, azi1, GEOD_LATITUDE | GEOD_LONGITUDE);
546  for (i = 0; i < 101; ++i) {
547  geod_genposition(&l, 1, i * a12 * 0.01,
548  lat + i, lon + i, 0, 0, 0, 0, 0, 0);
549  printf("%.5f %.5f\n", lat[i], lon[i]);
550  }
551  @endcode
552  **********************************************************************/
553  double geod_genposition(const struct geod_geodesicline* l,
554  unsigned flags, double s12_a12,
555  double* plat2, double* plon2, double* pazi2,
556  double* ps12, double* pm12,
557  double* pM12, double* pM21,
558  double* pS12);
559 
560  /**
561  * Initialize a geod_polygon object.
562  *
563  * @param[out] p a pointer to the object to be initialized.
564  * @param[in] polylinep non-zero if a polyline instead of a polygon.
565  *
566  * If \e polylinep is zero, then the sequence of vertices and edges added by
567  * geod_polygon_addpoint() and geod_polygon_addedge() define a polygon and
568  * the perimeter and area are returned by geod_polygon_compute(). If \e
569  * polylinep is non-zero, then the vertices and edges define a polyline and
570  * only the perimeter is returned by geod_polygon_compute().
571  *
572  * The area and perimeter are accumulated at two times the standard floating
573  * point precision to guard against the loss of accuracy with many-sided
574  * polygons. At any point you can ask for the perimeter and area so far.
575  *
576  * An example of the use of this function is given in the documentation for
577  * geod_polygon_compute().
578  **********************************************************************/
579  void geod_polygon_init(struct geod_polygon* p, int polylinep);
580 
581  /**
582  * Add a point to the polygon or polyline.
583  *
584  * @param[in] g a pointer to the geod_geodesic object specifying the
585  * ellipsoid.
586  * @param[in,out] p a pointer to the geod_polygon object specifying the
587  * polygon.
588  * @param[in] lat the latitude of the point (degrees).
589  * @param[in] lon the longitude of the point (degrees).
590  *
591  * \e g and \e p must have been initialized with calls to geod_init() and
592  * geod_polygon_init(), respectively. The same \e g must be used for all the
593  * points and edges in a polygon. \e lat should be in the range
594  * [&minus;90&deg;, 90&deg;] and \e lon should be in the range
595  * [&minus;540&deg;, 540&deg;).
596  *
597  * An example of the use of this function is given in the documentation for
598  * geod_polygon_compute().
599  **********************************************************************/
600  void geod_polygon_addpoint(const struct geod_geodesic* g,
601  struct geod_polygon* p,
602  double lat, double lon);
603 
604  /**
605  * Add an edge to the polygon or polyline.
606  *
607  * @param[in] g a pointer to the geod_geodesic object specifying the
608  * ellipsoid.
609  * @param[in,out] p a pointer to the geod_polygon object specifying the
610  * polygon.
611  * @param[in] azi azimuth at current point (degrees).
612  * @param[in] s distance from current point to next point (meters).
613  *
614  * \e g and \e p must have been initialized with calls to geod_init() and
615  * geod_polygon_init(), respectively. The same \e g must be used for all the
616  * points and edges in a polygon. \e azi should be in the range
617  * [&minus;540&deg;, 540&deg;). This does nothing if no points have been
618  * added yet. The \e lat and \e lon fields of \e p give the location of
619  * the new vertex.
620  **********************************************************************/
621  void geod_polygon_addedge(const struct geod_geodesic* g,
622  struct geod_polygon* p,
623  double azi, double s);
624 
625  /**
626  * Return the results for a polygon.
627  *
628  * @param[in] g a pointer to the geod_geodesic object specifying the
629  * ellipsoid.
630  * @param[in] p a pointer to the geod_polygon object specifying the polygon.
631  * @param[in] reverse if non-zero then clockwise (instead of
632  * counter-clockwise) traversal counts as a positive area.
633  * @param[in] sign if non-zero then return a signed result for the area if
634  * the polygon is traversed in the "wrong" direction instead of returning
635  * the area for the rest of the earth.
636  * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
637  * only set if \e polyline is non-zero in the call to geod_polygon_init().
638  * @param[out] pP pointer to the perimeter of the polygon or length of the
639  * polyline (meters).
640  * @return the number of points.
641  *
642  * The area and perimeter are accumulated at two times the standard floating
643  * point precision to guard against the loss of accuracy with many-sided
644  * polygons. Only simple polygons (which are not self-intersecting) are
645  * allowed. There's no need to "close" the polygon by repeating the first
646  * vertex. Set \e pA or \e pP to zero, if you do not want the corresponding
647  * quantity returned.
648  *
649  * Example, compute the perimeter and area of the geodesic triangle with
650  * vertices (0&deg;N,0&deg;E), (0&deg;N,90&deg;E), (90&deg;N,0&deg;E).
651  @code
652  double A, P;
653  int n;
654  struct geod_geodesic g;
655  struct geod_polygon p;
656  geod_init(&g, 6378137, 1/298.257223563);
657  geod_polygon_init(&p, 0);
658 
659  geod_polygon_addpoint(&g, &p, 0, 0);
660  geod_polygon_addpoint(&g, &p, 0, 90);
661  geod_polygon_addpoint(&g, &p, 90, 0);
662  n = geod_polygon_compute(&g, &p, 0, 1, &A, &P);
663  printf("%d %.8f %.3f\n", n, P, A);
664  @endcode
665  **********************************************************************/
666  unsigned geod_polygon_compute(const struct geod_geodesic* g,
667  const struct geod_polygon* p,
668  int reverse, int sign,
669  double* pA, double* pP);
670 
671  /**
672  * Return the results assuming a tentative final test point is added;
673  * however, the data for the test point is not saved. This lets you report a
674  * running result for the perimeter and area as the user moves the mouse
675  * cursor. Ordinary floating point arithmetic is used to accumulate the data
676  * for the test point; thus the area and perimeter returned are less accurate
677  * than if geod_polygon_addpoint() and geod_polygon_compute() are used.
678  *
679  * @param[in] g a pointer to the geod_geodesic object specifying the
680  * ellipsoid.
681  * @param[in] p a pointer to the geod_polygon object specifying the polygon.
682  * @param[in] lat the latitude of the test point (degrees).
683  * @param[in] lon the longitude of the test point (degrees).
684  * @param[in] reverse if non-zero then clockwise (instead of
685  * counter-clockwise) traversal counts as a positive area.
686  * @param[in] sign if non-zero then return a signed result for the area if
687  * the polygon is traversed in the "wrong" direction instead of returning
688  * the area for the rest of the earth.
689  * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
690  * only set if \e polyline is non-zero in the call to geod_polygon_init().
691  * @param[out] pP pointer to the perimeter of the polygon or length of the
692  * polyline (meters).
693  * @return the number of points.
694  *
695  * \e lat should be in the range [&minus;90&deg;, 90&deg;] and \e
696  * lon should be in the range [&minus;540&deg;, 540&deg;).
697  **********************************************************************/
698  unsigned geod_polygon_testpoint(const struct geod_geodesic* g,
699  const struct geod_polygon* p,
700  double lat, double lon,
701  int reverse, int sign,
702  double* pA, double* pP);
703 
704  /**
705  * Return the results assuming a tentative final test point is added via an
706  * azimuth and distance; however, the data for the test point is not saved.
707  * This lets you report a running result for the perimeter and area as the
708  * user moves the mouse cursor. Ordinary floating point arithmetic is used
709  * to accumulate the data for the test point; thus the area and perimeter
710  * returned are less accurate than if geod_polygon_addedge() and
711  * geod_polygon_compute() are used.
712  *
713  * @param[in] g a pointer to the geod_geodesic object specifying the
714  * ellipsoid.
715  * @param[in] p a pointer to the geod_polygon object specifying the polygon.
716  * @param[in] azi azimuth at current point (degrees).
717  * @param[in] s distance from current point to final test point (meters).
718  * @param[in] reverse if non-zero then clockwise (instead of
719  * counter-clockwise) traversal counts as a positive area.
720  * @param[in] sign if non-zero then return a signed result for the area if
721  * the polygon is traversed in the "wrong" direction instead of returning
722  * the area for the rest of the earth.
723  * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
724  * only set if \e polyline is non-zero in the call to geod_polygon_init().
725  * @param[out] pP pointer to the perimeter of the polygon or length of the
726  * polyline (meters).
727  * @return the number of points.
728  *
729  * \e azi should be in the range [&minus;540&deg;, 540&deg;).
730  **********************************************************************/
731  unsigned geod_polygon_testedge(const struct geod_geodesic* g,
732  const struct geod_polygon* p,
733  double azi, double s,
734  int reverse, int sign,
735  double* pA, double* pP);
736 
737  /**
738  * A simple interface for computing the area of a geodesic polygon.
739  *
740  * @param[in] g a pointer to the geod_geodesic object specifying the
741  * ellipsoid.
742  * @param[in] lats an array of latitudes of the polygon vertices (degrees).
743  * @param[in] lons an array of longitudes of the polygon vertices (degrees).
744  * @param[in] n the number of vertices.
745  * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>).
746  * @param[out] pP pointer to the perimeter of the polygon (meters).
747  *
748  * \e lats should be in the range [&minus;90&deg;, 90&deg;]; \e lons should
749  * be in the range [&minus;540&deg;, 540&deg;).
750  *
751  * Only simple polygons (which are not self-intersecting) are allowed.
752  * There's no need to "close" the polygon by repeating the first vertex. The
753  * area returned is signed with counter-clockwise traversal being treated as
754  * positive.
755  *
756  * Example, compute the area of Antarctica:
757  @code
758  double
759  lats[] = {-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
760  -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7},
761  lons[] = {-74, -102, -102, -131, -163, 163, 172, 140, 113,
762  88, 59, 25, -4, -14, -33, -46, -61};
763  struct geod_geodesic g;
764  double A, P;
765  geod_init(&g, 6378137, 1/298.257223563);
766  geod_polygonarea(&g, lats, lons, (sizeof lats) / (sizeof lats[0]), &A, &P);
767  printf("%.0f %.2f\n", A, P);
768  @endcode
769  **********************************************************************/
770  void geod_polygonarea(const struct geod_geodesic* g,
771  double lats[], double lons[], int n,
772  double* pA, double* pP);
773 
774  /**
775  * mask values for the \e caps argument to geod_lineinit().
776  **********************************************************************/
777  enum geod_mask {
778  GEOD_NONE = 0U, /**< Calculate nothing */
779  GEOD_LATITUDE = 1U<<7 | 0U, /**< Calculate latitude */
780  GEOD_LONGITUDE = 1U<<8 | 1U<<3, /**< Calculate longitude */
781  GEOD_AZIMUTH = 1U<<9 | 0U, /**< Calculate azimuth */
782  GEOD_DISTANCE = 1U<<10 | 1U<<0, /**< Calculate distance */
783  GEOD_DISTANCE_IN = 1U<<11 | 1U<<0 | 1U<<1, /**< Allow distance as input */
784  GEOD_REDUCEDLENGTH= 1U<<12 | 1U<<0 | 1U<<2, /**< Calculate reduced length */
785  GEOD_GEODESICSCALE= 1U<<13 | 1U<<0 | 1U<<2, /**< Calculate geodesic scale */
786  GEOD_AREA = 1U<<14 | 1U<<4, /**< Calculate reduced length */
787  GEOD_ALL = 0x7F80U| 0x1FU /**< Calculate everything */
788  };
789 
790  /**
791  * flag values for the \e flags argument to geod_gendirect() and
792  * geod_genposition()
793  **********************************************************************/
794  enum geod_flags {
795  GEOD_NOFLAGS = 0U, /**< No flags */
796  GEOD_ARCMODE = 1U<<0, /**< Position given in terms of arc distance */
797  GEOD_LONG_NOWRAP = 1U<<15 /**< Don't wrap longitude */
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:794
geod_mask
Definition: geodesic.h:777
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