diff options
| author | Charles Karney <charles@karney.com> | 2015-08-16 15:05:36 -0400 |
|---|---|---|
| committer | Charles Karney <charles@karney.com> | 2015-08-16 15:05:36 -0400 |
| commit | 8066dcd0e9ce33222b167dbb2a8ddab79465d299 (patch) | |
| tree | 617a14fc6f9e8301a75b8cb805d50b50c7948146 /src/geodesic.h | |
| parent | 6c16367e152747133bba7a8fbcdbabef1232cd93 (diff) | |
| download | PROJ-8066dcd0e9ce33222b167dbb2a8ddab79465d299.tar.gz PROJ-8066dcd0e9ce33222b167dbb2a8ddab79465d299.zip | |
Drop in the latest geodesic library from GeographicLib (version 1.44).
http://geographiclib.sourceforge.net/1.44/C/index.html
The changes are:
- Improve accuracy of calculations by evaluating trigonometric
functions more carefully and replacing the series for the reduced
length with one with a smaller truncation error.
- The allowed ranges for longitudes and azimuths is now unlimited; it
used to be [-540d, 540d).
- Enforce the restriction of latitude to [-90d, 90d] by returning NaNs
if the latitude is outside this range.
- The inverse calculation sets s12 to zero for coincident points at
pole (instead of returning a tiny quantity).
This commit also includes a work-around for an inaccurate value for
pi/180 in dmstor.c (see the definitions of DEG_IN and DEG_OUT in
geod_interface.c).
Diffstat (limited to 'src/geodesic.h')
| -rw-r--r-- | src/geodesic.h | 77 |
1 files changed, 30 insertions, 47 deletions
diff --git a/src/geodesic.h b/src/geodesic.h index 1a1892a7..7bd8270f 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -113,7 +113,7 @@ * http://geographiclib.sourceforge.net/ * * This library was distributed with - * <a href="../index.html">GeographicLib</a> 1.43. + * <a href="../index.html">GeographicLib</a> 1.44. **********************************************************************/ #if !defined(GEODESIC_H) @@ -128,7 +128,7 @@ * The minor version of the geodesic library. (This tracks the version of * GeographicLib.) **********************************************************************/ -#define GEODESIC_VERSION_MINOR 43 +#define GEODESIC_VERSION_MINOR 44 /** * The patch level of the geodesic library. (This tracks the version of * GeographicLib.) @@ -147,11 +147,11 @@ * where MM is the major version, mmmm is the minor version, and pp is the * patch level. Users should not rely on this particular packing of the * components of the version number. Instead they should use a test such as - * \code + * @code{.c} #if GEODESIC_VERSION >= GEODESIC_VERSION_NUM(1,40,0) ... #endif - * \endcode + * @endcode **********************************************************************/ #define GEODESIC_VERSION \ GEODESIC_VERSION_NUM(GEODESIC_VERSION_MAJOR, \ @@ -237,8 +237,7 @@ extern "C" { * geod_genposition(). * * \e g must have been initialized with a call to geod_init(). \e lat1 - * should be in the range [−90°, 90°]; \e lon1 and \e azi1 - * should be in the range [−540°, 540°). + * should be in the range [−90°, 90°]. * * The geod_mask values are [see geod_mask()]: * - \e caps |= GEOD_LATITUDE for the latitude \e lat2; this is @@ -278,8 +277,7 @@ extern "C" { * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). * * \e g must have been initialized with a call to geod_init(). \e lat1 - * should be in the range [−90°, 90°]; \e lon1 and \e azi1 - * should be in the range [−540°, 540°). The values of \e lon2 + * should be in the range [−90°, 90°]. The values of \e lon2 * and \e azi2 returned are in the range [−180°, 180°). Any of * the "return" arguments \e plat2, etc., may be replaced by 0, if you do not * need some quantities computed. @@ -292,7 +290,7 @@ extern "C" { * longitudinal extent must not exceed of 180°.) * * Example, determine the point 10000 km NE of JFK: - @code + @code{.c} struct geod_geodesic g; double lat, lon; geod_init(&g, 6378137, 1/298.257223563); @@ -318,11 +316,10 @@ extern "C" { * @param[out] pazi1 pointer to the azimuth at point 1 (degrees). * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees). * - * \e g must have been initialized with a call to geod_init(). \e lat1 - * and \e lat2 should be in the range [−90°, 90°]; \e lon1 and - * \e lon2 should be in the range [−540°, 540°). The values of + * \e g must have been initialized with a call to geod_init(). \e lat1 and + * \e lat2 should be in the range [−90°, 90°]. The values of * \e azi1 and \e azi2 returned are in the range [−180°, 180°). - * Any of the "return" arguments \e ps12, etc., may be replaced by 0, if you + * Any of the "return" arguments, \e ps12, etc., may be replaced by 0, if you * do not need some quantities computed. * * If either point is at a pole, the azimuth is defined by keeping the @@ -335,7 +332,7 @@ extern "C" { * is used to refine the solution. * * Example, determine the distance between JFK and Singapore Changi Airport: - @code + @code{.c} struct geod_geodesic g; double s12; geod_init(&g, 6378137, 1/298.257223563); @@ -367,7 +364,7 @@ extern "C" { * * Example, compute way points between JFK and Singapore Changi Airport * the "obvious" way using geod_direct(): - @code + @code{.c} struct geod_geodesic g; double s12, azi1, lat[101],lon[101]; int i; @@ -379,7 +376,7 @@ extern "C" { } @endcode * A faster way using geod_position(): - @code + @code{.c} struct geod_geodesic g; struct geod_geodesicline l; double s12, azi1, lat[101],lon[101]; @@ -425,18 +422,14 @@ extern "C" { * @return \e a12 arc length of between point 1 and point 2 (degrees). * * \e g must have been initialized with a call to geod_init(). \e lat1 - * should be in the range [−90°, 90°]; \e lon1 and \e azi1 - * should be in the range [−540°, 540°). The function - * value \e a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE. Any of the - * "return" arguments \e plat2, etc., may be replaced by 0, if you do not - * need some quantities computed. + * should be in the range [−90°, 90°]. The function value \e + * a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE. Any of the "return" + * arguments, \e plat2, etc., may be replaced by 0, if you do not need some + * quantities computed. * * With \e flags & GEOD_LONG_UNROLL bit set, the longitude is "unrolled" so * that the quantity \e lon2 − \e lon1 indicates how many times and in - * what sense the geodesic encircles the ellipsoid. Because \e lon2 might be - * outside the normal allowed range for longitudes, [−540°, - * 540°), be sure to normalize it, e.g., with fmod(\e lon2, 360.0) before - * using it in subsequent calculations + * what sense the geodesic encircles the ellipsoid. **********************************************************************/ double geod_gendirect(const struct geod_geodesic* g, double lat1, double lon1, double azi1, @@ -467,9 +460,8 @@ extern "C" { * (meters<sup>2</sup>). * @return \e a12 arc length of between point 1 and point 2 (degrees). * - * \e g must have been initialized with a call to geod_init(). \e lat1 - * and \e lat2 should be in the range [−90°, 90°]; \e lon1 and - * \e lon2 should be in the range [−540°, 540°). Any of the + * \e g must have been initialized with a call to geod_init(). \e lat1 and + * \e lat2 should be in the range [−90°, 90°]. Any of the * "return" arguments \e ps12, etc., may be replaced by 0, if you do not need * some quantities computed. **********************************************************************/ @@ -521,17 +513,14 @@ extern "C" { * * With \e flags & GEOD_LONG_UNROLL bit set, the longitude is "unrolled" so * that the quantity \e lon2 − \e lon1 indicates how many times and in - * what sense the geodesic encircles the ellipsoid. Because \e lon2 might be - * outside the normal allowed range for longitudes, [−540°, - * 540°), be sure to normalize it, e.g., with fmod(\e lon2, 360.0) before - * using it in subsequent calculations + * what sense the geodesic encircles the ellipsoid. * * Example, compute way points between JFK and Singapore Changi Airport * using geod_genposition(). In this example, the points are evenly space in * arc length (and so only approximately equally space in distance). This is * faster than using geod_position() would be appropriate if drawing the path * on a map. - @code + @code{.c} struct geod_geodesic g; struct geod_geodesicline l; double a12, azi1, lat[101], lon[101]; @@ -588,8 +577,7 @@ extern "C" { * \e g and \e p must have been initialized with calls to geod_init() and * geod_polygon_init(), respectively. The same \e g must be used for all the * points and edges in a polygon. \e lat should be in the range - * [−90°, 90°] and \e lon should be in the range - * [−540°, 540°). + * [−90°, 90°]. * * An example of the use of this function is given in the documentation for * geod_polygon_compute(). @@ -610,10 +598,9 @@ extern "C" { * * \e g and \e p must have been initialized with calls to geod_init() and * geod_polygon_init(), respectively. The same \e g must be used for all the - * points and edges in a polygon. \e azi should be in the range - * [−540°, 540°). This does nothing if no points have been - * added yet. The \e lat and \e lon fields of \e p give the location of - * the new vertex. + * points and edges in a polygon. This does nothing if no points have been + * added yet. The \e lat and \e lon fields of \e p give the location of the + * new vertex. **********************************************************************/ void geod_polygon_addedge(const struct geod_geodesic* g, struct geod_polygon* p, @@ -645,7 +632,7 @@ extern "C" { * * Example, compute the perimeter and area of the geodesic triangle with * vertices (0°N,0°E), (0°N,90°E), (90°N,0°E). - @code + @code{.c} double A, P; int n; struct geod_geodesic g; @@ -689,8 +676,7 @@ extern "C" { * polyline (meters). * @return the number of points. * - * \e lat should be in the range [−90°, 90°] and \e - * lon should be in the range [−540°, 540°). + * \e lat should be in the range [−90°, 90°]. **********************************************************************/ unsigned geod_polygon_testpoint(const struct geod_geodesic* g, const struct geod_polygon* p, @@ -722,8 +708,6 @@ extern "C" { * @param[out] pP pointer to the perimeter of the polygon or length of the * polyline (meters). * @return the number of points. - * - * \e azi should be in the range [−540°, 540°). **********************************************************************/ unsigned geod_polygon_testedge(const struct geod_geodesic* g, const struct geod_polygon* p, @@ -742,8 +726,7 @@ extern "C" { * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>). * @param[out] pP pointer to the perimeter of the polygon (meters). * - * \e lats should be in the range [−90°, 90°]; \e lons should - * be in the range [−540°, 540°). + * \e lats should be in the range [−90°, 90°]. * * Only simple polygons (which are not self-intersecting) are allowed. * There's no need to "close" the polygon by repeating the first vertex. The @@ -751,7 +734,7 @@ extern "C" { * positive. * * Example, compute the area of Antarctica: - @code + @code{.c} double lats[] = {-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7, -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7}, |