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