From 08c071938d16f99ae7510ebfefbb00b5ee008e96 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Bojan=20=C5=A0avri=C4=8D?= <bsavric@esri.com>
Date: Wed, 22 Aug 2018 14:18:04 -0700
Subject: Style and comment improvements to eqearth (#1102)

---
 src/PJ_eqearth.c | 75 +++++++++++++++++++++++++++-----------------------------
 1 file changed, 36 insertions(+), 39 deletions(-)

(limited to 'src')

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;
 }
-- 
cgit v1.2.3