From 48236a6b737223cbf43a334236ded47ff4c8c9be Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Sun, 9 Apr 2017 10:18:07 -0400 Subject: Merge is geodesic routines from GeographicLib 1.48. Changes: - http://geographiclib.sf.net -> http://geographiclib.sourceforge.io - backport fixes for warnings messages from some compilers - change default range for longitude and azimuth to (-180d, 180d] (instead of [-180d, 180d)) --- man/man1/geod.1 | 8 ++++---- man/man3/geodesic.3 | 12 ++++++------ src/geodesic.c | 39 ++++++++++++++++++++------------------- src/geodesic.h | 20 ++++++++++---------- src/geodtest.c | 24 ++++++++++++------------ 5 files changed, 52 insertions(+), 51 deletions(-) diff --git a/man/man1/geod.1 b/man/man1/geod.1 index 00400a5c..74587612 100644 --- a/man/man1/geod.1 +++ b/man/man1/geod.1 @@ -2,7 +2,7 @@ .nr LL 7.0i .ad b .hy 1 -.TH GEOD 1 "2016/02/16 Rel. 4.9.3" +.TH GEOD 1 "2017/04/09 Rel. 4.9.3" .SH NAME geod \- direct geodesic computations .br @@ -207,7 +207,7 @@ the precision of the Portland location. .SH SEE ALSO .B geodesic(3) .PP -\fBGeographicLib\fR, http://geographiclib.sf.net +\fBGeographicLib\fR, https://geographiclib.sourceforge.io .PP The \fBGeodSolve\fR utility in GeographicLib. With the \fB-E\fR option, this solves the geodesic problems in terms of elliptic integrals; the @@ -219,11 +219,11 @@ J. Geodesy \fB87\fR, 43-55 (2013); .br DOI: https://doi.org/10.1007/s00190-012-0578-z .br -http://geographiclib.sf.net/geod-addenda.html +https://geographiclib.sourceforge.io/geod-addenda.html .PP The \fIonline geodesic bibliography\fR, .br -http://geographiclib.sf.net/geodesic-papers/biblio.html +https://geographiclib.sourceforge.io/geodesic-papers/biblio.html .SH BUGS A list of known bugs can found at https://github.com/OSGeo/proj.4/issues where new bug reports can be submitted too. diff --git a/man/man3/geodesic.3 b/man/man3/geodesic.3 index 78f6c0f3..310fc3d5 100644 --- a/man/man3/geodesic.3 +++ b/man/man3/geodesic.3 @@ -1,6 +1,6 @@ .\" @(#)geodesic.3 .nr LL 7.0i -.TH GEODESIC 3 "2016/02/16 Rel. 4.9.3" +.TH GEODESIC 3 "2017/04/09 Rel. 4.9.3" .ad b .hy 1 .SH NAME @@ -53,7 +53,7 @@ measure angles (latitudes, longitudes, and azimuths) in degrees, unlike the rest of the \fBproj\fR library, which uses radians. The documentation for this library is included in geodesic.h. A formatted version of the documentation is available at -http://geographiclib.sf.net/1.46/C +https://geographiclib.sourceforge.io/1.48/C .SH EXAMPLE The following program reads in lines with the coordinates for two points in decimal degrees (\fIlat1\fR, \fIlon1\fR, \fIlat2\fR, \fIlon2\fR) and @@ -87,11 +87,11 @@ libproj.a \- library of projections and support procedures .SH SEE ALSO Full online documentation for \fBgeodesic(3)\fR, .br -http://geographiclib.sf.net/1.46/C +https://geographiclib.sourceforge.io/1.48/C .PP .B geod(1) .PP -\fBGeographicLib\fR, http://geographiclib.sf.net +\fBGeographicLib\fR, https://geographiclib.sourceforge.io .PP The \fBGeodesicExact\fR class in GeographicLib solves the geodesic problems in terms of elliptic integrals; the results are accurate for @@ -103,11 +103,11 @@ J. Geodesy \fB87\fR, 43-55 (2013); .br DOI: https://doi.org/10.1007/s00190-012-0578-z .br -http://geographiclib.sf.net/geod-addenda.html +https://geographiclib.sourceforge.io/geod-addenda.html .PP \fIA geodesic bibliography\fR, .br -http://geographiclib.sf.net/geodesic-papers/biblio.html +https://geographiclib.sourceforge.io/geodesic-papers/biblio.html .PP The Wikipedia page, \fIGeodesics on an ellipsoid\fR, .br diff --git a/src/geodesic.c b/src/geodesic.c index d93268c6..927e59eb 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -14,13 +14,13 @@ * Algorithms for geodesics, * J. Geodesy 87, 43--55 (2013); * https://doi.org/10.1007/s00190-012-0578-z - * Addenda: http://geographiclib.sourceforge.net/geod-addenda.html + * Addenda: https://geographiclib.sourceforge.io/geod-addenda.html * * See the comments in geodesic.h for documentation. * * Copyright (c) Charles Karney (2012-2017) and licensed * under the MIT/X11 License. For more information, see - * http://geographiclib.sourceforge.net/ + * https://geographiclib.sourceforge.io/ */ #include "geodesic.h" @@ -86,8 +86,8 @@ static void Init() { xthresh = 1000 * tol2; degree = pi/180; { - double minus1 = -1.0; - NaN = sqrt(minus1); + real minus1 = -1; + NaN = sqrt(minus1); } init = 1; } @@ -171,21 +171,21 @@ static void norm2(real* sinx, real* cosx) { static real AngNormalize(real x) { x = fmod(x, (real)(360)); - return x < -180 ? x + 360 : (x < 180 ? x : x - 360); + return x <= -180 ? x + 360 : (x <= 180 ? x : x - 360); } static real LatFix(real x) { return fabs(x) > 90 ? NaN : x; } static real AngDiff(real x, real y, real* e) { - real t, d = - AngNormalize(sumx(AngNormalize(x), AngNormalize(-y), &t)); - /* Here y - x = d - t (mod 360), exactly, where d is in (-180,180] and + real t, d = AngNormalize(sumx(AngNormalize(-x), AngNormalize(y), &t)); + /* Here y - x = d + t (mod 360), exactly, where d is in (-180,180] and * abs(t) <= eps (eps = 2^-45 for doubles). The only case where the * addition of t takes the result outside the range (-180,180] is d = 180 - * and t < 0. The case, d = -180 + eps, t = eps, can't happen, since + * and t > 0. The case, d = -180 + eps, t = -eps, can't happen, since * sum would have returned the exact result in such a case (i.e., given t * = 0). */ - return sumx(d == 180 && t < 0 ? -180 : d, -t, e); + return sumx(d == 180 && t > 0 ? -180 : d, t, e); } static real AngRound(real x) { @@ -210,11 +210,12 @@ static void sincosdx(real x, real* sinx, real* cosx) { /* Possibly could call the gnu extension sincos */ s = sin(r); c = cos(r); switch ((unsigned)q & 3U) { - case 0U: *sinx = s; *cosx = c; break; - case 1U: *sinx = c; *cosx = 0 - s; break; - case 2U: *sinx = 0 - s; *cosx = 0 - c; break; - default: *sinx = 0 - c; *cosx = s; break; /* case 3U */ + case 0U: *sinx = s; *cosx = c; break; + case 1U: *sinx = c; *cosx = -s; break; + case 2U: *sinx = -s; *cosx = -c; break; + default: *sinx = -c; *cosx = s; break; /* case 3U */ } + if (x) { *sinx += (real)(0); *cosx += (real)(0); } } static real atan2dx(real y, real x) { @@ -233,7 +234,7 @@ static real atan2dx(real y, real x) { * * case 0: ang = 0 + ang; break; */ - case 1: ang = (y > 0 ? 180 : -180) - ang; break; + case 1: ang = (y >= 0 ? 180 : -180) - ang; break; case 2: ang = 90 - ang; break; case 3: ang = -90 + ang; break; } @@ -1752,8 +1753,8 @@ int transit(real lon1, real lon2) { lon1 = AngNormalize(lon1); lon2 = AngNormalize(lon2); lon12 = AngDiff(lon1, lon2, 0); - return lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 : - (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0); + return lon1 <= 0 && lon2 > 0 && lon12 > 0 ? 1 : + (lon2 <= 0 && lon1 > 0 && lon12 < 0 ? -1 : 0); } int transitdirect(real lon1, real lon2) { @@ -1816,7 +1817,7 @@ void geod_polygon_addpoint(const struct geod_geodesic* g, p->lat0 = p->lat = lat; p->lon0 = p->lon = lon; } else { - real s12, S12; + real s12, S12 = 0; /* Initialize S12 to stop Visual Studio warning */ geod_geninverse(g, p->lat, p->lon, lat, lon, &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12); accadd(p->P, s12); @@ -1833,7 +1834,7 @@ void geod_polygon_addedge(const struct geod_geodesic* g, struct geod_polygon* p, real azi, real s) { if (p->num) { /* Do nothing is num is zero */ - real lat, lon, S12; + real lat, lon, S12 = 0; /* Initialize S12 to stop Visual Studio warning */ geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s, &lat, &lon, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12); @@ -1908,7 +1909,7 @@ unsigned geod_polygon_testpoint(const struct geod_geodesic* g, tempsum = p->polyline ? 0 : p->A[0]; crossings = p->crossings; for (i = 0; i < (p->polyline ? 1 : 2); ++i) { - real s12, S12; + real s12, S12 = 0; /* Initialize S12 to stop Visual Studio warning */ geod_geninverse(g, i == 0 ? p->lat : lat, i == 0 ? p->lon : lon, i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon, diff --git a/src/geodesic.h b/src/geodesic.h index cd8989db..b864ccf4 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -9,7 +9,7 @@ * J. Geodesy 87, 43--55 (2013); * DOI: * 10.1007/s00190-012-0578-z; - * addenda: + * addenda: * geod-addenda.html. * . * The principal advantages of these algorithms over previous ones (e.g., @@ -96,7 +96,7 @@ * [\e d, \e d], for arbitrary \e d. * * These routines are a simple transcription of the corresponding C++ classes - * in GeographicLib. The + * in GeographicLib. The * "class data" is represented by the structs geod_geodesic, geod_geodesicline, * geod_polygon and pointers to these objects are passed as initial arguments * to the member functions. Most of the internal comments have been retained. @@ -107,12 +107,12 @@ * twice about restructuring the internals of the C code since this may make * porting fixes from the C++ code more difficult. * - * Copyright (c) Charles Karney (2012-2016) and licensed + * Copyright (c) Charles Karney (2012-2017) and licensed * under the MIT/X11 License. For more information, see - * http://geographiclib.sourceforge.net/ + * https://geographiclib.sourceforge.io/ * * This library was distributed with - * GeographicLib 1.47. + * GeographicLib 1.48. **********************************************************************/ #if !defined(GEODESIC_H) @@ -127,7 +127,7 @@ * The minor version of the geodesic library. (This tracks the version of * GeographicLib.) **********************************************************************/ -#define GEODESIC_VERSION_MINOR 47 +#define GEODESIC_VERSION_MINOR 48 /** * The patch level of the geodesic library. (This tracks the version of * GeographicLib.) @@ -242,7 +242,7 @@ extern "C" { * * \e g must have been initialized with a call to geod_init(). \e lat1 * should be in the range [−90°, 90°]. The values of \e lon2 - * and \e azi2 returned are in the range [−180°, 180°). Any of + * 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. * @@ -327,7 +327,7 @@ extern "C" { * * \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°). + * \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 * do not need some quantities computed. * @@ -525,7 +525,7 @@ extern "C" { * * \e l must have been initialized with a call, e.g., to geod_lineinit(), * with \e caps |= GEOD_DISTANCE_IN. The values of \e lon2 and \e azi2 - * returned are in the range [−180°, 180°). Any of the + * 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. * @@ -594,7 +594,7 @@ extern "C" { * * \e l must have been initialized with a call to geod_lineinit() with \e * caps |= GEOD_DISTANCE_IN. The value \e azi2 returned is in the range - * [−180°, 180°). Any of the "return" arguments \e plat2, + * [−180°, 180°]. Any of the "return" arguments \e plat2, * etc., may be replaced by 0, if you do not need some quantities * computed. Requesting a value which \e l is not capable of computing * is not an error; the corresponding argument will not be altered. diff --git a/src/geodtest.c b/src/geodtest.c index 7b36f64b..5ca741b1 100644 --- a/src/geodtest.c +++ b/src/geodtest.c @@ -6,7 +6,7 @@ * * Copyright (c) Charles Karney (2015-2017) and licensed * under the MIT/X11 License. For more information, see - * http://geographiclib.sourceforge.net/ + * https://geographiclib.sourceforge.io/ **********************************************************************/ /** @cond SKIP */ @@ -254,7 +254,7 @@ static int GeodSolve5() { result += assertEquals(lat2, 90, 0.5e-5); if (lon2 < 0) { result += assertEquals(lon2, -150, 0.5e-5); - result += assertEquals(azi2, -180, 0.5e-5); + result += assertEquals(fabs(azi2), 180, 0.5e-5); } else { result += assertEquals(lon2, 30, 0.5e-5); result += assertEquals(azi2, 0, 0.5e-5); @@ -340,7 +340,7 @@ static int GeodSolve14() { struct geod_geodesic g; int result = 0; { - double minus1 = -1.0; + double minus1 = -1; nan = sqrt(minus1); } geod_init(&g, wgs84_a, wgs84_f); @@ -434,11 +434,11 @@ static int GeodSolve33() { result += assertEquals(s12, 19980862, 0.5); geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2); result += assertEquals(azi1, 0.00000, 0.5e-5); - result += assertEquals(azi2, -180.00000, 0.5e-5); + result += assertEquals(fabs(azi2), 180.00000, 0.5e-5); result += assertEquals(s12, 20003931, 0.5); geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2); result += assertEquals(azi1, 0.00000, 0.5e-5); - result += assertEquals(azi2, -180.00000, 0.5e-5); + result += assertEquals(fabs(azi2), 180.00000, 0.5e-5); result += assertEquals(s12, 19893357, 0.5); geod_init(&g, 6.4e6, 0); geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2); @@ -447,11 +447,11 @@ static int GeodSolve33() { result += assertEquals(s12, 19994492, 0.5); geod_inverse(&g, 0, 0, 0, 180, &s12, &azi1, &azi2); result += assertEquals(azi1, 0.00000, 0.5e-5); - result += assertEquals(azi2, -180.00000, 0.5e-5); + result += assertEquals(fabs(azi2), 180.00000, 0.5e-5); result += assertEquals(s12, 20106193, 0.5); geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2); result += assertEquals(azi1, 0.00000, 0.5e-5); - result += assertEquals(azi2, -180.00000, 0.5e-5); + result += assertEquals(fabs(azi2), 180.00000, 0.5e-5); result += assertEquals(s12, 19994492, 0.5); geod_init(&g, 6.4e6, -1/300.0); geod_inverse(&g, 0, 0, 0, 179, &s12, &azi1, &azi2); @@ -468,7 +468,7 @@ static int GeodSolve33() { result += assertEquals(s12, 20082617, 0.5); geod_inverse(&g, 0, 0, 1, 180, &s12, &azi1, &azi2); result += assertEquals(azi1, 0.00000, 0.5e-5); - result += assertEquals(azi2, -180.00000, 0.5e-5); + result += assertEquals(fabs(azi2), 180.00000, 0.5e-5); result += assertEquals(s12, 20027270, 0.5); return result; @@ -481,7 +481,7 @@ static int GeodSolve55() { struct geod_geodesic g; int result = 0; { - double minus1 = -1.0; + double minus1 = -1; nan = sqrt(minus1); } geod_init(&g, wgs84_a, wgs84_f); @@ -521,12 +521,12 @@ static int GeodSolve61() { &lat2, &lon2, &azi2, 0, 0, 0, 0, 0); result += assertEquals(lat2, 45.30632, 0.5e-5); result += assertEquals(lon2, -180, 0.5e-5); - result += assertEquals(azi2, -180, 0.5e-5); + result += assertEquals(fabs(azi2), 180, 0.5e-5); geod_inverseline(&l, &g, 45, 0, 80, -0.000000000000000003, 0); geod_genposition(&l, flags, 1e7, &lat2, &lon2, &azi2, 0, 0, 0, 0, 0); result += assertEquals(lat2, 45.30632, 0.5e-5); result += assertEquals(lon2, -180, 0.5e-5); - result += assertEquals(azi2, -180, 0.5e-5); + result += assertEquals(fabs(azi2), 180, 0.5e-5); return result; } @@ -545,7 +545,7 @@ static int GeodSolve65() { &lat2, &lon2, &azi2, &s12, &m12, &M12, &M21, &S12); result += assertEquals(lat2, -60.23169, 0.5e-5); result += assertEquals(lon2, -0.00000, 0.5e-5); - result += assertEquals(azi2, -180.00000, 0.5e-5); + result += assertEquals(fabs(azi2), 180.00000, 0.5e-5); result += assertEquals(s12, 10000000, 0.5); result += assertEquals(a12, 90.06544, 0.5e-5); result += assertEquals(m12, 6363636, 0.5); -- cgit v1.2.3