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