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