Merge pull request #527 from kbevers/issue-526

Correct calculation of meridian convergence for non-conformal projections

2 files changed, 138 insertions, 107 deletions

diff --git a/src/pj_deriv.c b/src/pj_deriv.c
index 2f886476..1c969a9d 100644
--- a/src/pj_deriv.c
+++ b/src/pj_deriv.c
@@ -1,33 +1,58 @@
 /* dervative of (*P->fwd) projection */
 #define PJ_LIB__
 #include "projects.h"
-	int
-pj_deriv(LP lp, double h, PJ *P, struct DERIVS *der) {
-	XY t;
-
-	lp.lam += h;
-	lp.phi += h;
-	if (fabs(lp.phi) > M_HALFPI) return 1;
-	h += h;
-	t = (*P->fwd)(lp, P);
-	if (t.x == HUGE_VAL) return 1;
-	der->x_l = t.x; der->y_p = t.y; der->x_p = -t.x; der->y_l = -t.y;
-	lp.phi -= h;
-	if (fabs(lp.phi) > M_HALFPI) return 1;
-	t = (*P->fwd)(lp, P);
-	if (t.x == HUGE_VAL) return 1;
-	der->x_l += t.x; der->y_p -= t.y; der->x_p += t.x; der->y_l -= t.y;
-	lp.lam -= h;
-	t = (*P->fwd)(lp, P);
-	if (t.x == HUGE_VAL) return 1;
-	der->x_l -= t.x; der->y_p -= t.y; der->x_p += t.x; der->y_l += t.y;
-	lp.phi += h;
-	t = (*P->fwd)(lp, P);
-	if (t.x == HUGE_VAL) return 1;
-	der->x_l -= t.x; der->y_p += t.y; der->x_p -= t.x; der->y_l += t.y;
-	der->x_l /= (h += h);
-	der->y_p /= h;
-	der->x_p /= h;
-	der->y_l /= h;
-	return 0;
+
+int pj_deriv(LP lp, double h, PJ *P, struct DERIVS *der) {
+    XY t;
+
+    lp.lam += h;
+    lp.phi += h;
+    if (fabs(lp.phi) > M_HALFPI)
+        return 1;
+
+    h += h;
+    t = (*P->fwd)(lp, P);
+    if (t.x == HUGE_VAL)
+        return 1;
+
+    der->x_l = t.x;
+    der->y_p = t.y;
+    der->x_p = t.x;
+    der->y_l = t.y;
+    lp.phi -= h;
+    if (fabs(lp.phi) > M_HALFPI)
+        return 1;
+
+    t = (*P->fwd)(lp, P);
+    if (t.x == HUGE_VAL)
+        return 1;
+
+    der->x_l += t.x;
+    der->y_p -= t.y;
+    der->x_p -= t.x;
+    der->y_l += t.y;
+    lp.lam -= h;
+    t = (*P->fwd)(lp, P);
+    if (t.x == HUGE_VAL)
+        return 1;
+
+    der->x_l -= t.x;
+    der->y_p -= t.y;
+    der->x_p -= t.x;
+    der->y_l -= t.y;
+    lp.phi += h;
+    t = (*P->fwd)(lp, P);
+    if (t.x == HUGE_VAL)
+        return 1;
+
+    der->x_l -= t.x;
+    der->y_p += t.y;
+    der->x_p += t.x;
+    der->y_l -= t.y;
+    der->x_l /= (h += h);
+    der->y_p /= h;
+    der->x_p /= h;
+    der->y_l /= h;
+
+    return 0;
 }
diff --git a/src/pj_factors.c b/src/pj_factors.c
index b87b49ff..6c1b4978 100644
--- a/src/pj_factors.c
+++ b/src/pj_factors.c
@@ -2,90 +2,96 @@
 #define PJ_LIB__
 #include <projects.h>
 #include <errno.h>
+
 #ifndef DEFAULT_H
 #define DEFAULT_H   1e-5    /* radian default for numeric h */
 #endif
+
 #define EPS 1.0e-12
-	int
-pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
-	struct DERIVS der;
-	double cosphi, t, n, r;
 
-        der.x_l = 0.0;
-        der.x_p = 0.0;
-        der.y_l = 0.0;
-        der.y_p = 0.0;
+int pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
+    struct DERIVS der;
+    double cosphi, t, n, r;
+
+    der.x_l = 0.0;
+    der.x_p = 0.0;
+    der.y_l = 0.0;
+    der.y_p = 0.0;
 
-	/* check for forward and latitude or longitude overange */
-	if ((t = fabs(lp.phi)-M_HALFPI) > EPS || fabs(lp.lam) > 10.) {
-                pj_ctx_set_errno( P->ctx, -14);
-		return 1;
-	} else { /* proceed */
-		errno = pj_errno = 0;
+    /* check for forward and latitude or longitude overange */
+    if ((t = fabs(lp.phi)-M_HALFPI) > EPS || fabs(lp.lam) > 10.) {
+        pj_ctx_set_errno( P->ctx, -14);
+        return 1;
+    } else { /* proceed */
+        errno = pj_errno = 0;
                 P->ctx->last_errno = 0;
 
-		if (h < EPS)
-			h = DEFAULT_H;
-		if (fabs(lp.phi) > (M_HALFPI - h))
-                /* adjust to value around pi/2 where derived still exists*/
-		        lp.phi = lp.phi < 0. ? (-M_HALFPI+h) : (M_HALFPI-h);
-		else if (P->geoc)
-			lp.phi = atan(P->rone_es * tan(lp.phi));
-		lp.lam -= P->lam0;	/* compute del lp.lam */
-		if (!P->over)
-			lp.lam = adjlon(lp.lam); /* adjust del longitude */
-		if (P->spc)	/* get what projection analytic values */
-			P->spc(lp, P, fac);
-		if (((fac->code & (IS_ANAL_XL_YL+IS_ANAL_XP_YP)) !=
-			  (IS_ANAL_XL_YL+IS_ANAL_XP_YP)) &&
-			  pj_deriv(lp, h, P, &der))
-			return 1;
-		if (!(fac->code & IS_ANAL_XL_YL)) {
-			fac->der.x_l = der.x_l;
-			fac->der.y_l = der.y_l;
-		}
-		if (!(fac->code & IS_ANAL_XP_YP)) {
-			fac->der.x_p = der.x_p;
-			fac->der.y_p = der.y_p;
-		}
-		cosphi = cos(lp.phi);
-		if (!(fac->code & IS_ANAL_HK)) {
-			fac->h = hypot(fac->der.x_p, fac->der.y_p);
-			fac->k = hypot(fac->der.x_l, fac->der.y_l) / cosphi;
-			if (P->es != 0.0) {
-				t = sin(lp.phi);
-				t = 1. - P->es * t * t;
-				n = sqrt(t);
-				fac->h *= t * n / P->one_es;
-				fac->k *= n;
-				r = t * t / P->one_es;
-			} else
-				r = 1.;
-		} else if (P->es != 0.0) {
-			r = sin(lp.phi);
-			r = 1. - P->es * r * r;
-			r = r * r / P->one_es;
-		} else
-			r = 1.;
-		/* convergence */
-		if (!(fac->code & IS_ANAL_CONV)) {
-			fac->conv = - atan2(fac->der.y_l, fac->der.x_l);
-			if (fac->code & IS_ANAL_XL_YL)
-				fac->code |= IS_ANAL_CONV;
-		}
-		/* areal scale factor */
-		fac->s = (fac->der.y_p * fac->der.x_l - fac->der.x_p * fac->der.y_l) *
-			r / cosphi;
-		/* meridian-parallel angle theta prime */
-		fac->thetap = aasin(P->ctx,fac->s / (fac->h * fac->k));
-		/* Tissot ellips axis */
-		t = fac->k * fac->k + fac->h * fac->h;
-		fac->a = sqrt(t + 2. * fac->s);
-		t = (t = t - 2. * fac->s) <= 0. ? 0. : sqrt(t);
-		fac->b = 0.5 * (fac->a - t);
-		fac->a = 0.5 * (fac->a + t);
-		/* omega */
-		fac->omega = 2. * aasin(P->ctx,(fac->a - fac->b)/(fac->a + fac->b));
-	}
-	return 0;
+        if (h < EPS)
+            h = DEFAULT_H;
+        if (fabs(lp.phi) > (M_HALFPI - h))
+            /* adjust to value around pi/2 where derived still exists*/
+            lp.phi = lp.phi < 0. ? (-M_HALFPI+h) : (M_HALFPI-h);
+        else if (P->geoc)
+            lp.phi = atan(P->rone_es * tan(lp.phi));
+        lp.lam -= P->lam0;  /* compute del lp.lam */
+        if (!P->over)
+            lp.lam = adjlon(lp.lam); /* adjust del longitude */
+        if (P->spc) /* get what projection analytic values */
+            P->spc(lp, P, fac);
+        if (((fac->code & (IS_ANAL_XL_YL+IS_ANAL_XP_YP)) !=
+              (IS_ANAL_XL_YL+IS_ANAL_XP_YP)) && pj_deriv(lp, h, P, &der))
+            return 1;
+
+        if (!(fac->code & IS_ANAL_XL_YL)) {
+            fac->der.x_l = der.x_l;
+            fac->der.y_l = der.y_l;
+        }
+        if (!(fac->code & IS_ANAL_XP_YP)) {
+            fac->der.x_p = der.x_p;
+            fac->der.y_p = der.y_p;
+        }
+        cosphi = cos(lp.phi);
+        if (!(fac->code & IS_ANAL_HK)) {
+            fac->h = hypot(fac->der.x_p, fac->der.y_p);
+            fac->k = hypot(fac->der.x_l, fac->der.y_l) / cosphi;
+            if (P->es != 0.0) {
+                t = sin(lp.phi);
+                t = 1. - P->es * t * t;
+                n = sqrt(t);
+                fac->h *= t * n / P->one_es;
+                fac->k *= n;
+                r = t * t / P->one_es;
+            } else
+                r = 1.;
+        } else if (P->es != 0.0) {
+            r = sin(lp.phi);
+            r = 1. - P->es * r * r;
+            r = r * r / P->one_es;
+        } else
+            r = 1.;
+
+        /* convergence */
+        if (!(fac->code & IS_ANAL_CONV)) {
+            fac->conv = - atan2(fac->der.x_p, fac->der.y_p);
+            if (fac->code & IS_ANAL_XL_YL)
+                fac->code |= IS_ANAL_CONV;
+        }
+        /* areal scale factor */
+        fac->s = (fac->der.y_p * fac->der.x_l - fac->der.x_p * fac->der.y_l) * r / cosphi;
+
+        /* meridian-parallel angle theta prime */
+        fac->thetap = aasin(P->ctx,fac->s / (fac->h * fac->k));
+
+        /* Tissot ellips axis */
+        t = fac->k * fac->k + fac->h * fac->h;
+        fac->a = sqrt(t + 2. * fac->s);
+        t = (t = t - 2. * fac->s) <= 0. ? 0. : sqrt(t);
+        fac->b = 0.5 * (fac->a - t);
+        fac->a = 0.5 * (fac->a + t);
+
+        /* omega */
+        fac->omega = 2. * aasin(P->ctx,(fac->a - fac->b)/(fac->a + fac->b));
+    }
+
+    return 0;
 }