aboutsummaryrefslogtreecommitdiff
path: root/src/geodesic.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/geodesic.c')
-rw-r--r--src/geodesic.c380
1 files changed, 239 insertions, 141 deletions
diff --git a/src/geodesic.c b/src/geodesic.c
index fd0214c7..9a8f043c 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -18,7 +18,7 @@
*
* See the comments in geodesic.h for documentation.
*
- * Copyright (c) Charles Karney (2012-2014) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2015) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/
*/
@@ -27,8 +27,10 @@
#include <math.h>
#define GEOGRAPHICLIB_GEODESIC_ORDER 6
+#define nA1 GEOGRAPHICLIB_GEODESIC_ORDER
#define nC1 GEOGRAPHICLIB_GEODESIC_ORDER
#define nC1p GEOGRAPHICLIB_GEODESIC_ORDER
+#define nA2 GEOGRAPHICLIB_GEODESIC_ORDER
#define nC2 GEOGRAPHICLIB_GEODESIC_ORDER
#define nA3 GEOGRAPHICLIB_GEODESIC_ORDER
#define nA3x nA3
@@ -137,6 +139,12 @@ static real sumx(real u, real v, real* t) {
return s;
}
+static real polyval(int N, const real p[], real x) {
+ real y = N < 0 ? 0 : *p++;
+ while (--N >= 0) y = y * x + *p++;
+ return y;
+}
+
static real minx(real x, real y)
{ return x < y ? x : y; }
@@ -146,7 +154,7 @@ static real maxx(real x, real y)
static void swapx(real* x, real* y)
{ real t = *x; *x = *y; *y = t; }
-static void SinCosNorm(real* sinx, real* cosx) {
+static void norm2(real* sinx, real* cosx) {
real r = hypotx(*sinx, *cosx);
*sinx /= r;
*cosx /= r;
@@ -171,7 +179,7 @@ static real AngRound(real x) {
volatile real y = fabs(x);
/* The compiler mustn't "simplify" z - (z - y) to y */
y = y < z ? z - (z - y) : y;
- return x < 0 ? -y : y;
+ return x < 0 ? 0 - y : y;
}
static void A3coeff(struct geod_geodesic* g);
@@ -270,7 +278,8 @@ void geod_lineinit(struct geod_geodesicline* l,
l->f1 = g->f1;
/* If caps is 0 assume the standard direct calculation */
l->caps = (caps ? caps : GEOD_DISTANCE_IN | GEOD_LONGITUDE) |
- GEOD_LATITUDE | GEOD_AZIMUTH; /* Always allow latitude and azimuth */
+ /* always allow latitude and azimuth and unrolling of longitude */
+ GEOD_LATITUDE | GEOD_AZIMUTH | GEOD_LONG_UNROLL;
l->lat1 = lat1;
l->lon1 = lon1;
@@ -286,7 +295,7 @@ void geod_lineinit(struct geod_geodesicline* l,
/* Ensure cbet1 = +epsilon at poles */
sbet1 = l->f1 * sin(phi);
cbet1 = fabs(lat1) == 90 ? tiny : cos(phi);
- SinCosNorm(&sbet1, &cbet1);
+ norm2(&sbet1, &cbet1);
l->dn1 = sqrt(1 + g->ep2 * sq(sbet1));
/* Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0), */
@@ -305,8 +314,8 @@ void geod_lineinit(struct geod_geodesicline* l,
* With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi. */
l->ssig1 = sbet1; l->somg1 = l->salp0 * sbet1;
l->csig1 = l->comg1 = sbet1 != 0 || l->calp1 != 0 ? cbet1 * l->calp1 : 1;
- SinCosNorm(&l->ssig1, &l->csig1); /* sig1 in (-pi, pi] */
- /* SinCosNorm(somg1, comg1); -- don't need to normalize! */
+ norm2(&l->ssig1, &l->csig1); /* sig1 in (-pi, pi] */
+ /* norm2(somg1, comg1); -- don't need to normalize! */
l->k2 = sq(l->calp0) * g->ep2;
eps = l->k2 / (2 * (1 + sqrt(1 + l->k2)) + l->k2);
@@ -452,12 +461,14 @@ real geod_genposition(const struct geod_geodesicline* l,
s12 = flags & GEOD_ARCMODE ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12;
if (outmask & GEOD_LONGITUDE) {
+ int E = l->salp0 < 0 ? -1 : 1; /* east or west going? */
/* tan(omg2) = sin(alp0) * tan(sig2) */
somg2 = l->salp0 * ssig2; comg2 = csig2; /* No need to normalize */
/* omg12 = omg2 - omg1 */
- omg12 = flags & GEOD_LONG_NOWRAP ? sig12
- - (atan2(ssig2, csig2) - atan2(l->ssig1, l->csig1))
- + (atan2(somg2, comg2) - atan2(l->somg1, l->comg1))
+ omg12 = flags & GEOD_LONG_UNROLL
+ ? E * (sig12
+ - (atan2( ssig2, csig2) - atan2( l->ssig1, l->csig1))
+ + (atan2(E * somg2, comg2) - atan2(E * l->somg1, l->comg1)))
: atan2(somg2 * l->comg1 - comg2 * l->somg1,
comg2 * l->comg1 + somg2 * l->somg1);
lam12 = omg12 + l->A3c *
@@ -466,7 +477,7 @@ real geod_genposition(const struct geod_geodesicline* l,
lon12 = lam12 / degree;
/* Use AngNormalize2 because longitude might have wrapped multiple
* times. */
- lon2 = flags & GEOD_LONG_NOWRAP ? l->lon1 + lon12 :
+ lon2 = flags & GEOD_LONG_UNROLL ? l->lon1 + lon12 :
AngNormalize(AngNormalize(l->lon1) + AngNormalize2(lon12));
}
@@ -499,7 +510,7 @@ real geod_genposition(const struct geod_geodesicline* l,
B42 = SinCosSeries(FALSE, ssig2, csig2, l->C4a, nC4);
real salp12, calp12;
if (l->calp0 == 0 || l->salp0 == 0) {
- /* alp12 = alp2 - alp1, used in atan2 so no need to normalized */
+ /* alp12 = alp2 - alp1, used in atan2 so no need to normalize */
salp12 = salp2 * l->calp1 - calp2 * l->salp1;
calp12 = calp2 * l->calp1 + salp2 * l->salp1;
/* The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz
@@ -645,13 +656,13 @@ real geod_geninverse(const struct geod_geodesic* g,
/* Ensure cbet1 = +epsilon at poles */
sbet1 = g->f1 * sin(phi);
cbet1 = lat1 == -90 ? tiny : cos(phi);
- SinCosNorm(&sbet1, &cbet1);
+ norm2(&sbet1, &cbet1);
phi = lat2 * degree;
/* Ensure cbet2 = +epsilon at poles */
sbet2 = g->f1 * sin(phi);
cbet2 = fabs(lat2) == 90 ? tiny : cos(phi);
- SinCosNorm(&sbet2, &cbet2);
+ norm2(&sbet2, &cbet2);
/* If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the
* |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is
@@ -771,7 +782,7 @@ real geod_geninverse(const struct geod_geodesic* g,
for (tripn = FALSE, tripb = FALSE; numit < maxit2; ++numit) {
/* the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
* WGS84 and random input: mean = 2.85, sd = 0.60 */
- real dv,
+ real dv = 0,
v = (Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
&salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2,
&eps, &omg12, numit < maxit1, &dv, C1a, C2a, C3a)
@@ -793,7 +804,7 @@ real geod_geninverse(const struct geod_geodesic* g,
if (nsalp1 > 0 && fabs(dalp1) < pi) {
calp1 = calp1 * cdalp1 - salp1 * sdalp1;
salp1 = nsalp1;
- SinCosNorm(&salp1, &calp1);
+ 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). */
@@ -811,7 +822,7 @@ real geod_geninverse(const struct geod_geodesic* g,
* WGS84 and random input: mean = 4.74, sd = 0.99 */
salp1 = (salp1a + salp1b)/2;
calp1 = (calp1a + calp1b)/2;
- SinCosNorm(&salp1, &calp1);
+ norm2(&salp1, &calp1);
tripn = FALSE;
tripb = (fabs(salp1a - salp1) + (calp1a - calp1) < tolb ||
fabs(salp1 - salp1b) + (calp1 - calp1b) < tolb);
@@ -852,8 +863,8 @@ real geod_geninverse(const struct geod_geodesic* g,
A4 = sq(g->a) * calp0 * salp0 * g->e2;
real C4a[nC4];
real B41, B42;
- SinCosNorm(&ssig1, &csig1);
- SinCosNorm(&ssig2, &csig2);
+ norm2(&ssig1, &csig1);
+ norm2(&ssig2, &csig2);
C4f(g, eps, C4a);
B41 = SinCosSeries(FALSE, ssig1, csig1, C4a, nC4);
B42 = SinCosSeries(FALSE, ssig2, csig2, C4a, nC4);
@@ -1119,7 +1130,7 @@ real InverseStart(const struct geod_geodesic* g,
salp2 = cbet1 * somg12;
calp2 = sbet12 - cbet1 * sbet2 *
(comg12 >= 0 ? sq(somg12) / (1 + comg12) : 1 - comg12);
- SinCosNorm(&salp2, &calp2);
+ norm2(&salp2, &calp2);
/* Set return value */
sig12 = atan2(ssig12, csig12);
} else if (fabs(g->n) > (real)(0.1) || /* No astroid calc if too eccentric */
@@ -1219,7 +1230,7 @@ real InverseStart(const struct geod_geodesic* g,
}
/* Sanity check on starting guess. Backwards check allows NaN through. */
if (!(salp1 <= 0))
- SinCosNorm(&salp1, &calp1);
+ norm2(&salp1, &calp1);
else {
salp1 = 1; calp1 = 0;
}
@@ -1266,8 +1277,8 @@ real Lambda12(const struct geod_geodesic* g,
* tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1) */
ssig1 = sbet1; somg1 = salp0 * sbet1;
csig1 = comg1 = calp1 * cbet1;
- SinCosNorm(&ssig1, &csig1);
- /* SinCosNorm(&somg1, &comg1); -- don't need to normalize! */
+ norm2(&ssig1, &csig1);
+ /* norm2(&somg1, &comg1); -- don't need to normalize! */
/* Enforce symmetries in the case abs(bet2) = -bet1. Need to be careful
* about this case, since this can yield singularities in the Newton
@@ -1288,8 +1299,8 @@ real Lambda12(const struct geod_geodesic* g,
* tan(omg2) = sin(alp0) * tan(sig2). */
ssig2 = sbet2; somg2 = salp0 * sbet2;
csig2 = comg2 = calp2 * cbet2;
- SinCosNorm(&ssig2, &csig2);
- /* SinCosNorm(&somg2, &comg2); -- don't need to normalize! */
+ norm2(&ssig2, &csig2);
+ /* norm2(&somg2, &comg2); -- don't need to normalize! */
/* sig12 = sig2 - sig1, limit to [0, pi] */
sig12 = atan2(maxx(csig1 * ssig2 - ssig1 * csig2, (real)(0)),
@@ -1335,177 +1346,264 @@ real Lambda12(const struct geod_geodesic* g,
}
real A3f(const struct geod_geodesic* g, real eps) {
- /* Evaluate sum(A3x[k] * eps^k, k, 0, nA3x-1) by Horner's method */
- real v = 0;
- int i;
- for (i = nA3x; i; )
- v = eps * v + g->A3x[--i];
- return v;
+ /* Evaluate A3 */
+ return polyval(nA3 - 1, g->A3x, eps);
}
void C3f(const struct geod_geodesic* g, real eps, real c[]) {
- /* Evaluate C3 coeffs by Horner's method
+ /* Evaluate C3 coeffs
* Elements c[1] thru c[nC3 - 1] are set */
- int i, j, k;
real mult = 1;
- for (j = nC3x, k = nC3 - 1; k; ) {
- real t = 0;
- for (i = nC3 - k; i; --i)
- t = eps * t + g->C3x[--j];
- c[k--] = t;
- }
-
- for (k = 1; k < nC3; ) {
+ int o = 0, l;
+ 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[k++] *= mult;
+ c[l] = mult * polyval(m, g->C3x + o, eps);
+ o += m + 1;
}
}
void C4f(const struct geod_geodesic* g, real eps, real c[]) {
- /* Evaluate C4 coeffs by Horner's method
+ /* Evaluate C4 coeffs
* Elements c[0] thru c[nC4 - 1] are set */
- int i, j, k;
real mult = 1;
- for (j = nC4x, k = nC4; k; ) {
- real t = 0;
- for (i = nC4 - k + 1; i; --i)
- t = eps * t + g->C4x[--j];
- c[--k] = t;
- }
-
- for (k = 1; k < nC4; ) {
+ 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);
+ o += m + 1;
mult *= eps;
- c[k++] *= mult;
}
}
-/* Generated by Maxima on 2010-09-04 10:26:17-04:00 */
-
/* The scale factor A1-1 = mean value of (d/dsigma)I1 - 1 */
real A1m1f(real eps) {
- real
- eps2 = sq(eps),
- t = eps2*(eps2*(eps2+4)+64)/256;
+ static const real coeff[] = {
+ /* (1-eps)*A1-1, polynomial in eps2 of order 3 */
+ 1, 4, 64, 0, 256,
+ };
+ int m = nA1/2;
+ real t = polyval(m, coeff, sq(eps)) / coeff[m + 1];
return (t + eps) / (1 - eps);
}
/* The coefficients C1[l] in the Fourier expansion of B1 */
void C1f(real eps, real c[]) {
+ static const real coeff[] = {
+ /* C1[1]/eps^1, polynomial in eps2 of order 2 */
+ -1, 6, -16, 32,
+ /* C1[2]/eps^2, polynomial in eps2 of order 2 */
+ -9, 64, -128, 2048,
+ /* C1[3]/eps^3, polynomial in eps2 of order 1 */
+ 9, -16, 768,
+ /* C1[4]/eps^4, polynomial in eps2 of order 1 */
+ 3, -5, 512,
+ /* C1[5]/eps^5, polynomial in eps2 of order 0 */
+ -7, 1280,
+ /* C1[6]/eps^6, polynomial in eps2 of order 0 */
+ -7, 2048,
+ };
real
eps2 = sq(eps),
d = eps;
- c[1] = d*((6-eps2)*eps2-16)/32;
- d *= eps;
- c[2] = d*((64-9*eps2)*eps2-128)/2048;
- d *= eps;
- c[3] = d*(9*eps2-16)/768;
- d *= eps;
- c[4] = d*(3*eps2-5)/512;
- d *= eps;
- c[5] = -7*d/1280;
- d *= eps;
- c[6] = -7*d/2048;
+ 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];
+ o += m + 2;
+ d *= eps;
+ }
}
/* The coefficients C1p[l] in the Fourier expansion of B1p */
void C1pf(real eps, real c[]) {
+ static const real coeff[] = {
+ /* C1p[1]/eps^1, polynomial in eps2 of order 2 */
+ 205, -432, 768, 1536,
+ /* C1p[2]/eps^2, polynomial in eps2 of order 2 */
+ 4005, -4736, 3840, 12288,
+ /* C1p[3]/eps^3, polynomial in eps2 of order 1 */
+ -225, 116, 384,
+ /* C1p[4]/eps^4, polynomial in eps2 of order 1 */
+ -7173, 2695, 7680,
+ /* C1p[5]/eps^5, polynomial in eps2 of order 0 */
+ 3467, 7680,
+ /* C1p[6]/eps^6, polynomial in eps2 of order 0 */
+ 38081, 61440,
+ };
real
eps2 = sq(eps),
d = eps;
- c[1] = d*(eps2*(205*eps2-432)+768)/1536;
- d *= eps;
- c[2] = d*(eps2*(4005*eps2-4736)+3840)/12288;
- d *= eps;
- c[3] = d*(116-225*eps2)/384;
- d *= eps;
- c[4] = d*(2695-7173*eps2)/7680;
- d *= eps;
- c[5] = 3467*d/7680;
- d *= eps;
- c[6] = 38081*d/61440;
+ 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];
+ o += m + 2;
+ d *= eps;
+ }
}
/* The scale factor A2-1 = mean value of (d/dsigma)I2 - 1 */
real A2m1f(real eps) {
- real
- eps2 = sq(eps),
- t = eps2*(eps2*(25*eps2+36)+64)/256;
+ static const real coeff[] = {
+ /* A2/(1-eps)-1, polynomial in eps2 of order 3 */
+ 25, 36, 64, 0, 256,
+ };
+ int m = nA2/2;
+ real t = polyval(m, coeff, sq(eps)) / coeff[m + 1];
return t * (1 - eps) - eps;
}
/* The coefficients C2[l] in the Fourier expansion of B2 */
void C2f(real eps, real c[]) {
+ static const real coeff[] = {
+ /* C2[1]/eps^1, polynomial in eps2 of order 2 */
+ 1, 2, 16, 32,
+ /* C2[2]/eps^2, polynomial in eps2 of order 2 */
+ 35, 64, 384, 2048,
+ /* C2[3]/eps^3, polynomial in eps2 of order 1 */
+ 15, 80, 768,
+ /* C2[4]/eps^4, polynomial in eps2 of order 1 */
+ 7, 35, 512,
+ /* C2[5]/eps^5, polynomial in eps2 of order 0 */
+ 63, 1280,
+ /* C2[6]/eps^6, polynomial in eps2 of order 0 */
+ 77, 2048,
+ };
real
eps2 = sq(eps),
d = eps;
- c[1] = d*(eps2*(eps2+2)+16)/32;
- d *= eps;
- c[2] = d*(eps2*(35*eps2+64)+384)/2048;
- d *= eps;
- c[3] = d*(15*eps2+80)/768;
- d *= eps;
- c[4] = d*(7*eps2+35)/512;
- d *= eps;
- c[5] = 63*d/1280;
- d *= eps;
- c[6] = 77*d/2048;
+ 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];
+ o += m + 2;
+ d *= eps;
+ }
}
/* The scale factor A3 = mean value of (d/dsigma)I3 */
void A3coeff(struct geod_geodesic* g) {
- g->A3x[0] = 1;
- g->A3x[1] = (g->n-1)/2;
- g->A3x[2] = (g->n*(3*g->n-1)-2)/8;
- g->A3x[3] = ((-g->n-3)*g->n-1)/16;
- g->A3x[4] = (-2*g->n-3)/64;
- g->A3x[5] = -3/(real)(128);
+ static const real coeff[] = {
+ /* A3, coeff of eps^5, polynomial in n of order 0 */
+ -3, 128,
+ /* A3, coeff of eps^4, polynomial in n of order 1 */
+ -2, -3, 64,
+ /* A3, coeff of eps^3, polynomial in n of order 2 */
+ -1, -3, -1, 16,
+ /* A3, coeff of eps^2, polynomial in n of order 2 */
+ 3, -1, -2, 8,
+ /* A3, coeff of eps^1, polynomial in n of order 1 */
+ 1, -1, 2,
+ /* A3, coeff of eps^0, polynomial in n of order 0 */
+ 1, 1,
+ };
+ 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];
+ o += m + 2;
+ }
}
/* The coefficients C3[l] in the Fourier expansion of B3 */
void C3coeff(struct geod_geodesic* g) {
- g->C3x[0] = (1-g->n)/4;
- g->C3x[1] = (1-g->n*g->n)/8;
- g->C3x[2] = ((3-g->n)*g->n+3)/64;
- g->C3x[3] = (2*g->n+5)/128;
- g->C3x[4] = 3/(real)(128);
- g->C3x[5] = ((g->n-3)*g->n+2)/32;
- g->C3x[6] = ((-3*g->n-2)*g->n+3)/64;
- g->C3x[7] = (g->n+3)/128;
- g->C3x[8] = 5/(real)(256);
- g->C3x[9] = (g->n*(5*g->n-9)+5)/192;
- g->C3x[10] = (9-10*g->n)/384;
- g->C3x[11] = 7/(real)(512);
- g->C3x[12] = (7-14*g->n)/512;
- g->C3x[13] = 7/(real)(512);
- g->C3x[14] = 21/(real)(2560);
+ static const real coeff[] = {
+ /* C3[1], coeff of eps^5, polynomial in n of order 0 */
+ 3, 128,
+ /* C3[1], coeff of eps^4, polynomial in n of order 1 */
+ 2, 5, 128,
+ /* C3[1], coeff of eps^3, polynomial in n of order 2 */
+ -1, 3, 3, 64,
+ /* C3[1], coeff of eps^2, polynomial in n of order 2 */
+ -1, 0, 1, 8,
+ /* C3[1], coeff of eps^1, polynomial in n of order 1 */
+ -1, 1, 4,
+ /* C3[2], coeff of eps^5, polynomial in n of order 0 */
+ 5, 256,
+ /* C3[2], coeff of eps^4, polynomial in n of order 1 */
+ 1, 3, 128,
+ /* C3[2], coeff of eps^3, polynomial in n of order 2 */
+ -3, -2, 3, 64,
+ /* C3[2], coeff of eps^2, polynomial in n of order 2 */
+ 1, -3, 2, 32,
+ /* C3[3], coeff of eps^5, polynomial in n of order 0 */
+ 7, 512,
+ /* C3[3], coeff of eps^4, polynomial in n of order 1 */
+ -10, 9, 384,
+ /* C3[3], coeff of eps^3, polynomial in n of order 2 */
+ 5, -9, 5, 192,
+ /* C3[4], coeff of eps^5, polynomial in n of order 0 */
+ 7, 512,
+ /* C3[4], coeff of eps^4, polynomial in n of order 1 */
+ -14, 7, 512,
+ /* C3[5], coeff of eps^5, polynomial in n of order 0 */
+ 21, 2560,
+ };
+ int o = 0, k = 0, l, j;
+ 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];
+ o += m + 2;
+ }
+ }
}
-/* Generated by Maxima on 2012-10-19 08:02:34-04:00 */
-
/* The coefficients C4[l] in the Fourier expansion of I4 */
void C4coeff(struct geod_geodesic* g) {
- g->C4x[0] = (g->n*(g->n*(g->n*(g->n*(100*g->n+208)+572)+3432)-12012)+30030)/
- 45045;
- g->C4x[1] = (g->n*(g->n*(g->n*(64*g->n+624)-4576)+6864)-3003)/15015;
- g->C4x[2] = (g->n*((14144-10656*g->n)*g->n-4576)-858)/45045;
- g->C4x[3] = ((-224*g->n-4784)*g->n+1573)/45045;
- g->C4x[4] = (1088*g->n+156)/45045;
- g->C4x[5] = 97/(real)(15015);
- g->C4x[6] = (g->n*(g->n*((-64*g->n-624)*g->n+4576)-6864)+3003)/135135;
- g->C4x[7] = (g->n*(g->n*(5952*g->n-11648)+9152)-2574)/135135;
- g->C4x[8] = (g->n*(5792*g->n+1040)-1287)/135135;
- g->C4x[9] = (468-2944*g->n)/135135;
- g->C4x[10] = 1/(real)(9009);
- g->C4x[11] = (g->n*((4160-1440*g->n)*g->n-4576)+1716)/225225;
- g->C4x[12] = ((4992-8448*g->n)*g->n-1144)/225225;
- g->C4x[13] = (1856*g->n-936)/225225;
- g->C4x[14] = 8/(real)(10725);
- g->C4x[15] = (g->n*(3584*g->n-3328)+1144)/315315;
- g->C4x[16] = (1024*g->n-208)/105105;
- g->C4x[17] = -136/(real)(63063);
- g->C4x[18] = (832-2560*g->n)/405405;
- g->C4x[19] = -128/(real)(135135);
- g->C4x[20] = 128/(real)(99099);
+ static const real coeff[] = {
+ /* C4[0], coeff of eps^5, polynomial in n of order 0 */
+ 97, 15015,
+ /* C4[0], coeff of eps^4, polynomial in n of order 1 */
+ 1088, 156, 45045,
+ /* C4[0], coeff of eps^3, polynomial in n of order 2 */
+ -224, -4784, 1573, 45045,
+ /* C4[0], coeff of eps^2, polynomial in n of order 3 */
+ -10656, 14144, -4576, -858, 45045,
+ /* C4[0], coeff of eps^1, polynomial in n of order 4 */
+ 64, 624, -4576, 6864, -3003, 15015,
+ /* C4[0], coeff of eps^0, polynomial in n of order 5 */
+ 100, 208, 572, 3432, -12012, 30030, 45045,
+ /* C4[1], coeff of eps^5, polynomial in n of order 0 */
+ 1, 9009,
+ /* C4[1], coeff of eps^4, polynomial in n of order 1 */
+ -2944, 468, 135135,
+ /* C4[1], coeff of eps^3, polynomial in n of order 2 */
+ 5792, 1040, -1287, 135135,
+ /* C4[1], coeff of eps^2, polynomial in n of order 3 */
+ 5952, -11648, 9152, -2574, 135135,
+ /* C4[1], coeff of eps^1, polynomial in n of order 4 */
+ -64, -624, 4576, -6864, 3003, 135135,
+ /* C4[2], coeff of eps^5, polynomial in n of order 0 */
+ 8, 10725,
+ /* C4[2], coeff of eps^4, polynomial in n of order 1 */
+ 1856, -936, 225225,
+ /* C4[2], coeff of eps^3, polynomial in n of order 2 */
+ -8448, 4992, -1144, 225225,
+ /* C4[2], coeff of eps^2, polynomial in n of order 3 */
+ -1440, 4160, -4576, 1716, 225225,
+ /* C4[3], coeff of eps^5, polynomial in n of order 0 */
+ -136, 63063,
+ /* C4[3], coeff of eps^4, polynomial in n of order 1 */
+ 1024, -208, 105105,
+ /* C4[3], coeff of eps^3, polynomial in n of order 2 */
+ 3584, -3328, 1144, 315315,
+ /* C4[4], coeff of eps^5, polynomial in n of order 0 */
+ -128, 135135,
+ /* C4[4], coeff of eps^4, polynomial in n of order 1 */
+ -2560, 832, 405405,
+ /* C4[5], coeff of eps^5, polynomial in n of order 0 */
+ 128, 99099,
+ };
+ int o = 0, k = 0, l, j;
+ 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];
+ o += m + 2;
+ }
+ }
}
int transit(real lon1, real lon2) {
@@ -1594,7 +1692,7 @@ void geod_polygon_addedge(const struct geod_geodesic* g,
real azi, real s) {
if (p->num) { /* Do nothing is num is zero */
real lat, lon, S12;
- geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_NOWRAP, s,
+ geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s,
&lat, &lon, 0,
0, 0, 0, 0, p->polyline ? 0 : &S12);
accadd(p->P, s);
@@ -1731,7 +1829,7 @@ unsigned geod_polygon_testedge(const struct geod_geodesic* g,
crossings = p->crossings;
{
real lat, lon, s12, S12;
- geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_NOWRAP, s,
+ geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s,
&lat, &lon, 0,
0, 0, 0, 0, &S12);
tempsum += S12;
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-16 03:29:53 +0000