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