diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/PJ_eqearth.c | 75 |
1 files changed, 36 insertions, 39 deletions
diff --git a/src/PJ_eqearth.c b/src/PJ_eqearth.c index 2633cbc8..5a2e74fd 100644 --- a/src/PJ_eqearth.c +++ b/src/PJ_eqearth.c @@ -20,6 +20,7 @@ Added ellipsoidal equations by Bojan Savric, 22 August 2018 PROJ_HEAD(eqearth, "Equal Earth") "\n\tPCyl., Sph&Ell"; +/* A1..A4, polynomial coefficients */ #define A1 1.340264 #define A2 -0.081106 #define A3 0.000893 @@ -36,28 +37,27 @@ struct pj_opaque { double *apa; }; -static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal/Spheroidal, forward */ +static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal/spheroidal, forward */ XY xy = {0.0,0.0}; struct pj_opaque *Q = P->opaque; double sbeta; double psi, psi2, psi6; + /* Spheroidal case, using sine latitude */ + sbeta = sin(lp.phi); + + /* In the ellipsoidal case, we convert sbeta to sine of authalic latitude */ if (P->es != 0.0) { - /* Ellipsoidal case, converting to authalic latitude */ - sbeta = pj_qsfn(sin(lp.phi), P->e, 1.0 - P->es) / Q->qp; - if (fabs(sbeta) > 1) { - /* Rounding error. */ + sbeta = pj_qsfn(sbeta, P->e, 1.0 - P->es) / Q->qp; + + /* Rounding error. */ + if (fabs(sbeta) > 1) sbeta = sbeta > 0 ? 1 : -1; - } - } - else { - /* Spheroidal case, using latitude */ - sbeta = sin(lp.phi); } /* Equal Earth projection */ psi = asin(M * sbeta); - psi2 = psi * psi; + psi2 = psi * psi; psi6 = psi2 * psi2 * psi2; xy.x = lp.lam * cos(psi) / (M * (A1 + 3 * A2 * psi2 + psi6 * (7 * A3 + 9 * A4 * psi2))); @@ -71,11 +71,10 @@ static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal/Spheroidal, forwar } -static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal/Spheroidal, inverse */ +static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal/spheroidal, inverse */ LP lp = {0.0,0.0}; struct pj_opaque *Q = P->opaque; - double yc, tol, y2, y6, f, fder; - double beta; + double yc, y2, y6; int i; /* Adjusting x and y for authalic radius */ @@ -83,24 +82,28 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal/Spheroidal, invers xy.y /= Q->rqda; /* Make sure y is inside valid range */ - if (xy.y > MAX_Y) { + if (xy.y > MAX_Y) xy.y = MAX_Y; - } else if (xy.y < -MAX_Y) { + else if (xy.y < -MAX_Y) xy.y = -MAX_Y; - } yc = xy.y; - for (i = MAX_ITER; i ; --i) { /* Newton-Raphson */ + /* Newton-Raphson */ + for (i = MAX_ITER; i ; --i) { + double f, fder, tol; + y2 = yc * yc; y6 = y2 * y2 * y2; - f = yc * (A1 + A2 * y2 + y6 * (A3 + A4 * y2)) - xy.y; + + f = yc * (A1 + A2 * y2 + y6 * (A3 + A4 * y2)) - xy.y; fder = A1 + 3 * A2 * y2 + y6 * (7 * A3 + 9 * A4 * y2); - tol = f / fder; - yc -= tol; - if (fabs(tol) < EPS) { + + tol = f / fder; + yc -= tol; + + if (fabs(tol) < EPS) break; - } } if( i == 0 ) { @@ -114,16 +117,12 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal/Spheroidal, invers lp.lam = M * xy.x * (A1 + 3 * A2 * y2 + y6 * (7 * A3 + 9 * A4 * y2)) / cos(yc); - /* Latitude */ - beta = asin(sin(yc) / M); - if (P->es != 0.0) { - /* Ellipsoidal case, converting authalic latitude */ - lp.phi = pj_authlat(beta, Q->apa); - } - else { - /* Spheroidal case, using latitude */ - lp.phi = beta; - } + /* Latitude (for spheroidal case, this is latitude */ + lp.phi = asin(sin(yc) / M); + + /* Ellipsoidal case, converting auth. latitude */ + if (P->es != 0.0) + lp.phi = pj_authlat(lp.phi, Q->apa); return lp; } @@ -146,7 +145,11 @@ PJ *PROJECTION(eqearth) { return pj_default_destructor (P, ENOMEM); P->opaque = Q; P->destructor = destructor; + P->fwd = e_forward; + P->inv = e_inverse; + Q->rqda = 1.0; + /* Ellipsoidal case */ if (P->es != 0.0) { Q->apa = pj_authset(P->es); /* For auth_lat(). */ if (0 == Q->apa) @@ -154,12 +157,6 @@ PJ *PROJECTION(eqearth) { Q->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */ Q->rqda = sqrt(0.5*Q->qp); /* Authalic radius devided by major axis */ } - else { - Q->rqda = 1.0; - } - - P->fwd = e_forward; - P->inv = e_inverse; return P; } |