From 2dcf018318e9513af7c96dad194780ec66fc42a0 Mon Sep 17 00:00:00 2001
From: Kristian Evers <kristianevers@gmail.com>
Date: Mon, 7 May 2018 11:38:23 +0200
Subject: Refactor calcofi so that global parameters are handled in the
 projection setup

---
 src/PJ_calcofi.c | 37 +++++++++++++++++++------------------
 1 file changed, 19 insertions(+), 18 deletions(-)

(limited to 'src')

diff --git a/src/PJ_calcofi.c b/src/PJ_calcofi.c
index b51b409b..f322d2bd 100644
--- a/src/PJ_calcofi.c
+++ b/src/PJ_calcofi.c
@@ -44,13 +44,12 @@ static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
     double l2;
     double ry; /* r is the point on the same station as o (60) and the same
                line as xy xy, r, o form a right triangle */
-    /* if the user has specified +lon_0 or +k0 for some reason,
-    we're going to ignore it so that xy is consistent with point O */
-    lp.lam = lp.lam + P->lam0;
+
     if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
         proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
         return xy;
     }
+
     xy.x = lp.lam;
     xy.y = -log(pj_tsfn(lp.phi, sin(lp.phi), P->e)); /* Mercator transform xy*/
     oy = -log(pj_tsfn(PT_O_PHI, sin(PT_O_PHI), P->e));
@@ -64,9 +63,7 @@ static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
         (ry - lp.phi) * DEG_TO_STATION / sin(ROTATION_ANGLE);
     /* set a = 1, x0 = 0, and y0 = 0 so that no further unit adjustments
     are done */
-    P->a = 1;
-    P->x0 = 0;
-    P->y0 = 0;
+
     return xy;
 }
 
@@ -77,7 +74,6 @@ static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
     double l1;
     double l2;
     double ry;
-    lp.lam = lp.lam + P->lam0;
     if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
         proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
         return xy;
@@ -93,24 +89,20 @@ static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
         (ry - PT_O_PHI) * DEG_TO_LINE / cos(ROTATION_ANGLE);
     xy.y = PT_O_STATION + RAD_TO_DEG *
         (ry - lp.phi) * DEG_TO_STATION / sin(ROTATION_ANGLE);
-    P->a = 1;
-    P->x0 = 0;
-    P->y0 = 0;
+
     return xy;
 }
 
 
 static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
     LP lp = {0.0,0.0};
-    double ry; /* y value of point r */
+    double ry;     /* y value of point r */
     double oymctr; /* Mercator-transformed y value of point O */
     double rymctr; /* Mercator-transformed ry */
     double xymctr; /* Mercator-transformed xy.y */
     double l1;
     double l2;
-    /* turn x and y back into Line/Station */
-    xy.x /= P->ra;
-    xy.y /= P->ra;
+
     ry = PT_O_PHI - LINE_TO_RAD * (xy.x - PT_O_LINE) *
         cos(ROTATION_ANGLE);
     lp.phi = ry - STATION_TO_RAD * (xy.y - PT_O_STATION) * sin(ROTATION_ANGLE);
@@ -120,7 +112,7 @@ static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
     l1 = (xymctr - oymctr) * tan(ROTATION_ANGLE);
     l2 = (rymctr - xymctr) / (cos(ROTATION_ANGLE) * sin(ROTATION_ANGLE));
     lp.lam = PT_O_LAMBDA - (l1 + l2);
-    P->over = 1;
+
     return lp;
 }
 
@@ -133,8 +125,8 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
     double xymctr;
     double l1;
     double l2;
-    xy.x /= P->ra;
-    xy.y /= P->ra;
+    (void) P;
+
     ry = PT_O_PHI - LINE_TO_RAD * (xy.x - PT_O_LINE) *
         cos(ROTATION_ANGLE);
     lp.phi = ry - STATION_TO_RAD * (xy.y - PT_O_STATION) * sin(ROTATION_ANGLE);
@@ -144,7 +136,7 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
     l1 = (xymctr - oymctr) * tan(ROTATION_ANGLE);
     l2 = (rymctr - xymctr) / (cos(ROTATION_ANGLE) * sin(ROTATION_ANGLE));
     lp.lam = PT_O_LAMBDA - (l1 + l2);
-    P->over = 1;
+
     return lp;
 }
 
@@ -152,6 +144,15 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
 PJ *PROJECTION(calcofi) {
     P->opaque = 0;
 
+    /* if the user has specified +lon_0 or +k0 for some reason,
+    we're going to ignore it so that xy is consistent with point O */
+    P->lam0 = 0;
+    P->ra = 1;
+    P->a = 1;
+    P->x0 = 0;
+    P->y0 = 0;
+    P->over = 1;
+
     if (P->es != 0.0) { /* ellipsoid */
         P->inv = e_inverse;
         P->fwd = e_forward;
-- 
cgit v1.2.3