gooderp18绿色标准版
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

925 lignes
43KB

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