preliminary addition of Icosahedral Snyder Equal Area projection (#111)

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@2165 4e78687f-474d-0410-85f9-8d5e500ac6b2

5 files changed, 1280 insertions, 31 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 17ad7ffb..c57342ae 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -25,7 +25,7 @@ libproj_la_SOURCES = \
 	projects.h pj_list.h \
 	PJ_aeqd.c PJ_gnom.c PJ_laea.c PJ_mod_ster.c \
 	PJ_nsper.c PJ_nzmg.c PJ_ortho.c PJ_stere.c PJ_sterea.c \
-	PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c \
+	PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c PJ_isea.c \
 	PJ_imw_p.c PJ_krovak.c PJ_lcc.c PJ_poly.c \
 	PJ_rpoly.c PJ_sconics.c proj_rouss.c \
 	PJ_cass.c PJ_cc.c PJ_cea.c PJ_eqc.c \
diff --git a/src/Makefile.in b/src/Makefile.in
index 4c7a0dfb..014aa181 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -77,34 +77,35 @@ libproj_la_LIBADD =
 am_libproj_la_OBJECTS = PJ_aeqd.lo PJ_gnom.lo PJ_laea.lo \
 	PJ_mod_ster.lo PJ_nsper.lo PJ_nzmg.lo PJ_ortho.lo PJ_stere.lo \
 	PJ_sterea.lo PJ_aea.lo PJ_bipc.lo PJ_bonne.lo PJ_eqdc.lo \
-	PJ_imw_p.lo PJ_krovak.lo PJ_lcc.lo PJ_poly.lo PJ_rpoly.lo \
-	PJ_sconics.lo proj_rouss.lo PJ_cass.lo PJ_cc.lo PJ_cea.lo \
-	PJ_eqc.lo PJ_gall.lo PJ_labrd.lo PJ_lsat.lo PJ_merc.lo \
-	PJ_mill.lo PJ_ocea.lo PJ_omerc.lo PJ_somerc.lo PJ_tcc.lo \
-	PJ_tcea.lo PJ_tmerc.lo PJ_airy.lo PJ_aitoff.lo PJ_august.lo \
-	PJ_bacon.lo PJ_chamb.lo PJ_hammer.lo PJ_lagrng.lo PJ_larr.lo \
-	PJ_lask.lo PJ_nocol.lo PJ_ob_tran.lo PJ_oea.lo PJ_tpeqd.lo \
-	PJ_vandg.lo PJ_vandg2.lo PJ_vandg4.lo PJ_wag7.lo PJ_lcca.lo \
-	PJ_geos.lo proj_etmerc.lo PJ_boggs.lo PJ_collg.lo PJ_crast.lo \
-	PJ_denoy.lo PJ_eck1.lo PJ_eck2.lo PJ_eck3.lo PJ_eck4.lo \
-	PJ_eck5.lo PJ_fahey.lo PJ_fouc_s.lo PJ_gins8.lo PJ_gstmerc.lo \
-	PJ_gn_sinu.lo PJ_goode.lo PJ_igh.lo PJ_hatano.lo PJ_loxim.lo \
-	PJ_mbt_fps.lo PJ_mbtfpp.lo PJ_mbtfpq.lo PJ_moll.lo PJ_nell.lo \
-	PJ_nell_h.lo PJ_putp2.lo PJ_putp3.lo PJ_putp4p.lo PJ_putp5.lo \
-	PJ_putp6.lo PJ_robin.lo PJ_sts.lo PJ_urm5.lo PJ_urmfps.lo \
-	PJ_wag2.lo PJ_wag3.lo PJ_wink1.lo PJ_wink2.lo pj_latlong.lo \
-	pj_geocent.lo aasincos.lo adjlon.lo bch2bps.lo bchgen.lo \
-	biveval.lo dmstor.lo mk_cheby.lo pj_auth.lo pj_deriv.lo \
-	pj_ell_set.lo pj_ellps.lo pj_errno.lo pj_factors.lo pj_fwd.lo \
-	pj_init.lo pj_inv.lo pj_list.lo pj_malloc.lo pj_mlfn.lo \
-	pj_msfn.lo proj_mdist.lo pj_open_lib.lo pj_param.lo pj_phi2.lo \
-	pj_pr_list.lo pj_qsfn.lo pj_strerrno.lo pj_tsfn.lo pj_units.lo \
-	pj_ctx.lo pj_log.lo pj_zpoly1.lo rtodms.lo vector1.lo \
-	pj_release.lo pj_gauss.lo PJ_healpix.lo nad_cvt.lo nad_init.lo \
-	nad_intr.lo emess.lo pj_apply_gridshift.lo pj_datums.lo \
-	pj_datum_set.lo pj_transform.lo geocent.lo pj_utils.lo \
-	pj_gridinfo.lo pj_gridlist.lo jniproj.lo pj_mutex.lo \
-	pj_initcache.lo pj_apply_vgridshift.lo
+	PJ_isea.lo PJ_imw_p.lo PJ_krovak.lo PJ_lcc.lo PJ_poly.lo \
+	PJ_rpoly.lo PJ_sconics.lo proj_rouss.lo PJ_cass.lo PJ_cc.lo \
+	PJ_cea.lo PJ_eqc.lo PJ_gall.lo PJ_labrd.lo PJ_lsat.lo \
+	PJ_merc.lo PJ_mill.lo PJ_ocea.lo PJ_omerc.lo PJ_somerc.lo \
+	PJ_tcc.lo PJ_tcea.lo PJ_tmerc.lo PJ_airy.lo PJ_aitoff.lo \
+	PJ_august.lo PJ_bacon.lo PJ_chamb.lo PJ_hammer.lo PJ_lagrng.lo \
+	PJ_larr.lo PJ_lask.lo PJ_nocol.lo PJ_ob_tran.lo PJ_oea.lo \
+	PJ_tpeqd.lo PJ_vandg.lo PJ_vandg2.lo PJ_vandg4.lo PJ_wag7.lo \
+	PJ_lcca.lo PJ_geos.lo proj_etmerc.lo PJ_boggs.lo PJ_collg.lo \
+	PJ_crast.lo PJ_denoy.lo PJ_eck1.lo PJ_eck2.lo PJ_eck3.lo \
+	PJ_eck4.lo PJ_eck5.lo PJ_fahey.lo PJ_fouc_s.lo PJ_gins8.lo \
+	PJ_gstmerc.lo PJ_gn_sinu.lo PJ_goode.lo PJ_igh.lo PJ_hatano.lo \
+	PJ_loxim.lo PJ_mbt_fps.lo PJ_mbtfpp.lo PJ_mbtfpq.lo PJ_moll.lo \
+	PJ_nell.lo PJ_nell_h.lo PJ_putp2.lo PJ_putp3.lo PJ_putp4p.lo \
+	PJ_putp5.lo PJ_putp6.lo PJ_robin.lo PJ_sts.lo PJ_urm5.lo \
+	PJ_urmfps.lo PJ_wag2.lo PJ_wag3.lo PJ_wink1.lo PJ_wink2.lo \
+	pj_latlong.lo pj_geocent.lo aasincos.lo adjlon.lo bch2bps.lo \
+	bchgen.lo biveval.lo dmstor.lo mk_cheby.lo pj_auth.lo \
+	pj_deriv.lo pj_ell_set.lo pj_ellps.lo pj_errno.lo \
+	pj_factors.lo pj_fwd.lo pj_init.lo pj_inv.lo pj_list.lo \
+	pj_malloc.lo pj_mlfn.lo pj_msfn.lo proj_mdist.lo \
+	pj_open_lib.lo pj_param.lo pj_phi2.lo pj_pr_list.lo pj_qsfn.lo \
+	pj_strerrno.lo pj_tsfn.lo pj_units.lo pj_ctx.lo pj_log.lo \
+	pj_zpoly1.lo rtodms.lo vector1.lo pj_release.lo pj_gauss.lo \
+	PJ_healpix.lo nad_cvt.lo nad_init.lo nad_intr.lo emess.lo \
+	pj_apply_gridshift.lo pj_datums.lo pj_datum_set.lo \
+	pj_transform.lo geocent.lo pj_utils.lo pj_gridinfo.lo \
+	pj_gridlist.lo jniproj.lo pj_mutex.lo pj_initcache.lo \
+	pj_apply_vgridshift.lo
 libproj_la_OBJECTS = $(am_libproj_la_OBJECTS)
 libproj_la_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
 	$(LIBTOOLFLAGS) --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \
@@ -280,7 +281,7 @@ libproj_la_SOURCES = \
 	projects.h pj_list.h \
 	PJ_aeqd.c PJ_gnom.c PJ_laea.c PJ_mod_ster.c \
 	PJ_nsper.c PJ_nzmg.c PJ_ortho.c PJ_stere.c PJ_sterea.c \
-	PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c \
+	PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c PJ_isea.c \
 	PJ_imw_p.c PJ_krovak.c PJ_lcc.c PJ_poly.c \
 	PJ_rpoly.c PJ_sconics.c proj_rouss.c \
 	PJ_cass.c PJ_cc.c PJ_cea.c PJ_eqc.c \
@@ -500,6 +501,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_healpix.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_igh.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_imw_p.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_isea.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_krovak.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_labrd.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_laea.Plo@am__quote@
diff --git a/src/PJ_isea.c b/src/PJ_isea.c
new file mode 100644
index 00000000..14404cb2
--- /dev/null
+++ b/src/PJ_isea.c
@@ -0,0 +1,1246 @@
+/*
+ * This code was entirely written by Nathan Wagner
+ * and is in the public domain.
+ */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <float.h>
+
+/*
+ * 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
+
+struct hex {
+        int iso;
+        int x, y, z;
+};
+
+/* y *must* be positive down as the xy /iso conversion assumes this */
+ISEA_STATIC
+int hex_xy(struct hex *h) {
+	if (!h->iso) return 1;
+	if (h->x >= 0) {
+		h->y = -h->y - (h->x+1)/2;
+	} else {
+		/* need to round toward -inf, not toward zero, so x-1 */
+		h->y = -h->y - h->x/2;
+	}
+	h->iso = 0;
+
+	return 1;
+}
+
+ISEA_STATIC
+int hex_iso(struct hex *h) {
+	if (h->iso) return 1;
+
+	if (h->x >= 0) {
+		h->y = (-h->y - (h->x+1)/2);
+	} else {
+		/* need to round toward -inf, not toward zero, so x-1 */
+		h->y = (-h->y - (h->x)/2);
+	}
+
+	h->z = -h->x - h->y;
+	h->iso = 1;
+	return 1;
+}
+
+ISEA_STATIC
+int hexbin2(int horizontal, double width, double x, double y,
+                int *i, int *j) {
+	double z, rx, ry, rz;
+	double abs_dx, abs_dy, abs_dz;
+	int ix, iy, iz, s;
+	struct hex h;
+
+	x = x / cos(30 * M_PI / 180.0); /* rotated X coord */
+	y = y - x / 2.0; /* adjustment for rotated X */
+
+	/* adjust for actual hexwidth */
+	x /= width;
+	y /= width;
+
+	z = -x - y;
+
+	ix = rx = floor(x + 0.5);
+	iy = ry = floor(y + 0.5);
+	iz = rz = floor(z + 0.5);
+
+	s = ix + iy + iz;
+
+	if (s) {
+		abs_dx = fabs(rx - x);
+		abs_dy = fabs(ry - y);
+		abs_dz = fabs(rz - z);
+
+		if (abs_dx >= abs_dy && abs_dx >= abs_dz) {
+			ix -= s;
+		} else if (abs_dy >= abs_dx && abs_dy >= abs_dz) {
+			iy -= s;
+		} else {
+			iz -= s;
+		}
+	}
+	h.x = ix;
+	h.y = iy;
+	h.z = iz;
+	h.iso = 1;
+
+	hex_xy(&h);
+	*i = h.x;
+	*j = h.y;
+        return ix * 100 + iy;
+}
+#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 };
+enum isea_address_form { ISEA_GEO, ISEA_Q2DI, ISEA_SEQNUM, ISEA_INTERLEAVE,
+	ISEA_PLANE, ISEA_Q2DD, ISEA_PROJTRI, ISEA_VERTEX2DD, ISEA_HEX
+};
+
+struct isea_dgg {
+	int	polyhedron; /* ignored, icosahedron */
+	double	o_lat, o_lon, o_az; /* orientation, radians */
+	int	pole; /* true if standard snyder */
+	int	topology; /* ignored, hexagon */
+	int	aperture; /* valid values depend on partitioning method */
+	int	resolution;
+	double	radius; /* radius of the earth in meters, ignored 1.0 */
+	int	output; /* an isea_address_form */
+	int	triangle; /* triangle of last transformed point */
+	int	quad; /* quad of last transformed point */
+	unsigned long serial;
+};
+
+struct isea_pt {
+	double x, y;
+};
+
+struct isea_geo {
+	double lon, lat;
+};
+
+struct isea_address {
+	int	type; /* enum isea_address_form */
+	int	number;
+	double	x,y; /* or i,j or lon,lat depending on type */
+};
+
+/* ENDINC */
+
+enum snyder_polyhedron {
+	SNYDER_POLY_HEXAGON, SNYDER_POLY_PENTAGON,
+	SNYDER_POLY_TETRAHEDRON, SNYDER_POLY_CUBE,
+	SNYDER_POLY_OCTAHEDRON, SNYDER_POLY_DODECAHEDRON,
+	SNYDER_POLY_ICOSAHEDRON
+};
+
+struct snyder_constants {
+	double          g, G, theta, ea_w, ea_a, ea_b, g_w, g_a, g_b;
+};
+
+/* TODO put these in radians to avoid a later conversion */
+ISEA_STATIC
+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},
+	{37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0},
+};
+
+#define E 52.62263186
+#define F 10.81231696
+
+#define DEG60 1.04719755119659774614
+#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[] = {
+	{0.0, DEG90},
+	{DEG180, V_LAT},
+	{-DEG108, V_LAT},
+	{-DEG36, V_LAT},
+	{DEG36, V_LAT},
+	{DEG108, V_LAT},
+	{-DEG144, -V_LAT},
+	{-DEG72, -V_LAT},
+	{0.0, -V_LAT},
+	{DEG72, -V_LAT},
+	{DEG144, -V_LAT},
+	{0.0, -DEG90}
+};
+
+/* 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
+
+/* triangle Centers */
+struct isea_geo icostriangles[] = {
+	{0.0, 0.0},
+	{-DEG144, E_RAD},
+	{-DEG72, E_RAD},
+	{0.0, E_RAD},
+	{DEG72, E_RAD},
+	{DEG144, E_RAD},
+	{-DEG144, F_RAD},
+	{-DEG72, F_RAD},
+	{0.0, F_RAD},
+	{DEG72, F_RAD},
+	{DEG144, F_RAD},
+	{-DEG108, -F_RAD},
+	{-DEG36, -F_RAD},
+	{DEG36, -F_RAD},
+	{DEG108, -F_RAD},
+	{DEG180, -F_RAD},
+	{-DEG108, -E_RAD},
+	{-DEG36, -E_RAD},
+	{DEG36, -E_RAD},
+	{DEG108, -E_RAD},
+	{DEG180, -E_RAD},
+};
+
+static double
+az_adjustment(int triangle)
+{
+	double          adj;
+
+	struct isea_geo v;
+	struct isea_geo c;
+
+	v = vertex[tri_v1[triangle]];
+	c = icostriangles[triangle];
+
+	/* TODO looks like the adjustment is always either 0 or 180 */
+	/* at least if you pick your vertex carefully */
+	adj = atan2(cos(v.lat) * sin(v.lon - c.lon),
+		    cos(c.lat) * sin(v.lat)
+		    - sin(c.lat) * cos(v.lat) * cos(v.lon - c.lon));
+	return adj;
+}
+
+/* R tan(g) sin(60) */
+#define TABLE_G 0.6615845383
+
+/* H = 0.25 R tan g = */
+#define TABLE_H 0.1909830056
+
+#define RPRIME 0.91038328153090290025
+
+ISEA_STATIC
+struct isea_pt
+isea_triangle_xy(int triangle)
+{
+	struct isea_pt  c;
+	double Rprime = 0.91038328153090290025;
+
+	triangle = (triangle - 1) % 20;
+
+	c.x = TABLE_G * ((triangle % 5) - 2) * 2.0;
+	if (triangle > 9) {
+		c.x += TABLE_G;
+	}
+	switch (triangle / 5) {
+	case 0:
+		c.y = 5.0 * TABLE_H;
+		break;
+	case 1:
+		c.y = TABLE_H;
+		break;
+	case 2:
+		c.y = -TABLE_H;
+		break;
+	case 3:
+		c.y = -5.0 * TABLE_H;
+		break;
+	default:
+		/* should be impossible */
+		exit(EXIT_FAILURE);
+	};
+	c.x *= Rprime;
+	c.y *= Rprime;
+
+	return c;
+}
+
+/* snyder eq 14 */
+static double
+sph_azimuth(double f_lon, double f_lat, double t_lon, double t_lat)
+{
+	double          az;
+
+	az = atan2(cos(t_lat) * sin(t_lon - f_lon),
+		   cos(f_lat) * sin(t_lat)
+		   - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon)
+		);
+	return az;
+}
+
+/* coord needs to be in radians */
+ISEA_STATIC
+int
+isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)
+{
+	int             i;
+
+	/*
+	 * spherical distance from center of polygon face to any of its
+	 * vertexes on the globe
+	 */
+	double          g;
+
+	/*
+	 * spherical angle between radius vector to center and adjacent edge
+	 * of spherical polygon on the globe
+	 */
+	double          G;
+
+	/*
+	 * plane angle between radius vector to center and adjacent edge of
+	 * plane polygon
+	 */
+	double          theta;
+
+	/* additional variables from snyder */
+	double          q, Rprime, H, Ag, Azprime, Az, dprime, f, rho,
+	                x, y;
+
+	/* variables used to store intermediate results */
+	double          cot_theta, tan_g, az_offset;
+
+	/* how many multiples of 60 degrees we adjust the azimuth */
+	int             Az_adjust_multiples;
+
+	struct snyder_constants c;
+
+	/*
+	 * TODO by locality of reference, start by trying the same triangle
+	 * as last time
+	 */
+
+	/* 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;
+
+	for (i = 1; i <= 20; i++) {
+		double          z;
+		struct isea_geo center;
+
+		center = icostriangles[i];
+
+		/* step 1 */
+#if 0
+		z = sph_distance(center.lon, center.lat, ll->lon, ll->lat);
+#else
+		z = acos(sin(center.lat) * sin(ll->lat)
+			 + cos(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon));
+#endif
+
+		/* not on this triangle */
+		if (z > g + 0.000005) { /* TODO DBL_EPSILON */
+			continue;
+		}
+		Az = sph_azimuth(ll->lon, ll->lat, center.lon, center.lat);
+
+		Az = atan2(cos(ll->lat) * sin(ll->lon - center.lon),
+			   cos(center.lat) * sin(ll->lat)
+		- sin(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon)
+			);
+
+		/* step 2 */
+
+		/* This calculates "some" vertex coordinate */
+		az_offset = az_adjustment(i);
+
+		Az -= az_offset;
+
+		/* TODO I don't know why we do this.  It's not in snyder */
+		/* maybe because we should have picked a better vertex */
+		if (Az < 0.0) {
+			Az += 2.0 * M_PI;
+		}
+		/*
+		 * adjust Az for the point to fall within the range of 0 to
+		 * 2(90 - theta) or 60 degrees for the hexagon, by
+		 * and therefore 120 degrees for the triangle
+		 * of the icosahedron
+		 * subtracting or adding multiples of 60 degrees to Az and
+		 * recording the amount of adjustment
+		 */
+
+		Az_adjust_multiples = 0;
+		while (Az < 0.0) {
+			Az += DEG120;
+			Az_adjust_multiples--;
+		}
+		while (Az > DEG120 + DBL_EPSILON) {
+			Az -= DEG120;
+			Az_adjust_multiples++;
+		}
+
+		/* step 3 */
+		cot_theta = 1.0 / tan(theta);
+		tan_g = tan(g);	/* TODO this is a constant */
+
+		/* Calculate q from eq 9. */
+		/* TODO cot_theta is cot(30) */
+		q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta);
+
+		/* not in this triangle */
+		if (z > q + 0.000005) {
+			continue;
+		}
+		/* step 4 */
+
+		/* Apply equations 5-8 and 10-12 in order */
+
+		/* eq 5 */
+		/* Rprime = 0.9449322893 * R; */
+		/* R' in the paper is for the truncated */
+		Rprime = 0.91038328153090290025;
+
+		/* eq 6 */
+		H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G));
+
+		/* eq 7 */
+		/* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */
+		Ag = Az + G + H - DEG180;
+
+		/* eq 8 */
+		Azprime = atan2(2.0 * Ag, Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta);
+
+		/* eq 10 */
+		/* cot(theta) = 1.73205080756887729355 */
+		dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta);
+
+		/* eq 11 */
+		f = dprime / (2.0 * Rprime * sin(q / 2.0));
+
+		/* eq 12 */
+		rho = 2.0 * Rprime * f * sin(z / 2.0);
+
+		/*
+		 * add back the same 60 degree multiple adjustment from step
+		 * 2 to Azprime
+		 */
+
+		Azprime += DEG120 * Az_adjust_multiples;
+
+		/* calculate rectangular coordinates */
+
+		x = rho * sin(Azprime);
+		y = rho * cos(Azprime);
+
+		/*
+		 * TODO
+		 * translate coordinates to the origin for the particular
+		 * hexagon on the flattened polyhedral map plot
+		 */
+
+		out->x = x;
+		out->y = y;
+
+		return i;
+	}
+
+	/*
+	 * should be impossible, this implies that the coordinate is not on
+	 * any triangle
+	 */
+
+	fprintf(stderr, "impossible transform: %f %f is not on any triangle\n",
+		ll->lon * RAD2DEG, ll->lat * RAD2DEG);
+
+	exit(EXIT_FAILURE);
+
+	/* not reached */
+	return 0;		/* supresses a warning */
+}
+
+/*
+ * return the new coordinates of any point in orginal coordinate system.
+ * Define a point (newNPold) in orginal coordinate system as the North Pole in
+ * new coordinate system, and the great circle connect the original and new
+ * North Pole as the lon0 longitude in new coordinate system, given any point
+ * in orginal coordinate system, this function return the new coordinates.
+ */
+
+#define PRECISION 0.0000000000005
+
+/*
+ * TODO rename to isea_coordtrans and take an struct pointer to the result
+ * struct TODO remove function entirely.  we use the synder one. the only
+ * difference is a 180 degree longitude rotation
+ */
+ISEA_STATIC
+struct isea_geo
+coordtrans(struct isea_geo * newNPold, struct isea_geo * ptold, double lon0)
+{
+	double          cosptnewlat, cosptnewlon;
+	struct isea_geo ptnew;
+
+	cosptnewlat = sin(newNPold->lat) * sin(ptold->lat) +
+		cos(newNPold->lat) * cos(ptold->lat) * cos(newNPold->lon - ptold->lon);
+
+	if (cosptnewlat > 1.)
+		cosptnewlat = 1.0;
+	if (cosptnewlat < -1.)
+		cosptnewlat = -1.0;
+
+	ptnew.lat = acos(cosptnewlat);
+
+	if (fabs(ptnew.lat - 0.0) < PRECISION * 100000)
+		ptnew.lon = 0.0;
+	else if (fabs(ptnew.lat - M_PI) < PRECISION * 100000)
+		ptnew.lon = 0.0;
+	else {
+		cosptnewlon = (sin(ptold->lat) * cos(newNPold->lat) - cos(ptold->lat) *
+			       sin(newNPold->lat) * cos(newNPold->lon - ptold->lon)) / sin(ptnew.lat);
+		if (cosptnewlon > 1.)
+			cosptnewlon = 1.0;
+		if (cosptnewlon < -1.)
+			cosptnewlon = -1.0;
+
+		ptnew.lon = acos(cosptnewlon);
+
+		if ((ptold->lon - newNPold->lon) >= 0 && (ptold->lon - newNPold->lon) < 180)
+			ptnew.lon = -ptnew.lon + lon0;
+
+		else
+			ptnew.lon = ptnew.lon + lon0;
+
+		if (ptnew.lon > M_PI)
+			ptnew.lon -= 2 * M_PI;
+		if (ptnew.lon < -M_PI)
+			ptnew.lon += 2 * M_PI;
+	}
+	ptnew.lat = M_PI / 2 - ptnew.lat;
+	return ptnew;
+}
+
+/* formula from Snyder, Map Projections: A working manual, p31 */
+/*
+ * old north pole at np in new coordinates
+ * could be simplified a bit with fewer intermediates
+ *
+ * TODO take a result pointer
+ */
+ISEA_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;
+	double          sin_phip;
+	double          lp_b;	/* lambda prime minus beta */
+	double          cos_p, sin_a;
+
+	phi = pt->lat;
+	lambda = pt->lon;
+	alpha = np->lat;
+	beta = np->lon;
+	lambda0 = beta;
+
+	cos_p = cos(phi);
+	sin_a = sin(alpha);
+
+	/* mpawm 5-7 */
+	sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0);
+
+	/* mpawm 5-8b */
+
+	/* use the two argument form so we end up in the right quadrant */
+	lp_b = atan2(cos_p * sin(lambda - lambda0),
+	   (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));
+
+	lambdap = lp_b + beta;
+
+	/* normalize longitude */
+	/* TODO can we just do a modulus ? */
+	lambdap = fmod(lambdap, 2 * M_PI);
+	while (lambdap > M_PI)
+		lambdap -= 2 * M_PI;
+	while (lambdap < -M_PI)
+		lambdap += 2 * M_PI;
+
+	phip = asin(sin_phip);
+
+	npt.lat = phip;
+	npt.lon = lambdap;
+
+	return npt;
+}
+
+ISEA_STATIC
+struct isea_geo
+isea_ctran(struct isea_geo * np, struct isea_geo * pt, double lon0)
+{
+	struct isea_geo npt;
+
+	np->lon += M_PI;
+	npt = snyder_ctran(np, pt);
+	np->lon -= M_PI;
+
+	npt.lon -= (M_PI - lon0 + np->lon);
+
+	/*
+	 * snyder is down tri 3, isea is along side of tri1 from vertex 0 to
+	 * vertex 1 these are 180 degrees apart
+	 */
+	npt.lon += M_PI;
+	/* normalize longitude */
+	npt.lon = fmod(npt.lon, 2 * M_PI);
+	while (npt.lon > M_PI)
+		npt.lon -= 2 * M_PI;
+	while (npt.lon < -M_PI)
+		npt.lon += 2 * M_PI;
+
+	return npt;
+}
+
+/* in radians */
+#define ISEA_STD_LAT 1.01722196792335072101
+#define ISEA_STD_LON .19634954084936207740
+
+ISEA_STATIC
+struct isea_dgg isea_standard_dgg = {
+	20, 58.28252559, 11.25, 0.0,
+	6, 4, 0, 6, 1.0
+};
+
+/* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */
+
+ISEA_STATIC
+int
+isea_grid_init(struct isea_dgg * g)
+{
+	if (!g)
+		return 0;
+
+	g->polyhedron = 20;
+	g->o_lat = ISEA_STD_LAT;
+	g->o_lon = ISEA_STD_LON;
+	g->o_az = 0.0;
+	g->aperture = 4;
+	g->resolution = 6;
+	g->radius = 1.0;
+	g->topology = 6;
+
+	return 1;
+}
+
+ISEA_STATIC
+int
+isea_orient_isea(struct isea_dgg * g)
+{
+	if (!g)
+		return 0;
+	g->o_lat = ISEA_STD_LAT;
+	g->o_lon = ISEA_STD_LON;
+	g->o_az = 0.0;
+	return 1;
+}
+
+ISEA_STATIC
+int
+isea_orient_pole(struct isea_dgg * g)
+{
+	if (!g)
+		return 0;
+	g->o_lat = M_PI / 2.0;
+	g->o_lon = 0.0;
+	g->o_az = 0;
+	return 1;
+}
+
+ISEA_STATIC
+int
+isea_orient_dymax(struct isea_dgg * g)
+{
+	if (!g)
+		return 0;
+	g->o_lat = 2.300882 * M_PI / 180.0;
+	g->o_lon = -5.24539 * M_PI / 180.0;
+	g->o_az = 7.46658 * M_PI / 180.0;
+	return 1;
+}
+
+ISEA_STATIC
+int
+isea_transform(struct isea_dgg * g, struct isea_geo * in,
+	       struct isea_pt * out)
+{
+	struct isea_geo i, pole;
+	int             tri;
+
+	pole.lat = g->o_lat;
+	pole.lon = g->o_lon;
+
+	i = isea_ctran(&pole, in, g->o_az);
+
+	tri = isea_snyder_forward(&i, out);
+	out->x *= g->radius;
+	out->y *= g->radius;
+	g->triangle = tri;
+
+	return tri;
+}
+
+#define DOWNTRI(tri) (((tri - 1) / 5) % 2 == 1)
+
+/* get which triangle a point is in, assumes radius = 1.0 */
+ISEA_STATIC
+int isea_xy_triangle(struct isea_pt *p, double radius) {
+	double dist, test;
+	int tri, closest;
+	struct isea_pt tc, pt;
+	pt.x = p->x / radius;
+	pt.y = p->y / radius;
+
+	tc = isea_triangle_xy(1);
+	dist = (pt.x - tc.x) * (pt.x - tc.x) + (pt.y - tc.y) * (pt.y - tc.y);
+	closest = 1;
+
+	/* TODO just calculate it directly, no need to actually do all this
+	 * work
+	 */
+	for (tri = 2; tri <= 20; tri++) {
+		tc = isea_triangle_xy(tri);
+		test = (pt.x - tc.x) * (pt.x - tc.x) + (pt.y - tc.y) * (pt.y - tc.y);
+		if (test < dist) {
+			dist = test;
+			closest = tri;
+		}
+	}
+	return closest;
+}
+
+ISEA_STATIC
+int isea_projtri(struct isea_pt *xy, double radius) {
+	int tri;
+
+	tri = isea_xy_triangle(xy, radius);
+
+	xy->x /= radius;
+	xy->y /= radius;
+	xy->x += 0.5;
+	xy->y += 2.0 * .14433756729740644112;
+
+	return tri;
+}
+
+ISEA_STATIC
+void
+isea_rotate(struct isea_pt * pt, double degrees)
+{
+	double          rad;
+
+	double          x, y;
+
+	rad = -degrees * M_PI / 180.0;
+	while (rad >= 2.0 * M_PI) rad -= 2.0 * M_PI;
+	while (rad <= -2.0 * M_PI) rad += 2.0 * M_PI;
+
+	x = pt->x * cos(rad) + pt->y * sin(rad);
+	y = -pt->x * sin(rad) + pt->y * cos(rad);
+
+	pt->x = x;
+	pt->y = y;
+}
+
+ISEA_STATIC
+int isea_tri_plane(int tri, struct isea_pt *pt, double radius) {
+	struct isea_pt tc; /* center of triangle */
+
+	if (DOWNTRI(tri)) {
+		isea_rotate(pt, 180.0);
+	}
+	tc = isea_triangle_xy(tri);
+	tc.x *= radius;
+	tc.y *= radius;
+	pt->x += tc.x;
+	pt->y += tc.y;
+
+	return tri;
+}
+
+ISEA_STATIC
+int isea_projection(struct isea_dgg *g, double xin, double yin, double *xout,
+		double *yout) {
+	return 1;
+}
+
+/* convert projected triangle coords to quad xy coords, return quad number */
+ISEA_STATIC
+int
+isea_ptdd(int tri, struct isea_pt *pt) {
+	int             downtri, quad;
+
+	downtri = (((tri - 1) / 5) % 2 == 1);
+	quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
+
+	isea_rotate(pt, downtri ? 240.0 : 60.0);
+	if (downtri) {
+		pt->x += 0.5;
+		/* pt->y += cos(30.0 * M_PI / 180.0); */
+		pt->y += .86602540378443864672;
+	}
+	return quad;
+}
+
+ISEA_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;
+	double          sidelength;	/* in hexes */
+	int             d, i;
+	int             maxcoord;
+	struct hex      h;
+
+	/* This is the number of hexes from apex to base of a triangle */
+	sidelength = (pow(2.0, g->resolution) + 1.0) / 2.0;
+
+	/* apex to base is cos(30deg) */
+	hexwidth = cos(M_PI / 6.0) / sidelength;
+
+	/* TODO I think sidelength is always x.5, so
+	 * (int)sidelength * 2 + 1 might be just as good
+	 */
+	maxcoord = (int) (sidelength * 2.0 + 0.5);
+
+	v = *pt;
+	hexbin2(0, hexwidth, v.x, v.y, &h.x, &h.y);
+	h.iso = 0;
+	hex_iso(&h);
+
+	d = h.x - h.z;
+	i = h.x + h.y + h.y;
+	
+	/*
+	 * you want to test for max coords for the next quad in the same
+	 * "row" first to get the case where both are max
+	 */
+	if (quad <= 5) {
+		if (d == 0 && i == maxcoord) {
+			/* north pole */
+			quad = 0;
+			d = 0;
+			i = 0;
+		} else if (i == maxcoord) {
+			/* upper right in next quad */
+			quad += 1;
+			if (quad == 6)
+				quad = 1;
+			i = maxcoord - d;
+			d = 0;
+		} else if (d == maxcoord) {
+			/* lower right in quad to lower right */
+			quad += 5;
+			d = 0;
+		}
+	} else if (quad >= 6) {
+		if (i == 0 && d == maxcoord) {
+			/* south pole */
+			quad = 11;
+			d = 0;
+			i = 0;
+		} else if (d == maxcoord) {
+			/* lower right in next quad */
+			quad += 1;
+			if (quad == 11)
+				quad = 6;
+			d = maxcoord - i;
+			i = 0;
+		} else if (i == maxcoord) {
+			/* upper right in quad to upper right */
+			quad = (quad - 4) % 5;
+			i = 0;
+		}
+	}
+
+	di->x = d;
+	di->y = i;
+
+	g->quad = quad;
+	return quad;
+}
+
+ISEA_STATIC
+int
+isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di) {
+	struct isea_pt  v;
+	double          hexwidth;
+	int             sidelength;	/* in hexes */
+	struct hex      h;
+
+	if (g->aperture == 3 && g->resolution % 2 != 0) {
+		return isea_dddi_ap3odd(g, quad, pt, di);
+	}
+	/* todo might want to do this as an iterated loop */
+	if (g->aperture >0) {
+		sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
+	} else {
+		sidelength = g->resolution;
+	}
+
+	hexwidth = 1.0 / sidelength;
+
+	v = *pt;
+	isea_rotate(&v, -30.0);
+	hexbin2(0, hexwidth, v.x, v.y, &h.x, &h.y);
+	h.iso = 0;
+	hex_iso(&h);
+
+	/* we may actually be on another quad */
+	if (quad <= 5) {
+		if (h.x == 0 && h.z == -sidelength) {
+			/* north pole */
+			quad = 0;
+			h.z = 0;
+			h.y = 0;
+			h.x = 0;
+		} else if (h.z == -sidelength) {
+			quad = quad + 1;
+			if (quad == 6)
+				quad = 1;
+			h.y = sidelength - h.x;
+			h.z = h.x - sidelength;
+			h.x = 0;
+		} else if (h.x == sidelength) {
+			quad += 5;
+			h.y = -h.z;
+			h.x = 0;
+		}
+	} else if (quad >= 6) {
+		if (h.z == 0 && h.x == sidelength) {
+			/* south pole */
+			quad = 11;
+			h.x = 0;
+			h.y = 0;
+			h.z = 0;
+		} else if (h.x == sidelength) {
+			quad = quad + 1;
+			if (quad == 11)
+				quad = 6;
+			h.x = h.y + sidelength;
+			h.y = 0;
+			h.z = -h.x;
+		} else if (h.y == -sidelength) {
+			quad -= 4;
+			h.y = 0;
+			h.z = -h.x;
+		}
+	}
+	di->x = h.x;
+	di->y = -h.z;
+
+	g->quad = quad;
+	return quad;
+}
+
+ISEA_STATIC
+int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,
+	       	struct isea_pt *di) {
+	struct isea_pt  v;
+	int             quad;
+
+	v = *pt;
+	quad = isea_ptdd(tri, &v);
+	quad = isea_dddi(g, quad, &v, di);
+	return quad;
+}
+
+/* q2di to seqnum */
+ISEA_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;
+		return g->serial;
+	}
+	/* hexes in a quad */
+	hexes = (int) (pow(g->aperture, g->resolution) + 0.5);
+	if (quad == 11) {
+		g->serial = 1 + 10 * hexes + 1;
+		return g->serial;
+	}
+	if (g->aperture == 3 && g->resolution % 2 == 1) {
+		height = (int) (pow(g->aperture, (g->resolution - 1) / 2.0));
+		sn = ((int) di->x) * height;
+		sn += ((int) di->y) / height;
+		sn += (quad - 1) * hexes;
+		sn += 2;
+	} else {
+		sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
+		sn = (quad - 1) * hexes + sidelength * di->x + di->y + 2;
+	}
+
+	g->serial = sn;
+	return sn;
+}
+
+/* TODO just encode the quad in the d or i coordinate
+ * quad is 0-11, which can be four bits.
+ * 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) {
+	struct isea_pt v;
+	int sidelength;
+	int d, i, x, y, quad;
+	double oddsl;
+
+	quad = isea_ptdi(g, tri, pt, &v);
+
+	hex->x = ((int)v.x << 4) + quad;
+	hex->y = v.y;
+
+	return 1;
+
+	d = v.x;
+	i = v.y;
+
+	/* Aperture 3 odd resolutions */
+	if (g->aperture == 3 && g->resolution % 2 != 0) {
+		int offset = (int)(pow(3.0, g->resolution - 1) + 0.5);
+
+		oddsl = (pow(2.0, g->resolution) + 1.0) / 2.0;
+
+		d += offset * ((g->quad-1) % 5);
+		i += offset * ((g->quad-1) % 5);
+
+		if (quad == 0) {
+			d = 0;
+			i = offset;
+		} else if (quad == 11) {
+			d = 2 * offset;
+			i = 0;
+		} else if (quad > 5) {
+			d += offset;
+		}
+
+		x = (2*d - i) /3;
+		y = (2*i - d) /3;
+
+		hex->x = x + offset / 3;
+		hex->y = y + 2 * offset / 3;
+		return 1;
+	}
+
+	/* aperture 3 even resolutions and aperture 4 */
+	sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
+	if (g->quad == 0) {
+		hex->x = 0;
+		hex->y = sidelength;
+	} else if (g->quad == 11) {
+		hex->x = sidelength * 2;
+		hex->y = 0;
+	} else {
+		hex->x = d + sidelength * ((g->quad-1) % 5);
+		if (g->quad > 5) hex->x += sidelength;
+		hex->y = i + sidelength * ((g->quad-1) % 5);
+	}
+
+	return 1;
+}
+
+ISEA_STATIC
+struct isea_pt
+isea_forward(struct isea_dgg *g, struct isea_geo *in)
+{
+	int             tri, downtri, quad;
+	struct isea_pt  out, coord;
+
+	tri = isea_transform(g, in, &out);
+
+	downtri = (((tri - 1) / 5) % 2 == 1);
+	quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
+
+	if (g->output == ISEA_PLANE) {
+		isea_tri_plane(tri, &out, g->radius);
+		return out;
+	}
+
+	/* convert to isea standard triangle size */
+	out.x = out.x / g->radius * ISEA_SCALE;
+	out.y = out.y / g->radius * ISEA_SCALE;
+	out.x += 0.5;
+	out.y += 2.0 * .14433756729740644112;
+
+	switch (g->output) {
+	case ISEA_PROJTRI:
+		/* nothing to do, already in projected triangle */
+		break;
+	case ISEA_VERTEX2DD:
+		g->quad = isea_ptdd(tri, &out);
+		break;
+	case ISEA_Q2DD:
+		/* Same as above, we just don't print as much */
+		g->quad = isea_ptdd(tri, &out);
+		break;
+	case ISEA_Q2DI:
+		g->quad = isea_ptdi(g, tri, &out, &coord);
+		return coord;
+		break;
+	case ISEA_SEQNUM:
+		isea_ptdi(g, tri, &out, &coord);
+		/* disn will set g->serial */
+		isea_disn(g, g->quad, &coord);
+		return coord;
+		break;
+	case ISEA_HEX:
+		isea_hex(g, tri, &out, &coord);
+		return coord;
+		break;
+	}
+
+	return out;
+}
+/*
+ * Proj 4 integration code follows
+ */
+
+#define PROJ_PARMS__ \
+	struct isea_dgg dgg;
+
+#define PJ_LIB__
+#include <projects.h>
+
+PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph";
+
+FORWARD(s_forward);
+	struct isea_pt out;
+	struct isea_geo in;
+
+	in.lon = lp.lam;
+	in.lat = lp.phi;
+
+	out = isea_forward(&P->dgg, &in);
+	
+	xy.x = out.x;
+	xy.y = out.y;
+
+	return xy;
+}
+FREEUP; if (P) pj_dalloc(P); }
+
+ENTRY0(isea)
+	char *opt;
+
+        P->fwd = s_forward;
+        isea_grid_init(&P->dgg);
+
+        P->dgg.output = ISEA_PLANE;
+/*		P->dgg.radius = P->a; /* otherwise defaults to 1 */
+	/* calling library will scale, I think */
+
+	opt = pj_param(P->ctx,P->params, "sorient").s;
+	if (opt) {
+		if (!strcmp(opt, "isea")) {
+			isea_orient_isea(&P->dgg);
+		} else if (!strcmp(opt, "pole")) {
+			isea_orient_pole(&P->dgg);
+		} else {
+			E_ERROR(-34);
+		}
+	}
+
+	if (pj_param(P->ctx,P->params, "tazi").i) {
+		P->dgg.o_az = pj_param(P->ctx,P->params, "razi").f;
+	}
+
+	if (pj_param(P->ctx,P->params, "tlon_0").i) {
+		P->dgg.o_lon = pj_param(P->ctx,P->params, "rlon_0").f;
+	}
+
+	if (pj_param(P->ctx,P->params, "tlat_0").i) {
+		P->dgg.o_lat = pj_param(P->ctx,P->params, "rlat_0").f;
+	}
+
+	if (pj_param(P->ctx,P->params, "taperture").i) {
+		P->dgg.aperture = pj_param(P->ctx,P->params, "iaperture").i;
+	}
+
+	if (pj_param(P->ctx,P->params, "tresolution").i) {
+		P->dgg.resolution = pj_param(P->ctx,P->params, "iresolution").i;
+	}
+
+	opt = pj_param(P->ctx,P->params, "smode").s;
+	if (opt) {
+		if (!strcmp(opt, "plane")) {
+			P->dgg.output = ISEA_PLANE;
+		} else if (!strcmp(opt, "di")) {
+			P->dgg.output = ISEA_Q2DI;
+		}
+		else if (!strcmp(opt, "dd")) {
+			P->dgg.output = ISEA_Q2DD;
+		}
+		else if (!strcmp(opt, "hex")) {
+			P->dgg.output = ISEA_HEX;
+		}
+		else {
+			/* TODO verify error code.  Possibly eliminate magic */
+			E_ERROR(-34);
+		}
+	}
+
+	if (pj_param(P->ctx,P->params, "trescale").i) {
+		P->dgg.radius = ISEA_SCALE;
+	}
+
+	if (pj_param(P->ctx,P->params, "tresolution").i) {
+		P->dgg.resolution = pj_param(P->ctx,P->params, "iresolution").i;
+	} else {
+		P->dgg.resolution = 4;
+	}
+
+	if (pj_param(P->ctx,P->params, "taperture").i) {
+		P->dgg.aperture = pj_param(P->ctx,P->params, "iaperture").i;
+	} else {
+		P->dgg.aperture = 3;
+	}
+
+ENDENTRY(P)
diff --git a/src/makefile.vc b/src/makefile.vc
index 94a4d484..220468ab 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -38,7 +38,7 @@ pseudo = \
 	PJ_nell.obj PJ_nell_h.obj PJ_putp2.obj PJ_putp3.obj \
 	PJ_putp4p.obj PJ_putp5.obj PJ_putp6.obj PJ_robin.obj \
 	PJ_sts.obj PJ_urm5.obj PJ_urmfps.obj PJ_wag2.obj \
-	PJ_wag3.obj PJ_wink1.obj PJ_wink2.obj
+	PJ_wag3.obj PJ_wink1.obj PJ_wink2.obj PJ_isea.obj 
 
 support = \
 	aasincos.obj adjlon.obj bch2bps.obj bchgen.obj pj_gauss.obj \
diff --git a/src/pj_list.h b/src/pj_list.h
index 89e9e7aa..ac88f9d8 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -51,6 +51,7 @@ PROJ_HEAD(healpix, "HEALPix")
 PROJ_HEAD(rhealpix, "rHEALPix")
 PROJ_HEAD(igh,  "Interrupted Goode Homolosine")
 PROJ_HEAD(imw_p, "Internation Map of the World Polyconic")
+PROJ_HEAD(isea, "Icosahedral Snyder Equal Area")
 PROJ_HEAD(kav5, "Kavraisky V")
 PROJ_HEAD(kav7, "Kavraisky VII")
 PROJ_HEAD(krovak, "Krovak")