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