Merge branch 'master' into return-nans-quickly

143 files changed, 749 insertions, 517 deletions

diff --git a/src/PJ_aitoff.c b/src/PJ_aitoff.c
index 89035f3a..82c981d5 100644
--- a/src/PJ_aitoff.c
+++ b/src/PJ_aitoff.c
@@ -29,8 +29,11 @@
  *****************************************************************************/
 
 #define PJ_LIB__
-#include "proj.h"
+
 #include <errno.h>
+#include <math.h>
+
+#include "proj.h"
 #include "projects.h"
 
 
@@ -192,5 +195,3 @@ PJ *PROJECTION(wintri) {
         Q->cosphi1 = 0.636619772367581343;
     return setup(P);
 }
-
-
diff --git a/src/PJ_august.c b/src/PJ_august.c
index ba9ea5cd..4c197114 100644
--- a/src/PJ_august.c
+++ b/src/PJ_august.c
@@ -1,6 +1,8 @@
 #define PJ_LIB__
-#include	"projects.h"
 
+#include <math.h>
+
+#include "projects.h"
 
 PROJ_HEAD(august, "August Epicycloidal") "\n\tMisc Sph, no inv.";
 #define M 1.333333333333333
@@ -30,4 +32,3 @@ PJ *PROJECTION(august) {
     P->es = 0.;
     return P;
 }
-
diff --git a/src/PJ_axisswap.c b/src/PJ_axisswap.c
index 0b81a733..d31e927e 100644
--- a/src/PJ_axisswap.c
+++ b/src/PJ_axisswap.c
@@ -54,6 +54,9 @@ operation:
 
 #define PJ_LIB__
 #include <errno.h>
+#include <stdlib.h>
+#include <string.h>
+
 #include "proj_internal.h"
 #include "projects.h"
 
diff --git a/src/PJ_bacon.c b/src/PJ_bacon.c
index 4ba70db4..ae25bce1 100644
--- a/src/PJ_bacon.c
+++ b/src/PJ_bacon.c
@@ -1,8 +1,10 @@
 # define HLFPI2	2.46740110027233965467      /* (pi/2)^2 */
 # define EPS	1e-10
 #define PJ_LIB__
-#include	<errno.h>
-#include	"projects.h"
+#include <errno.h>
+#include <math.h>
+
+#include "projects.h"
 
 
 struct pj_opaque {
@@ -75,5 +77,3 @@ PJ *PROJECTION(ortel) {
     P->fwd = s_forward;
     return P;
 }
-
-
diff --git a/src/PJ_bipc.c b/src/PJ_bipc.c
index 2f60808d..ea7883a5 100644
--- a/src/PJ_bipc.c
+++ b/src/PJ_bipc.c
@@ -1,6 +1,8 @@
 #define PJ_LIB__
-#include "proj.h"
 #include <errno.h>
+#include <math.h>
+
+#include "proj.h"
 #include "projects.h"
 #include "proj_math.h"
 
@@ -170,5 +172,3 @@ PJ *PROJECTION(bipc) {
     P->es = 0.;
     return P;
 }
-
-
diff --git a/src/PJ_boggs.c b/src/PJ_boggs.c
index ed753505..83bafc18 100644
--- a/src/PJ_boggs.c
+++ b/src/PJ_boggs.c
@@ -1,6 +1,8 @@
 #define PJ_LIB__
-# include	<errno.h>
-# include	"projects.h"
+#include <math.h>
+
+#include "projects.h"
+
 PROJ_HEAD(boggs, "Boggs Eumorphic") "\n\tPCyl., no inv., Sph.";
 # define NITER	20
 # define EPS	1e-7
@@ -39,5 +41,3 @@ PJ *PROJECTION(boggs) {
     P->fwd = s_forward;
     return P;
 }
-
-
diff --git a/src/PJ_calcofi.c b/src/PJ_calcofi.c
index b51b409b..ed4cfe86 100644
--- a/src/PJ_calcofi.c
+++ b/src/PJ_calcofi.c
@@ -1,15 +1,14 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
+#include "proj_api.h"
 
 PROJ_HEAD(calcofi,
     "Cal Coop Ocean Fish Invest Lines/Stations") "\n\tCyl, Sph&Ell";
 
-#include <string.h>
-#include <stdio.h>
-#include <math.h>
-#include "proj_api.h"
-#include <errno.h>
 
 /* Conversions for the California Cooperative Oceanic Fisheries Investigations
 Line/Station coordinate system following the algorithm of:
@@ -44,13 +43,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 +62,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 +73,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 +88,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 +111,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 +124,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 +135,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 +143,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;
@@ -161,5 +161,3 @@ PJ *PROJECTION(calcofi) {
     }
     return P;
 }
-
-
diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index a4fd3254..b4fc60b2 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -41,11 +41,11 @@
  *****************************************************************************/
 
 #define PJ_LIB__
+
 #include "proj_internal.h"
 #include "projects.h"
 #include "proj_math.h"
-#include <stddef.h>
-#include <errno.h>
+
 PROJ_HEAD(cart,    "Geodetic/cartesian conversions");
 
 
diff --git a/src/PJ_cass.c b/src/PJ_cass.c
index acf779a1..4e3b4251 100644
--- a/src/PJ_cass.c
+++ b/src/PJ_cass.c
@@ -1,6 +1,10 @@
 #define PJ_LIB__
-# include   <errno.h>
-# include   "projects.h"
+
+#include <errno.h>
+#include <math.h>
+
+#include "projects.h"
+
 PROJ_HEAD(cass, "Cassini") "\n\tCyl, Sph&Ell";
 
 
@@ -115,5 +119,3 @@ PJ *PROJECTION(cass) {
 
     return P;
 }
-
-
diff --git a/src/PJ_cc.c b/src/PJ_cc.c
index 4f41e3e1..152e6e4a 100644
--- a/src/PJ_cc.c
+++ b/src/PJ_cc.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -36,5 +39,3 @@ PJ *PROJECTION(cc) {
 
     return P;
 }
-
-
diff --git a/src/PJ_cea.c b/src/PJ_cea.c
index e7a33899..e05e524b 100644
--- a/src/PJ_cea.c
+++ b/src/PJ_cea.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -96,5 +99,3 @@ PJ *PROJECTION(cea) {
 
     return P;
 }
-
-
diff --git a/src/PJ_chamb.c b/src/PJ_chamb.c
index 571bdf10..0d14b95e 100644
--- a/src/PJ_chamb.c
+++ b/src/PJ_chamb.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -134,5 +137,3 @@ PJ *PROJECTION(chamb) {
 
     return P;
 }
-
-
diff --git a/src/PJ_collg.c b/src/PJ_collg.c
index 3d7f0d80..ef712905 100644
--- a/src/PJ_collg.c
+++ b/src/PJ_collg.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -48,4 +51,3 @@ PJ *PROJECTION(collg) {
 
     return P;
 }
-
diff --git a/src/PJ_comill.c b/src/PJ_comill.c
index ba0744c6..ce5bb392 100644
--- a/src/PJ_comill.c
+++ b/src/PJ_comill.c
@@ -7,7 +7,10 @@ Port to PROJ.4 by Bojan Savric, 4 April 2016
 */
 
 #define PJ_LIB__
-#include    "projects.h"
+
+#include <math.h>
+
+#include "projects.h"
 
 PROJ_HEAD(comill, "Compact Miller") "\n\tCyl., Sph.";
 
@@ -79,5 +82,3 @@ PJ *PROJECTION(comill) {
 
     return P;
 }
-
-
diff --git a/src/PJ_crast.c b/src/PJ_crast.c
index dec22faf..6db617d0 100644
--- a/src/PJ_crast.c
+++ b/src/PJ_crast.c
@@ -1,5 +1,7 @@
 #define PJ_LIB__
-# include   "projects.h"
+#include <math.h>
+
+#include "projects.h"
 
 PROJ_HEAD(crast, "Craster Parabolic (Putnins P4)") "\n\tPCyl., Sph.";
 
@@ -36,4 +38,3 @@ PJ *PROJECTION(crast) {
 
     return P;
 }
-
diff --git a/src/PJ_denoy.c b/src/PJ_denoy.c
index a5fc2fe7..64de1f9f 100644
--- a/src/PJ_denoy.c
+++ b/src/PJ_denoy.c
@@ -1,5 +1,7 @@
 #define PJ_LIB__
-#include    "projects.h"
+#include <math.h>
+
+#include "projects.h"
 
 PROJ_HEAD(denoy, "Denoyer Semi-Elliptical") "\n\tPCyl., no inv., Sph.";
 
@@ -28,4 +30,3 @@ PJ *PROJECTION(denoy) {
 
     return P;
 }
-
diff --git a/src/PJ_eck1.c b/src/PJ_eck1.c
index 4a9ac06b..069aa185 100644
--- a/src/PJ_eck1.c
+++ b/src/PJ_eck1.c
@@ -1,4 +1,6 @@
 #define PJ_LIB__
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(eck1, "Eckert I") "\n\tPCyl., Sph.";
@@ -37,5 +39,3 @@ PJ *PROJECTION(eck1) {
     return P
 ;
 }
-
-
diff --git a/src/PJ_eck2.c b/src/PJ_eck2.c
index b87f3d97..e9ef1229 100644
--- a/src/PJ_eck2.c
+++ b/src/PJ_eck2.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -51,5 +54,3 @@ PJ *PROJECTION(eck2) {
 
     return P;
 }
-
-
diff --git a/src/PJ_eck3.c b/src/PJ_eck3.c
index 21353be2..674f6919 100644
--- a/src/PJ_eck3.c
+++ b/src/PJ_eck3.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(eck3, "Eckert III") "\n\tPCyl, Sph.";
@@ -105,4 +108,3 @@ PJ *PROJECTION(putp1) {
 
     return setup(P);
 }
-
diff --git a/src/PJ_eck4.c b/src/PJ_eck4.c
index 7cd24720..367d9cc9 100644
--- a/src/PJ_eck4.c
+++ b/src/PJ_eck4.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
-#include    "projects.h"
+
+#include <math.h>
+
+#include "projects.h"
 
 PROJ_HEAD(eck4, "Eckert IV") "\n\tPCyl, Sph.";
 
@@ -58,4 +61,3 @@ PJ *PROJECTION(eck4) {
 
     return P;
 }
-
diff --git a/src/PJ_eck5.c b/src/PJ_eck5.c
index 34a258ee..275430d8 100644
--- a/src/PJ_eck5.c
+++ b/src/PJ_eck5.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(eck5, "Eckert V") "\n\tPCyl, Sph.";
@@ -35,4 +38,3 @@ PJ *PROJECTION(eck5) {
 
     return P;
 }
-
diff --git a/src/PJ_eqc.c b/src/PJ_eqc.c
index 9c80f82f..86845aa3 100644
--- a/src/PJ_eqc.c
+++ b/src/PJ_eqc.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -47,5 +50,3 @@ PJ *PROJECTION(eqc) {
 
     return P;
 }
-
-
diff --git a/src/PJ_eqdc.c b/src/PJ_eqdc.c
index 5f65e7d0..8caca87a 100644
--- a/src/PJ_eqdc.c
+++ b/src/PJ_eqdc.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 #include "proj_math.h"
@@ -115,5 +118,3 @@ PJ *PROJECTION(eqdc) {
 
     return P;
 }
-
-
diff --git a/src/PJ_fahey.c b/src/PJ_fahey.c
index 90ed691f..18be5ea0 100644
--- a/src/PJ_fahey.c
+++ b/src/PJ_fahey.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(fahey, "Fahey") "\n\tPcyl, Sph.";
@@ -36,4 +39,3 @@ PJ *PROJECTION(fahey) {
 
     return P;
 }
-
diff --git a/src/PJ_fouc_s.c b/src/PJ_fouc_s.c
index 87526265..10bd1765 100644
--- a/src/PJ_fouc_s.c
+++ b/src/PJ_fouc_s.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -65,5 +68,3 @@ PJ *PROJECTION(fouc_s) {
     P->fwd = s_forward;
     return P;
 }
-
-
diff --git a/src/PJ_gall.c b/src/PJ_gall.c
index ca8a1bcb..a8697482 100644
--- a/src/PJ_gall.c
+++ b/src/PJ_gall.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(gall, "Gall (Gall Stereographic)") "\n\tCyl, Sph";
@@ -39,5 +42,3 @@ PJ *PROJECTION(gall) {
 
     return P;
 }
-
-
diff --git a/src/PJ_geoc.c b/src/PJ_geoc.c
index e7a7abf1..bed0db49 100644
--- a/src/PJ_geoc.c
+++ b/src/PJ_geoc.c
@@ -27,8 +27,10 @@
  *****************************************************************************/
 
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "proj.h"
-#include <errno.h>
 #include "projects.h"
 
 PROJ_HEAD(geoc, "Geocentric Latitude");
diff --git a/src/PJ_geos.c b/src/PJ_geos.c
index 2988e62e..0eb25610 100644
--- a/src/PJ_geos.c
+++ b/src/PJ_geos.c
@@ -29,6 +29,9 @@
 
 #define PJ_LIB__
 #include <errno.h>
+#include <math.h>
+#include <stddef.h>
+
 #include "proj.h"
 #include "projects.h"
 #include "proj_math.h"
@@ -200,9 +203,6 @@ PJ *PROJECTION(geos) {
     if ((Q->h = pj_param(P->ctx, P->params, "dh").f) <= 0.)
         return pj_default_destructor (P, PJD_ERR_H_LESS_THAN_ZERO);
 
-    if (P->phi0 != 0.0)
-        return pj_default_destructor (P, PJD_ERR_UNKNOWN_PRIME_MERIDIAN);
-
     sweep_axis = pj_param(P->ctx, P->params, "ssweep").s;
     if (sweep_axis == NULL)
       Q->flip_axis = 0;
@@ -234,5 +234,3 @@ PJ *PROJECTION(geos) {
 
     return P;
 }
-
-
diff --git a/src/PJ_gn_sinu.c b/src/PJ_gn_sinu.c
index 36f438e7..21cba308 100644
--- a/src/PJ_gn_sinu.c
+++ b/src/PJ_gn_sinu.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -182,4 +185,3 @@ PJ *PROJECTION(gn_sinu) {
 
     return P;
 }
-
diff --git a/src/PJ_gnom.c b/src/PJ_gnom.c
index b4bab0e2..9c2c7106 100644
--- a/src/PJ_gnom.c
+++ b/src/PJ_gnom.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 #include "proj_math.h"
@@ -138,4 +141,3 @@ PJ *PROJECTION(gnom) {
 
     return P;
 }
-
diff --git a/src/PJ_goode.c b/src/PJ_goode.c
index 65761c4b..4f299b54 100644
--- a/src/PJ_goode.c
+++ b/src/PJ_goode.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -77,5 +80,3 @@ PJ *PROJECTION(goode) {
 
     return P;
 }
-
-
diff --git a/src/PJ_gstmerc.c b/src/PJ_gstmerc.c
index 81ec4e17..5b5109b0 100644
--- a/src/PJ_gstmerc.c
+++ b/src/PJ_gstmerc.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(gstmerc, "Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reunion)")
@@ -67,4 +70,3 @@ PJ *PROJECTION(gstmerc) {
 
     return P;
 }
-
diff --git a/src/PJ_hammer.c b/src/PJ_hammer.c
index e7546ac0..f3e0d64e 100644
--- a/src/PJ_hammer.c
+++ b/src/PJ_hammer.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -70,4 +73,3 @@ PJ *PROJECTION(hammer) {
 
     return P;
 }
-
diff --git a/src/PJ_hatano.c b/src/PJ_hatano.c
index 985bf0cc..4e852888 100644
--- a/src/PJ_hatano.c
+++ b/src/PJ_hatano.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -78,4 +81,3 @@ PJ *PROJECTION(hatano) {
 
     return P;
 }
-
diff --git a/src/PJ_healpix.c b/src/PJ_healpix.c
index c3eae14f..d76ca33d 100644
--- a/src/PJ_healpix.c
+++ b/src/PJ_healpix.c
@@ -28,16 +28,18 @@
  * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  * SOFTWARE.
  *****************************************************************************/
-# define PJ_LIB__
-# include <errno.h>
-# include "proj_internal.h"
-# include "proj.h"
-# include "projects.h"
+#define PJ_LIB__
+
+#include <errno.h>
+#include <math.h>
+
+#include "proj_internal.h"
+#include "proj.h"
+#include "projects.h"
 
 PROJ_HEAD(healpix, "HEALPix") "\n\tSph., Ellps.";
 PROJ_HEAD(rhealpix, "rHEALPix") "\n\tSph., Ellps.\n\tnorth_square= south_square=";
 
-# include <stdio.h>
 /* Matrix for counterclockwise rotation by pi/2: */
 # define R1 {{ 0,-1},{ 1, 0}}
 /* Matrix for counterclockwise rotation by pi: */
@@ -668,4 +670,3 @@ PJ *PROJECTION(rhealpix) {
 
     return P;
 }
-
diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c
index 5ff8c5d2..8ba955bb 100644
--- a/src/PJ_helmert.c
+++ b/src/PJ_helmert.c
@@ -43,11 +43,14 @@ Last update: 2017-05-15
 ***********************************************************************/
 
 #define PJ_LIB__
-#include <stddef.h>
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj_internal.h"
 #include "projects.h"
 #include "geocent.h"
+
 PROJ_HEAD(helmert, "3(6)-, 4(8)- and 7(14)-parameter Helmert shift");
 
 static XYZ helmert_forward_3d (LPZ lpz, PJ *P);
@@ -145,7 +148,6 @@ static void update_parameters(PJ *P) {
         proj_log_trace(P, "rz: %g", Q->opk.k);
         proj_log_trace(P, "theta: %g", Q->theta);
     }
-    return;
 }
 
 /**************************************************************************/
@@ -302,8 +304,6 @@ static void build_rot_matrix(PJ *P) {
         proj_log_trace(P, "  | % 6.6g  % 6.6g  % 6.6g |", R10, R11, R12);
         proj_log_trace(P, "  | % 6.6g  % 6.6g  % 6.6g |", R20, R21, R22);
     }
-
-    return;
 }
 
 
@@ -610,4 +610,3 @@ PJ *TRANSFORMATION(helmert, 0) {
 
     return P;
 }
-
diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c
index 54440822..3932ef51 100644
--- a/src/PJ_hgridshift.c
+++ b/src/PJ_hgridshift.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <stddef.h>
+
 #include "proj_internal.h"
 #include "projects.h"
 
@@ -72,4 +75,3 @@ PJ *TRANSFORMATION(hgridshift,0) {
 
     return P;
 }
-
diff --git a/src/PJ_horner.c b/src/PJ_horner.c
index 24e1cbe9..09554a7f 100644
--- a/src/PJ_horner.c
+++ b/src/PJ_horner.c
@@ -76,13 +76,17 @@
  *****************************************************************************/
 
 #define PJ_LIB__
+
+#include <errno.h>
+#include <math.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <string.h>
+
 #include "proj_internal.h"
 #include "projects.h"
-#include <stddef.h>
-#include <math.h>
-#include <errno.h>
 
-PROJ_HEAD(horner,    "Horner polynomial evaluation");
+PROJ_HEAD(horner, "Horner polynomial evaluation");
 
 /* make horner.h interface with proj's memory management */
 #define horner_dealloc(x) pj_dealloc(x)
@@ -450,9 +454,9 @@ PJ *PROJECTION(horner) {
     /* Polynomial degree specified? */
     if (pj_param (P->ctx, P->params, "tdeg").i) { /* degree specified? */
         degree = pj_param(P->ctx, P->params, "ideg").i;
-        if (degree > 10000) {
-            /* What is a reasonable maximum for the degree? */
-            proj_log_debug (P, "Horner: Degree too large: %d", degree);
+        if (degree < 0 || degree > 10000) {
+            /* What are reasonable minimum and maximums for degree? */
+            proj_log_debug (P, "Horner: Degree is unreasonable: %d", degree);
             return horner_freeup (P, PJD_ERR_INVALID_ARG);
         }
     } else {
@@ -503,4 +507,3 @@ PJ *PROJECTION(horner) {
 
     return P;
 }
-
diff --git a/src/PJ_igh.c b/src/PJ_igh.c
index bb97a8b5..23614ac5 100644
--- a/src/PJ_igh.c
+++ b/src/PJ_igh.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(igh, "Interrupted Goode Homolosine") "\n\tPCyl, Sph.";
@@ -220,5 +223,3 @@ PJ *PROJECTION(igh) {
 
     return P;
 }
-
-
diff --git a/src/PJ_imw_p.c b/src/PJ_imw_p.c
index c7939d3e..4c1b5e45 100644
--- a/src/PJ_imw_p.c
+++ b/src/PJ_imw_p.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -208,4 +211,3 @@ PJ *PROJECTION(imw_p) {
 
     return P;
 }
-
diff --git a/src/PJ_isea.c b/src/PJ_isea.c
index 14aec843..c1ac2f00 100644
--- a/src/PJ_isea.c
+++ b/src/PJ_isea.c
@@ -2,31 +2,49 @@
  * This code was entirely written by Nathan Wagner
  * and is in the public domain.
  */
-#define PJ_LIB__
 
 #include <errno.h>
 #include <math.h>
+#include <float.h>
 #include <stdio.h>
 #include <stdlib.h>
-#include <float.h>
+#include <string.h>
 
+#define PJ_LIB__
+#include "proj_internal.h"
+#include "proj_math.h"
 #include "proj.h"
 #include "projects.h"
-#include "proj_math.h"
 
-#ifndef M_PI
-#  define M_PI 3.14159265358979323846
-#endif
+#define DEG36 0.62831853071795864768
+#define DEG72 1.25663706143591729537
+#define DEG90 M_PI_2
+#define DEG108 1.88495559215387594306
+#define DEG120 2.09439510239319549229
+#define DEG144 2.51327412287183459075
+#define DEG180 M_PI
 
-/*
- * Proj 4 provides its own entry points into
- * the code, so none of the library functions
- * need to be global
- */
-#define ISEA_STATIC static
-#ifndef ISEA_STATIC
-#define ISEA_STATIC
-#endif
+/* sqrt(5)/M_PI */
+#define ISEA_SCALE 0.8301572857837594396028083
+
+/* 26.565051177 degrees */
+#define V_LAT 0.46364760899944494524
+
+/* 52.62263186 */
+#define E_RAD 0.91843818702186776133
+
+/* 10.81231696 */
+#define F_RAD 0.18871053072122403508
+
+/* R tan(g) sin(60) */
+#define TABLE_G 0.6615845383
+
+/* H = 0.25 R tan g = */
+#define TABLE_H 0.1909830056
+
+/* in radians */
+#define ISEA_STD_LAT 1.01722196792335072101
+#define ISEA_STD_LON .19634954084936207740
 
 struct hex {
         int iso;
@@ -34,8 +52,7 @@ struct hex {
 };
 
 /* y *must* be positive down as the xy /iso conversion assumes this */
-ISEA_STATIC
-void hex_xy(struct hex *h) {
+static void hex_xy(struct hex *h) {
     if (!h->iso) return;
     if (h->x >= 0) {
         h->y = -h->y - (h->x+1)/2;
@@ -46,8 +63,7 @@ void hex_xy(struct hex *h) {
     h->iso = 0;
 }
 
-ISEA_STATIC
-void hex_iso(struct hex *h) {
+static void hex_iso(struct hex *h) {
     if (h->iso) return;
 
     if (h->x >= 0) {
@@ -61,8 +77,7 @@ void hex_iso(struct hex *h) {
     h->iso = 1;
 }
 
-ISEA_STATIC
-void hexbin2(double width, double x, double y, int *i, int *j) {
+static void hexbin2(double width, double x, double y, int *i, int *j) {
     double z, rx, ry, rz;
     double abs_dx, abs_dy, abs_dz;
     long ix, iy, iz, s;
@@ -108,9 +123,6 @@ void hexbin2(double width, double x, double y, int *i, int *j) {
     *i = h.x;
     *j = h.y;
 }
-#ifndef ISEA_STATIC
-#define ISEA_STATIC
-#endif
 
 enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = 20 };
 enum isea_topology { ISEA_HEXAGON=6, ISEA_TRIANGLE=3, ISEA_DIAMOND=4 };
@@ -156,8 +168,7 @@ struct snyder_constants {
 };
 
 /* TODO put these in radians to avoid a later conversion */
-ISEA_STATIC const
-struct snyder_constants constants[] = {
+static const struct snyder_constants constants[] = {
     {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0},
     {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027},
     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
@@ -167,24 +178,7 @@ struct snyder_constants constants[] = {
     {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0},
 };
 
-#define DEG120 2.09439510239319549229
-#define DEG72 1.25663706143591729537
-#define DEG90 1.57079632679489661922
-#define DEG144 2.51327412287183459075
-#define DEG36 0.62831853071795864768
-#define DEG108 1.88495559215387594306
-#define DEG180 M_PI
-/* sqrt(5)/M_PI */
-#define ISEA_SCALE 0.8301572857837594396028083
-
-/* 26.565051177 degrees */
-#define V_LAT 0.46364760899944494524
-
-#define RAD2DEG (180.0/M_PI)
-#define DEG2RAD (M_PI/180.0)
-
-ISEA_STATIC
-struct isea_geo vertex[] = {
+static struct isea_geo vertex[] = {
     {0.0, DEG90},
     {DEG180, V_LAT},
     {-DEG108, V_LAT},
@@ -201,13 +195,7 @@ struct isea_geo vertex[] = {
 
 /* TODO make an isea_pt array of the vertices as well */
 
-static int      tri_v1[] = {0, 0, 0, 0, 0, 0, 6, 7, 8, 9, 10, 2, 3, 4, 5, 1, 11, 11, 11, 11, 11};
-
-/* 52.62263186 */
-#define E_RAD 0.91843818702186776133
-
-/* 10.81231696 */
-#define F_RAD 0.18871053072122403508
+static int tri_v1[] = {0, 0, 0, 0, 0, 0, 6, 7, 8, 9, 10, 2, 3, 4, 5, 1, 11, 11, 11, 11, 11};
 
 /* triangle Centers */
 static const struct isea_geo icostriangles[] = {
@@ -234,8 +222,7 @@ static const struct isea_geo icostriangles[] = {
     {DEG180, -E_RAD},
 };
 
-static double
-az_adjustment(int triangle)
+static double az_adjustment(int triangle)
 {
     double          adj;
 
@@ -253,15 +240,7 @@ az_adjustment(int triangle)
     return adj;
 }
 
-/* R tan(g) sin(60) */
-#define TABLE_G 0.6615845383
-
-/* H = 0.25 R tan g = */
-#define TABLE_H 0.1909830056
-
-ISEA_STATIC
-struct isea_pt
-isea_triangle_xy(int triangle)
+static struct isea_pt isea_triangle_xy(int triangle)
 {
     struct isea_pt  c;
     const double Rprime = 0.91038328153090290025;
@@ -296,8 +275,8 @@ isea_triangle_xy(int triangle)
 }
 
 /* snyder eq 14 */
-static double
-sph_azimuth(double f_lon, double f_lat, double t_lon, double t_lat)
+static double sph_azimuth(double f_lon, double f_lat,
+                          double t_lon, double t_lat)
 {
     double          az;
 
@@ -315,9 +294,7 @@ sph_azimuth(double f_lon, double f_lat, double t_lon, double t_lat)
 #endif
 
 /* coord needs to be in radians */
-ISEA_STATIC
-int
-isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)
+static int isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)
 {
     int             i;
 
@@ -358,9 +335,9 @@ isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)
 
     /* TODO put these constants in as radians to begin with */
     c = constants[SNYDER_POLY_ICOSAHEDRON];
-    theta = c.theta * DEG2RAD;
-    g = c.g * DEG2RAD;
-    G = c.G * DEG2RAD;
+    theta = PJ_TORAD(c.theta);
+    g = PJ_TORAD(c.g);
+    G = PJ_TORAD(c.G);
 
     for (i = 1; i <= 20; i++) {
         double          z;
@@ -480,7 +457,7 @@ isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)
      */
 
     fprintf(stderr, "impossible transform: %f %f is not on any triangle\n",
-        ll->lon * RAD2DEG, ll->lat * RAD2DEG);
+            PJ_TODEG(ll->lon), PJ_TODEG(ll->lat));
 
     exit(EXIT_FAILURE);
 
@@ -507,9 +484,7 @@ isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)
  *
  * TODO take a result pointer
  */
-ISEA_STATIC
-struct isea_geo
-snyder_ctran(struct isea_geo * np, struct isea_geo * pt)
+static struct isea_geo snyder_ctran(struct isea_geo * np, struct isea_geo * pt)
 {
     struct isea_geo npt;
     double          alpha, phi, lambda, lambda0, beta, lambdap, phip;
@@ -553,9 +528,8 @@ snyder_ctran(struct isea_geo * np, struct isea_geo * pt)
     return npt;
 }
 
-ISEA_STATIC
-struct isea_geo
-isea_ctran(struct isea_geo * np, struct isea_geo * pt, double lon0)
+static struct isea_geo isea_ctran(struct isea_geo * np, struct isea_geo * pt,
+                                  double lon0)
 {
     struct isea_geo npt;
 
@@ -580,15 +554,9 @@ isea_ctran(struct isea_geo * np, struct isea_geo * pt, double lon0)
     return npt;
 }
 
-/* in radians */
-#define ISEA_STD_LAT 1.01722196792335072101
-#define ISEA_STD_LON .19634954084936207740
-
 /* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */
 
-ISEA_STATIC
-int
-isea_grid_init(struct isea_dgg * g)
+static int isea_grid_init(struct isea_dgg * g)
 {
     if (!g)
         return 0;
@@ -605,34 +573,26 @@ isea_grid_init(struct isea_dgg * g)
     return 1;
 }
 
-ISEA_STATIC
-void
-isea_orient_isea(struct isea_dgg * g)
+static void isea_orient_isea(struct isea_dgg * g)
 {
     if (!g)
         return;
     g->o_lat = ISEA_STD_LAT;
     g->o_lon = ISEA_STD_LON;
     g->o_az = 0.0;
-    return;
 }
 
-ISEA_STATIC
-void
-isea_orient_pole(struct isea_dgg * g)
+static void isea_orient_pole(struct isea_dgg * g)
 {
     if (!g)
         return;
     g->o_lat = M_PI / 2.0;
     g->o_lon = 0.0;
     g->o_az = 0;
-    return;
 }
 
-ISEA_STATIC
-int
-isea_transform(struct isea_dgg * g, struct isea_geo * in,
-           struct isea_pt * out)
+static int isea_transform(struct isea_dgg * g, struct isea_geo * in,
+                          struct isea_pt * out)
 {
     struct isea_geo i, pole;
     int             tri;
@@ -652,9 +612,7 @@ isea_transform(struct isea_dgg * g, struct isea_geo * in,
 
 #define DOWNTRI(tri) (((tri - 1) / 5) % 2 == 1)
 
-ISEA_STATIC
-void
-isea_rotate(struct isea_pt * pt, double degrees)
+static void isea_rotate(struct isea_pt * pt, double degrees)
 {
     double          rad;
 
@@ -671,8 +629,7 @@ isea_rotate(struct isea_pt * pt, double degrees)
     pt->y = y;
 }
 
-ISEA_STATIC
-int isea_tri_plane(int tri, struct isea_pt *pt, double radius) {
+static int isea_tri_plane(int tri, struct isea_pt *pt, double radius) {
     struct isea_pt tc; /* center of triangle */
 
     if (DOWNTRI(tri)) {
@@ -688,9 +645,7 @@ int isea_tri_plane(int tri, struct isea_pt *pt, double radius) {
 }
 
 /* convert projected triangle coords to quad xy coords, return quad number */
-ISEA_STATIC
-int
-isea_ptdd(int tri, struct isea_pt *pt) {
+static int isea_ptdd(int tri, struct isea_pt *pt) {
     int             downtri, quad;
 
     downtri = (((tri - 1) / 5) % 2 == 1);
@@ -705,9 +660,8 @@ isea_ptdd(int tri, struct isea_pt *pt) {
     return quad;
 }
 
-ISEA_STATIC
-int
-isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di)
+static int isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt,
+                            struct isea_pt *di)
 {
     struct isea_pt  v;
     double          hexwidth;
@@ -784,9 +738,8 @@ isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_p
     return quad;
 }
 
-ISEA_STATIC
-int
-isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di) {
+static int isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt,
+                     struct isea_pt *di) {
     struct isea_pt  v;
     double          hexwidth;
     long            sidelength; /* in hexes */
@@ -857,9 +810,8 @@ isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di)
     return quad;
 }
 
-ISEA_STATIC
-int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,
-            struct isea_pt *di) {
+static int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,
+                     struct isea_pt *di) {
     struct isea_pt  v;
     int             quad;
 
@@ -870,11 +822,11 @@ int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,
 }
 
 /* q2di to seqnum */
-ISEA_STATIC
-int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
-    long            sidelength;
-    long            sn, height;
-    long            hexes;
+
+static int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
+    int             sidelength;
+    int             sn, height;
+    int             hexes;
 
     if (quad == 0) {
         g->serial = 1;
@@ -906,9 +858,8 @@ int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
  * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf
  */
 /* convert a q2di to global hex coord */
-ISEA_STATIC
-int isea_hex(struct isea_dgg *g, int tri,
-        struct isea_pt *pt, struct isea_pt *hex) {
+static int isea_hex(struct isea_dgg *g, int tri,
+                    struct isea_pt *pt, struct isea_pt *hex) {
     struct isea_pt v;
 #ifdef FIXME
     long sidelength;
@@ -969,9 +920,7 @@ int isea_hex(struct isea_dgg *g, int tri,
 #endif
 }
 
-ISEA_STATIC
-struct isea_pt
-isea_forward(struct isea_dgg *g, struct isea_geo *in)
+static struct isea_pt isea_forward(struct isea_dgg *g, struct isea_geo *in)
 {
     int             tri;
     struct isea_pt  out, coord;
@@ -1018,11 +967,11 @@ isea_forward(struct isea_dgg *g, struct isea_geo *in)
 
     return out;
 }
+
 /*
  * Proj 4 integration code follows
  */
 
-
 PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph";
 
 struct pj_opaque {
diff --git a/src/PJ_krovak.c b/src/PJ_krovak.c
index ef00d715..9b269c88 100644
--- a/src/PJ_krovak.c
+++ b/src/PJ_krovak.c
@@ -75,9 +75,11 @@
  *
  *****************************************************************************/
 
-
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(krovak, "Krovak") "\n\tPCyl., Ellps.";
@@ -218,4 +220,3 @@ PJ *PROJECTION(krovak) {
 
     return P;
 }
-
diff --git a/src/PJ_labrd.c b/src/PJ_labrd.c
index 16c45a0d..502adf0f 100644
--- a/src/PJ_labrd.c
+++ b/src/PJ_labrd.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(labrd, "Laborde") "\n\tCyl, Sph\n\tSpecial for Madagascar";
@@ -127,4 +130,3 @@ PJ *PROJECTION(labrd) {
 
     return P;
 }
-
diff --git a/src/PJ_lagrng.c b/src/PJ_lagrng.c
index cf6d97d6..d2b53b09 100644
--- a/src/PJ_lagrng.c
+++ b/src/PJ_lagrng.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -60,4 +63,3 @@ PJ *PROJECTION(lagrng) {
 
     return P;
 }
-
diff --git a/src/PJ_larr.c b/src/PJ_larr.c
index a8d1af56..e4e39dc9 100644
--- a/src/PJ_larr.c
+++ b/src/PJ_larr.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(larr, "Larrivee") "\n\tMisc Sph, no inv.";
@@ -23,4 +26,3 @@ PJ *PROJECTION(larr) {
 
     return P;
 }
-
diff --git a/src/PJ_lcca.c b/src/PJ_lcca.c
index efae5c7c..f0f256b1 100644
--- a/src/PJ_lcca.c
+++ b/src/PJ_lcca.c
@@ -46,7 +46,10 @@
 *****************************************************************************/
 
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -157,4 +160,3 @@ PJ *PROJECTION(lcca) {
 
     return P;
 }
-
diff --git a/src/PJ_loxim.c b/src/PJ_loxim.c
index 36b05b86..512bfc23 100644
--- a/src/PJ_loxim.c
+++ b/src/PJ_loxim.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -70,4 +73,3 @@ PJ *PROJECTION(loxim) {
 
    return P;
 }
-
diff --git a/src/PJ_lsat.c b/src/PJ_lsat.c
index 7e1567e1..810a1cba 100644
--- a/src/PJ_lsat.c
+++ b/src/PJ_lsat.c
@@ -1,6 +1,9 @@
 /* based upon Snyder and Linck, USGS-NMD */
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -205,4 +208,3 @@ PJ *PROJECTION(lsat) {
 
     return P;
 }
-
diff --git a/src/PJ_mbt_fps.c b/src/PJ_mbt_fps.c
index 2fccf3aa..b689014f 100644
--- a/src/PJ_mbt_fps.c
+++ b/src/PJ_mbt_fps.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(mbt_fps, "McBryde-Thomas Flat-Pole Sine (No. 2)") "\n\tCyl., Sph.";
@@ -52,4 +55,3 @@ PJ *PROJECTION(mbt_fps) {
 
     return P;
 }
-
diff --git a/src/PJ_mbtfpp.c b/src/PJ_mbtfpp.c
index 9b90444b..5e5b2d50 100644
--- a/src/PJ_mbtfpp.c
+++ b/src/PJ_mbtfpp.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -60,4 +63,3 @@ PJ *PROJECTION(mbtfpp) {
 
     return P;
 }
-
diff --git a/src/PJ_mbtfpq.c b/src/PJ_mbtfpq.c
index cff7fbdd..75c71708 100644
--- a/src/PJ_mbtfpq.c
+++ b/src/PJ_mbtfpq.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -69,4 +72,3 @@ PJ *PROJECTION(mbtfpq) {
 
     return P;
 }
-
diff --git a/src/PJ_merc.c b/src/PJ_merc.c
index 9123eae1..b75f9b11 100644
--- a/src/PJ_merc.c
+++ b/src/PJ_merc.c
@@ -1,9 +1,12 @@
 #define PJ_LIB__
+
+#include <float.h>
+#include <math.h>
+
 #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";
diff --git a/src/PJ_mill.c b/src/PJ_mill.c
index fdb0b2ad..3ea9636f 100644
--- a/src/PJ_mill.c
+++ b/src/PJ_mill.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(mill, "Miller Cylindrical") "\n\tCyl, Sph";
@@ -32,4 +35,3 @@ PJ *PROJECTION(mill) {
 
     return P;
 }
-
diff --git a/src/PJ_misrsom.c b/src/PJ_misrsom.c
index b9b97117..0308fc42 100644
--- a/src/PJ_misrsom.c
+++ b/src/PJ_misrsom.c
@@ -21,7 +21,10 @@
  *****************************************************************************/
 /* based upon Snyder and Linck, USGS-NMD */
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -212,4 +215,3 @@ PJ *PROJECTION(misrsom) {
 
    return P;
 }
-
diff --git a/src/PJ_moll.c b/src/PJ_moll.c
index 66e6315a..7c0e6f85 100644
--- a/src/PJ_moll.c
+++ b/src/PJ_moll.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(moll, "Mollweide") "\n\tPCyl., Sph.";
@@ -105,4 +108,3 @@ PJ *PROJECTION(wag5) {
 
     return P;
 }
-
diff --git a/src/PJ_molodensky.c b/src/PJ_molodensky.c
index 92969d67..2db51038 100644
--- a/src/PJ_molodensky.c
+++ b/src/PJ_molodensky.c
@@ -43,7 +43,10 @@
 *
 ***********************************************************************/
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "proj_internal.h"
 #include "projects.h"
@@ -311,4 +314,3 @@ PJ *TRANSFORMATION(molodensky,1) {
 
     return P;
 }
-
diff --git a/src/PJ_natearth.c b/src/PJ_natearth.c
index 10abb23c..df015222 100644
--- a/src/PJ_natearth.c
+++ b/src/PJ_natearth.c
@@ -13,6 +13,9 @@ and designed in collaboration with Tom Patterson.
 Port to PROJ.4 by Bernhard Jenny, 6 June 2011
 */
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(natearth, "Natural Earth") "\n\tPCyl., Sph.";
@@ -95,4 +98,3 @@ PJ *PROJECTION(natearth) {
 
     return P;
 }
-
diff --git a/src/PJ_natearth2.c b/src/PJ_natearth2.c
index dbdc8fed..0498f186 100644
--- a/src/PJ_natearth2.c
+++ b/src/PJ_natearth2.c
@@ -6,6 +6,9 @@ and Atmospheric Sciences, Oregon State University.
 Port to PROJ.4 by Bojan Savric, 4 April 2016
 */
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(natearth2, "Natural Earth 2") "\n\tPCyl., Sph.";
@@ -92,4 +95,3 @@ PJ *PROJECTION(natearth2) {
 
     return P;
 }
-
diff --git a/src/PJ_nell.c b/src/PJ_nell.c
index 51f9f99a..751c01c8 100644
--- a/src/PJ_nell.c
+++ b/src/PJ_nell.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(nell, "Nell") "\n\tPCyl., Sph.";
@@ -46,4 +49,3 @@ PJ *PROJECTION(nell) {
 
     return P;
 }
-
diff --git a/src/PJ_nell_h.c b/src/PJ_nell_h.c
index 24957786..db9e9de0 100644
--- a/src/PJ_nell_h.c
+++ b/src/PJ_nell_h.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(nell_h, "Nell-Hammer") "\n\tPCyl., Sph.";
@@ -48,4 +51,3 @@ PJ *PROJECTION(nell_h) {
 
     return P;
 }
-
diff --git a/src/PJ_nocol.c b/src/PJ_nocol.c
index 13688a9f..93e128ea 100644
--- a/src/PJ_nocol.c
+++ b/src/PJ_nocol.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(nicol, "Nicolosi Globular") "\n\tMisc Sph, no inv.";
@@ -49,4 +52,3 @@ PJ *PROJECTION(nicol) {
 
     return P;
 }
-
diff --git a/src/PJ_nzmg.c b/src/PJ_nzmg.c
index b30deabf..bf0862fb 100644
--- a/src/PJ_nzmg.c
+++ b/src/PJ_nzmg.c
@@ -26,6 +26,9 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(nzmg, "New Zealand Map Grid") "\n\tfixed Earth";
@@ -118,4 +121,3 @@ PJ *PROJECTION(nzmg) {
 
     return P;
 }
-
diff --git a/src/PJ_ob_tran.c b/src/PJ_ob_tran.c
index 8b13eb65..f5a05cf9 100644
--- a/src/PJ_ob_tran.c
+++ b/src/PJ_ob_tran.c
@@ -1,8 +1,11 @@
 #define PJ_LIB__
 #include <errno.h>
+#include <math.h>
+#include <stddef.h>
+#include <string.h>
+
 #include "proj.h"
 #include "projects.h"
-#include <string.h>
 
 struct pj_opaque {
     struct PJconsts *link;
@@ -238,4 +241,3 @@ PJ *PROJECTION(ob_tran) {
 
     return P;
 }
-
diff --git a/src/PJ_ocea.c b/src/PJ_ocea.c
index 4bd120e1..7a9353a6 100644
--- a/src/PJ_ocea.c
+++ b/src/PJ_ocea.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(ocea, "Oblique Cylindrical Equal Area") "\n\tCyl, Sph"
@@ -95,4 +98,3 @@ PJ *PROJECTION(ocea) {
 
     return P;
 }
-
diff --git a/src/PJ_omerc.c b/src/PJ_omerc.c
index 1e80b4ca..70c12e27 100644
--- a/src/PJ_omerc.c
+++ b/src/PJ_omerc.c
@@ -22,7 +22,10 @@
 ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -124,9 +127,9 @@ PJ *PROJECTION(omerc) {
 
     Q->no_rot = pj_param(P->ctx, P->params, "bno_rot").i;
         if ((alp = pj_param(P->ctx, P->params, "talpha").i) != 0)
-        alpha_c = pj_param(P->ctx, P->params, "ralpha").f;
+            alpha_c = pj_param(P->ctx, P->params, "ralpha").f;
         if ((gam = pj_param(P->ctx, P->params, "tgamma").i) != 0)
-        gamma = pj_param(P->ctx, P->params, "rgamma").f;
+            gamma = pj_param(P->ctx, P->params, "rgamma").f;
     if (alp || gam) {
         lamc    = pj_param(P->ctx, P->params, "rlonc").f;
         no_off =
@@ -222,4 +225,3 @@ PJ *PROJECTION(omerc) {
 
     return P;
 }
-
diff --git a/src/PJ_patterson.c b/src/PJ_patterson.c
index e58b49bd..a15c730f 100644
--- a/src/PJ_patterson.c
+++ b/src/PJ_patterson.c
@@ -39,6 +39,9 @@
  */
 
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(patterson, "Patterson Cylindrical") "\n\tCyl.";
@@ -112,4 +115,3 @@ PJ *PROJECTION(patterson) {
 
     return P;
 }
-
diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
index aa7d76f8..d623423b 100644
--- a/src/PJ_pipeline.c
+++ b/src/PJ_pipeline.c
@@ -95,13 +95,17 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-20
 ********************************************************************************/
 
 #define PJ_LIB__
+
+#include <errno.h>
+#include <math.h>
+#include <stddef.h>
+#include <string.h>
+
 #include "geodesic.h"
 #include "proj.h"
 #include "proj_internal.h"
 #include "projects.h"
 
-#include <stddef.h>
-#include <errno.h>
 PROJ_HEAD(pipeline,         "Transformation pipeline manager");
 
 /* Projection specific elements for the PJ object */
@@ -114,12 +118,12 @@ struct pj_opaque {
 
 
 
-static PJ_COORD pipeline_forward_4d (PJ_COORD, PJ *P);
-static PJ_COORD pipeline_reverse_4d (PJ_COORD, PJ *P);
+static PJ_COORD pipeline_forward_4d (PJ_COORD point, PJ *P);
+static PJ_COORD pipeline_reverse_4d (PJ_COORD point, PJ *P);
 static XYZ    pipeline_forward_3d (LPZ lpz, PJ *P);
 static LPZ    pipeline_reverse_3d (XYZ xyz, PJ *P);
-static XY     pipeline_forward (LP lpz, PJ *P);
-static LP     pipeline_reverse (XY xyz, PJ *P);
+static XY     pipeline_forward (LP lp, PJ *P);
+static LP     pipeline_reverse (XY xy, PJ *P);
 
 
 
diff --git a/src/PJ_poly.c b/src/PJ_poly.c
index 7cf8bce1..03fa84a4 100644
--- a/src/PJ_poly.c
+++ b/src/PJ_poly.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -163,4 +166,3 @@ PJ *PROJECTION(poly) {
 
     return P;
 }
-
diff --git a/src/PJ_putp2.c b/src/PJ_putp2.c
index 51cf263d..0b4397e5 100644
--- a/src/PJ_putp2.c
+++ b/src/PJ_putp2.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(putp2, "Putnins P2") "\n\tPCyl., Sph.";
@@ -56,4 +59,3 @@ PJ *PROJECTION(putp2) {
 
     return P;
 }
-
diff --git a/src/PJ_putp4p.c b/src/PJ_putp4p.c
index 958f9681..b716060a 100644
--- a/src/PJ_putp4p.c
+++ b/src/PJ_putp4p.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 struct pj_opaque {
@@ -69,4 +72,3 @@ PJ *PROJECTION(weren) {
 
     return P;
 }
-
diff --git a/src/PJ_putp5.c b/src/PJ_putp5.c
index 14bfea23..6836f71d 100644
--- a/src/PJ_putp5.c
+++ b/src/PJ_putp5.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 struct pj_opaque {
@@ -68,4 +71,3 @@ PJ *PROJECTION(putp5p) {
 
     return P;
 }
-
diff --git a/src/PJ_putp6.c b/src/PJ_putp6.c
index 938cef09..68c10701 100644
--- a/src/PJ_putp6.c
+++ b/src/PJ_putp6.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 struct pj_opaque {
@@ -90,4 +93,3 @@ PJ *PROJECTION(putp6p) {
 
     return P;
 }
-
diff --git a/src/PJ_qsc.c b/src/PJ_qsc.c
index 0f5d2299..4b463b46 100644
--- a/src/PJ_qsc.c
+++ b/src/PJ_qsc.c
@@ -39,7 +39,10 @@
  */
 
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 /* The six cube faces. */
@@ -397,4 +400,3 @@ PJ *PROJECTION(qsc) {
 
     return P;
 }
-
diff --git a/src/PJ_rpoly.c b/src/PJ_rpoly.c
index d7241719..57c4fc92 100644
--- a/src/PJ_rpoly.c
+++ b/src/PJ_rpoly.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 struct pj_opaque {
@@ -51,4 +54,3 @@ PJ *PROJECTION(rpoly) {
 
     return P;
 }
-
diff --git a/src/PJ_sch.c b/src/PJ_sch.c
index e9e67ff8..64e5cdb8 100644
--- a/src/PJ_sch.c
+++ b/src/PJ_sch.c
@@ -31,7 +31,10 @@
  ****************************************************************************/
 
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 #include "geocent.h"
@@ -225,4 +228,3 @@ PJ *PROJECTION(sch) {
 
     return setup(P);
 }
-
diff --git a/src/PJ_somerc.c b/src/PJ_somerc.c
index 02df5868..c6c3ff21 100644
--- a/src/PJ_somerc.c
+++ b/src/PJ_somerc.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -87,4 +90,3 @@ PJ *PROJECTION(somerc) {
     P->fwd = e_forward;
     return P;
 }
-
diff --git a/src/PJ_sts.c b/src/PJ_sts.c
index d8866b66..8409d86a 100644
--- a/src/PJ_sts.c
+++ b/src/PJ_sts.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(kav5,    "Kavraisky V")         "\n\tPCyl., Sph.";
@@ -102,4 +105,3 @@ PJ *PROJECTION(mbt_s) {
     P->opaque = Q;
     return setup(P, 1.48875, 1.36509, 0);
 }
-
diff --git a/src/PJ_tcc.c b/src/PJ_tcc.c
index c355a602..14255325 100644
--- a/src/PJ_tcc.c
+++ b/src/PJ_tcc.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -29,4 +32,3 @@ PJ *PROJECTION(tcc) {
 
     return P;
 }
-
diff --git a/src/PJ_tcea.c b/src/PJ_tcea.c
index 415b42de..d30f3df0 100644
--- a/src/PJ_tcea.c
+++ b/src/PJ_tcea.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
-#include    "projects.h"
+
+#include <math.h>
+
+#include "projects.h"
 
 PROJ_HEAD(tcea, "Transverse Cylindrical Equal Area") "\n\tCyl, Sph";
 
@@ -31,4 +34,3 @@ PJ *PROJECTION(tcea) {
     P->es = 0.;
     return P;
 }
-
diff --git a/src/PJ_times.c b/src/PJ_times.c
index 18fe7233..e8b4499f 100644
--- a/src/PJ_times.c
+++ b/src/PJ_times.c
@@ -30,6 +30,9 @@
  *****************************************************************************/
 
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(times, "Times") "\n\tCyl, Sph";
@@ -74,4 +77,3 @@ PJ *PROJECTION(times) {
 
     return P;
 }
-
diff --git a/src/PJ_tmerc.c b/src/PJ_tmerc.c
index f0f34c93..069cdc2c 100644
--- a/src/PJ_tmerc.c
+++ b/src/PJ_tmerc.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -203,4 +206,3 @@ PJ *PROJECTION(tmerc) {
 
     return setup(P);
 }
-
diff --git a/src/PJ_unitconvert.c b/src/PJ_unitconvert.c
index 53bf5e81..a3051ad7 100644
--- a/src/PJ_unitconvert.c
+++ b/src/PJ_unitconvert.c
@@ -64,11 +64,14 @@ Last update: 2017-05-16
 ***********************************************************************/
 
 #define PJ_LIB__
-#include <time.h>
+
 #include <errno.h>
+#include <math.h>
+#include <string.h>
+#include <time.h>
 
-#include "proj_math.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 PROJ_HEAD(unitconvert, "Unit conversion");
diff --git a/src/PJ_urm5.c b/src/PJ_urm5.c
index 091fe9d6..416cb866 100644
--- a/src/PJ_urm5.c
+++ b/src/PJ_urm5.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -49,4 +52,3 @@ PJ *PROJECTION(urm5) {
 
     return P;
 }
-
diff --git a/src/PJ_urmfps.c b/src/PJ_urmfps.c
index bbceab51..5749059b 100644
--- a/src/PJ_urmfps.c
+++ b/src/PJ_urmfps.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -69,4 +72,3 @@ PJ *PROJECTION(wag1) {
     P->opaque->n = 0.8660254037844386467637231707;
     return setup(P);
 }
-
diff --git a/src/PJ_vandg2.c b/src/PJ_vandg2.c
index 81f28a74..20b7601a 100644
--- a/src/PJ_vandg2.c
+++ b/src/PJ_vandg2.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "projects.h"
 
 struct pj_opaque {
@@ -69,4 +72,3 @@ PJ *PROJECTION(vandg3) {
 
     return P;
 }
-
diff --git a/src/PJ_vandg4.c b/src/PJ_vandg4.c
index afffe41f..770594f1 100644
--- a/src/PJ_vandg4.c
+++ b/src/PJ_vandg4.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 PROJ_HEAD(vandg4, "van der Grinten IV") "\n\tMisc Sph, no inv.";
@@ -50,4 +53,3 @@ PJ *PROJECTION(vandg4) {
 
     return P;
 }
-
diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c
index bb8b4d4d..205806b2 100644
--- a/src/PJ_vgridshift.c
+++ b/src/PJ_vgridshift.c
@@ -1,4 +1,7 @@
 #define PJ_LIB__
+
+#include <stddef.h>
+
 #include "proj_internal.h"
 #include "projects.h"
 
@@ -74,4 +77,3 @@ PJ *TRANSFORMATION(vgridshift,0) {
 
     return P;
 }
-
diff --git a/src/PJ_wag2.c b/src/PJ_wag2.c
index 059494d7..d12cade1 100644
--- a/src/PJ_wag2.c
+++ b/src/PJ_wag2.c
@@ -1,6 +1,11 @@
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
+
 PROJ_HEAD(wag2, "Wagner II") "\n\tPCyl., Sph.";
+
 #define C_x 0.92483
 #define C_y 1.38725
 #define C_p1 0.88022
@@ -31,4 +36,3 @@ PJ *PROJECTION(wag2) {
     P->fwd = s_forward;
     return P;
 }
-
diff --git a/src/PJ_wag3.c b/src/PJ_wag3.c
index 5db1e50f..6ae9ad28 100644
--- a/src/PJ_wag3.c
+++ b/src/PJ_wag3.c
@@ -1,6 +1,9 @@
 #define PJ_LIB__
-#include "projects.h"
+
 #include <errno.h>
+#include <math.h>
+
+#include "projects.h"
 
 PROJ_HEAD(wag3, "Wagner III") "\n\tPCyl., Sph.\n\tlat_ts=";
 
@@ -43,4 +46,3 @@ PJ *PROJECTION(wag3) {
 
     return P;
 }
-
diff --git a/src/PJ_wag7.c b/src/PJ_wag7.c
index a820b2cd..3cd332ad 100644
--- a/src/PJ_wag7.c
+++ b/src/PJ_wag7.c
@@ -1,5 +1,8 @@
 #define PJ_LIB__
-#include	"projects.h"
+
+#include <math.h>
+
+#include "projects.h"
 
 PROJ_HEAD(wag7, "Wagner VII") "\n\tMisc Sph, no inv.";
 
@@ -25,4 +28,3 @@ PJ *PROJECTION(wag7) {
     P->es = 0.;
     return P;
 }
-
diff --git a/src/PJ_wink1.c b/src/PJ_wink1.c
index f64f97d2..10392fe9 100644
--- a/src/PJ_wink1.c
+++ b/src/PJ_wink1.c
@@ -1,6 +1,9 @@
 #define PJ_LIB__
-#include "projects.h"
+
 #include <errno.h>
+#include <math.h>
+
+#include "projects.h"
 
 PROJ_HEAD(wink1, "Winkel I") "\n\tPCyl., Sph.\n\tlat_ts=";
 
@@ -39,4 +42,3 @@ PJ *PROJECTION(wink1) {
 
     return P;
 }
-
diff --git a/src/PJ_wink2.c b/src/PJ_wink2.c
index d715074e..d02ca567 100644
--- a/src/PJ_wink2.c
+++ b/src/PJ_wink2.c
@@ -1,6 +1,9 @@
 #define PJ_LIB__
-#include "projects.h"
+
 #include <errno.h>
+#include <math.h>
+
+#include "projects.h"
 
 PROJ_HEAD(wink2, "Winkel II") "\n\tPCyl., Sph., no inv.\n\tlat_1=";
 
@@ -47,4 +50,3 @@ PJ *PROJECTION(wink2) {
 
     return P;
 }
-
diff --git a/src/aasincos.c b/src/aasincos.c
index c3915613..4f9ea25f 100644
--- a/src/aasincos.c
+++ b/src/aasincos.c
@@ -1,34 +1,38 @@
 /* arc sin, cosine, tan2 and sqrt that will NOT fail */
+
+#include <math.h>
+
 #include "projects.h"
+
 #define ONE_TOL	 1.00000000000001
 #define ATOL 1e-50
 
-	double
+        double
 aasin(projCtx ctx,double v) {
-	double av;
-        
-	if ((av = fabs(v)) >= 1.) {
-		if (av > ONE_TOL)
+        double av;
+
+        if ((av = fabs(v)) >= 1.) {
+                if (av > ONE_TOL)
                         pj_ctx_set_errno( ctx, PJD_ERR_ACOS_ASIN_ARG_TOO_LARGE );
-		return (v < 0. ? -M_HALFPI : M_HALFPI);
-	}
-	return asin(v);
+                return (v < 0. ? -M_HALFPI : M_HALFPI);
+        }
+        return asin(v);
 }
 
-	double
+        double
 aacos(projCtx ctx, double v) {
-	double av;
+        double av;
 
-	if ((av = fabs(v)) >= 1.) {
-		if (av > ONE_TOL)
+        if ((av = fabs(v)) >= 1.) {
+                if (av > ONE_TOL)
                         pj_ctx_set_errno( ctx, PJD_ERR_ACOS_ASIN_ARG_TOO_LARGE );
-		return (v < 0. ? M_PI : 0.);
-	}
-	return acos(v);
+                return (v < 0. ? M_PI : 0.);
+        }
+        return acos(v);
 }
-	double
+        double
 asqrt(double v) { return ((v <= 0) ? 0. : sqrt(v)); }
-	double
+        double
 aatan2(double n, double d) {
-	return ((fabs(n) < ATOL && fabs(d) < ATOL) ? 0. : atan2(n,d));
+        return ((fabs(n) < ATOL && fabs(d) < ATOL) ? 0. : atan2(n,d));
 }
diff --git a/src/cct.c b/src/cct.c
index d9d25329..8dd1e0ad 100644
--- a/src/cct.c
+++ b/src/cct.c
@@ -88,7 +88,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-10-26
 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, ...);
+static void print(PJ_LOG_LEVEL log_level, const char *fmt, ...);
 
 /* Prototypes from functions in this file */
 char *column (char *buf, int n);
diff --git a/src/cs2cs.c b/src/cs2cs.c
index de3259d2..167eafa3 100644
--- a/src/cs2cs.c
+++ b/src/cs2cs.c
@@ -26,15 +26,16 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
-#include "proj.h"
-#include "projects.h"
+#include <ctype.h>
+#include <locale.h>
+#include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
-#include <ctype.h>
 #include <string.h>
-#include <math.h>
+
+#include "proj.h"
+#include "projects.h"
 #include "emess.h"
-#include <locale.h>
 
 #define MAX_LINE 1000
 #define MAX_PARGS 100
diff --git a/src/dmstor.c b/src/dmstor.c
index 23871a85..d3496cd9 100644
--- a/src/dmstor.c
+++ b/src/dmstor.c
@@ -1,7 +1,11 @@
 /* Convert DMS string to radians */
-#include "projects.h"
-#include <string.h>
+
 #include <ctype.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "projects.h"
 
 static double proj_strtod(char *nptr, char **endptr);
 
diff --git a/src/emess.c b/src/emess.c
index a645bac9..330b8412 100644
--- a/src/emess.c
+++ b/src/emess.c
@@ -9,14 +9,16 @@
 #  endif
 #endif
 
+#include <errno.h>
+#include <stdarg.h>
 #include <stdio.h>
 #include <stdlib.h>
-#include <stdarg.h>
-#include <errno.h>
 #include <string.h>
+
 #include "proj_api.h"
 #define EMESS_ROUTINE
 #include "emess.h"
+
 	void
 emess(int code, const char *fmt, ...) {
 	va_list args;
diff --git a/src/geocent.c b/src/geocent.c
index a8cec043..c023bdd3 100644
--- a/src/geocent.c
+++ b/src/geocent.c
@@ -432,6 +432,5 @@ void pj_Convert_Geocentric_To_Geodetic (GeocentricInfo *gi,
 /*	ellipsoidal (geodetic) latitude */
     *Latitude=atan2(SPHI, fabs(CPHI));
 
-    return;
 #endif /* defined(USE_ITERATIVE_METHOD) */
 } /* END OF Convert_Geocentric_To_Geodetic */
diff --git a/src/geod_set.c b/src/geod_set.c
index 571a7d20..479f18da 100644
--- a/src/geod_set.c
+++ b/src/geod_set.c
@@ -1,11 +1,14 @@
-
 #define _IN_GEOD_SET
 
+#include <math.h>
+#include <stdlib.h>
 #include <string.h>
+
 #include "proj.h"
 #include "projects.h"
 #include "geod_interface.h"
 #include "emess.h"
+
 	void
 geod_set(int argc, char **argv) {
 	paralist *start = 0, *curr;
diff --git a/src/geodesic.c b/src/geodesic.c
index 9904c7fa..3fcfd1c9 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -331,7 +331,7 @@ static real Lambda12(const struct geod_geodesic* g,
                      real* pssig1, real* pcsig1,
                      real* pssig2, real* pcsig2,
                      real* peps,
-                     real* pgomg12,
+                     real* pdomg12,
                      boolx diffp, real* pdlam12,
                      /* Scratch area of the right size */
                      real Ca[]);
@@ -481,10 +481,10 @@ void geod_lineinit(struct geod_geodesicline* l,
 void geod_gendirectline(struct geod_geodesicline* l,
                         const struct geod_geodesic* g,
                         real lat1, real lon1, real azi1,
-                        unsigned flags, real a12_s12,
+                        unsigned flags, real s12_a12,
                         unsigned caps) {
   geod_lineinit(l, g, lat1, lon1, azi1, caps);
-  geod_gensetdistance(l, flags, a12_s12);
+  geod_gensetdistance(l, flags, s12_a12);
 }
 
 void geod_directline(struct geod_geodesicline* l,
diff --git a/src/gie.c b/src/gie.c
index f9f76452..66da8d85 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -108,6 +108,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-10-01/2017-10-08
 #include <errno.h>
 #include <math.h>
 #include <stdarg.h>
+#include <stddef.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -134,15 +135,15 @@ typedef struct ffio {
     size_t level;
 }  ffio;
 
-static int get_inp (ffio *F);
-static int skip_to_next_tag (ffio *F);
-static int step_into_gie_block (ffio *F);
-static int locate_tag (ffio *F, const char *tag);
-static int nextline (ffio *F);
-static int at_end_delimiter (ffio *F);
-static const char *at_tag (ffio *F);
-static int at_decorative_element (ffio *F);
-static ffio *ffio_destroy (ffio *F);
+static int get_inp (ffio *G);
+static int skip_to_next_tag (ffio *G);
+static int step_into_gie_block (ffio *G);
+static int locate_tag (ffio *G, const char *tag);
+static int nextline (ffio *G);
+static int at_end_delimiter (ffio *G);
+static const char *at_tag (ffio *G);
+static int at_decorative_element (ffio *G);
+static ffio *ffio_destroy (ffio *G);
 static ffio *ffio_create (const char **tags, size_t n_tags, size_t max_record_size);
 
 static const char *gie_tags[] = {
@@ -303,9 +304,11 @@ int main (int argc, char **argv) {
         process_file (o->fargv[i]);
 
     if (T.verbosity > 0) {
-        if (o->fargc > 1)
-        fprintf (T.fout, "%sGrand total: %d. Success: %d, Skipped: %d, Failure: %d\n",
-                 delim, T.grand_ok+T.grand_ko+T.grand_skip, T.grand_ok, T.grand_skip, T.grand_ko);
+        if (o->fargc > 1) {
+            fprintf (T.fout, "%sGrand total: %d. Success: %d, Skipped: %d, Failure: %d\n",
+                     delim, T.grand_ok+T.grand_ko+T.grand_skip, T.grand_ok, T.grand_skip,
+                     T.grand_ko);
+        }
         fprintf (T.fout, "%s", delim);
         if (T.verbosity > 1) {
             fprintf (T.fout, "Failing roundtrips: %4d,    Succeeding roundtrips: %4d\n", fail_rtps, succ_rtps);
@@ -414,10 +417,11 @@ static int process_file (const char *fname) {
     T.grand_ok   += T.total_ok;
     T.grand_ko   += T.total_ko;
     T.grand_skip += T.grand_skip;
-    if (T.verbosity > 0)
-    fprintf (T.fout, "%stotal: %2d tests succeeded, %2d tests skipped, %2d tests %s\n",
-             delim, T.total_ok, T.total_skip, T.total_ko, T.total_ko? "FAILED!": "failed.");
-
+    if (T.verbosity > 0) {
+        fprintf (T.fout, "%stotal: %2d tests succeeded, %2d tests skipped, %2d tests %s\n",
+                 delim, T.total_ok, T.total_skip, T.total_ko,
+                 T.total_ko? "FAILED!": "failed.");
+    }
     if (F->level==0)
         return errmsg (-3, "File '%s':Missing '<gie>' cmnd - bye!\n", fname);
     if (F->level && F->level%2)
@@ -1160,30 +1164,6 @@ See the PROJ ".gie" test suites for examples of supported formatting.
 
 ****************************************************************************************/
 
-#include <stdio.h>
-#include <stdlib.h>
-#include <stdarg.h>
-
-#include <string.h>
-#include <ctype.h>
-
-#include <math.h>
-#include <errno.h>
-
-
-
-static int get_inp (ffio *F);
-static int skip_to_next_tag (ffio *F);
-static int step_into_gie_block (ffio *F);
-static int locate_tag (ffio *F, const char *tag);
-static int nextline (ffio *F);
-static int at_end_delimiter (ffio *F);
-static const char *at_tag (ffio *F);
-static int at_decorative_element (ffio *F);
-static ffio *ffio_destroy (ffio *F);
-static ffio *ffio_create (const char **tags, size_t n_tags, size_t max_record_size);
-
-
 
 /***************************************************************************************/
 static ffio *ffio_create (const char **tags, size_t n_tags, size_t max_record_size) {
diff --git a/src/nad_cvt.c b/src/nad_cvt.c
index d4330d54..ec4a2b47 100644
--- a/src/nad_cvt.c
+++ b/src/nad_cvt.c
@@ -1,6 +1,11 @@
 #define PJ_LIB__
+
+#include <stdio.h>
+#include <stdlib.h>
+
 #include "projects.h"
 #include "proj_math.h"
+
 #define MAX_ITERATIONS 10
 #define TOL 1e-12
 
diff --git a/src/nad_init.c b/src/nad_init.c
index 844dcad2..8b024ac7 100644
--- a/src/nad_init.c
+++ b/src/nad_init.c
@@ -27,11 +27,14 @@
 
 #define PJ_LIB__
 
-#include "projects.h"
-#include <stdio.h>
 #include <errno.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
 #include <string.h>
 
+#include "projects.h"
+
 /************************************************************************/
 /*                             swap_words()                             */
 /*                                                                      */
diff --git a/src/pj_apply_gridshift.c b/src/pj_apply_gridshift.c
index 1ef39b20..31e0124e 100644
--- a/src/pj_apply_gridshift.c
+++ b/src/pj_apply_gridshift.c
@@ -30,8 +30,10 @@
 
 #define PJ_LIB__
 
-#include <string.h>
 #include <math.h>
+#include <stdio.h>
+#include <string.h>
+
 #include "proj_internal.h"
 #include "projects.h"
 
@@ -174,7 +176,7 @@ static struct CTABLE* find_ctable(projCtx ctx, LP input, int grid_count, PJ_GRID
 /*      This is the real workhorse, given a gridlist.                   */
 /************************************************************************/
 
-int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count,
+int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **gridlist, int gridlist_count,
                           int inverse, long point_count, int point_offset,
                           double *x, double *y, double *z )
 {
@@ -183,7 +185,7 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count,
     static int debug_count = 0;
     (void) z;
 
-    if( tables == NULL || grid_count == 0 )
+    if( gridlist== NULL || gridlist_count == 0 )
     {
         pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
         return PJD_ERR_FAILED_TO_LOAD_GRID;
@@ -202,7 +204,7 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count,
         output.phi = HUGE_VAL;
         output.lam = HUGE_VAL;
 
-        ct = find_ctable(ctx, input, grid_count, tables);
+        ct = find_ctable(ctx, input, gridlist_count, gridlist);
         if( ct != NULL )
         {
             output = nad_cvt( input, inverse, ct );
@@ -220,9 +222,9 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count,
                     "                      location (%.7fdW,%.7fdN)",
                     x[io] * RAD_TO_DEG,
                     y[io] * RAD_TO_DEG );
-                for( itable = 0; itable < grid_count; itable++ )
+                for( itable = 0; itable < gridlist_count; itable++ )
                 {
-                    PJ_GRIDINFO *gi = tables[itable];
+                    PJ_GRIDINFO *gi = gridlist[itable];
                     if( itable == 0 )
                         pj_log( ctx, PJ_LOG_DEBUG_MAJOR, "   tried: %s", gi->gridname );
                     else
diff --git a/src/pj_apply_vgridshift.c b/src/pj_apply_vgridshift.c
index b5b35640..c1344951 100644
--- a/src/pj_apply_vgridshift.c
+++ b/src/pj_apply_vgridshift.c
@@ -28,11 +28,22 @@
 
 #define PJ_LIB__
 
+#include <stdio.h>
 #include <string.h>
+
 #include "proj_math.h"
 #include "proj_internal.h"
 #include "projects.h"
 
+static int is_nodata(float value)
+{
+    /* nodata?  */
+    /* GTX official nodata value if  -88.88880f, but some grids also */
+    /* use other  big values for nodata (e.g naptrans2008.gtx has */
+    /* nodata values like -2147479936), so test them too */
+    return value > 1000 || value < -1000 || value == -88.88880f;
+}
+
 static double read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GRIDINFO **tables, struct CTABLE *ct) {
     int  itable = 0;
     double value = HUGE_VAL;
@@ -110,23 +121,49 @@ static double read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GR
             grid_iy2 = ct->lim.phi - 1;
 
         cvs = (float *) ct->cvs;
-        value = cvs[grid_ix + grid_iy * ct->lim.lam]
-            * (1.0-grid_x) * (1.0-grid_y)
-            + cvs[grid_ix2 + grid_iy * ct->lim.lam]
-            * (grid_x) * (1.0-grid_y)
-            + cvs[grid_ix + grid_iy2 * ct->lim.lam]
-            * (1.0-grid_x) * (grid_y)
-            + cvs[grid_ix2 + grid_iy2 * ct->lim.lam]
-            * (grid_x) * (grid_y);
+        {
+            float value_a = cvs[grid_ix + grid_iy * ct->lim.lam];
+            float value_b = cvs[grid_ix2 + grid_iy * ct->lim.lam];
+            float value_c = cvs[grid_ix + grid_iy2 * ct->lim.lam];
+            float value_d = cvs[grid_ix2 + grid_iy2 * ct->lim.lam];
+            double total_weight = 0.0;
+            int n_weights = 0;
+            value = 0.0f;
+            if( !is_nodata(value_a) )
+            {
+                double weight = (1.0-grid_x) * (1.0-grid_y);
+                value += value_a * weight;
+                total_weight += weight;
+                n_weights ++;
+            }
+            if( !is_nodata(value_b) )
+            {
+                double weight = (grid_x) * (1.0-grid_y);
+                value += value_b * weight;
+                total_weight += weight;
+                n_weights ++;
+            }
+            if( !is_nodata(value_c) )
+            {
+                double weight = (1.0-grid_x) * (grid_y);
+                value += value_c * weight;
+                total_weight += weight;
+                n_weights ++;
+            }
+            if( !is_nodata(value_d) )
+            {
+                double weight = (grid_x) * (grid_y);
+                value += value_d * weight;
+                total_weight += weight;
+                n_weights ++;
+            }
+            if( n_weights == 0 )
+                value = HUGE_VAL;
+            else if( n_weights != 4 )
+                value /= total_weight;
+        }
 
     }
-    /* nodata?  */
-    /* GTX official nodata value if  -88.88880f, but some grids also */
-    /* use other  big values for nodata (e.g naptrans2008.gtx has */
-    /* nodata values like -2147479936), so test them too */
-    if( value > 1000 || value < -1000 || value == -88.88880f )
-        value = HUGE_VAL;
-
 
     return value;
 }
@@ -292,4 +329,3 @@ double proj_vgrid_value(PJ *P, LP lp){
 
     return value;
 }
-
diff --git a/src/pj_auth.c b/src/pj_auth.c
index ed716fde..d6024671 100644
--- a/src/pj_auth.c
+++ b/src/pj_auth.c
@@ -1,5 +1,10 @@
 /* determine latitude from authalic latitude */
+
+#include <math.h>
+#include <stddef.h>
+
 #include "projects.h"
+
 # define P00 .33333333333333333333 /*   1 /     3 */
 # define P01 .17222222222222222222 /*  31 /   180 */
 # define P02 .10257936507936507937 /* 517 /  5040 */
@@ -7,6 +12,7 @@
 # define P11 .06640211640211640212 /* 251 /  3780 */
 # define P20 .01677689594356261023 /* 761 / 45360 */
 #define APA_SIZE 3
+
 	double *
 pj_authset(double es) {
 	double t, *APA;
diff --git a/src/pj_ctx.c b/src/pj_ctx.c
index 1b2aea27..d77bc88a 100644
--- a/src/pj_ctx.c
+++ b/src/pj_ctx.c
@@ -25,9 +25,11 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
-#include "projects.h"
-#include <string.h>
 #include <errno.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "projects.h"
 
 static projCtx_t default_context;
 static volatile int       default_context_initialized = 0;
diff --git a/src/pj_datum_set.c b/src/pj_datum_set.c
index e527ff84..466b56c5 100644
--- a/src/pj_datum_set.c
+++ b/src/pj_datum_set.c
@@ -26,9 +26,10 @@
  *****************************************************************************/
 
 #include <errno.h>
-#include "projects.h"
 #include <string.h>
 
+#include "projects.h"
+
 /* SEC_TO_RAD = Pi/180/3600 */
 #define SEC_TO_RAD 4.84813681109535993589914102357e-6
 
diff --git a/src/pj_datums.c b/src/pj_datums.c
index 4bb23f24..2011d1bb 100644
--- a/src/pj_datums.c
+++ b/src/pj_datums.c
@@ -25,6 +25,8 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
+#include <stddef.h>
+
 #include "proj.h"
 
 #define PJ_DATUMS__
diff --git a/src/pj_deriv.c b/src/pj_deriv.c
index 7c1fcded..2c91e0cc 100644
--- a/src/pj_deriv.c
+++ b/src/pj_deriv.c
@@ -1,5 +1,8 @@
 /* dervative of (*P->fwd) projection */
 #define PJ_LIB__
+
+#include <math.h>
+
 #include "projects.h"
 
 int pj_deriv(LP lp, double h, const PJ *P, struct DERIVS *der) {
diff --git a/src/pj_ell_set.c b/src/pj_ell_set.c
index fc655289..fd27bd9d 100644
--- a/src/pj_ell_set.c
+++ b/src/pj_ell_set.c
@@ -1,6 +1,9 @@
 /* set ellipsoid parameters a and es */
+
+#include <math.h>
+#include <stddef.h>
 #include <string.h>
-#include <errno.h>
+
 #include "proj.h"
 #include "proj_internal.h"
 #include "projects.h"
diff --git a/src/pj_ellps.c b/src/pj_ellps.c
index 4005d1ce..2bc1e2bb 100644
--- a/src/pj_ellps.c
+++ b/src/pj_ellps.c
@@ -1,5 +1,7 @@
 /* definition of standard geoids */
 
+#include <stddef.h>
+
 #include "proj.h"
 
 #define PJ_ELLPS__
diff --git a/src/pj_fileapi.c b/src/pj_fileapi.c
index d7d3f00b..b08a8b95 100644
--- a/src/pj_fileapi.c
+++ b/src/pj_fileapi.c
@@ -27,9 +27,13 @@
  *****************************************************************************/
 
 #include <errno.h>
-#include "projects.h"
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
 #include <string.h>
 
+#include "projects.h"
+
 static PAFile stdio_fopen(projCtx ctx, const char *filename, 
                              const char *access);
 static size_t stdio_fread(void *buffer, size_t size, size_t nmemb, 
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index 4b6cfc58..b5f1b36e 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -27,7 +27,9 @@
  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
+
 #include <errno.h>
+#include <math.h>
 
 #include "proj_internal.h"
 #include "proj_math.h"
diff --git a/src/pj_gauss.c b/src/pj_gauss.c
index 108ac059..4520bb39 100644
--- a/src/pj_gauss.c
+++ b/src/pj_gauss.c
@@ -24,6 +24,10 @@
 ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */
 #define PJ_LIB__
+
+#include <math.h>
+#include <stdlib.h>
+
 #include "projects.h"
 
 #define MAX_ITER 20
diff --git a/src/pj_gc_reader.c b/src/pj_gc_reader.c
index c38cbf17..493fc075 100644
--- a/src/pj_gc_reader.c
+++ b/src/pj_gc_reader.c
@@ -27,10 +27,12 @@
 
 #define PJ_LIB__
 
+#include <ctype.h>
 #include <errno.h>
-#include "projects.h"
+#include <stdlib.h>
 #include <string.h>
-#include <ctype.h>
+
+#include "projects.h"
 
 static int gc_readentry(projCtx ctx, PAFile fid, PJ_GridCatalogEntry *entry);
 
@@ -46,7 +48,7 @@ PJ_GridCatalog *pj_gc_readcatalog( projCtx ctx, const char *catalog_name )
     PJ_GridCatalog *catalog;
     int entry_max;
     char line[302];
-    
+
     fid = pj_open_lib( ctx, catalog_name, "r" );
     if (fid == NULL) 
         return NULL;
diff --git a/src/pj_gridcatalog.c b/src/pj_gridcatalog.c
index 0498d5b0..ea9c4aa1 100644
--- a/src/pj_gridcatalog.c
+++ b/src/pj_gridcatalog.c
@@ -27,9 +27,12 @@
 
 #define PJ_LIB__
 
-#include "projects.h"
-#include <string.h>
 #include <assert.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "projects.h"
 
 static PJ_GridCatalog *grid_catalog_list = NULL;
 
@@ -234,7 +237,7 @@ int pj_gc_apply_gridshift( PJ *defn, int inverse,
 
 PJ_GRIDINFO *pj_gc_findgrid( projCtx ctx, PJ_GridCatalog *catalog, int after,
                              LP location, double date,
-                             PJ_Region *optimal_region,
+                             PJ_Region *optional_region,
                              double *grid_date ) 
 {
     int iEntry;
@@ -264,15 +267,15 @@ PJ_GRIDINFO *pj_gc_findgrid( projCtx ctx, PJ_GridCatalog *catalog, int after,
     {
         if( grid_date )
             *grid_date = 0.0;
-        if( optimal_region != NULL )
-            memset( optimal_region, 0, sizeof(PJ_Region));
+        if( optional_region != NULL )
+            memset( optional_region, 0, sizeof(PJ_Region));
         return NULL;
     }
 
     if( grid_date )
         *grid_date = entry->date;
 
-    if( optimal_region )
+    if( optional_region )
     {
         
     }
diff --git a/src/pj_gridinfo.c b/src/pj_gridinfo.c
index 442b3fe6..f201f39e 100644
--- a/src/pj_gridinfo.c
+++ b/src/pj_gridinfo.c
@@ -28,11 +28,15 @@
 
 #define PJ_LIB__
 
+#include <errno.h>
+#include <math.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
 #include "proj_internal.h"
 #include "projects.h"
-#include <string.h>
-#include <math.h>
-#include <errno.h>
 
 /************************************************************************/
 /*                             swap_words()                             */
diff --git a/src/pj_gridlist.c b/src/pj_gridlist.c
index 9daabbe1..332de8dd 100644
--- a/src/pj_gridlist.c
+++ b/src/pj_gridlist.c
@@ -29,9 +29,10 @@
 #define PJ_LIB__
 
 #include <errno.h>
-#include "projects.h"
+#include <stddef.h>
 #include <string.h>
-#include <math.h>
+
+#include "projects.h"
 
 static PJ_GRIDINFO *grid_list = NULL;
 #define PJ_MAX_PATH_LENGTH 1024
diff --git a/src/pj_init.c b/src/pj_init.c
index bdaf64f7..c5e92105 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -27,14 +27,16 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
-
-
 #define PJ_LIB__
-#include "geodesic.h"
+
+#include <ctype.h>
+#include <errno.h>
+#include <math.h>
+#include <stddef.h>
 #include <stdio.h>
 #include <string.h>
-#include <errno.h>
-#include <ctype.h>
+
+#include "geodesic.h"
 #include "proj.h"
 #include "proj_internal.h"
 #include "proj_math.h"
diff --git a/src/pj_initcache.c b/src/pj_initcache.c
index a0593138..3c347e4b 100644
--- a/src/pj_initcache.c
+++ b/src/pj_initcache.c
@@ -25,9 +25,10 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
-#include "projects.h"
 #include <string.h>
 
+#include "projects.h"
+
 static int cache_count = 0;
 static int cache_alloc = 0;
 static char **cache_key = NULL;
@@ -51,7 +52,7 @@ paralist *pj_clone_paralist( const paralist *list)
       newitem->used = 0;
       newitem->next = 0;
       strcpy( newitem->param, list->param );
-      
+
       if( list_copy == NULL )
 	list_copy = newitem;
       else
diff --git a/src/pj_internal.c b/src/pj_internal.c
index dc528649..58ea9b48 100644
--- a/src/pj_internal.c
+++ b/src/pj_internal.c
@@ -28,10 +28,14 @@
  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
+
 #include <ctype.h>
-#include <stddef.h>
-#include <stdarg.h>
 #include <errno.h>
+#include <math.h>
+#include <stdarg.h>
+#include <stddef.h>
+#include <stdlib.h>
+#include <string.h>
 
 #include "geodesic.h"
 #include "proj_internal.h"
@@ -134,7 +138,6 @@ void proj_context_set (PJ *P, PJ_CONTEXT *ctx) {
     if (0==ctx)
         ctx = pj_get_default_ctx ();
     pj_set_ctx (P, ctx);
-    return;
 }
 
 
@@ -143,7 +146,6 @@ void proj_context_inherit (PJ *parent, PJ *child) {
         pj_set_ctx (child, pj_get_default_ctx());
     else
         pj_set_ctx (child, pj_get_ctx(parent));
-    return;
 }
 
 
@@ -361,7 +363,6 @@ to that context.
     if (0==ctx)
         ctx = pj_get_default_ctx();
     pj_ctx_set_errno (ctx, err);
-    return;
 }
 
 /*  logging  */
diff --git a/src/pj_inv.c b/src/pj_inv.c
index 25d416c0..ca149674 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -28,6 +28,7 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 #include <errno.h>
+#include <math.h>
 
 #include "proj_internal.h"
 #include "proj_math.h"
diff --git a/src/pj_log.c b/src/pj_log.c
index 19151800..6654691c 100644
--- a/src/pj_log.c
+++ b/src/pj_log.c
@@ -25,10 +25,13 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
+#include <stdarg.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
 #include "proj.h"
 #include "projects.h"
-#include <string.h>
-#include <stdarg.h>
 
 /************************************************************************/
 /*                          pj_stderr_logger()                          */
diff --git a/src/pj_malloc.c b/src/pj_malloc.c
index 2d670ae2..c45da85d 100644
--- a/src/pj_malloc.c
+++ b/src/pj_malloc.c
@@ -34,15 +34,18 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
-
 /* allocate and deallocate memory */
 /* These routines are used so that applications can readily replace
 ** projection system memory allocation/deallocation call with custom
 ** application procedures.  */
 
+#include <errno.h>
+#include <stddef.h>
+#include <stdlib.h>
+#include <string.h>
+
 #include "proj.h"
 #include "projects.h"
-#include <errno.h>
 
 /**********************************************************************/
 void *pj_malloc(size_t size) {
diff --git a/src/pj_mlfn.c b/src/pj_mlfn.c
index 24e7b63a..02e05c3a 100644
--- a/src/pj_mlfn.c
+++ b/src/pj_mlfn.c
@@ -1,4 +1,7 @@
+#include <math.h>
+
 #include "projects.h"
+
 /* meridional distance for ellipsoid and inverse
 **	8th degree - accurate to < 1e-5 meters when used in conjunction
 **		with typical major axis values.
diff --git a/src/pj_open_lib.c b/src/pj_open_lib.c
index 8ecbbd79..6b908360 100644
--- a/src/pj_open_lib.c
+++ b/src/pj_open_lib.c
@@ -29,11 +29,15 @@
  *****************************************************************************/
 
 #define PJ_LIB__
-#include "proj_internal.h"
-#include "projects.h"
+
+#include <errno.h>
+#include <stddef.h>
 #include <stdio.h>
+#include <stdlib.h>
 #include <string.h>
-#include <errno.h>
+
+#include "proj_internal.h"
+#include "projects.h"
 
 static const char *(*pj_finder)(const char *) = NULL;
 static int path_count = 0;
diff --git a/src/pj_param.c b/src/pj_param.c
index b52ac940..6cee4d1f 100644
--- a/src/pj_param.c
+++ b/src/pj_param.c
@@ -1,6 +1,9 @@
 /* put parameters in linked list and retrieve */
+
 #include <ctype.h>
+#include <stddef.h>
 #include <stdio.h>
+#include <stdlib.h>
 #include <string.h>
 
 #include "projects.h"
diff --git a/src/pj_pr_list.c b/src/pj_pr_list.c
index baf8125f..4e71e471 100644
--- a/src/pj_pr_list.c
+++ b/src/pj_pr_list.c
@@ -1,7 +1,11 @@
 /* print projection's list of parameters */
-#include "projects.h"
+
+#include <stddef.h>
 #include <stdio.h>
 #include <string.h>
+
+#include "projects.h"
+
 #define LINE_LEN 72
 	static int
 pr_list(PJ *P, int not_used) {
diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c
index d881ae7b..e06ea58b 100644
--- a/src/pj_strerrno.c
+++ b/src/pj_strerrno.c
@@ -1,7 +1,9 @@
 /* list of projection system pj_errno values */
+
+#include <stddef.h>
 #include <stdio.h>
-#include <errno.h>
 #include <string.h>
+
 #include "proj.h"
 #include "projects.h"
 
diff --git a/src/pj_strtod.c b/src/pj_strtod.c
index 1cf94c9b..45159d11 100644
--- a/src/pj_strtod.c
+++ b/src/pj_strtod.c
@@ -27,6 +27,11 @@
  * DEALINGS IN THE SOFTWARE.
  ****************************************************************************/
 
+#include <errno.h>
+#include <locale.h>
+#include <stdlib.h>
+#include <string.h>
+
 #include "projects.h"
 
 /* Windows nmake build doesn't have a proj_config.h, but HAVE_LOCALECONV */
@@ -35,10 +40,6 @@
 #include "proj_config.h"
 #endif
 
-#include <stdlib.h>
-#include <locale.h>
-#include <errno.h>
-
 #define PJ_STRTOD_WORK_BUFFER_SIZE 64
 
 /************************************************************************/
@@ -193,4 +194,3 @@ double pj_strtod( const char *nptr, char **endptr )
     errno = nError;
     return dfValue;
 }
-
diff --git a/src/pj_transform.c b/src/pj_transform.c
index 168711ba..5863b39e 100644
--- a/src/pj_transform.c
+++ b/src/pj_transform.c
@@ -27,9 +27,10 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
-#include "projects.h"
-#include <string.h>
 #include <math.h>
+#include <string.h>
+
+#include "projects.h"
 #include "geocent.h"
 
 
@@ -495,66 +496,66 @@ static int long_wrap (PJ *P, long n, int dist, double *x) {
 /************************************************************************/
 
 int pj_transform(
-    PJ *srcdefn, PJ *dstdefn,
+    PJ *src, PJ *dst,
     long point_count, int point_offset,
     double *x, double *y, double *z
 ){
     int       err;
 
-    srcdefn->ctx->last_errno = 0;
-    dstdefn->ctx->last_errno = 0;
+    src->ctx->last_errno = 0;
+    dst->ctx->last_errno = 0;
 
     if( point_offset == 0 )
         point_offset = 1;
 
     /* Bring input to "normal form": longitude, latitude, ellipsoidal height */
 
-    err = adjust_axes (srcdefn, PJ_INV, point_count, point_offset, x, y, z);
+    err = adjust_axes (src, PJ_INV, point_count, point_offset, x, y, z);
     if (err)
         return err;
-    err = geographic_to_cartesian (srcdefn, PJ_INV, point_count, point_offset, x, y, z);
+    err = geographic_to_cartesian (src, PJ_INV, point_count, point_offset, x, y, z);
     if (err)
         return err;
-    err = projected_to_geographic (srcdefn, point_count, point_offset, x, y, z);
+    err = projected_to_geographic (src, point_count, point_offset, x, y, z);
     if (err)
         return err;
-    err = prime_meridian (srcdefn, PJ_INV, point_count, point_offset, x);
+    err = prime_meridian (src, PJ_INV, point_count, point_offset, x);
     if (err)
         return err;
-    err = height_unit (srcdefn, PJ_INV, point_count, point_offset, z);
+    err = height_unit (src, PJ_INV, point_count, point_offset, z);
     if (err)
         return err;
-    err = geometric_to_orthometric (srcdefn, PJ_INV, point_count, point_offset, x, y, z);
+    err = geometric_to_orthometric (src, PJ_INV, point_count, point_offset, x, y, z);
     if (err)
         return err;
 
     /* At the center of the process we do the datum shift (if needed) */
 
-    err = datum_transform(srcdefn, dstdefn, point_count, point_offset, x, y, z );
+    err = datum_transform(src, dst, point_count, point_offset, x, y, z );
     if (err)
         return err;
 
     /* Now get out on the other side: Bring "normal form" to output form */
 
-    err = geometric_to_orthometric (dstdefn, PJ_FWD, point_count, point_offset, x, y, z);
+    err = geometric_to_orthometric (dst, PJ_FWD, point_count, point_offset, x, y, z);
     if (err)
         return err;
-    err = height_unit (dstdefn, PJ_FWD, point_count, point_offset, z);
+    err = height_unit (dst, PJ_FWD, point_count, point_offset, z);
     if (err)
         return err;
-    err = prime_meridian (dstdefn, PJ_FWD, point_count, point_offset, x);
+    err = prime_meridian (dst, PJ_FWD, point_count, point_offset, x);
     if (err)
         return err;
-    err = geographic_to_cartesian (dstdefn, PJ_FWD, point_count, point_offset, x, y, z);
+    err = geographic_to_cartesian (dst, PJ_FWD, point_count, point_offset, x, y, z);
     if (err)
         return err;
-    err = geographic_to_projected (dstdefn, point_count, point_offset, x, y, z);
+    err = geographic_to_projected (dst, point_count, point_offset, x, y, z);
     if (err)
         return err;
-    err = long_wrap (dstdefn, point_count, point_offset, x);
+    err = long_wrap (dst, point_count, point_offset, x);
     if (err)
         return err;
-    err = adjust_axes (dstdefn, PJ_FWD, point_count, point_offset, x, y, z);
+    err = adjust_axes (dst, PJ_FWD, point_count, point_offset, x, y, z);
     if (err)
         return err;
 
@@ -799,7 +800,7 @@ int pj_geocentric_from_wgs84( PJ *defn,
 /*      coordinates in radians in the destination datum.                */
 /************************************************************************/
 
-int pj_datum_transform( PJ *srcdefn, PJ *dstdefn,
+int pj_datum_transform( PJ *src, PJ *dst,
                         long point_count, int point_offset,
                         double *x, double *y, double *z )
 
@@ -813,21 +814,21 @@ int pj_datum_transform( PJ *srcdefn, PJ *dstdefn,
 /*      (ie. only a +ellps declaration, no +datum).  This is new        */
 /*      behavior for PROJ 4.6.0.                                        */
 /* -------------------------------------------------------------------- */
-    if( srcdefn->datum_type == PJD_UNKNOWN
-        || dstdefn->datum_type == PJD_UNKNOWN )
+    if( src->datum_type == PJD_UNKNOWN
+        || dst->datum_type == PJD_UNKNOWN )
         return 0;
 
 /* -------------------------------------------------------------------- */
 /*      Short cut if the datums are identical.                          */
 /* -------------------------------------------------------------------- */
-    if( pj_compare_datums( srcdefn, dstdefn ) )
+    if( pj_compare_datums( src, dst ) )
         return 0;
 
-    src_a = srcdefn->a_orig;
-    src_es = srcdefn->es_orig;
+    src_a = src->a_orig;
+    src_es = src->es_orig;
 
-    dst_a = dstdefn->a_orig;
-    dst_es = dstdefn->es_orig;
+    dst_a = dst->a_orig;
+    dst_es = dst->es_orig;
 
 /* -------------------------------------------------------------------- */
 /*      Create a temporary Z array if one is not provided.              */
@@ -846,16 +847,16 @@ int pj_datum_transform( PJ *srcdefn, PJ *dstdefn,
 /*	If this datum requires grid shifts, then apply it to geodetic   */
 /*      coordinates.                                                    */
 /* -------------------------------------------------------------------- */
-    if( srcdefn->datum_type == PJD_GRIDSHIFT )
+    if( src->datum_type == PJD_GRIDSHIFT )
     {
-        pj_apply_gridshift_2( srcdefn, 0, point_count, point_offset, x, y, z );
-        CHECK_RETURN(srcdefn);
+        pj_apply_gridshift_2( src, 0, point_count, point_offset, x, y, z );
+        CHECK_RETURN(src);
 
         src_a = SRS_WGS84_SEMIMAJOR;
         src_es = SRS_WGS84_ESQUARED;
     }
 
-    if( dstdefn->datum_type == PJD_GRIDSHIFT )
+    if( dst->datum_type == PJD_GRIDSHIFT )
     {
         dst_a = SRS_WGS84_SEMIMAJOR;
         dst_es = SRS_WGS84_ESQUARED;
@@ -865,52 +866,52 @@ int pj_datum_transform( PJ *srcdefn, PJ *dstdefn,
 /*      Do we need to go through geocentric coordinates?                */
 /* ==================================================================== */
     if( src_es != dst_es || src_a != dst_a
-        || srcdefn->datum_type == PJD_3PARAM
-        || srcdefn->datum_type == PJD_7PARAM
-        || dstdefn->datum_type == PJD_3PARAM
-        || dstdefn->datum_type == PJD_7PARAM)
+        || src->datum_type == PJD_3PARAM
+        || src->datum_type == PJD_7PARAM
+        || dst->datum_type == PJD_3PARAM
+        || dst->datum_type == PJD_7PARAM)
     {
 /* -------------------------------------------------------------------- */
 /*      Convert to geocentric coordinates.                              */
 /* -------------------------------------------------------------------- */
-        srcdefn->ctx->last_errno =
+        src->ctx->last_errno =
             pj_geodetic_to_geocentric( src_a, src_es,
                                        point_count, point_offset, x, y, z );
-        CHECK_RETURN(srcdefn);
+        CHECK_RETURN(src);
 
 /* -------------------------------------------------------------------- */
 /*      Convert between datums.                                         */
 /* -------------------------------------------------------------------- */
-        if( srcdefn->datum_type == PJD_3PARAM
-            || srcdefn->datum_type == PJD_7PARAM )
+        if( src->datum_type == PJD_3PARAM
+            || src->datum_type == PJD_7PARAM )
         {
-            pj_geocentric_to_wgs84( srcdefn, point_count, point_offset,x,y,z);
-            CHECK_RETURN(srcdefn);
+            pj_geocentric_to_wgs84( src, point_count, point_offset,x,y,z);
+            CHECK_RETURN(src);
         }
 
-        if( dstdefn->datum_type == PJD_3PARAM
-            || dstdefn->datum_type == PJD_7PARAM )
+        if( dst->datum_type == PJD_3PARAM
+            || dst->datum_type == PJD_7PARAM )
         {
-            pj_geocentric_from_wgs84( dstdefn, point_count,point_offset,x,y,z);
-            CHECK_RETURN(dstdefn);
+            pj_geocentric_from_wgs84( dst, point_count,point_offset,x,y,z);
+            CHECK_RETURN(dst);
         }
 
 /* -------------------------------------------------------------------- */
 /*      Convert back to geodetic coordinates.                           */
 /* -------------------------------------------------------------------- */
-        dstdefn->ctx->last_errno =
+        dst->ctx->last_errno =
             pj_geocentric_to_geodetic( dst_a, dst_es,
                                        point_count, point_offset, x, y, z );
-        CHECK_RETURN(dstdefn);
+        CHECK_RETURN(dst);
     }
 
 /* -------------------------------------------------------------------- */
 /*      Apply grid shift to destination if required.                    */
 /* -------------------------------------------------------------------- */
-    if( dstdefn->datum_type == PJD_GRIDSHIFT )
+    if( dst->datum_type == PJD_GRIDSHIFT )
     {
-        pj_apply_gridshift_2( dstdefn, 1, point_count, point_offset, x, y, z );
-        CHECK_RETURN(dstdefn);
+        pj_apply_gridshift_2( dst, 1, point_count, point_offset, x, y, z );
+        CHECK_RETURN(dst);
     }
 
     if( z_is_temp )
diff --git a/src/pj_units.c b/src/pj_units.c
index 84b8894e..59ff1ebe 100644
--- a/src/pj_units.c
+++ b/src/pj_units.c
@@ -1,5 +1,7 @@
 /* definition of standard cartesian units */
 
+#include <stddef.h>
+
 #include "proj.h"
 
 #define PJ_UNITS__
diff --git a/src/pj_utils.c b/src/pj_utils.c
index 5a1faf72..3f18cc7e 100644
--- a/src/pj_utils.c
+++ b/src/pj_utils.c
@@ -28,9 +28,10 @@
 
 #define PJ_LIB__
 
-#include "projects.h"
+#include <stdio.h>
 #include <string.h>
-#include <math.h>
+
+#include "projects.h"
 
 /************************************************************************/
 /*                           pj_is_latlong()                            */
diff --git a/src/proj.c b/src/proj.c
index 01b32065..acd34da6 100644
--- a/src/proj.c
+++ b/src/proj.c
@@ -200,6 +200,7 @@ static void vprocess(FILE *fid) {
         emess(1,"binary I/O not available in -V option");
 
     for (;;) {
+        proj_errno_reset(Proj);
         ++emess_dat.File_line;
 
         if (!(s = fgets(line, MAX_LINE, fid)))
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 70110746..bd27cc17 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -26,8 +26,13 @@
  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
-#include <stddef.h>
+
 #include <errno.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
 #include "proj.h"
 #include "proj_internal.h"
 #include "proj_math.h"
@@ -1080,4 +1085,3 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
 
     return factors;
 }
-
diff --git a/src/proj_api.h b/src/proj_api.h
index 5bba5887..36187441 100644
--- a/src/proj_api.h
+++ b/src/proj_api.h
@@ -25,16 +25,16 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
- /*
+/*
   * This version number should be updated with every release!
   *
   * This file is expected to be removed from the PROJ distribution
   * when a few minor-version releases has been made.
   *
   */
- #ifndef PJ_VERSION
- #define PJ_VERSION 500
- #endif
+#ifndef PJ_VERSION
+#define PJ_VERSION 500
+#endif
 
 
 /* If we're not asked for PJ_VERSION only, give them everything */
@@ -43,16 +43,17 @@
 #ifndef PROJ_API_H
 #define PROJ_API_H
 
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-
 /* standard inclusions */
 #include <math.h>
+#include <stddef.h>
 #include <stdlib.h>
 
 
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
 /* pj_init() and similar functions can be used with a non-C locale */
 /* Can be detected too at runtime if the symbol pj_atof exists */
 #define PJ_LOCALE_SAFE 1
diff --git a/src/proj_etmerc.c b/src/proj_etmerc.c
index a4088ad9..457edcb8 100644
--- a/src/proj_etmerc.c
+++ b/src/proj_etmerc.c
@@ -38,10 +38,10 @@
  *
 */
 
-
 #define PJ_LIB__
 
 #include <errno.h>
+
 #include "proj.h"
 #include "projects.h"
 #include "proj_math.h"
@@ -349,4 +349,3 @@ PJ *PROJECTION(utm) {
 
     return setup (P);
 }
-
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 30a4a89e..bd289419 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -31,7 +31,9 @@
 #define _USE_MATH_DEFINES
 #endif
 #endif
+
 #include <math.h>   /* For M_PI */
+#include <stddef.h>
 
 #include "proj.h"
 
diff --git a/src/proj_math.h b/src/proj_math.h
index ca3ab401..f6dd7e09 100644
--- a/src/proj_math.h
+++ b/src/proj_math.h
@@ -58,10 +58,13 @@ int     pj_isnan(double x);
 #define round   pj_round
 #define lround  pj_lround
 
-#ifndef isnan
-#define isnan   pj_isnan
+
+#ifdef isnan
+#undef isnan
 #endif
 
+#define isnan   pj_isnan
+
 #endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */
 
 #ifdef __cplusplus
diff --git a/src/proj_mdist.c b/src/proj_mdist.c
index 02d9ebee..777f704d 100644
--- a/src/proj_mdist.c
+++ b/src/proj_mdist.c
@@ -28,7 +28,12 @@
 ** Precision commensurate with double precision.
 */
 #define PJ_LIB__
+
+#include <math.h>
+#include <stdlib.h>
+
 #include "projects.h"
+
 #define MAX_ITER 20
 #define TOL 1e-14
 
diff --git a/src/proj_rouss.c b/src/proj_rouss.c
index bbf6ff0f..ca4476bb 100644
--- a/src/proj_rouss.c
+++ b/src/proj_rouss.c
@@ -24,7 +24,10 @@
 ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */
 #define PJ_LIB__
+
 #include <errno.h>
+#include <math.h>
+
 #include "proj.h"
 #include "projects.h"
 
@@ -151,4 +154,3 @@ PJ *PROJECTION(rouss) {
 
     return P;
 }
-
diff --git a/src/rtodms.c b/src/rtodms.c
index 77d438b7..f0e2f675 100644
--- a/src/rtodms.c
+++ b/src/rtodms.c
@@ -1,7 +1,12 @@
 /* Convert radian argument to DMS ascii format */
-#include "projects.h"
+
+#include <math.h>
+#include <stddef.h>
 #include <stdio.h>
 #include <string.h>
+
+#include "projects.h"
+
 /*
 ** RES is fractional second figures
 ** RES60 = 60 * RES