Support numerical factors only (#664)

* Support numerical factors only
* Make sure h positive. Improve some comments
* Let pole overshoot check have effect even for geocentric latitudes
* Factor-typological constants, now all returning false, for backwards compatibility

7 files changed, 92 insertions, 138 deletions

diff --git a/src/PJ_eqdc.c b/src/PJ_eqdc.c
index e1f8ea16..6846abaa 100644
--- a/src/PJ_eqdc.c
+++ b/src/PJ_eqdc.c
@@ -54,19 +54,6 @@ static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
 }
 
 
-static void special(LP lp, PJ *P, struct FACTORS *fac) {
-    struct pj_opaque *Q = P->opaque;
-    double sinphi, cosphi;
-
-    sinphi = sin(lp.phi);
-    cosphi = cos(lp.phi);
-    fac->code |= IS_ANAL_HK;
-    fac->h = 1.;
-    fac->k = Q->n * (Q->c - (Q->ellips ? pj_mlfn(lp.phi, sinphi,
-        cosphi, Q->en) : lp.phi)) / pj_msfn(sinphi, cosphi, P->es);
-}
-
-
 static void *destructor (PJ *P, int errlev) {                        /* Destructor */
     if (0==P)
         return 0;
@@ -124,7 +111,6 @@ PJ *PROJECTION(eqdc) {
 
     P->inv = e_inverse;
     P->fwd = e_forward;
-    P->spc = special;
 
     return P;
 }
diff --git a/src/PJ_lcc.c b/src/PJ_lcc.c
index bc02667a..a24d5ac0 100644
--- a/src/PJ_lcc.c
+++ b/src/PJ_lcc.c
@@ -73,21 +73,6 @@ static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
     return lp;
 }
 
-static void special(LP lp, PJ *P, struct FACTORS *fac) {
-    struct pj_opaque *Q = P->opaque;
-    double rho;
-    if (fabs(fabs(lp.phi) - M_HALFPI) < EPS10) {
-        if ((lp.phi * Q->n) <= 0.) return;
-        rho = 0.;
-    } else
-        rho = Q->c * (Q->ellips ? pow(pj_tsfn(lp.phi, sin(lp.phi),
-            P->e), Q->n) : pow(tan(M_FORTPI + .5 * lp.phi), -Q->n));
-    fac->code |= IS_ANAL_HK + IS_ANAL_CONV;
-    fac->k = fac->h = P->k0 * Q->n * rho /
-        pj_msfn(sin(lp.phi), cos(lp.phi), P->es);
-    fac->conv = Q->n * lp.lam;
-}
-
 
 PJ *PROJECTION(lcc) {
     double cosphi, sinphi;
@@ -139,7 +124,6 @@ PJ *PROJECTION(lcc) {
 
     P->inv = e_inverse;
     P->fwd = e_forward;
-    P->spc = special;
 
     return P;
 }
diff --git a/src/pj_factors.c b/src/pj_factors.c
index 6c1b4978..69d59e49 100644
--- a/src/pj_factors.c
+++ b/src/pj_factors.c
@@ -1,6 +1,8 @@
 /* projection scale factors */
 #define PJ_LIB__
-#include <projects.h>
+#include <proj.h>
+#include "projects.h"
+
 #include <errno.h>
 
 #ifndef DEFAULT_H
@@ -10,88 +12,85 @@
 #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;
+    int err;
+
+    err = proj_errno_reset (P);
+
+    if (0==fac) {
+        proj_errno_set (P, ENOMEM);
+        return 1;
+    }
 
-    der.x_l = 0.0;
-    der.x_p = 0.0;
-    der.y_l = 0.0;
-    der.y_p = 0.0;
+    /* Indicate that all factors are numerical approximations */
+    fac->code = 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);
+    /* Check for latitude or longitude overange */
+    if ((fabs (lp.phi)-M_HALFPI) > EPS || fabs (lp.lam) > 10.) {
+        proj_errno_set (P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
         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.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));
     }
 
+    /* Set a reasonable step size for the numerical derivatives */
+    h = fabs (h);
+    if (h < EPS)
+        h = DEFAULT_H;
+
+    /* If input latitudes are geocentric, convert to geographic */
+    if (P->geoc)
+        lp.phi = atan(P->rone_es * tan(lp.phi));
+
+    /* If latitude + one step overshoots the pole, move it slightly inside, */
+    /* so the numerical derivative still exists */
+    if (fabs (lp.phi) > (M_HALFPI - h))
+        lp.phi = lp.phi < 0. ? -(M_HALFPI-h) : (M_HALFPI-h);
+
+    /* Longitudinal distance from central meridian */
+    lp.lam -= P->lam0;
+    if (!P->over)
+        lp.lam = adjlon(lp.lam);
+
+    /* Derivatives */
+    if (pj_deriv (lp, h, P, &(fac->der))) {
+        proj_errno_set (P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
+        return 1;
+    }
+
+    /* Scale factors */
+    cosphi = cos (lp.phi);
+    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.;
+
+    /* Convergence */
+    fac->conv = -atan2 (fac->der.x_p, fac->der.y_p);
+
+    /* 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 ellipse axis */
+    t = fac->k * fac->k + fac->h * fac->h;
+    fac->a = sqrt(t + 2. * fac->s);
+    t = t - 2. * fac->s;
+    t = t > 0? sqrt(t): 0;
+    fac->b = 0.5 * (fac->a - t);
+    fac->a = 0.5 * (fac->a + t);
+
+    /* Angular distortion */
+    fac->omega = 2. * aasin(P->ctx, (fac->a - fac->b) / (fac->a + fac->b) );
+
+    proj_errno_restore (P, err);
     return 0;
 }
diff --git a/src/proj.c b/src/proj.c
index 57d33a73..77d6b28b 100644
--- a/src/proj.c
+++ b/src/proj.c
@@ -287,21 +287,15 @@ static void vprocess(FILE *fid) {
         (void)printf(oform, dat_xy.u); putchar('\n');
         (void)fputs("Northing (y):  ", stdout);
         (void)printf(oform, dat_xy.v); putchar('\n');
-        (void)printf("Meridian scale (h)%c: %.8f  ( %.4g %% error )\n",
-            facs.code & IS_ANAL_HK ? '*' : ' ', facs.h, (facs.h-1.)*100.);
-        (void)printf("Parallel scale (k)%c: %.8f  ( %.4g %% error )\n",
-            facs.code & IS_ANAL_HK ? '*' : ' ', facs.k, (facs.k-1.)*100.);
-        (void)printf("Areal scale (s):     %.8f  ( %.4g %% error )\n",
-            facs.s, (facs.s-1.)*100.);
-        (void)printf("Angular distortion (w): %.3f\n", facs.omega *
-            RAD_TO_DEG);
-        (void)printf("Meridian/Parallel angle: %.5f\n",
-            facs.thetap * RAD_TO_DEG);
-        (void)printf("Convergence%c: ",facs.code & IS_ANAL_CONV ? '*' : ' ');
+        (void)printf("Meridian scale (h) : %.8f  ( %.4g %% error )\n", facs.h, (facs.h-1.)*100.);
+        (void)printf("Parallel scale (k) : %.8f  ( %.4g %% error )\n", facs.k, (facs.k-1.)*100.);
+        (void)printf("Areal scale (s):     %.8f  ( %.4g %% error )\n", facs.s, (facs.s-1.)*100.);
+        (void)printf("Angular distortion (w): %.3f\n", facs.omega * RAD_TO_DEG);
+        (void)printf("Meridian/Parallel angle: %.5f\n", facs.thetap * RAD_TO_DEG);
+        (void)printf("Convergence : ");
         (void)fputs(rtodms(pline, facs.conv, 0, 0), stdout);
         (void)printf(" [ %.8f ]\n", facs.conv * RAD_TO_DEG);
-        (void)printf("Max-min (Tissot axis a-b) scale error: %.5f %.5f\n\n",
-            facs.a, facs.b);
+        (void)printf("Max-min (Tissot axis a-b) scale error: %.5f %.5f\n\n", facs.a, facs.b);
     }
 }
 
diff --git a/src/proj.h b/src/proj.h
index 19bb2192..cb689621 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -289,12 +289,6 @@ union PJ_PAIR {
     double v[2]; /* Yes - It's really just a vector! */
 };
 
-
-#define PJ_IS_ANAL_XL_YL 01    /* derivatives of lon analytic */
-#define PJ_IS_ANAL_XP_YP 02    /* derivatives of lat analytic */
-#define PJ_IS_ANAL_HK    04    /* h and k analytic */
-#define PJ_IS_ANAL_CONV 010    /* convergence analytic */
-
 struct PJ_INFO {
     char        release[64];        /* Release info. Version + date         */
     char        version[64];        /* Full version number                  */
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index ee4cb20f..6e89fb60 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -816,9 +816,6 @@ PJ_FACTORS proj_factors(PJ *P, const LP lp) {
 ******************************************************************************/
     PJ_FACTORS factors;
 
-    /* pj_factors rely code being zero */
-    factors.code = 0;
-
     if (pj_factors(lp, P, 0.0, &factors)) {
         /* errno set in pj_factors */
         memset(&factors, 0, sizeof(PJ_FACTORS));
diff --git a/src/projects.h b/src/projects.h
index 1d4fff50..32b74b08 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -257,9 +257,6 @@ struct PJconsts {
     PJ_COORD (*fwd4d)(PJ_COORD, PJ *);
     PJ_COORD (*inv4d)(PJ_COORD, PJ *);
 
-
-    void (*spc)(LP, PJ *, struct FACTORS *);
-
     void *(*destructor)(PJ *, int);
 
 
@@ -460,14 +457,17 @@ struct FACTORS {
     double conv;           /* convergence */
     double s;              /* areal scale factor */
     double a, b;           /* max-min scale error */
-    int code;              /* info as to analytics, see following */
+    int    code;           /* always 0 */
+};
+
+enum deprecated_constants_for_now_dropped_analytical_factors {
+    IS_ANAL_XL_YL =  01,   /* derivatives of lon analytic */
+    IS_ANAL_XP_YP =  02,   /* derivatives of lat analytic */
+    IS_ANAL_HK    =  04,   /* h and k analytic */
+    IS_ANAL_CONV  = 010    /* convergence analytic */
 };
 
 
-#define IS_ANAL_XL_YL 01    /* derivatives of lon analytic */
-#define IS_ANAL_XP_YP 02    /* derivatives of lat analytic */
-#define IS_ANAL_HK    04    /* h and k analytic */
-#define IS_ANAL_CONV 010    /* convergence analytic */
 
 /* datum_type values */
 #define PJD_UNKNOWN   0