00001 /** 00002 * \file geodesic.h 00003 * \brief Header for the geodesic routines in C 00004 * 00005 * This an implementation in C of the geodesic algorithms described in 00006 * - C. F. F. Karney, 00007 * <a href="http://dx.doi.org/10.1007/s00190-012-0578-z"> 00008 * Algorithms for geodesics</a>, 00009 * J. Geodesy <b>87</b>, 43--55 (2013); 00010 * DOI: <a href="http://dx.doi.org/10.1007/s00190-012-0578-z"> 00011 * 10.1007/s00190-012-0578-z</a>; 00012 * addenda: <a href="http://geographiclib.sf.net/geod-addenda.html"> 00013 * geod-addenda.html</a>. 00014 * . 00015 * The principal advantages of these algorithms over previous ones (e.g., 00016 * Vincenty, 1975) are 00017 * - accurate to round off for |<i>f</i>| < 1/50; 00018 * - the solution of the inverse problem is always found; 00019 * - differential and integral properties of geodesics are computed. 00020 * 00021 * The shortest path between two points on the ellipsoid at (\e lat1, \e 00022 * lon1) and (\e lat2, \e lon2) is called the geodesic. Its length is 00023 * \e s12 and the geodesic from point 1 to point 2 has forward azimuths 00024 * \e azi1 and \e azi2 at the two end points. 00025 * 00026 * Traditionally two geodesic problems are considered: 00027 * - the direct problem -- given \e lat1, \e lon1, \e s12, and \e azi1, 00028 * determine \e lat2, \e lon2, and \e azi2. This is solved by the function 00029 * geod_direct(). 00030 * - the inverse problem -- given \e lat1, \e lon1, and \e lat2, \e lon2, 00031 * determine \e s12, \e azi1, and \e azi2. This is solved by the function 00032 * geod_inverse(). 00033 * 00034 * The ellipsoid is specified by its equatorial radius \e a (typically in 00035 * meters) and flattening \e f. The routines are accurate to round off with 00036 * double precision arithmetic provided that |<i>f</i>| < 1/50; for the 00037 * WGS84 ellipsoid, the errors are less than 15 nanometers. (Reasonably 00038 * accurate results are obtained for |<i>f</i>| < 1/5.) For a prolate 00039 * ellipsoid, specify \e f < 0. 00040 * 00041 * The routines also calculate several other quantities of interest 00042 * - \e S12 is the area between the geodesic from point 1 to point 2 and the 00043 * equator; i.e., it is the area, measured counter-clockwise, of the 00044 * quadrilateral with corners (\e lat1,\e lon1), (0,\e lon1), (0,\e lon2), 00045 * and (\e lat2,\e lon2). 00046 * - \e m12, the reduced length of the geodesic is defined such that if 00047 * the initial azimuth is perturbed by \e dazi1 (radians) then the 00048 * second point is displaced by \e m12 \e dazi1 in the direction 00049 * perpendicular to the geodesic. On a curved surface the reduced 00050 * length obeys a symmetry relation, \e m12 + \e m21 = 0. On a flat 00051 * surface, we have \e m12 = \e s12. 00052 * - \e M12 and \e M21 are geodesic scales. If two geodesics are 00053 * parallel at point 1 and separated by a small distance \e dt, then 00054 * they are separated by a distance \e M12 \e dt at point 2. \e M21 00055 * is defined similarly (with the geodesics being parallel to one 00056 * another at point 2). On a flat surface, we have \e M12 = \e M21 00057 * = 1. 00058 * - \e a12 is the arc length on the auxiliary sphere. This is a 00059 * construct for converting the problem to one in spherical 00060 * trigonometry. \e a12 is measured in degrees. The spherical arc 00061 * length from one equator crossing to the next is always 180°. 00062 * 00063 * If points 1, 2, and 3 lie on a single geodesic, then the following 00064 * addition rules hold: 00065 * - \e s13 = \e s12 + \e s23 00066 * - \e a13 = \e a12 + \e a23 00067 * - \e S13 = \e S12 + \e S23 00068 * - \e m13 = \e m12 \e M23 + \e m23 \e M21 00069 * - \e M13 = \e M12 \e M23 − (1 − \e M12 \e M21) \e 00070 * m23 / \e m12 00071 * - \e M31 = \e M32 \e M21 − (1 − \e M23 \e M32) \e 00072 * m12 / \e m23 00073 * 00074 * The shortest distance returned by the solution of the inverse problem is 00075 * (obviously) uniquely defined. However, in a few special cases there are 00076 * multiple azimuths which yield the same shortest distance. Here is a 00077 * catalog of those cases: 00078 * - \e lat1 = −\e lat2 (with neither point at a pole). If \e azi1 = \e 00079 * azi2, the geodesic is unique. Otherwise there are two geodesics 00080 * and the second one is obtained by setting [\e azi1, \e azi2] = [\e 00081 * azi2, \e azi1], [\e M12, \e M21] = [\e M21, \e M12], \e S12 = 00082 * −\e S12. (This occurs when the longitude difference is near 00083 * ±180° for oblate ellipsoids.) 00084 * - \e lon2 = \e lon1 ± 180° (with neither point at a pole). If 00085 * \e azi1 = 0° or ±180°, the geodesic is unique. 00086 * Otherwise there are two geodesics and the second one is obtained by 00087 * setting [\e azi1, \e azi2] = [−\e azi1, −\e azi2], \e 00088 * S12 = −\e S12. (This occurs when \e lat2 is near 00089 * −\e lat1 for prolate ellipsoids.) 00090 * - Points 1 and 2 at opposite poles. There are infinitely many 00091 * geodesics which can be generated by setting [\e azi1, \e azi2] = 00092 * [\e azi1, \e azi2] + [\e d, −\e d], for arbitrary \e d. (For 00093 * spheres, this prescription applies when points 1 and 2 are 00094 * antipodal.) 00095 * - \e s12 = 0 (coincident points). There are infinitely many geodesics 00096 * which can be generated by setting [\e azi1, \e azi2] = [\e azi1, \e 00097 * azi2] + [\e d, \e d], for arbitrary \e d. 00098 * 00099 * These routines are a simple transcription of the corresponding C++ classes 00100 * in <a href="http://geographiclib.sf.net"> GeographicLib</a>. The "class 00101 * data" is represented by the structs geod_geodesic, geod_geodesicline, 00102 * geod_polygon and pointers to these objects are passed as initial arguments 00103 * to the member functions. Most of the internal comments have been retained. 00104 * However, in the process of transcription some documentation has been lost 00105 * and the documentation for the C++ classes, GeographicLib::Geodesic, 00106 * GeographicLib::GeodesicLine, and GeographicLib::PolygonArea, should be 00107 * consulted. The C++ code remains the "reference implementation". Think 00108 * twice about restructuring the internals of the C code since this may make 00109 * porting fixes from the C++ code more difficult. 00110 * 00111 * Copyright (c) Charles Karney (2012-2013) <charles@karney.com> and licensed 00112 * under the MIT/X11 License. For more information, see 00113 * http://geographiclib.sourceforge.net/ 00114 * 00115 * This library was distributed with 00116 * <a href="../index.html">GeographicLib</a> 1.32. 00117 **********************************************************************/ 00118 00119 #if !defined(GEODESIC_H) 00120 #define GEODESIC_H 1 00121 00122 /** 00123 * The major version of the geodesic library. (This tracks the version of 00124 * GeographicLib.) 00125 **********************************************************************/ 00126 #define GEODESIC_VERSION_MAJOR 1 00127 /** 00128 * The minor version of the geodesic library. (This tracks the version of 00129 * GeographicLib.) 00130 **********************************************************************/ 00131 #define GEODESIC_VERSION_MINOR 32 00132 /** 00133 * The patch level of the geodesic library. (This tracks the version of 00134 * GeographicLib.) 00135 **********************************************************************/ 00136 #define GEODESIC_VERSION_PATCH 0 00137 00138 #if defined(__cplusplus) 00139 extern "C" { 00140 #endif 00141 00142 /** 00143 * The struct containing information about the ellipsoid. This must be 00144 * initialized by geod_init() before use. 00145 **********************************************************************/ 00146 struct geod_geodesic { 00147 double a; /**< the equatorial radius */ 00148 double f; /**< the flattening */ 00149 /**< @cond SKIP */ 00150 double f1, e2, ep2, n, b, c2, etol2; 00151 double A3x[6], C3x[15], C4x[21]; 00152 /**< @endcond */ 00153 }; 00154 00155 /** 00156 * The struct containing information about a single geodesic. This must be 00157 * initialized by geod_lineinit() before use. 00158 **********************************************************************/ 00159 struct geod_geodesicline { 00160 double lat1; /**< the starting latitude */ 00161 double lon1; /**< the starting longitude */ 00162 double azi1; /**< the starting azimuth */ 00163 double a; /**< the equatorial radius */ 00164 double f; /**< the flattening */ 00165 /**< @cond SKIP */ 00166 double b, c2, f1, salp0, calp0, k2, 00167 salp1, calp1, ssig1, csig1, dn1, stau1, ctau1, somg1, comg1, 00168 A1m1, A2m1, A3c, B11, B21, B31, A4, B41; 00169 double C1a[6+1], C1pa[6+1], C2a[6+1], C3a[6], C4a[6]; 00170 /**< @endcond */ 00171 unsigned caps; /**< the capabilities */ 00172 }; 00173 00174 /** 00175 * The struct for accumulating information about a geodesic polygon. This is 00176 * used for computing the perimeter and area of a polygon. This must be 00177 * initialized by geod_polygon_init() before use. 00178 **********************************************************************/ 00179 struct geod_polygon { 00180 double lat; /**< the current latitude */ 00181 double lon; /**< the current longitude */ 00182 /**< @cond SKIP */ 00183 double lat0; 00184 double lon0; 00185 double A[2]; 00186 double P[2]; 00187 int polyline; 00188 int crossings; 00189 /**< @endcond */ 00190 unsigned num; /**< the number of points so far */ 00191 }; 00192 00193 /** 00194 * Initialize a geod_geodesic object. 00195 * 00196 * @param[out] g a pointer to the object to be initialized. 00197 * @param[in] a the equatorial radius (meters). 00198 * @param[in] f the flattening. 00199 **********************************************************************/ 00200 void geod_init(struct geod_geodesic* g, double a, double f); 00201 00202 /** 00203 * Initialize a geod_geodesicline object. 00204 * 00205 * @param[out] l a pointer to the object to be initialized. 00206 * @param[in] g a pointer to the geod_geodesic object specifying the 00207 * ellipsoid. 00208 * @param[in] lat1 latitude of point 1 (degrees). 00209 * @param[in] lon1 longitude of point 1 (degrees). 00210 * @param[in] azi1 azimuth at point 1 (degrees). 00211 * @param[in] caps bitor'ed combination of geod_mask() values specifying the 00212 * capabilities the geod_geodesicline object should possess, i.e., which 00213 * quantities can be returned in calls to geod_position() and 00214 * geod_genposition(). 00215 * 00216 * \e g must have been initialized with a call to geod_init(). \e lat1 00217 * should be in the range [−90°, 90°]; \e lon1 and \e azi1 00218 * should be in the range [−540°, 540°). 00219 * 00220 * The geod_mask values are [see geod_mask()]: 00221 * - \e caps |= GEOD_LATITUDE for the latitude \e lat2; this is 00222 * added automatically, 00223 * - \e caps |= GEOD_LONGITUDE for the latitude \e lon2, 00224 * - \e caps |= GEOD_AZIMUTH for the latitude \e azi2; this is 00225 * added automatically, 00226 * - \e caps |= GEOD_DISTANCE for the distance \e s12, 00227 * - \e caps |= GEOD_REDUCEDLENGTH for the reduced length \e m12, 00228 * - \e caps |= GEOD_GEODESICSCALE for the geodesic scales \e M12 00229 * and \e M21, 00230 * - \e caps |= GEOD_AREA for the area \e S12, 00231 * - \e caps |= GEOD_DISTANCE_IN permits the length of the 00232 * geodesic to be given in terms of \e s12; without this capability the 00233 * length can only be specified in terms of arc length. 00234 * . 00235 * A value of \e caps = 0 is treated as GEOD_LATITUDE | GEOD_LONGITUDE | 00236 * GEOD_AZIMUTH | GEOD_DISTANCE_IN (to support the solution of the "standard" 00237 * direct problem). 00238 **********************************************************************/ 00239 void geod_lineinit(struct geod_geodesicline* l, 00240 const struct geod_geodesic* g, 00241 double lat1, double lon1, double azi1, unsigned caps); 00242 00243 /** 00244 * Solve the direct geodesic problem. 00245 * 00246 * @param[in] g a pointer to the geod_geodesic object specifying the 00247 * ellipsoid. 00248 * @param[in] lat1 latitude of point 1 (degrees). 00249 * @param[in] lon1 longitude of point 1 (degrees). 00250 * @param[in] azi1 azimuth at point 1 (degrees). 00251 * @param[in] s12 distance between point 1 and point 2 (meters); it can be 00252 * negative. 00253 * @param[out] plat2 pointer to the latitude of point 2 (degrees). 00254 * @param[out] plon2 pointer to the longitude of point 2 (degrees). 00255 * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). 00256 * 00257 * \e g must have been initialized with a call to geod_init(). \e lat1 00258 * should be in the range [−90°, 90°]; \e lon1 and \e azi1 00259 * should be in the range [−540°, 540°). The values of \e lon2 00260 * and \e azi2 returned are in the range [−180°, 180°). Any of 00261 * the "return" arguments \e plat2, etc., may be replaced by 0, if you do not 00262 * need some quantities computed. 00263 * 00264 * If either point is at a pole, the azimuth is defined by keeping the 00265 * longitude fixed, writing \e lat = ±(90° − ε), and 00266 * taking the limit ε → 0+. An arc length greater that 180° 00267 * signifies a geodesic which is not a shortest path. (For a prolate 00268 * ellipsoid, an additional condition is necessary for a shortest path: the 00269 * longitudinal extent must not exceed of 180°.) 00270 * 00271 * Example, determine the point 10000 km NE of JFK: 00272 @code 00273 struct geod_geodesic g; 00274 double lat, lon; 00275 geod_init(&g, 6378137, 1/298.257223563); 00276 geod_direct(&g, 40.64, -73.78, 45.0, 10e6, &lat, &lon, 0); 00277 printf("%.5f %.5f\n", lat, lon); 00278 @endcode 00279 **********************************************************************/ 00280 void geod_direct(const struct geod_geodesic* g, 00281 double lat1, double lon1, double azi1, double s12, 00282 double* plat2, double* plon2, double* pazi2); 00283 00284 /** 00285 * Solve the inverse geodesic problem. 00286 * 00287 * @param[in] g a pointer to the geod_geodesic object specifying the 00288 * ellipsoid. 00289 * @param[in] lat1 latitude of point 1 (degrees). 00290 * @param[in] lon1 longitude of point 1 (degrees). 00291 * @param[in] lat2 latitude of point 2 (degrees). 00292 * @param[in] lon2 longitude of point 2 (degrees). 00293 * @param[out] ps12 pointer to the distance between point 1 and point 2 00294 * (meters). 00295 * @param[out] pazi1 pointer to the azimuth at point 1 (degrees). 00296 * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). 00297 * 00298 * \e g must have been initialized with a call to geod_init(). \e lat1 00299 * and \e lat2 should be in the range [−90°, 90°]; \e lon1 and 00300 * \e lon2 should be in the range [−540°, 540°). The values of 00301 * \e azi1 and \e azi2 returned are in the range [−180°, 180°). 00302 * Any of the "return" arguments \e ps12, etc., may be replaced by 0, if you 00303 * do not need some quantities computed. 00304 * 00305 * If either point is at a pole, the azimuth is defined by keeping the 00306 * longitude fixed, writing \e lat = ±(90° − ε), and 00307 * taking the limit ε → 0+. 00308 * 00309 * The solution to the inverse problem is found using Newton's method. If 00310 * this fails to converge (this is very unlikely in geodetic applications 00311 * but does occur for very eccentric ellipsoids), then the bisection method 00312 * is used to refine the solution. 00313 * 00314 * Example, determine the distance between JFK and Singapore Changi Airport: 00315 @code 00316 struct geod_geodesic g; 00317 double s12; 00318 geod_init(&g, 6378137, 1/298.257223563); 00319 geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, 0, 0); 00320 printf("%.3f\n", s12); 00321 @endcode 00322 **********************************************************************/ 00323 void geod_inverse(const struct geod_geodesic* g, 00324 double lat1, double lon1, double lat2, double lon2, 00325 double* ps12, double* pazi1, double* pazi2); 00326 00327 /** 00328 * Compute the position along a geod_geodesicline. 00329 * 00330 * @param[in] l a pointer to the geod_geodesicline object specifying the 00331 * geodesic line. 00332 * @param[in] s12 distance between point 1 and point 2 (meters); it can be 00333 * negative. 00334 * @param[out] plat2 pointer to the latitude of point 2 (degrees). 00335 * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires 00336 * that \e l was initialized with \e caps |= GEOD_LONGITUDE. 00337 * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). 00338 * 00339 * \e l must have been initialized with a call to geod_lineinit() with \e 00340 * caps |= GEOD_DISTANCE_IN. The values of \e lon2 and \e azi2 returned are 00341 * in the range [−180°, 180°). Any of the "return" arguments 00342 * \e plat2, etc., may be replaced by 0, if you do not need some quantities 00343 * computed. 00344 * 00345 * Example, compute way points between JFK and Singapore Changi Airport 00346 * the "obvious" way using geod_direct(): 00347 @code 00348 struct geod_geodesic g; 00349 double s12, azi1, lat[101],lon[101]; 00350 int i; 00351 geod_init(&g, 6378137, 1/298.257223563); 00352 geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0); 00353 for (i = 0; i < 101; ++i) { 00354 geod_direct(&g, 40.64, -73.78, azi1, i * s12 * 0.01, lat + i, lon + i, 0); 00355 printf("%.5f %.5f\n", lat[i], lon[i]); 00356 } 00357 @endcode 00358 * A faster way using geod_position(): 00359 @code 00360 struct geod_geodesic g; 00361 struct geod_geodesicline l; 00362 double s12, azi1, lat[101],lon[101]; 00363 int i; 00364 geod_init(&g, 6378137, 1/298.257223563); 00365 geod_inverse(&g, 40.64, -73.78, 1.36, 103.99, &s12, &azi1, 0); 00366 geod_lineinit(&l, &g, 40.64, -73.78, azi1, 0); 00367 for (i = 0; i < 101; ++i) { 00368 geod_position(&l, i * s12 * 0.01, lat + i, lon + i, 0); 00369 printf("%.5f %.5f\n", lat[i], lon[i]); 00370 } 00371 @endcode 00372 **********************************************************************/ 00373 void geod_position(const struct geod_geodesicline* l, double s12, 00374 double* plat2, double* plon2, double* pazi2); 00375 00376 /** 00377 * The general direct geodesic problem. 00378 * 00379 * @param[in] g a pointer to the geod_geodesic object specifying the 00380 * ellipsoid. 00381 * @param[in] lat1 latitude of point 1 (degrees). 00382 * @param[in] lon1 longitude of point 1 (degrees). 00383 * @param[in] azi1 azimuth at point 1 (degrees). 00384 * @param[in] arcmode flag determining the meaning of the \e 00385 * s12_a12. 00386 * @param[in] s12_a12 if \e arcmode is 0, this is the distance between 00387 * point 1 and point 2 (meters); otherwise it is the arc length between 00388 * point 1 and point 2 (degrees); it can be negative. 00389 * @param[out] plat2 pointer to the latitude of point 2 (degrees). 00390 * @param[out] plon2 pointer to the longitude of point 2 (degrees). 00391 * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). 00392 * @param[out] ps12 pointer to the distance between point 1 and point 2 00393 * (meters). 00394 * @param[out] pm12 pointer to the reduced length of geodesic (meters). 00395 * @param[out] pM12 pointer to the geodesic scale of point 2 relative to 00396 * point 1 (dimensionless). 00397 * @param[out] pM21 pointer to the geodesic scale of point 1 relative to 00398 * point 2 (dimensionless). 00399 * @param[out] pS12 pointer to the area under the geodesic 00400 * (meters<sup>2</sup>). 00401 * @return \e a12 arc length of between point 1 and point 2 (degrees). 00402 * 00403 * \e g must have been initialized with a call to geod_init(). \e lat1 00404 * should be in the range [−90°, 90°]; \e lon1 and \e azi1 00405 * should be in the range [−540°, 540°). The function value \e 00406 * a12 equals \e s12_a12 is \e arcmode is non-zero. Any of the "return" 00407 * arguments \e plat2, etc., may be replaced by 0, if you do not need some 00408 * quantities computed. 00409 **********************************************************************/ 00410 double geod_gendirect(const struct geod_geodesic* g, 00411 double lat1, double lon1, double azi1, 00412 int arcmode, double s12_a12, 00413 double* plat2, double* plon2, double* pazi2, 00414 double* ps12, double* pm12, double* pM12, double* pM21, 00415 double* pS12); 00416 00417 /** 00418 * The general inverse geodesic calculation. 00419 * 00420 * @param[in] g a pointer to the geod_geodesic object specifying the 00421 * ellipsoid. 00422 * @param[in] lat1 latitude of point 1 (degrees). 00423 * @param[in] lon1 longitude of point 1 (degrees). 00424 * @param[in] lat2 latitude of point 2 (degrees). 00425 * @param[in] lon2 longitude of point 2 (degrees). 00426 * @param[out] ps12 pointer to the distance between point 1 and point 2 00427 * (meters). 00428 * @param[out] pazi1 pointer to the azimuth at point 1 (degrees). 00429 * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). 00430 * @param[out] pm12 pointer to the reduced length of geodesic (meters). 00431 * @param[out] pM12 pointer to the geodesic scale of point 2 relative to 00432 * point 1 (dimensionless). 00433 * @param[out] pM21 pointer to the geodesic scale of point 1 relative to 00434 * point 2 (dimensionless). 00435 * @param[out] pS12 pointer to the area under the geodesic 00436 * (meters<sup>2</sup>). 00437 * @return \e a12 arc length of between point 1 and point 2 (degrees). 00438 * 00439 * \e g must have been initialized with a call to geod_init(). \e lat1 00440 * and \e lat2 should be in the range [−90°, 90°]; \e lon1 and 00441 * \e lon2 should be in the range [−540°, 540°). Any of the 00442 * "return" arguments \e ps12, etc., may be replaced by 0, if you do not need 00443 * some quantities computed. 00444 **********************************************************************/ 00445 double geod_geninverse(const struct geod_geodesic* g, 00446 double lat1, double lon1, double lat2, double lon2, 00447 double* ps12, double* pazi1, double* pazi2, 00448 double* pm12, double* pM12, double* pM21, 00449 double* pS12); 00450 00451 /** 00452 * The general position function. 00453 * 00454 * @param[in] l a pointer to the geod_geodesicline object specifying the 00455 * geodesic line. 00456 * @param[in] arcmode flag determining the meaning of the second parameter; 00457 * if arcmode is 0, then \e l must have been initialized with \e caps |= 00458 * GEOD_DISTANCE_IN. 00459 * @param[in] s12_a12 if \e arcmode is 0, this is the distance between 00460 * point 1 and point 2 (meters); otherwise it is the arc length between 00461 * point 1 and point 2 (degrees); it can be negative. 00462 * @param[out] plat2 pointer to the latitude of point 2 (degrees). 00463 * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires 00464 * that \e l was initialized with \e caps |= GEOD_LONGITUDE. 00465 * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). 00466 * @param[out] ps12 pointer to the distance between point 1 and point 2 00467 * (meters); requires that \e l was initialized with \e caps |= 00468 * GEOD_DISTANCE. 00469 * @param[out] pm12 pointer to the reduced length of geodesic (meters); 00470 * requires that \e l was initialized with \e caps |= GEOD_REDUCEDLENGTH. 00471 * @param[out] pM12 pointer to the geodesic scale of point 2 relative to 00472 * point 1 (dimensionless); requires that \e l was initialized with \e caps 00473 * |= GEOD_GEODESICSCALE. 00474 * @param[out] pM21 pointer to the geodesic scale of point 1 relative to 00475 * point 2 (dimensionless); requires that \e l was initialized with \e caps 00476 * |= GEOD_GEODESICSCALE. 00477 * @param[out] pS12 pointer to the area under the geodesic 00478 * (meters<sup>2</sup>); requires that \e l was initialized with \e caps |= 00479 * GEOD_AREA. 00480 * @return \e a12 arc length of between point 1 and point 2 (degrees). 00481 * 00482 * \e l must have been initialized with a call to geod_lineinit() with \e 00483 * caps |= GEOD_DISTANCE_IN. The values of \e lon2 and \e azi2 returned are 00484 * in the range [−180°, 180°). Any of the "return" arguments 00485 * \e plat2, etc., may be replaced by 0, if you do not need some quantities 00486 * computed. Requesting a value which \e l is not capable of computing is 00487 * not an error; the corresponding argument will not be altered. 00488 * 00489 * Example, compute way points between JFK and Singapore Changi Airport 00490 * using geod_genposition(). In this example, the points are evenly space in 00491 * arc length (and so only approximately equally space in distance). This is 00492 * faster than using geod_position() would be appropriate if drawing the path 00493 * on a map. 00494 @code 00495 struct geod_geodesic g; 00496 struct geod_geodesicline l; 00497 double a12, azi1, lat[101],lon[101]; 00498 int i; 00499 geod_init(&g, 6378137, 1/298.257223563); 00500 a12 = geod_geninverse(&g, 40.64, -73.78, 1.36, 103.99, 00501 0, &azi1, 0, 0, 0, 0, 0); 00502 geod_lineinit(&l, &g, 40.64, -73.78, azi1, GEOD_LATITUDE | GEOD_LONGITUDE); 00503 for (i = 0; i < 101; ++i) { 00504 geod_genposition(&l, 1, i * a12 * 0.01, 00505 lat + i, lon + i, 0, 0, 0, 0, 0, 0); 00506 printf("%.5f %.5f\n", lat[i], lon[i]); 00507 } 00508 @endcode 00509 **********************************************************************/ 00510 double geod_genposition(const struct geod_geodesicline* l, 00511 int arcmode, double s12_a12, 00512 double* plat2, double* plon2, double* pazi2, 00513 double* ps12, double* pm12, 00514 double* pM12, double* pM21, 00515 double* pS12); 00516 00517 /** 00518 * Initialize a geod_polygon object. 00519 * 00520 * @param[out] p a pointer to the object to be initialized. 00521 * @param[in] polylinep non-zero if a polyline instead of a polygon. 00522 * 00523 * If \e polylinep is zero, then the sequence of vertices and edges added by 00524 * geod_polygon_addpoint() and geod_polygon_addedge() define a polygon and 00525 * the perimeter and area are returned by geod_polygon_compute(). If \e 00526 * polylinep is non-zero, then the vertices and edges define a polyline and 00527 * only the perimeter is returned by geod_polygon_compute(). 00528 * 00529 * An example of the use of this function is given in the documentation for 00530 * geod_polygon_compute(). 00531 **********************************************************************/ 00532 void geod_polygon_init(struct geod_polygon* p, int polylinep); 00533 00534 /** 00535 * Add a point to the polygon or polyline. 00536 * 00537 * @param[in] g a pointer to the geod_geodesic object specifying the 00538 * ellipsoid. 00539 * @param[in,out] p a pointer to the geod_polygon object specifying the 00540 * polygon. 00541 * @param[in] lat the latitude of the point (degrees). 00542 * @param[in] lon the longitude of the point (degrees). 00543 * 00544 * \e g and \e p must have been initialized with calls to geod_init() and 00545 * geod_polygon_init(), respectively. The same \e g must be used for all the 00546 * points and edges in a polygon. \e lat should be in the range 00547 * [−90°, 90°] and \e lon should be in the range 00548 * [−540°, 540°). 00549 * 00550 * An example of the use of this function is given in the documentation for 00551 * geod_polygon_compute(). 00552 **********************************************************************/ 00553 void geod_polygon_addpoint(const struct geod_geodesic* g, 00554 struct geod_polygon* p, 00555 double lat, double lon); 00556 00557 /** 00558 * Add an edge to the polygon or polyline. 00559 * 00560 * @param[in] g a pointer to the geod_geodesic object specifying the 00561 * ellipsoid. 00562 * @param[in,out] p a pointer to the geod_polygon object specifying the 00563 * polygon. 00564 * @param[in] azi azimuth at current point (degrees). 00565 * @param[in] s distance from current point to next point (meters). 00566 * 00567 * \e g and \e p must have been initialized with calls to geod_init() and 00568 * geod_polygon_init(), respectively. The same \e g must be used for all the 00569 * points and edges in a polygon. \e azi should be in the range 00570 * [−540°, 540°). This does nothing if no points have been 00571 * added yet. The \e lat and \e lon fields of \e p give the location of 00572 * the new vertex. 00573 **********************************************************************/ 00574 void geod_polygon_addedge(const struct geod_geodesic* g, 00575 struct geod_polygon* p, 00576 double azi, double s); 00577 00578 /** 00579 * Return the results for a polygon. 00580 * 00581 * @param[in] g a pointer to the geod_geodesic object specifying the 00582 * ellipsoid. 00583 * @param[in] p a pointer to the geod_polygon object specifying the polygon. 00584 * @param[in] reverse if non-zero then clockwise (instead of 00585 * counter-clockwise) traversal counts as a positive area. 00586 * @param[in] sign if non-zero then return a signed result for the area if 00587 * the polygon is traversed in the "wrong" direction instead of returning 00588 * the area for the rest of the earth. 00589 * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>); 00590 * only set if \e polyline is non-zero in the call to geod_polygon_init(). 00591 * @param[out] pP pointer to the perimeter of the polygon or length of the 00592 * polyline (meters). 00593 * @return the number of points. 00594 * 00595 * Only simple polygons (which are not self-intersecting) are allowed. 00596 * There's no need to "close" the polygon by repeating the first vertex. Set 00597 * \e pA or \e pP to zero, if you do not want the corresponding quantity 00598 * returned. 00599 * 00600 * Example, compute the perimeter and area of the geodesic triangle with 00601 * vertices (0°N,0°E), (0°N,90°E), (90°N,0°E). 00602 @code 00603 double A, P; 00604 int n; 00605 struct geod_geodesic g; 00606 struct geod_polygon p; 00607 geod_init(&g, 6378137, 1/298.257223563); 00608 geod_polygon_init(&p, 0); 00609 00610 geod_polygon_addpoint(&g, &p, 0, 0); 00611 geod_polygon_addpoint(&g, &p, 0, 90); 00612 geod_polygon_addpoint(&g, &p, 90, 0); 00613 n = geod_polygon_compute(&g, &p, 0, 1, &A, &P); 00614 printf("%d %.8f %.3f\n", n, P, A); 00615 @endcode 00616 **********************************************************************/ 00617 unsigned geod_polygon_compute(const struct geod_geodesic* g, 00618 const struct geod_polygon* p, 00619 int reverse, int sign, 00620 double* pA, double* pP); 00621 00622 /** 00623 * Return the results assuming a tentative final test point is added; 00624 * however, the data for the test point is not saved. This lets you report a 00625 * running result for the perimeter and area as the user moves the mouse 00626 * cursor. Ordinary floating point arithmetic is used to accumulate the data 00627 * for the test point; thus the area and perimeter returned are less accurate 00628 * than if geod_polygon_addpoint() and geod_polygon_compute() are used. 00629 * 00630 * @param[in] g a pointer to the geod_geodesic object specifying the 00631 * ellipsoid. 00632 * @param[in] p a pointer to the geod_polygon object specifying the polygon. 00633 * @param[in] lat the latitude of the test point (degrees). 00634 * @param[in] lon the longitude of the test point (degrees). 00635 * @param[in] reverse if non-zero then clockwise (instead of 00636 * counter-clockwise) traversal counts as a positive area. 00637 * @param[in] sign if non-zero then return a signed result for the area if 00638 * the polygon is traversed in the "wrong" direction instead of returning 00639 * the area for the rest of the earth. 00640 * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>); 00641 * only set if \e polyline is non-zero in the call to geod_polygon_init(). 00642 * @param[out] pP pointer to the perimeter of the polygon or length of the 00643 * polyline (meters). 00644 * @return the number of points. 00645 * 00646 * \e lat should be in the range [−90°, 90°] and \e 00647 * lon should be in the range [−540°, 540°). 00648 **********************************************************************/ 00649 unsigned geod_polygon_testpoint(const struct geod_geodesic* g, 00650 const struct geod_polygon* p, 00651 double lat, double lon, 00652 int reverse, int sign, 00653 double* pA, double* pP); 00654 00655 /** 00656 * Return the results assuming a tentative final test point is added via an 00657 * azimuth and distance; however, the data for the test point is not saved. 00658 * This lets you report a running result for the perimeter and area as the 00659 * user moves the mouse cursor. Ordinary floating point arithmetic is used 00660 * to accumulate the data for the test point; thus the area and perimeter 00661 * returned are less accurate than if geod_polygon_addedge() and 00662 * geod_polygon_compute() are used. 00663 * 00664 * @param[in] g a pointer to the geod_geodesic object specifying the 00665 * ellipsoid. 00666 * @param[in] p a pointer to the geod_polygon object specifying the polygon. 00667 * @param[in] azi azimuth at current point (degrees). 00668 * @param[in] s distance from current point to final test point (meters). 00669 * @param[in] reverse if non-zero then clockwise (instead of 00670 * counter-clockwise) traversal counts as a positive area. 00671 * @param[in] sign if non-zero then return a signed result for the area if 00672 * the polygon is traversed in the "wrong" direction instead of returning 00673 * the area for the rest of the earth. 00674 * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>); 00675 * only set if \e polyline is non-zero in the call to geod_polygon_init(). 00676 * @param[out] pP pointer to the perimeter of the polygon or length of the 00677 * polyline (meters). 00678 * @return the number of points. 00679 * 00680 * \e azi should be in the range [−540°, 540°). 00681 **********************************************************************/ 00682 unsigned geod_polygon_testedge(const struct geod_geodesic* g, 00683 const struct geod_polygon* p, 00684 double azi, double s, 00685 int reverse, int sign, 00686 double* pA, double* pP); 00687 00688 /** 00689 * A simple interface for computing the area of a geodesic polygon. 00690 * 00691 * @param[in] g a pointer to the geod_geodesic object specifying the 00692 * ellipsoid. 00693 * @param[in] lats an array of latitudes of the polygon vertices (degrees). 00694 * @param[in] lons an array of longitudes of the polygon vertices (degrees). 00695 * @param[in] n the number of vertices. 00696 * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>). 00697 * @param[out] pP pointer to the perimeter of the polygon (meters). 00698 * 00699 * \e lats should be in the range [−90°, 90°]; \e lons should 00700 * be in the range [−540°, 540°). 00701 * 00702 * Only simple polygons (which are not self-intersecting) are allowed. 00703 * There's no need to "close" the polygon by repeating the first vertex. The 00704 * area returned is signed with counter-clockwise traversal being treated as 00705 * positive. 00706 * 00707 * Example, compute the area of Antarctic: 00708 @code 00709 double 00710 lats[] = {-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7, 00711 -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7}, 00712 lons[] = {-74, -102, -102, -131, -163, 163, 172, 140, 113, 00713 88, 59, 25, -4, -14, -33, -46, -61}; 00714 struct geod_geodesic g; 00715 double A, P; 00716 geod_init(&g, 6378137, 1/298.257223563); 00717 geod_polygonarea(&g, lats, lons, (sizeof lats) / (sizeof lats[0]), &A, &P); 00718 printf("%.0f %.2f\n", A, P); 00719 @endcode 00720 **********************************************************************/ 00721 void geod_polygonarea(const struct geod_geodesic* g, 00722 double lats[], double lons[], int n, 00723 double* pA, double* pP); 00724 00725 /** 00726 * mask values for the the \e caps argument to geod_lineinit(). 00727 **********************************************************************/ 00728 enum geod_mask { 00729 GEOD_NONE = 0U, /**< Calculate nothing */ 00730 GEOD_LATITUDE = 1U<<7 | 0U, /**< Calculate latitude */ 00731 GEOD_LONGITUDE = 1U<<8 | 1U<<3, /**< Calculate longitude */ 00732 GEOD_AZIMUTH = 1U<<9 | 0U, /**< Calculate azimuth */ 00733 GEOD_DISTANCE = 1U<<10 | 1U<<0, /**< Calculate distance */ 00734 GEOD_DISTANCE_IN = 1U<<11 | 1U<<0 | 1U<<1, /**< Allow distance as input */ 00735 GEOD_REDUCEDLENGTH= 1U<<12 | 1U<<0 | 1U<<2, /**< Calculate reduced length */ 00736 GEOD_GEODESICSCALE= 1U<<13 | 1U<<0 | 1U<<2, /**< Calculate geodesic scale */ 00737 GEOD_AREA = 1U<<14 | 1U<<4, /**< Calculate reduced length */ 00738 GEOD_ALL = 0x7F80U| 0x1FU /**< Calculate everything */ 00739 }; 00740 00741 #if defined(__cplusplus) 00742 } 00743 #endif 00744 00745 #endif