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