Merge branch 'master' into return-nans-quickly

1 files changed, 74 insertions, 125 deletions

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 {