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