C library for Geodesics 2.2
Loading...
Searching...
No Matches
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) <karney@alum.mit.edu> 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 2
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)
79extern "C" {
80#endif
81
82 /**
83 * The struct containing information about the ellipsoid. This must be
84 * initialized by geod_init() before use.
85 **********************************************************************/
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 **********************************************************************/
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 **********************************************************************/
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 **********************************************************************/
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 **********************************************************************/
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 **********************************************************************/
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 lon
Definition geodesic.h:126
unsigned num
Definition geodesic.h:135
double lat
Definition geodesic.h:125