Merge remote-tracking branch 'osgeo/master' into doc-improvements

52 files changed, 357 insertions, 236 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index f94e5f26..971ac78c 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -10,7 +10,7 @@ AM_CPPFLAGS =	-DPROJ_LIB=\"$(pkgdatadir)\" \
 		-DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@
 
 include_HEADERS = proj.h proj_api.h projects.h geodesic.h \
-	org_proj4_Projections.h org_proj4_PJ.h
+	org_proj4_PJ.h
 
 EXTRA_DIST = makefile.vc proj.def bin_cct.cmake bin_gie.cmake bin_cs2cs.cmake \
 			 bin_geod.cmake bin_nad2bin.cmake bin_proj.cmake \
@@ -43,7 +43,7 @@ lib_LTLIBRARIES = libproj.la
 libproj_la_LDFLAGS = -no-undefined -version-info 13:0:0
 
 libproj_la_SOURCES = \
-	pj_list.h proj_internal.h\
+	pj_list.h proj_internal.h proj_math.h\
 	PJ_aeqd.c PJ_gnom.c PJ_laea.c PJ_mod_ster.c \
 	PJ_nsper.c PJ_nzmg.c PJ_ortho.c PJ_stere.c PJ_sterea.c \
 	PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c PJ_isea.c PJ_ccon.c\
@@ -83,7 +83,7 @@ libproj_la_SOURCES = \
 	pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \
 	geocent.c geocent.h pj_utils.c pj_gridinfo.c pj_gridlist.c \
 	jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c \
-	pj_strtod.c \
+	pj_strtod.c pj_math.c\
 	\
 	proj_4D_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c \
 	PJ_vgridshift.c PJ_hgridshift.c PJ_unitconvert.c PJ_molodensky.c \
diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index 3fe90524..640f1013 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -31,6 +31,8 @@
 #include "proj.h"
 #include <errno.h>
 #include "projects.h"
+#include "proj_math.h"
+
 
 # define EPS10  1.e-10
 # define TOL7   1.e-7
diff --git a/src/PJ_aeqd.c b/src/PJ_aeqd.c
index 84554b91..5d8c3d38 100644
--- a/src/PJ_aeqd.c
+++ b/src/PJ_aeqd.c
@@ -30,6 +30,7 @@
 #include "proj.h"
 #include <errno.h>
 #include "projects.h"
+#include "proj_math.h"
 
 enum Mode {
     N_POLE = 0,
diff --git a/src/PJ_bipc.c b/src/PJ_bipc.c
index 48fc93d0..2f60808d 100644
--- a/src/PJ_bipc.c
+++ b/src/PJ_bipc.c
@@ -2,6 +2,7 @@
 #include "proj.h"
 #include <errno.h>
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(bipc, "Bipolar conic of western hemisphere") "\n\tConic Sph.";
 
diff --git a/src/PJ_bonne.c b/src/PJ_bonne.c
index 40bb28ee..d3d5b757 100644
--- a/src/PJ_bonne.c
+++ b/src/PJ_bonne.c
@@ -2,6 +2,8 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
+
 
 PROJ_HEAD(bonne, "Bonne (Werner lat_1=90)")
     "\n\tConic Sph&Ell\n\tlat_1=";
diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index ec273f1b..a4fd3254 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -43,8 +43,8 @@
 #define PJ_LIB__
 #include "proj_internal.h"
 #include "projects.h"
+#include "proj_math.h"
 #include <stddef.h>
-#include <math.h>
 #include <errno.h>
 PROJ_HEAD(cart,    "Geodetic/cartesian conversions");
 
diff --git a/src/PJ_ccon.c b/src/PJ_ccon.c
index 02b14388..aed59fb6 100644
--- a/src/PJ_ccon.c
+++ b/src/PJ_ccon.c
@@ -24,6 +24,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 #define EPS10   1e-10
 
diff --git a/src/PJ_deformation.c b/src/PJ_deformation.c
index d6376e40..5511eed4 100644
--- a/src/PJ_deformation.c
+++ b/src/PJ_deformation.c
@@ -55,6 +55,7 @@ grid-values in units of mm/year in ENU-space.
 #include <errno.h>
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 PROJ_HEAD(deformation, "Kinematic grid shift");
@@ -82,6 +83,7 @@ static XYZ get_grid_shift(PJ* P, XYZ cartesian) {
 ********************************************************************************/
     PJ_COORD geodetic, shift, temp;
     double sp, cp, sl, cl;
+    int previous_errno = proj_errno_reset(P);
 
     /* cartesian to geodetic */
     geodetic.lpz = pj_inv3d(cartesian, P->opaque->cart);
@@ -90,6 +92,10 @@ static XYZ get_grid_shift(PJ* P, XYZ cartesian) {
     shift.lp    = proj_hgrid_value(P, geodetic.lp);
     shift.enu.u = proj_vgrid_value(P, geodetic.lp);
 
+    if (proj_errno(P) == PJD_ERR_GRID_AREA)
+        proj_log_debug(P, "deformation: coordinate (%.3f, %.3f) outside deformation model",
+                       proj_todeg(geodetic.lp.lam), proj_todeg(geodetic.lp.phi));
+
     /* grid values are stored as mm/yr, we need m/yr */
     shift.xyz.x /= 1000;
     shift.xyz.y /= 1000;
@@ -108,6 +114,8 @@ static XYZ get_grid_shift(PJ* P, XYZ cartesian) {
 
     shift.xyz = temp.xyz;
 
+    proj_errno_restore(P, previous_errno);
+
     return shift.xyz;
 }
 
diff --git a/src/PJ_eqdc.c b/src/PJ_eqdc.c
index f36e0518..5f65e7d0 100644
--- a/src/PJ_eqdc.c
+++ b/src/PJ_eqdc.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 struct pj_opaque {
     double phi1;
diff --git a/src/PJ_geos.c b/src/PJ_geos.c
index f0a9bc24..2988e62e 100644
--- a/src/PJ_geos.c
+++ b/src/PJ_geos.c
@@ -31,6 +31,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 struct pj_opaque {
     double h;
diff --git a/src/PJ_gnom.c b/src/PJ_gnom.c
index 8bb24b6a..b4bab0e2 100644
--- a/src/PJ_gnom.c
+++ b/src/PJ_gnom.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(gnom, "Gnomonic") "\n\tAzi, Sph.";
 
diff --git a/src/PJ_isea.c b/src/PJ_isea.c
index 57412a38..3e941c77 100644
--- a/src/PJ_isea.c
+++ b/src/PJ_isea.c
@@ -29,8 +29,8 @@ struct hex {
 
 /* y *must* be positive down as the xy /iso conversion assumes this */
 ISEA_STATIC
-int hex_xy(struct hex *h) {
-    if (!h->iso) return 1;
+void hex_xy(struct hex *h) {
+    if (!h->iso) return;
     if (h->x >= 0) {
         h->y = -h->y - (h->x+1)/2;
     } else {
@@ -38,13 +38,11 @@ int hex_xy(struct hex *h) {
         h->y = -h->y - h->x/2;
     }
     h->iso = 0;
-
-    return 1;
 }
 
 ISEA_STATIC
-int hex_iso(struct hex *h) {
-    if (h->iso) return 1;
+void hex_iso(struct hex *h) {
+    if (h->iso) return;
 
     if (h->x >= 0) {
         h->y = (-h->y - (h->x+1)/2);
@@ -55,12 +53,10 @@ int hex_iso(struct hex *h) {
 
     h->z = -h->x - h->y;
     h->iso = 1;
-    return 1;
 }
 
 ISEA_STATIC
-int hexbin2(double width, double x, double y,
-                int *i, int *j) {
+void hexbin2(double width, double x, double y, int *i, int *j) {
     double z, rx, ry, rz;
     double abs_dx, abs_dy, abs_dz;
     int ix, iy, iz, s;
@@ -105,7 +101,6 @@ int hexbin2(double width, double x, double y,
     hex_xy(&h);
     *i = h.x;
     *j = h.y;
-        return ix * 100 + iy;
 }
 #ifndef ISEA_STATIC
 #define ISEA_STATIC
@@ -605,27 +600,27 @@ isea_grid_init(struct isea_dgg * g)
 }
 
 ISEA_STATIC
-int
+void
 isea_orient_isea(struct isea_dgg * g)
 {
     if (!g)
-        return 0;
+        return;
     g->o_lat = ISEA_STD_LAT;
     g->o_lon = ISEA_STD_LON;
     g->o_az = 0.0;
-    return 1;
+    return;
 }
 
 ISEA_STATIC
-int
+void
 isea_orient_pole(struct isea_dgg * g)
 {
     if (!g)
-        return 0;
+        return;
     g->o_lat = M_PI / 2.0;
     g->o_lon = 0.0;
     g->o_az = 0;
-    return 1;
+    return;
 }
 
 ISEA_STATIC
@@ -1134,4 +1129,3 @@ PJ *PROJECTION(isea) {
 
     return P;
 }
-
diff --git a/src/PJ_laea.c b/src/PJ_laea.c
index 0422784b..bcf9c44d 100644
--- a/src/PJ_laea.c
+++ b/src/PJ_laea.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(laea, "Lambert Azimuthal Equal Area") "\n\tAzi, Sph&Ell";
 
diff --git a/src/PJ_lcc.c b/src/PJ_lcc.c
index 863f6be9..34ed99cb 100644
--- a/src/PJ_lcc.c
+++ b/src/PJ_lcc.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(lcc, "Lambert Conformal Conic")
     "\n\tConic, Sph&Ell\n\tlat_1= and lat_2= or lat_0";
diff --git a/src/PJ_merc.c b/src/PJ_merc.c
index 18303c45..9123eae1 100644
--- a/src/PJ_merc.c
+++ b/src/PJ_merc.c
@@ -1,10 +1,21 @@
 #define PJ_LIB__
+#include "proj_internal.h"
 #include "proj.h"
+#include "proj_math.h"
 #include "projects.h"
+#include <float.h>
 
 PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts=";
+PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t";
 
 #define EPS10 1.e-10
+static double logtanpfpim1(double x) {       /* log(tan(x/2 + M_FORTPI)) */
+    if (fabs(x) <= DBL_EPSILON) {
+        /* tan(M_FORTPI + .5 * x) can be approximated by  1.0 + x */
+        return log1p(x);
+    }
+    return log(tan(M_FORTPI + .5 * x));
+}
 
 static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
     XY xy = {0.0,0.0};
@@ -25,7 +36,7 @@ static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
         return xy;
 }
     xy.x = P->k0 * lp.lam;
-    xy.y = P->k0 * log(tan(M_FORTPI + .5 * lp.phi));
+    xy.y = P->k0 * logtanpfpim1(lp.phi);
     return xy;
 }
 
@@ -76,3 +87,18 @@ PJ *PROJECTION(merc) {
     return P;
 }
 
+PJ *PROJECTION(webmerc) {
+
+    /* Overriding k_0, lat_0 and lon_0 with fixed parameters */
+    P->k0 = 1.0;
+    P->phi0 = 0.0;
+    P->lam0 = 0.0;
+
+    P->b = P->a;
+    /* Clean up the ellipsoidal parameters to reflect the sphere */
+    P->es = P->e = P->f = 0;
+    pj_calc_ellipsoid_params (P, P->a, 0);
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+    return P;
+}
diff --git a/src/PJ_mod_ster.c b/src/PJ_mod_ster.c
index d807660c..5e6ce136 100644
--- a/src/PJ_mod_ster.c
+++ b/src/PJ_mod_ster.c
@@ -2,6 +2,7 @@
 #define PJ_LIB__
 #include <errno.h>
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(mil_os, "Miller Oblated Stereographic") "\n\tAzi(mod)";
 PROJ_HEAD(lee_os, "Lee Oblated Stereographic") "\n\tAzi(mod)";
diff --git a/src/PJ_nsper.c b/src/PJ_nsper.c
index 3c0b01f0..223cd75b 100644
--- a/src/PJ_nsper.c
+++ b/src/PJ_nsper.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 enum Mode {
     N_POLE = 0,
diff --git a/src/PJ_oea.c b/src/PJ_oea.c
index 5112896f..2bfeffa8 100644
--- a/src/PJ_oea.c
+++ b/src/PJ_oea.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(oea, "Oblated Equal Area") "\n\tMisc Sph\n\tn= m= theta=";
 
diff --git a/src/PJ_ortho.c b/src/PJ_ortho.c
index 7c27cd56..2b037819 100644
--- a/src/PJ_ortho.c
+++ b/src/PJ_ortho.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph.";
diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
index b60e05a8..aa7d76f8 100644
--- a/src/PJ_pipeline.c
+++ b/src/PJ_pipeline.c
@@ -306,6 +306,8 @@ static void set_ellipsoid(PJ *P) {
         /* the PJ you provided".                                         */
         proj_errno_reset (P);
     }
+    P->a_orig = P->a;
+    P->es_orig = P->es;
 
     pj_calc_ellipsoid_params (P, P->a, P->es);
 
diff --git a/src/PJ_robin.c b/src/PJ_robin.c
index 92aebfd3..a1d5ad4d 100644
--- a/src/PJ_robin.c
+++ b/src/PJ_robin.c
@@ -1,4 +1,5 @@
 #define PJ_LIB__
+#include "proj_math.h"
 #include "proj_internal.h"
 #include "proj.h"
 #include "projects.h"
@@ -82,7 +83,7 @@ static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
     (void) P;
 
     dphi = fabs(lp.phi);
-    i = pj_is_nan(lp.phi) ? -1 : (int)floor(dphi * C1);
+    i = isnan(lp.phi) ? -1 : (int)floor(dphi * C1);
     if( i < 0 ){
         proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
         return xy;
@@ -117,7 +118,7 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
         }
     } else { /* general problem */
         /* in Y space, reduce to table interval */
-        i = pj_is_nan(lp.phi) ? -1 : (int)floor(lp.phi * NODES);
+        i = isnan(lp.phi) ? -1 : (int)floor(lp.phi * NODES);
         if( i < 0 || i >= NODES ) {
             proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
             return lp;
diff --git a/src/PJ_sconics.c b/src/PJ_sconics.c
index 34a0f1fd..ce044c24 100644
--- a/src/PJ_sconics.c
+++ b/src/PJ_sconics.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 
 enum Type {
diff --git a/src/PJ_stere.c b/src/PJ_stere.c
index 62e42460..82fd5c07 100644
--- a/src/PJ_stere.c
+++ b/src/PJ_stere.c
@@ -2,6 +2,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 PROJ_HEAD(stere, "Stereographic") "\n\tAzi, Sph&Ell\n\tlat_ts=";
 PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Sph&Ell\n\tsouth";
diff --git a/src/PJ_sterea.c b/src/PJ_sterea.c
index 96b35609..eb4c9f2c 100644
--- a/src/PJ_sterea.c
+++ b/src/PJ_sterea.c
@@ -26,6 +26,7 @@
 #define PJ_LIB__
 #include <errno.h>
 #include "projects.h"
+#include "proj_math.h"
 
 
 struct pj_opaque {
diff --git a/src/PJ_tpeqd.c b/src/PJ_tpeqd.c
index 21e3c5ed..87877ec1 100644
--- a/src/PJ_tpeqd.c
+++ b/src/PJ_tpeqd.c
@@ -1,6 +1,7 @@
 #define PJ_LIB__
 #include <errno.h>
 #include "proj.h"
+#include "proj_math.h"
 #include "projects.h"
 
 
diff --git a/src/cct.c b/src/cct.c
index 6b8f35f3..956f91b3 100644
--- a/src/cct.c
+++ b/src/cct.c
@@ -76,6 +76,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-10-26
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <stdarg.h>
 
 #include "proj.h"
 #include "proj_internal.h"
@@ -86,13 +87,14 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-10-26
 /* Prototypes for functions in proj_strtod.c */
 double proj_strtod(const char *str, char **endptr);
 double proj_atof(const char *str);
+static void logger(void *data, int level, const char *msg);
+static void print(PJ_LOG_LEVEL verbosity, const char *fmt, ...);
 
 /* Prototypes from functions in this file */
 char *column (char *buf, int n);
 PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double fixed_time);
 
 
-
 static const char usage[] = {
     "--------------------------------------------------------------------------------\n"
     "Usage: %s [-options]... [+operator_specs]... infile...\n"
@@ -105,6 +107,7 @@ static const char usage[] = {
     "    -z value          Provide a fixed z value for all input data (e.g. -z 0)\n"
     "    -t value          Provide a fixed t value for all input data (e.g. -t 0)\n"
     "    -I                Do the inverse transformation\n"
+    "    -s n              Skip n first lines of a infile\n"
     "    -v                Verbose: Provide non-essential informational output.\n"
     "                      Repeat -v for more verbosity (e.g. -vv)\n"
     "--------------------------------------------------------------------------------\n"
@@ -116,6 +119,7 @@ static const char usage[] = {
     "    --time            Alias for -t\n"
     "    --verbose         Alias for -v\n"
     "    --inverse         Alias for -I\n"
+    "    --skip-lines      Alias for -s\n"
     "    --help            Alias for -h\n"
     "    --version         Print version number\n"
     "--------------------------------------------------------------------------------\n"
@@ -149,20 +153,73 @@ static const char usage[] = {
     "--------------------------------------------------------------------------------\n"
 };
 
+
+static void logger(void *data, int level, const char *msg) {
+    FILE *stream;
+    int log_tell = proj_log_level(PJ_DEFAULT_CTX, PJ_LOG_TELL);
+
+    stream  = (FILE *) data;
+
+    /* if we use PJ_LOG_NONE we always want to print stuff to stream */
+    if (level == PJ_LOG_NONE) {
+        fprintf(stream, "%s", msg);
+        return;
+    }
+
+    /* should always print to stderr if level == PJ_LOG_ERROR */
+    if (level == PJ_LOG_ERROR) {
+        fprintf(stderr, "%s", msg);
+        return;
+    }
+
+    /* otherwise only print if log level set by user is high enough */
+    if (level <= log_tell)
+        fprintf(stream, "%s", msg);
+}
+
+FILE *fout;
+
+static void print(PJ_LOG_LEVEL log_level, const char *fmt, ...) {
+
+    va_list args;
+    char *msg_buf;
+
+    va_start( args, fmt );
+
+    msg_buf = (char *) malloc(100000);
+    if( msg_buf == NULL )
+        return;
+
+    vsprintf( msg_buf, fmt, args );
+
+    logger((void *) fout, log_level, msg_buf);
+
+    va_end( args );
+    free( msg_buf );
+}
+
+
 int main(int argc, char **argv) {
     PJ *P;
     PJ_COORD point;
     PJ_PROJ_INFO info;
     OPTARGS *o;
-    FILE *fout = stdout;
     char *buf;
-    int nfields = 4, direction = 1, verbose;
+    int nfields = 4, direction = 1, skip_lines = 0, verbose;
     double fixed_z = HUGE_VAL, fixed_time = HUGE_VAL;
     int columns_xyzt[] = {1, 2, 3, 4};
     const char *longflags[]  = {"v=verbose", "h=help", "I=inverse", "version", 0};
-    const char *longkeys[]   = {"o=output",  "c=columns", "z=height", "t=time", 0};
+    const char *longkeys[]   = {
+        "o=output",
+        "c=columns",
+        "z=height",
+        "t=time",
+        "s=skip-lines",
+        0};
+
+    fout = stdout;
 
-    o = opt_parse (argc, argv, "hvI", "cozt", longflags, longkeys);
+    o = opt_parse (argc, argv, "hvI", "cozts", longflags, longkeys);
     if (0==o)
         return 0;
 
@@ -171,26 +228,26 @@ int main(int argc, char **argv) {
         return 0;
     }
 
-
     direction = opt_given (o, "I")? -1: 1;
-    verbose   = opt_given (o, "v");
+
+    verbose   = MIN(opt_given (o, "v"), 3); /* log level can't be larger than 3 */
+    proj_log_level (PJ_DEFAULT_CTX, verbose);
+    proj_log_func  (PJ_DEFAULT_CTX, (void *) fout, logger);
 
     if (opt_given (o, "version")) {
-        fprintf (stdout, "%s: %s\n", o->progname, pj_get_release ());
+        print (PJ_LOG_NONE, "%s: %s\n", o->progname, pj_get_release ());
         return 0;
     }
 
     if (opt_given (o, "o"))
         fout = fopen (opt_arg (o, "output"), "wt");
     if (0==fout) {
-        fprintf (stderr, "%s: Cannot open '%s' for output\n", o->progname, opt_arg (o, "output"));
+        print (PJ_LOG_ERROR, "%s: Cannot open '%s' for output\n", o->progname, opt_arg (o, "output"));
         free (o);
         return 1;
     }
-    if (verbose > 3)
-        fprintf (fout, "%s: Running in very verbose mode\n", o->progname);
-
 
+    print (PJ_LOG_TRACE, "%s: Running in very verbose mode\n", o->progname);
 
     if (opt_given (o, "z")) {
         fixed_z = proj_atof (opt_arg (o, "z"));
@@ -202,11 +259,15 @@ int main(int argc, char **argv) {
         nfields--;
     }
 
+    if (opt_given (o, "s")) {
+        skip_lines = atoi (opt_arg(o, "s"));
+    }
+
     if (opt_given (o, "c")) {
         /* cppcheck-suppress invalidscanf */
         int ncols = sscanf (opt_arg (o, "c"), "%d,%d,%d,%d", columns_xyzt, columns_xyzt+1, columns_xyzt+2, columns_xyzt+3);
         if (ncols != nfields) {
-            fprintf (stderr, "%s: Too few input columns given: '%s'\n", o->progname, opt_arg (o, "c"));
+            print (PJ_LOG_ERROR, "%s: Too few input columns given: '%s'\n", o->progname, opt_arg (o, "c"));
             free (o);
             if (stdout != fout)
                 fclose (fout);
@@ -217,7 +278,7 @@ int main(int argc, char **argv) {
     /* Setup transformation */
     P = proj_create_argv (0, o->pargc, o->pargv);
     if ((0==P) || (0==o->pargc)) {
-        fprintf (stderr, "%s: Bad transformation arguments - (%s)\n    '%s -h' for help\n",
+        print (PJ_LOG_ERROR, "%s: Bad transformation arguments - (%s)\n    '%s -h' for help\n",
                  o->progname, pj_strerrno (proj_errno(P)), o->progname);
         free (o);
         if (stdout != fout)
@@ -225,15 +286,13 @@ int main(int argc, char **argv) {
         return 1;
     }
 
-    if (verbose > 4) {
-        info = proj_pj_info (P);
-        fprintf (stdout, "Final: %s argc=%d pargc=%d\n", info.definition, argc, o->pargc);
-    }
+    info = proj_pj_info (P);
+    print (PJ_LOG_TRACE, "Final: %s argc=%d pargc=%d\n", info.definition, argc, o->pargc);
 
     if (direction==-1) {
         /* fail if an inverse operation is not available */
-        if (!proj_pj_info(P).has_inverse) {
-            fprintf (stderr, "Inverse operation not available\n");
+        if (!info.has_inverse) {
+            print (PJ_LOG_ERROR, "Inverse operation not available\n");
             if (stdout != fout)
                 fclose (fout);
             return 1;
@@ -246,7 +305,7 @@ int main(int argc, char **argv) {
     /* Allocate input buffer */
     buf = calloc (1, 10000);
     if (0==buf) {
-        fprintf (stderr, "%s: Out of memory\n", o->progname);
+        print (PJ_LOG_ERROR, "%s: Out of memory\n", o->progname);
         pj_free (P);
         free (o);
         if (stdout != fout)
@@ -259,25 +318,28 @@ int main(int argc, char **argv) {
     while (opt_input_loop (o, optargs_file_format_text)) {
         int err;
         void *ret = fgets (buf, 10000, o->input);
+        char *c = column (buf, 1);
         opt_eof_handler (o);
         if (0==ret) {
-            fprintf (stderr, "Read error in record %d\n", (int) o->record_index);
+            print (PJ_LOG_ERROR, "Read error in record %d\n", (int) o->record_index);
             continue;
         }
         point = parse_input_line (buf, columns_xyzt, fixed_z, fixed_time);
-        if (HUGE_VAL==point.xyzt.x) {
-            char *c = column (buf, 1);
+        if (skip_lines > 0) {
+            skip_lines--;
+            continue;
+        }
 
-            /* if it's a comment or blank line, we reflect it */
-            if (c && ((*c=='\0') || (*c=='#'))) {
-                fprintf (fout, "%s\n", buf);
-                continue;
-            }
+        /* if it's a comment or blank line, we reflect it */
+        if (c && ((*c=='\0') || (*c=='#'))) {
+            fprintf (fout, "%s", buf);
+            continue;
+        }
 
+        if (HUGE_VAL==point.xyzt.x) {
             /* otherwise, it must be a syntax error */
-            fprintf (fout, "# Record %d UNREADABLE: %s", (int) o->record_index, buf);
-            if (verbose)
-                fprintf (stderr, "%s: Could not parse file '%s' line %d\n", o->progname, opt_filename (o), opt_record (o));
+            print (PJ_LOG_NONE, "# Record %d UNREADABLE: %s", (int) o->record_index, buf);
+            print (PJ_LOG_ERROR, "%s: Could not parse file '%s' line %d\n", o->progname, opt_filename (o), opt_record (o));
             continue;
         }
 
@@ -290,8 +352,8 @@ int main(int argc, char **argv) {
 
         if (HUGE_VAL==point.xyzt.x) {
             /* transformation error */
-            fprintf (fout, "# Record %d TRANSFORMATION ERROR: %s (%s)",
-                            (int) o->record_index, buf, pj_strerrno (proj_errno(P)));
+            print (PJ_LOG_NONE, "# Record %d TRANSFORMATION ERROR: %s (%s)",
+                   (int) o->record_index, buf, pj_strerrno (proj_errno(P)));
             proj_errno_restore (P, err);
             continue;
         }
@@ -301,10 +363,10 @@ int main(int argc, char **argv) {
         if (proj_angular_output (P, direction)) {
             point.lpzt.lam = proj_todeg (point.lpzt.lam);
             point.lpzt.phi = proj_todeg (point.lpzt.phi);
-            fprintf (fout, "%14.10f  %14.10f  %12.4f  %12.4f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t);
+            print (PJ_LOG_NONE, "%14.10f  %14.10f  %12.4f  %12.4f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t);
         }
         else
-            fprintf (fout, "%13.4f  %13.4f  %12.4f  %12.4f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t);
+            print (PJ_LOG_NONE, "%13.4f  %13.4f  %12.4f  %12.4f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t);
     }
 
     if (stdout != fout)
diff --git a/src/geodesic.c b/src/geodesic.c
index eb956c65..9904c7fa 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -24,7 +24,11 @@
  */
 
 #include "geodesic.h"
+#ifdef PJ_LIB__
+#include "proj_math.h"
+#else
 #include <math.h>
+#endif
 
 #if !defined(HAVE_C99_MATH)
 #define HAVE_C99_MATH 0
@@ -239,8 +243,7 @@ static void sincosdx(real x, real* sinx, real* cosx) {
   r = remquo(x, (real)(90), &q);
 #else
   r = fmod(x, (real)(360));
-  /* check for NaN -- do not use pj_is_nan, since we want geodesic.c not to
-   * depend on the rest of proj.4 */
+  /* check for NaN */
   q = r == r ? (int)(floor(r / 90 + (real)(0.5))) : 0;
   r -= 90 * q;
 #endif
diff --git a/src/gie.c b/src/gie.c
index 8c735b88..b81e6bdc 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -114,6 +114,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-10-01/2017-10-08
 
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 #include "optargpm.h"
@@ -964,6 +965,7 @@ Indicate that the remaining material should be skipped. Mostly for debugging.
 ******************************************************************************/
     T.skip = 1;
     (void) args;
+    F->level = 2; /* Silence complaints about missing </gie> element */
     return 0;
 }
 
diff --git a/src/hypot.c b/src/hypot.c
deleted file mode 100644
index 822c4595..00000000
--- a/src/hypot.c
+++ /dev/null
@@ -1,36 +0,0 @@
-/* hypot - sqrt(x * x + y * y)
-**
-**	Because this was omitted from the ANSI standards, this version
-** 	is included for those systems that do not include hypot as an
-**	extension to libm.a.  Note: GNU version was not used because it
-**	was not properly coded to minimize potential overflow.
-**
-**	The proper technique for determining hypot is to factor out the
-**	larger of the two terms,  thus leaving a possible case of float
-**	overflow when max(x,y)*sqrt(2) > max machine value.  This allows
-**	a wider range of numbers than the alternative of the sum of the
-**	squares < max machine value.  For an Intel x87 IEEE double of
-**	approximately 1.8e308, only argument values > 1.27e308 are at
-**	risk of causing overflow.  Whereas, not using this method limits
-**	the range to values less that 9.5e153 --- a considerable reduction
-**	in range!
-*/
-extern double sqrt(double);
-	double
-hypot(double x, double y) {
-	if ( x < 0.)
-		x = -x;
-	else if (x == 0.)
-		return (y < 0. ? -y : y);
-	if (y < 0.)
-		y = -y;
-	else if (y == 0.)
-		return (x);
-	if ( x < y ) {
-		x /= y;
-		return ( y * sqrt( 1. + x * x ) );
-	} else {
-		y /= x;
-		return ( x * sqrt( 1. + y * y ) );
-	}
-}
diff --git a/src/jniproj.c b/src/jniproj.c
index 85d984da..072fa04c 100644
--- a/src/jniproj.c
+++ b/src/jniproj.c
@@ -5,7 +5,9 @@
  *           Martin Desruisseaux
  *
  ******************************************************************************
- * Copyright (c) 2005, Antonello Andrea
+ * Copyright (c) 2005, Andrea Antonello
+ * Copyright (c) 2011, Martin Desruisseaux
+ * Copyright (c) 2018, Even Rouault
  *
  * Permission is hereby granted, free of charge, to any person obtaining a
  * copy of this software and associated documentation files (the "Software"),
@@ -351,11 +353,8 @@ void convertAngularOrdinates(PJ *pj, double* data, jint numPts, int dimension, d
     if (pj_is_latlong(pj)) {
         /* Convert only the 2 first ordinates and skip all the other dimensions. */
         dimToSkip = dimension - 2;
-    } else if (pj_is_geocent(pj)) {
-        /* Convert only the 3 first ordinates and skip all the other dimensions. */
-        dimToSkip = dimension - 3;
     } else {
-        /* Not a geographic or geocentric CRS: nothing to convert. */
+        /* Not a geographic CRS: nothing to convert. */
         return;
     }
     double *stop = data + dimension*numPts;
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index c9e4d9e6..82a8fc81 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -199,6 +199,7 @@ SET(SRC_LIBPROJ_CORE
         pj_list.h
         pj_log.c
         pj_malloc.c
+        pj_math.c
         pj_mlfn.c
         pj_msfn.c
         pj_mutex.c
@@ -218,6 +219,7 @@ SET(SRC_LIBPROJ_CORE
         pj_utils.c
         pj_zpoly1.c
         proj_mdist.c
+        proj_math.h
         proj_rouss.c
         rtodms.c
         vector1.c
@@ -256,8 +258,7 @@ if(JNI_SUPPORT)
   set(SRC_LIBPROJ_CORE ${SRC_LIBPROJ_CORE}
                        jniproj.c )
   set(HEADERS_LIBPROJ ${HEADERS_LIBPROJ}
-                        org_proj4_PJ.h
-                        org_proj4_Projections.h)
+                        org_proj4_PJ.h)
   source_group("Source Files\\JNI" FILES ${SRC_LIBPROJ_JNI})
   add_definitions(-DJNI_ENABLED)
   include_directories( ${JNI_INCLUDE_DIRS})
diff --git a/src/makefile.vc b/src/makefile.vc
index ee89beb3..ef39b084 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -55,7 +55,7 @@ support = \
 	pj_utils.obj pj_gridlist.obj pj_gridinfo.obj \
 	proj_mdist.obj pj_mutex.obj pj_initcache.obj \
 	pj_ctx.obj pj_fileapi.obj pj_log.obj pj_apply_vgridshift.obj \
-	pj_strtod.obj pj_internal.obj
+	pj_strtod.obj pj_internal.obj pj_math.obj
 
 pipeline = \
 	proj_4D_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \
diff --git a/src/nad_cvt.c b/src/nad_cvt.c
index c913511f..d4330d54 100644
--- a/src/nad_cvt.c
+++ b/src/nad_cvt.c
@@ -1,5 +1,6 @@
 #define PJ_LIB__
 #include "projects.h"
+#include "proj_math.h"
 #define MAX_ITERATIONS 10
 #define TOL 1e-12
 
diff --git a/src/nad_intr.c b/src/nad_intr.c
index dc245831..80e428ed 100644
--- a/src/nad_intr.c
+++ b/src/nad_intr.c
@@ -1,6 +1,7 @@
 /* Determine nad table correction value */
 #define PJ_LIB__
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 	LP
@@ -13,9 +14,9 @@ nad_intr(LP t, struct CTABLE *ct) {
 	int in;
 
 	t.lam /= ct->del.lam;
-	indx.lam = pj_is_nan(t.lam) ? 0 : (int)floor(t.lam);
+	indx.lam = isnan(t.lam) ? 0 : (int)floor(t.lam);
 	t.phi /= ct->del.phi;
-	indx.phi = pj_is_nan(t.phi) ? 0 : (int)floor(t.phi);
+	indx.phi = isnan(t.phi) ? 0 : (int)floor(t.phi);
 
 	frct.lam = t.lam - indx.lam;
 	frct.phi = t.phi - indx.phi;
diff --git a/src/org_proj4_PJ.h b/src/org_proj4_PJ.h
index 7e690e7a..be5d3f58 100644
--- a/src/org_proj4_PJ.h
+++ b/src/org_proj4_PJ.h
@@ -44,7 +44,7 @@ JNIEXPORT jstring JNICALL Java_org_proj4_PJ_getDefinition
 /*
  * Class:     org_proj4_PJ
  * Method:    getType
- * Signature: ()Lorg/proj4/PJ$Type;
+ * Signature: ()Lorg/proj4/PJ/Type;
  */
 JNIEXPORT jobject JNICALL Java_org_proj4_PJ_getType
   (JNIEnv *, jobject);
diff --git a/src/org_proj4_Projections.h b/src/org_proj4_Projections.h
deleted file mode 100644
index 3841e057..00000000
--- a/src/org_proj4_Projections.h
+++ /dev/null
@@ -1,37 +0,0 @@
-/* DO NOT EDIT THIS FILE - it is machine generated */
-#include <jni.h>
-/* Header for class org_proj4_Projections */
-
-#ifndef _Included_org_proj4_Projections
-#define _Included_org_proj4_Projections
-#ifdef __cplusplus
-extern "C" {
-#endif
-/*
- * Class:     org_proj4_Projections
- * Method:    getProjInfo
- * Signature: (Ljava/lang/String;)Ljava/lang/String;
- */
-JNIEXPORT jstring JNICALL Java_org_proj4_Projections_getProjInfo
-  (JNIEnv *, jobject, jstring);
-
-/*
- * Class:     org_proj4_Projections
- * Method:    getEllipsInfo
- * Signature: (Ljava/lang/String;)Ljava/lang/String;
- */
-JNIEXPORT jstring JNICALL Java_org_proj4_Projections_getEllipsInfo
-  (JNIEnv *, jobject, jstring);
-
-/*
- * Class:     org_proj4_Projections
- * Method:    transform
- * Signature: ([D[D[DLjava/lang/String;Ljava/lang/String;JI)V
- */
-JNIEXPORT void JNICALL Java_org_proj4_Projections_transform
-  (JNIEnv *, jobject, jdoubleArray, jdoubleArray, jdoubleArray, jstring, jstring, jlong, jint);
-
-#ifdef __cplusplus
-}
-#endif
-#endif
diff --git a/src/pj_apply_gridshift.c b/src/pj_apply_gridshift.c
index c503ec0b..1ef39b20 100644
--- a/src/pj_apply_gridshift.c
+++ b/src/pj_apply_gridshift.c
@@ -306,13 +306,18 @@ int proj_hgrid_init(PJ* P, const char *grids) {
 /********************************************/
 LP proj_hgrid_value(PJ *P, LP lp) {
     struct CTABLE *ct;
-    LP out;
+    LP out = proj_coord_error().lp;
 
     ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist);
+    if (ct == 0) {
+        pj_ctx_set_errno( P->ctx, PJD_ERR_GRID_AREA);
+        return out;
+    }
 
     /* normalize input to ll origin */
     lp.lam -= ct->ll.lam;
     lp.phi -= ct->ll.phi;
+
     lp.lam = adjlon(lp.lam - M_PI) + M_PI;
 
     out = nad_intr(lp, ct);
diff --git a/src/pj_apply_vgridshift.c b/src/pj_apply_vgridshift.c
index 05803bc8..e7106b25 100644
--- a/src/pj_apply_vgridshift.c
+++ b/src/pj_apply_vgridshift.c
@@ -29,7 +29,7 @@
 #define PJ_LIB__
 
 #include <string.h>
-#include <math.h>
+#include "proj_math.h"
 #include "proj_internal.h"
 #include "projects.h"
 
@@ -42,7 +42,7 @@ static double read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GR
     float  *cvs;
     /* do not deal with NaN coordinates */
     /* cppcheck-suppress duplicateExpression */
-    if( pj_is_nan(input.phi) || pj_is_nan(input.lam) )
+    if( isnan(input.phi) || isnan(input.lam) )
         itable = *gridlist_count_p;
 
     /* keep trying till we find a table that works */
diff --git a/src/pj_factors.c b/src/pj_factors.c
index e682a12d..373bd743 100644
--- a/src/pj_factors.c
+++ b/src/pj_factors.c
@@ -1,6 +1,7 @@
 /* projection scale factors */
 #define PJ_LIB__
 #include "proj.h"
+#include "proj_math.h"
 #include "projects.h"
 
 #include <errno.h>
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index 2a064e58..347f8334 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -112,9 +112,6 @@ static PJ_COORD fwd_finalize (PJ *P, PJ_COORD coo) {
             coo = proj_trans (P->cart, PJ_FWD, coo);
         }
 
-        coo.xyz.x = P->fr_meter * (coo.xyz.x + P->x0);
-        coo.xyz.y = P->fr_meter * (coo.xyz.y + P->y0);
-        coo.xyz.z = P->fr_meter * (coo.xyz.z + P->z0);
         break;
 
     /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
diff --git a/src/pj_init.c b/src/pj_init.c
index 9b06ff28..bdaf64f7 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -37,6 +37,7 @@
 #include <ctype.h>
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 
diff --git a/src/pj_internal.c b/src/pj_internal.c
index 4da47051..dc528649 100644
--- a/src/pj_internal.c
+++ b/src/pj_internal.c
@@ -364,11 +364,6 @@ to that context.
     return;
 }
 
-
-
-
-
-
 /*  logging  */
 
 /* pj_vlog resides in pj_log.c and relates to pj_log as vsprintf relates to sprintf */
@@ -376,11 +371,11 @@ void pj_vlog( projCtx ctx, int level, const char *fmt, va_list args );
 
 
 /***************************************************************************************/
-enum proj_log_level proj_log_level (PJ_CONTEXT *ctx, enum proj_log_level log_level) {
+PJ_LOG_LEVEL proj_log_level (PJ_CONTEXT *ctx, PJ_LOG_LEVEL log_level) {
 /****************************************************************************************
    Set logging level 0-3. Higher number means more debug info. 0 turns it off
 ****************************************************************************************/
-    enum proj_log_level previous;
+    PJ_LOG_LEVEL previous;
     if (0==ctx)
         ctx = pj_get_default_ctx();
     if (0==ctx)
@@ -443,19 +438,3 @@ void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION logf) {
     if (0!=logf)
         ctx->logger = logf;
 }
-
-
-#if HAVE_C99_MATH
-/* proj_internal.h defines pj_is_nan as isnan */
-#else
-/*****************************************************************************/
-int pj_is_nan (double val) {
-/******************************************************************************
-    Returns 0 if not a NaN and non-zero if val is a NaN.
-
-    Provides an equivalent to isnan().
-******************************************************************************/
-    /* cppcheck-suppress duplicateExpression */
-    return val != val;
-}
-#endif
diff --git a/src/pj_inv.c b/src/pj_inv.c
index f66fd7d1..327fef9f 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -98,10 +98,6 @@ static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) {
 
     /* de-scale and de-offset */
     case PJ_IO_UNITS_CARTESIAN:
-        coo.xyz.x = P->to_meter * coo.xyz.x - P->x0;
-        coo.xyz.y = P->to_meter * coo.xyz.y - P->y0;
-        coo.xyz.z = P->to_meter * coo.xyz.z - P->z0;
-
         if (P->is_geocent)
             coo = proj_trans (P->cart, PJ_INV, coo);
 
diff --git a/src/pj_list.h b/src/pj_list.h
index 09a66034..440f7972 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -156,6 +156,7 @@ PROJ_HEAD(wag4, "Wagner IV")
 PROJ_HEAD(wag5, "Wagner V")
 PROJ_HEAD(wag6, "Wagner VI")
 PROJ_HEAD(wag7, "Wagner VII")
+PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator")
 PROJ_HEAD(weren, "Werenskiold I")
 PROJ_HEAD(wink1, "Winkel I")
 PROJ_HEAD(wink2, "Winkel II")
diff --git a/src/pj_math.c b/src/pj_math.c
new file mode 100644
index 00000000..751c2750
--- /dev/null
+++ b/src/pj_math.c
@@ -0,0 +1,71 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Make C99 math functions available on C89 systems
+ * Author:   Kristian Evers
+ *
+ ******************************************************************************
+ * Copyright (c) 2018, Kristian Evers
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#include "proj_math.h"
+
+#if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH)
+
+/* Compute hypotenuse */
+double pj_hypot(double x, double y) {
+    x = fabs(x);
+    y = fabs(y);
+    if ( x < y ) {
+        x /= y;
+        return ( y * sqrt( 1. + x * x ) );
+    } else {
+        y /= (x != 0.0 ? x : 1.0);
+        return ( x * sqrt( 1. + y * y ) );
+    }
+}
+
+/* Compute log(1+x) accurately */
+double pj_log1p(double x) {
+    volatile double
+      y = 1 + x,
+      z = y - 1;
+    /* Here's the explanation for this magic: y = 1 + z, exactly, and z
+     * approx x, thus log(y)/z (which is nearly constant near z = 0) returns
+     * a good approximation to the true log(1 + x)/x.  The multiplication x *
+     * (log(y)/z) introduces little additional error. */
+    return z == 0 ? x : x * log(y) / z;
+}
+
+/* Compute asinh(x) accurately */
+double pj_asinh(double x) {
+    double y = fabs(x); /* Enforce odd parity */
+    y = log1p(y * (1 + y/(hypot(1.0, y) + 1)));
+    return x > 0 ? y : (x < 0 ? -y : x);
+}
+
+/* Returns 0 if not a NaN and non-zero if val is a NaN */
+int pj_isnan (double x) {
+    /* cppcheck-suppress duplicateExpression */
+    return x != x;
+}
+
+
+#endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */
diff --git a/src/pj_transform.c b/src/pj_transform.c
index 22b685b2..aaa8c1e4 100644
--- a/src/pj_transform.c
+++ b/src/pj_transform.c
@@ -42,6 +42,8 @@ enum PJ_DIRECTION {
 };
 typedef enum PJ_DIRECTION PJ_DIRECTION;
 
+/* Copied from proj.h FIXME */
+int  proj_errno_reset (const PJ *P);
 
 
 static int adjust_axis( projCtx ctx, const char *axis, int denormalize_flag,
@@ -199,6 +201,7 @@ static int geographic_to_projected (PJ *P, long n, int dist, double *x, double *
             if (geodetic_loc.u == HUGE_VAL)
                 continue;
 
+            proj_errno_reset( P );
             projected_loc = pj_fwd3d( geodetic_loc, P);
             if( P->ctx->last_errno != 0 )
             {
@@ -234,6 +237,7 @@ static int geographic_to_projected (PJ *P, long n, int dist, double *x, double *
         if( geodetic_loc.u == HUGE_VAL )
             continue;
 
+        proj_errno_reset( P );
         projected_loc = pj_fwd( geodetic_loc, P );
         if( P->ctx->last_errno != 0 )
         {
@@ -301,6 +305,7 @@ static int projected_to_geographic (PJ *P, long n, int dist, double *x, double *
             if (projected_loc.u == HUGE_VAL)
                 continue;
 
+            proj_errno_reset( P );
             geodetic_loc = pj_inv3d(projected_loc, P);
             if( P->ctx->last_errno != 0 )
             {
@@ -337,6 +342,7 @@ static int projected_to_geographic (PJ *P, long n, int dist, double *x, double *
         if( projected_loc.u == HUGE_VAL )
             continue;
 
+        proj_errno_reset( P );
         geodetic_loc = pj_inv( projected_loc, P );
         if( P->ctx->last_errno != 0 )
         {
diff --git a/src/proj.h b/src/proj.h
index 7bc9b10e..d254eade 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -268,6 +268,17 @@ struct PJ_INIT_INFO {
     char        lastupdate[16];     /* Date of last update in YYYY-MM-DD format */
 };
 
+typedef enum PJ_LOG_LEVEL {
+    PJ_LOG_NONE  = 0,
+    PJ_LOG_ERROR = 1,
+    PJ_LOG_DEBUG = 2,
+    PJ_LOG_TRACE = 3,
+    PJ_LOG_TELL  = 4,
+    PJ_LOG_DEBUG_MAJOR = 2, /* for proj_api.h compatibility */
+    PJ_LOG_DEBUG_MINOR = 3  /* for proj_api.h compatibility */
+} PJ_LOG_LEVEL;
+
+typedef void (*PJ_LOG_FUNCTION)(void *, int, const char *);
 
 
 /* The context type - properly namespaced synonym for projCtx */
@@ -342,7 +353,6 @@ double proj_xyz_dist (PJ_COORD a, PJ_COORD b);
 PJ_COORD proj_geod (const PJ *P, PJ_COORD a, PJ_COORD b);
 
 
-
 /* Set or read error level */
 int  proj_context_errno (PJ_CONTEXT *ctx);
 int  proj_errno (const PJ *P);
@@ -351,6 +361,9 @@ int  proj_errno_reset (const PJ *P);
 int  proj_errno_restore (const PJ *P, int err);
 const char* proj_errno_string (int err);
 
+PJ_LOG_LEVEL proj_log_level (PJ_CONTEXT *ctx, PJ_LOG_LEVEL log_level);
+void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION logf);
+
 /* Scaling and angular distortion factors */
 PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp);
 
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 11a56ac0..70110746 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -30,6 +30,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 #include "geodesic.h"
 
@@ -485,7 +486,7 @@ Returns 1 on success, 0 on failure
         if (0==d[0] && 0==d[1] && 0==d[2] && 0==d[3] && 0==d[4] && 0==d[5] && 0==d[6]) {
             /* If the current ellipsoid is not WGS84, then make sure the */
             /* change in ellipsoid is still done. */
-            if (!(fabs(P->a - 6378137.0) < 1e-8 && fabs(P->f - 1./ 298.257223563) < 1e-15)) {
+            if (!(fabs(P->a_orig - 6378137.0) < 1e-8 && fabs(P->es_orig - 0.0066943799901413) < 1e-15)) {
                 do_cart = 1;
             }
             break;
@@ -512,7 +513,7 @@ Returns 1 on success, 0 on failure
     /* geocentric/cartesian space or we need to do a Helmert transform.         */
     if (P->is_geocent || P->helmert || do_cart) {
         char def[150];
-        sprintf (def, "break_cs2cs_recursion     proj=cart   a=%40.20g  f=%40.20g", P->a, P->f);
+        sprintf (def, "break_cs2cs_recursion     proj=cart   a=%40.20g  es=%40.20g", P->a_orig, P->es_orig);
         Q = proj_create (P->ctx, def);
         if (0==Q)
             return 0;
diff --git a/src/proj_etmerc.c b/src/proj_etmerc.c
index 99646d14..9a89e2b0 100644
--- a/src/proj_etmerc.c
+++ b/src/proj_etmerc.c
@@ -44,6 +44,7 @@
 #include <errno.h>
 #include "proj.h"
 #include "projects.h"
+#include "proj_math.h"
 
 
 struct pj_opaque {
@@ -62,32 +63,6 @@ PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)")
 
 #define PROJ_ETMERC_ORDER 6
 
-
-#ifdef _GNU_SOURCE
-    inline
-#endif
-static double log1py(double x) {              /* Compute log(1+x) accurately */
-    volatile double
-      y = 1 + x,
-      z = y - 1;
-    /* Here's the explanation for this magic: y = 1 + z, exactly, and z
-     * approx x, thus log(y)/z (which is nearly constant near z = 0) returns
-     * a good approximation to the true log(1 + x)/x.  The multiplication x *
-     * (log(y)/z) introduces little additional error. */
-    return z == 0 ? x : x * log(y) / z;
-}
-
-
-#ifdef _GNU_SOURCE
-    inline
-#endif
-static double asinhy(double x) {              /* Compute asinh(x) accurately */
-    double y = fabs(x);         /* Enforce odd parity */
-    y = log1py(y * (1 + y/(hypot(1.0, y) + 1)));
-    return x < 0 ? -y : y;
-}
-
-
 #ifdef _GNU_SOURCE
     inline
 #endif
@@ -182,7 +157,7 @@ static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
     Ce     = atan2 (sin_Ce*cos_Cn,  hypot (sin_Cn, cos_Cn*cos_Ce));
 
     /* compl. sph. N, E -> ell. norm. N, E */
-    Ce  = asinhy ( tan (Ce) );     /* Replaces: Ce  = log(tan(FORTPI + Ce*0.5)); */
+    Ce  = asinh ( tan (Ce) );     /* Replaces: Ce  = log(tan(FORTPI + Ce*0.5)); */
     Cn += clenS (Q->gtu, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
     Ce += dCe;
     if (fabs (Ce) <= 2.623395162778) {
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 3f6ccde0..30a4a89e 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -97,25 +97,9 @@ double          proj_vgrid_value(PJ *P, PJ_LP lp);
 PJ_LP           proj_hgrid_value(PJ *P, PJ_LP lp);
 PJ_LP           proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction);
 
-/* High level functionality for handling thread contexts */
-enum proj_log_level {
-    PJ_LOG_NONE  = 0,
-    PJ_LOG_ERROR = 1,
-    PJ_LOG_DEBUG = 2,
-    PJ_LOG_TRACE = 3,
-    PJ_LOG_TELL  = 4,
-    PJ_LOG_DEBUG_MAJOR = 2, /* for proj_api.h compatibility */
-    PJ_LOG_DEBUG_MINOR = 3  /* for proj_api.h compatibility */
-};
-
-/* Set logging level 0-3. Higher number means more debug info. 0 turns it off */
-enum proj_log_level proj_log_level (PJ_CONTEXT *ctx, enum proj_log_level log_level);
-typedef void (*PJ_LOG_FUNCTION)(void *, int, const char *);
-
 void proj_log_error (PJ *P, const char *fmt, ...);
 void proj_log_debug (PJ *P, const char *fmt, ...);
 void proj_log_trace (PJ *P, const char *fmt, ...);
-void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION logf);
 
 int pj_ellipsoid (PJ *);
 void pj_inherit_ellipsoid_def (const PJ *src, PJ *dst);
@@ -134,12 +118,6 @@ void proj_fileapi_set (PJ *P, void *fileapi);
 const char * const *proj_get_searchpath(void);
 int    proj_get_path_count(void);
 
-#if HAVE_C99_MATH
-#define pj_is_nan isnan
-#else
-int pj_is_nan (double val);
-#endif
-
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/proj_math.h b/src/proj_math.h
new file mode 100644
index 00000000..0108a439
--- /dev/null
+++ b/src/proj_math.h
@@ -0,0 +1,55 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Make C99 math functions available on C89 systems
+ * Author:   Kristian Evers
+ *
+ ******************************************************************************
+ * Copyright (c) 2018, Kristian Evers
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#include <math.h>
+
+#ifndef PROJ_MATH_H
+#define PROJ_MATH_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH)
+
+
+double  pj_hypot(double x, double y);
+double  pj_log1p(double x);
+double  pj_asinh(double x);
+int     pj_isnan(double x);
+
+#define hypot   pj_hypot
+#define log1p   pj_log1p
+#define asinh   pj_asinh
+#define isnan   pj_isnan
+
+
+#endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */
+
+#ifdef __cplusplus
+}
+#endif
+#endif /*PROJ_MATH_H */
diff --git a/src/projects.h b/src/projects.h
index 5c8d1301..d6d6e23e 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -92,11 +92,6 @@ typedef long pj_int32;
 #define MAX_PATH_FILENAME 1024
 #endif
 
-/* prototype hypot for systems where absent */
-#if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH)
-extern double hypot(double, double);
-#endif
-
 /* If we still haven't got M_PI*, we rely on our own defines.
  * For example, this is necessary when compiling with gcc and
  * the -ansi flag.