/*
 * 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>

#ifndef M_PI
#  define M_PI 3.14159265358979323846
#endif

/*
 * 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(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, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 0.0, 0.0, 0.0, 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 */
		z = acos(sin(center.lat) * sin(ll->lat)
			 + cos(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon));
		/* not on this triangle */
		if (z > g + 0.000005) { /* TODO DBL_EPSILON */
			continue;
		}

		Az = sph_azimuth(center.lon, center.lat, ll->lon, ll->lat);

		/* 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

/* 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

/* 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_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)

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;
}

/* 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(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(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;

	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);

		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;
	struct isea_pt  out, coord;

	tri = isea_transform(g, in, &out);

	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)