Skip to content

Commit

Permalink
Merge pull request #3528 from cffk/geod-2.1
Browse files Browse the repository at this point in the history 
Synchronize with 2.1 for geographiclib-c for geodesic support

Nothing much to cause concern with the PR.  Plus it simplifies matters for maintainers of proj.4-feedstock.
  • Loading branch information
@cffk
cffk committed Jan 4, 2023
2 parents b23965e + fec14d6 commit 83a0bf0
Show file tree
Hide file tree
Showing 4 changed files with 55 additions and 30 deletions.
62 changes: 34 additions & 28 deletions src/geodesic.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ static void Init(void) {
tol1 = 200 * tol0;
tol2 = sqrt(tol0);
/* Check on bisection interval */
tolb = tol0 * tol2;
tolb = tol0;
xthresh = 1000 * tol2;
degree = pi/hd;
NaN = nan("0");
Expand Down Expand Up @@ -111,7 +111,7 @@ static double sumx(double u, double v, double* t) {
return s;
}

static double polyval(int N, const double p[], double x) {
static double polyvalx(int N, const double p[], double x) {
double y = N < 0 ? 0 : *p++;
while (--N >= 0) y = y * x + *p++;
return y;
Expand Down Expand Up @@ -879,16 +879,20 @@ static double geod_geninverse_int(const struct geod_geodesic* g,
double salp1a = tiny, calp1a = 1, salp1b = tiny, calp1b = -1;
boolx tripn = FALSE;
boolx tripb = FALSE;
for (; numit < maxit2; ++numit) {
for (;; ++numit) {
/* the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
* WGS84 and random input: mean = 2.85, sd = 0.60 */
double dv = 0,
v = Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
slam12, clam12,
&salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2,
&eps, &domg12, numit < maxit1, &dv, Ca);
/* Reversed test to allow escape with NaNs */
if (tripb || !(fabs(v) >= (tripn ? 8 : 1) * tol0)) break;
if (tripb ||
/* Reversed test to allow escape with NaNs */
!(fabs(v) >= (tripn ? 8 : 1) * tol0) ||
/* Enough bisections to get accurate result */
numit == maxit2)
break;
/* Update bracketing values */
if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b))
{ salp1b = salp1; calp1b = calp1; }
Expand All @@ -897,18 +901,20 @@ static double geod_geninverse_int(const struct geod_geodesic* g,
if (numit < maxit1 && dv > 0) {
double
dalp1 = -v/dv;
double
sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
if (nsalp1 > 0 && fabs(dalp1) < pi) {
calp1 = calp1 * cdalp1 - salp1 * sdalp1;
salp1 = nsalp1;
norm2(&salp1, &calp1);
/* In some regimes we don't get quadratic convergence because
* slope -> 0. So use convergence conditions based on epsilon
* instead of sqrt(epsilon). */
tripn = fabs(v) <= 16 * tol0;
continue;
if (fabs(dalp1) < pi) {
double
sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),
nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;
if (nsalp1 > 0) {
calp1 = calp1 * cdalp1 - salp1 * sdalp1;
salp1 = nsalp1;
norm2(&salp1, &calp1);
/* In some regimes we don't get quadratic convergence because
* slope -> 0. So use convergence conditions based on epsilon
* instead of sqrt(epsilon). */
tripn = fabs(v) <= 16 * tol0;
continue;
}
}
}
/* Either dv was not positive or updated value was outside legal
Expand Down Expand Up @@ -1480,7 +1486,7 @@ double Lambda12(const struct geod_geodesic* g,

double A3f(const struct geod_geodesic* g, double eps) {
/* Evaluate A3 */
return polyval(nA3 - 1, g->A3x, eps);
return polyvalx(nA3 - 1, g->A3x, eps);
}

void C3f(const struct geod_geodesic* g, double eps, double c[]) {
Expand All @@ -1491,7 +1497,7 @@ void C3f(const struct geod_geodesic* g, double eps, double c[]) {
for (l = 1; l < nC3; ++l) { /* l is index of C3[l] */
int m = nC3 - l - 1; /* order of polynomial in eps */
mult *= eps;
c[l] = mult * polyval(m, g->C3x + o, eps);
c[l] = mult * polyvalx(m, g->C3x + o, eps);
o += m + 1;
}
}
Expand All @@ -1503,7 +1509,7 @@ void C4f(const struct geod_geodesic* g, double eps, double c[]) {
int o = 0, l;
for (l = 0; l < nC4; ++l) { /* l is index of C4[l] */
int m = nC4 - l - 1; /* order of polynomial in eps */
c[l] = mult * polyval(m, g->C4x + o, eps);
c[l] = mult * polyvalx(m, g->C4x + o, eps);
o += m + 1;
mult *= eps;
}
Expand All @@ -1516,7 +1522,7 @@ double A1m1f(double eps) {
1, 4, 64, 0, 256,
};
int m = nA1/2;
double t = polyval(m, coeff, sq(eps)) / coeff[m + 1];
double t = polyvalx(m, coeff, sq(eps)) / coeff[m + 1];
return (t + eps) / (1 - eps);
}

Expand All @@ -1542,7 +1548,7 @@ void C1f(double eps, double c[]) {
int o = 0, l;
for (l = 1; l <= nC1; ++l) { /* l is index of C1p[l] */
int m = (nC1 - l) / 2; /* order of polynomial in eps^2 */
c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
c[l] = d * polyvalx(m, coeff + o, eps2) / coeff[o + m + 1];
o += m + 2;
d *= eps;
}
Expand Down Expand Up @@ -1570,7 +1576,7 @@ void C1pf(double eps, double c[]) {
int o = 0, l;
for (l = 1; l <= nC1p; ++l) { /* l is index of C1p[l] */
int m = (nC1p - l) / 2; /* order of polynomial in eps^2 */
c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
c[l] = d * polyvalx(m, coeff + o, eps2) / coeff[o + m + 1];
o += m + 2;
d *= eps;
}
Expand All @@ -1583,7 +1589,7 @@ double A2m1f(double eps) {
-11, -28, -192, 0, 256,
};
int m = nA2/2;
double t = polyval(m, coeff, sq(eps)) / coeff[m + 1];
double t = polyvalx(m, coeff, sq(eps)) / coeff[m + 1];
return (t - eps) / (1 + eps);
}

Expand All @@ -1609,7 +1615,7 @@ void C2f(double eps, double c[]) {
int o = 0, l;
for (l = 1; l <= nC2; ++l) { /* l is index of C2[l] */
int m = (nC2 - l) / 2; /* order of polynomial in eps^2 */
c[l] = d * polyval(m, coeff + o, eps2) / coeff[o + m + 1];
c[l] = d * polyvalx(m, coeff + o, eps2) / coeff[o + m + 1];
o += m + 2;
d *= eps;
}
Expand All @@ -1634,7 +1640,7 @@ void A3coeff(struct geod_geodesic* g) {
int o = 0, k = 0, j;
for (j = nA3 - 1; j >= 0; --j) { /* coeff of eps^j */
int m = nA3 - j - 1 < j ? nA3 - j - 1 : j; /* order of polynomial in n */
g->A3x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
g->A3x[k++] = polyvalx(m, coeff + o, g->n) / coeff[o + m + 1];
o += m + 2;
}
}
Expand Down Expand Up @@ -1677,7 +1683,7 @@ void C3coeff(struct geod_geodesic* g) {
for (l = 1; l < nC3; ++l) { /* l is index of C3[l] */
for (j = nC3 - 1; j >= l; --j) { /* coeff of eps^j */
int m = nC3 - j - 1 < j ? nC3 - j - 1 : j; /* order of polynomial in n */
g->C3x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
g->C3x[k++] = polyvalx(m, coeff + o, g->n) / coeff[o + m + 1];
o += m + 2;
}
}
Expand Down Expand Up @@ -1733,7 +1739,7 @@ void C4coeff(struct geod_geodesic* g) {
for (l = 0; l < nC4; ++l) { /* l is index of C4[l] */
for (j = nC4 - 1; j >= l; --j) { /* coeff of eps^j */
int m = nC4 - j - 1; /* order of polynomial in n */
g->C4x[k++] = polyval(m, coeff + o, g->n) / coeff[o + m + 1];
g->C4x[k++] = polyvalx(m, coeff + o, g->n) / coeff[o + m + 1];
o += m + 2;
}
}
Expand Down
2 changes: 1 addition & 1 deletion src/geodesic.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@
* The minor version of the geodesic library. (This tracks the version of
* GeographicLib.)
**********************************************************************/
#define GEODESIC_VERSION_MINOR 0
#define GEODESIC_VERSION_MINOR 1
/**
* The patch level of the geodesic library. (This tracks the version of
* GeographicLib.)
Expand Down
2 changes: 1 addition & 1 deletion src/init.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -772,7 +772,7 @@ pj_init_ctx_with_allow_init_epsg(PJ_CONTEXT *ctx, int argc, char **argv, int all
PIN->geod = static_cast<struct geod_geodesic*>(calloc (1, sizeof (struct geod_geodesic)));
if (nullptr==PIN->geod)
return pj_default_destructor (PIN, PROJ_ERR_OTHER /*ENOMEM*/);
geod_init(PIN->geod, PIN->a, PIN->es / (1 + sqrt (PIN->one_es)));
geod_init(PIN->geod, PIN->a, PIN->f);

/* Projection specific initialization */
err = proj_errno_reset (PIN);
Expand Down
19 changes: 19 additions & 0 deletions src/tests/geodsigntest.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,19 @@ typedef double T;
#define nullptr 0
#endif

#if !defined(OLD_BUGGY_REMQUO)
/*
* glibc prior to version 2.22 had a bug in remquo. This was reported in 2014
* and fixed in 2015. See
* https://sourceware.org/bugzilla/show_bug.cgi?id=17569
*
* The bug causes some of the tests here to fail. The failures aren't terribly
* serious (just a loss of accuracy). If you're still using the buggy glibc,
* then define OLD_BUGGY_REMQUO to be 1.
*/
#define OLD_BUGGY_REMQUO 0
#endif

static const T wgs84_a = 6378137, wgs84_f = 1/298.257223563; /* WGS84 */

static int equiv(T x, T y) {
Expand Down Expand Up @@ -123,7 +136,9 @@ int main() {
check( geod_AngRound( 90.0 ), 90 );

checksincosd(- inf, nan, nan);
#if !OLD_BUGGY_REMQUO
checksincosd(-810.0, -1.0, +0.0);
#endif
checksincosd(-720.0, -0.0, +1.0);
checksincosd(-630.0, +1.0, +0.0);
checksincosd(-540.0, -0.0, -1.0);
Expand All @@ -142,10 +157,13 @@ int main() {
checksincosd(+540.0, +0.0, -1.0);
checksincosd(+630.0, -1.0, +0.0);
checksincosd(+720.0, +0.0, +1.0);
#if !OLD_BUGGY_REMQUO
checksincosd(+810.0, +1.0, +0.0);
#endif
checksincosd(+ inf, nan, nan);
checksincosd( nan, nan, nan);

#if !OLD_BUGGY_REMQUO
{
T s1, c1, s2, c2, s3, c3;
geod_sincosd( 9.0, &s1, &c1);
Expand All @@ -156,6 +174,7 @@ int main() {
++n;
}
}
#endif

check( geod_atan2d(+0.0 , -0.0 ), +180 );
check( geod_atan2d(-0.0 , -0.0 ), -180 );
Expand Down

0 comments on commit 83a0bf0

Please sign in to comment.