Update geographiclib-c to version 2.1 (#3)
This commit is contained in:
Stephen Booth
GitHub
parent
7ffdf978da
commit
21e5af22a3
|
|
@ -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");
|
||||
|
|
@ -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;
|
||||
|
|
@ -879,7 +879,7 @@ 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,
|
||||
|
|
@ -887,8 +887,12 @@ static double geod_geninverse_int(const struct geod_geodesic* g,
|
|||
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; }
|
||||
|
|
@ -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
|
||||
|
|
@ -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[]) {
|
||||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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);
|
||||
}
|
||||
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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);
|
||||
}
|
||||
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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.)
|
||||
|
|
|
|||
Loading…
Reference in New Issue