Expanded tabs and removed trailing whitespace in an attempt at getting cleaner diffs for upcomming commits.

10 files changed, 1488 insertions, 1488 deletions

diff --git a/src/PJ_gn_sinu.c b/src/PJ_gn_sinu.c
index bfd8bc2d..1ab43826 100644
--- a/src/PJ_gn_sinu.c
+++ b/src/PJ_gn_sinu.c
@@ -1,98 +1,98 @@
 #define PROJ_PARMS__ \
-	double	*en; \
-	double	m, n, C_x, C_y;
+    double  *en; \
+    double  m, n, C_x, C_y;
 #define PJ_LIB__
-#include	<projects.h>
+#include    <projects.h>
 PROJ_HEAD(gn_sinu, "General Sinusoidal Series") "\n\tPCyl, Sph.\n\tm= n=";
 PROJ_HEAD(sinu, "Sinusoidal (Sanson-Flamsteed)") "\n\tPCyl, Sph&Ell";
 PROJ_HEAD(eck6, "Eckert VI") "\n\tPCyl, Sph.";
 PROJ_HEAD(mbtfps, "McBryde-Thomas Flat-Polar Sinusoidal") "\n\tPCyl, Sph.";
-#define EPS10	1e-10
+#define EPS10   1e-10
 #define MAX_ITER 8
 #define LOOP_TOL 1e-7
 /* Ellipsoidal Sinusoidal only */
 FORWARD(e_forward); /* ellipsoid */
-	double s, c;
+    double s, c;
 
-	xy.y = pj_mlfn(lp.phi, s = sin(lp.phi), c = cos(lp.phi), P->en);
-	xy.x = lp.lam * c / sqrt(1. - P->es * s * s);
-	return (xy);
+    xy.y = pj_mlfn(lp.phi, s = sin(lp.phi), c = cos(lp.phi), P->en);
+    xy.x = lp.lam * c / sqrt(1. - P->es * s * s);
+    return (xy);
 }
 INVERSE(e_inverse); /* ellipsoid */
-	double s;
+    double s;
 
-	if ((s = fabs(lp.phi = pj_inv_mlfn(P->ctx, xy.y, P->es, P->en))) < HALFPI) {
-		s = sin(lp.phi);
-		lp.lam = xy.x * sqrt(1. - P->es * s * s) / cos(lp.phi);
-	} else if ((s - EPS10) < HALFPI)
-		lp.lam = 0.;
-	else I_ERROR;
-	return (lp);
+    if ((s = fabs(lp.phi = pj_inv_mlfn(P->ctx, xy.y, P->es, P->en))) < HALFPI) {
+        s = sin(lp.phi);
+        lp.lam = xy.x * sqrt(1. - P->es * s * s) / cos(lp.phi);
+    } else if ((s - EPS10) < HALFPI)
+        lp.lam = 0.;
+    else I_ERROR;
+    return (lp);
 }
 /* General spherical sinusoidals */
 FORWARD(s_forward); /* sphere */
-	if (!P->m)
-		lp.phi = P->n != 1. ? aasin(P->ctx,P->n * sin(lp.phi)): lp.phi;
-	else {
-		double k, V;
-		int i;
+    if (!P->m)
+        lp.phi = P->n != 1. ? aasin(P->ctx,P->n * sin(lp.phi)): lp.phi;
+    else {
+        double k, V;
+        int i;
 
-		k = P->n * sin(lp.phi);
-		for (i = MAX_ITER; i ; --i) {
-			lp.phi -= V = (P->m * lp.phi + sin(lp.phi) - k) /
-				(P->m + cos(lp.phi));
-			if (fabs(V) < LOOP_TOL)
-				break;
-		}
-		if (!i)
-			F_ERROR
-	}
-	xy.x = P->C_x * lp.lam * (P->m + cos(lp.phi));
-	xy.y = P->C_y * lp.phi;
-	return (xy);
+        k = P->n * sin(lp.phi);
+        for (i = MAX_ITER; i ; --i) {
+            lp.phi -= V = (P->m * lp.phi + sin(lp.phi) - k) /
+                (P->m + cos(lp.phi));
+            if (fabs(V) < LOOP_TOL)
+                break;
+        }
+        if (!i)
+            F_ERROR
+    }
+    xy.x = P->C_x * lp.lam * (P->m + cos(lp.phi));
+    xy.y = P->C_y * lp.phi;
+    return (xy);
 }
 INVERSE(s_inverse); /* sphere */
-	xy.y /= P->C_y;
-	lp.phi = P->m ? aasin(P->ctx,(P->m * xy.y + sin(xy.y)) / P->n) :
-		( P->n != 1. ? aasin(P->ctx,sin(xy.y) / P->n) : xy.y );
-	lp.lam = xy.x / (P->C_x * (P->m + cos(xy.y)));
-	return (lp);
+    xy.y /= P->C_y;
+    lp.phi = P->m ? aasin(P->ctx,(P->m * xy.y + sin(xy.y)) / P->n) :
+        ( P->n != 1. ? aasin(P->ctx,sin(xy.y) / P->n) : xy.y );
+    lp.lam = xy.x / (P->C_x * (P->m + cos(xy.y)));
+    return (lp);
 }
 FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
-	static void /* for spheres, only */
+    static void /* for spheres, only */
 setup(PJ *P) {
-	P->es = 0;
-	P->C_x = (P->C_y = sqrt((P->m + 1.) / P->n))/(P->m + 1.);
-	P->inv = s_inverse;
-	P->fwd = s_forward;
+    P->es = 0;
+    P->C_x = (P->C_y = sqrt((P->m + 1.) / P->n))/(P->m + 1.);
+    P->inv = s_inverse;
+    P->fwd = s_forward;
 }
 ENTRY1(sinu, en)
-	if (!(P->en = pj_enfn(P->es)))
-		E_ERROR_0;
-	if (P->es) {
-		P->inv = e_inverse;
-		P->fwd = e_forward;
-	} else {
-		P->n = 1.;
-		P->m = 0.;
-		setup(P);
-	}
+    if (!(P->en = pj_enfn(P->es)))
+        E_ERROR_0;
+    if (P->es) {
+        P->inv = e_inverse;
+        P->fwd = e_forward;
+    } else {
+        P->n = 1.;
+        P->m = 0.;
+        setup(P);
+    }
 ENDENTRY(P)
 ENTRY1(eck6, en)
-	P->m = 1.;
-	P->n = 2.570796326794896619231321691;
-	setup(P);
+    P->m = 1.;
+    P->n = 2.570796326794896619231321691;
+    setup(P);
 ENDENTRY(P)
 ENTRY1(mbtfps, en)
-	P->m = 0.5;
-	P->n = 1.785398163397448309615660845;
-	setup(P);
+    P->m = 0.5;
+    P->n = 1.785398163397448309615660845;
+    setup(P);
 ENDENTRY(P)
 ENTRY1(gn_sinu, en)
-	if (pj_param(P->ctx, P->params, "tn").i && pj_param(P->ctx, P->params, "tm").i) {
-		P->n = pj_param(P->ctx, P->params, "dn").f;
-		P->m = pj_param(P->ctx, P->params, "dm").f;
-	} else
-		E_ERROR(-99)
-	setup(P);
+    if (pj_param(P->ctx, P->params, "tn").i && pj_param(P->ctx, P->params, "tm").i) {
+        P->n = pj_param(P->ctx, P->params, "dn").f;
+        P->m = pj_param(P->ctx, P->params, "dm").f;
+    } else
+        E_ERROR(-99)
+    setup(P);
 ENDENTRY(P)
diff --git a/src/PJ_goode.c b/src/PJ_goode.c
index 387557e6..b03e15d1 100644
--- a/src/PJ_goode.c
+++ b/src/PJ_goode.c
@@ -1,49 +1,49 @@
 #define PROJ_PARMS__ \
-	struct PJconsts	*sinu; \
-	struct PJconsts	*moll;
+    struct PJconsts *sinu; \
+    struct PJconsts *moll;
 #define PJ_LIB__
-#include	<projects.h>
+#include    <projects.h>
 PROJ_HEAD(goode, "Goode Homolosine") "\n\tPCyl, Sph.";
-	C_NAMESPACE PJ
+    C_NAMESPACE PJ
 *pj_sinu(PJ *), *pj_moll(PJ *);
-#define Y_COR		0.05280
-#define PHI_LIM	.71093078197902358062
+#define Y_COR       0.05280
+#define PHI_LIM .71093078197902358062
 FORWARD(s_forward); /* spheroid */
-	if (fabs(lp.phi) <= PHI_LIM)
-		xy = P->sinu->fwd(lp, P->sinu);
-	else {
-		xy = P->moll->fwd(lp, P->moll);
-		xy.y -= lp.phi >= 0.0 ? Y_COR : -Y_COR;
-	}
-	return (xy);
+    if (fabs(lp.phi) <= PHI_LIM)
+        xy = P->sinu->fwd(lp, P->sinu);
+    else {
+        xy = P->moll->fwd(lp, P->moll);
+        xy.y -= lp.phi >= 0.0 ? Y_COR : -Y_COR;
+    }
+    return (xy);
 }
 INVERSE(s_inverse); /* spheroid */
-	if (fabs(xy.y) <= PHI_LIM)
-		lp = P->sinu->inv(xy, P->sinu);
-	else {
-		xy.y += xy.y >= 0.0 ? Y_COR : -Y_COR;
-		lp = P->moll->inv(xy, P->moll);
-	}
-	return (lp);
+    if (fabs(xy.y) <= PHI_LIM)
+        lp = P->sinu->inv(xy, P->sinu);
+    else {
+        xy.y += xy.y >= 0.0 ? Y_COR : -Y_COR;
+        lp = P->moll->inv(xy, P->moll);
+    }
+    return (lp);
 }
 FREEUP;
-	if (P) {
-		if (P->sinu)
-			(*(P->sinu->pfree))(P->sinu);
-		if (P->moll)
-			(*(P->moll->pfree))(P->moll);
-		pj_dalloc(P);
-	}
+    if (P) {
+        if (P->sinu)
+            (*(P->sinu->pfree))(P->sinu);
+        if (P->moll)
+            (*(P->moll->pfree))(P->moll);
+        pj_dalloc(P);
+    }
 }
 ENTRY2(goode, sinu, moll)
-	P->es = 0.;
-	if (!(P->sinu = pj_sinu(0)) || !(P->moll = pj_moll(0)))
-		E_ERROR_0;
-	P->sinu->es = 0.;
+    P->es = 0.;
+    if (!(P->sinu = pj_sinu(0)) || !(P->moll = pj_moll(0)))
+        E_ERROR_0;
+    P->sinu->es = 0.;
         P->sinu->ctx = P->ctx;
         P->moll->ctx = P->ctx;
-	if (!(P->sinu = pj_sinu(P->sinu)) || !(P->moll = pj_moll(P->moll)))
-		E_ERROR_0;
-	P->fwd = s_forward;
-	P->inv = s_inverse;
+    if (!(P->sinu = pj_sinu(P->sinu)) || !(P->moll = pj_moll(P->moll)))
+        E_ERROR_0;
+    P->fwd = s_forward;
+    P->inv = s_inverse;
 ENDENTRY(P)
diff --git a/src/PJ_gstmerc.c b/src/PJ_gstmerc.c
index bffe0b26..da67389f 100644
--- a/src/PJ_gstmerc.c
+++ b/src/PJ_gstmerc.c
@@ -1,36 +1,36 @@
 #define PROJ_PARMS__ \
-	double lamc;\
-	double phic;\
-	double c;\
-	double n1;\
-	double n2;\
+    double lamc;\
+    double phic;\
+    double c;\
+    double n1;\
+    double n2;\
     double XS;\
     double YS;
 
 #define PJ_LIB__
-# include	<projects.h>
+# include   <projects.h>
 PROJ_HEAD(gstmerc, "Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reunion)")
-	"\n\tCyl, Sph&Ell\n\tlat_0= lon_0= k_0=";
+    "\n\tCyl, Sph&Ell\n\tlat_0= lon_0= k_0=";
 FORWARD(s_forward); /* spheroid */
-	double L, Ls, sinLs1, Ls1;
-	L= P->n1*lp.lam;
+    double L, Ls, sinLs1, Ls1;
+    L= P->n1*lp.lam;
     Ls= P->c+P->n1*log(pj_tsfn(-1.0*lp.phi,-1.0*sin(lp.phi),P->e));
     sinLs1= sin(L)/cosh(Ls);
     Ls1= log(pj_tsfn(-1.0*asin(sinLs1),0.0,0.0));
     xy.x= (P->XS + P->n2*Ls1)*P->ra;
     xy.y= (P->YS + P->n2*atan(sinh(Ls)/cos(L)))*P->ra;
     /*fprintf(stderr,"fwd:\nL      =%16.13f\nLs     =%16.13f\nLs1    =%16.13f\nLP(%16.13f,%16.13f)=XY(%16.4f,%16.4f)\n",L,Ls,Ls1,lp.lam+P->lam0,lp.phi,(xy.x*P->a + P->x0)*P->to_meter,(xy.y*P->a + P->y0)*P->to_meter);*/
-	return (xy);
+    return (xy);
 }
 INVERSE(s_inverse); /* spheroid */
-	double L, LC, sinC;
-	L= atan(sinh((xy.x*P->a - P->XS)/P->n2)/cos((xy.y*P->a - P->YS)/P->n2));
+    double L, LC, sinC;
+    L= atan(sinh((xy.x*P->a - P->XS)/P->n2)/cos((xy.y*P->a - P->YS)/P->n2));
     sinC= sin((xy.y*P->a - P->YS)/P->n2)/cosh((xy.x*P->a - P->XS)/P->n2);
     LC= log(pj_tsfn(-1.0*asin(sinC),0.0,0.0));
     lp.lam= L/P->n1;
     lp.phi= -1.0*pj_phi2(P->ctx, exp((LC-P->c)/P->n1),P->e);
     /*fprintf(stderr,"inv:\nL      =%16.13f\nsinC   =%16.13f\nLC     =%16.13f\nXY(%16.4f,%16.4f)=LP(%16.13f,%16.13f)\n",L,sinC,LC,((xy.x/P->ra)+P->x0)/P->to_meter,((xy.y/P->ra)+P->y0)/P->to_meter,lp.lam+P->lam0,lp.phi);*/
-	return (lp);
+    return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(gstmerc)
@@ -42,7 +42,7 @@ ENTRY0(gstmerc)
     P->n2= P->k0*P->a*sqrt(1.0-P->es)/(1.0-P->es*sin(P->phi0)*sin(P->phi0));
     P->XS= 0;/* -P->x0 */
     P->YS= -1.0*P->n2*P->phic;/* -P->y0 */
-	P->inv= s_inverse;
-	P->fwd= s_forward;
+    P->inv= s_inverse;
+    P->fwd= s_forward;
     /*fprintf(stderr,"a  (m) =%16.4f\ne      =%16.13f\nl0(rad)=%16.13f\np0(rad)=%16.13f\nk0     =%16.4f\nX0  (m)=%16.4f\nY0  (m)=%16.4f\n\nlC(rad)=%16.13f\npC(rad)=%16.13f\nc      =%16.13f\nn1     =%16.13f\nn2 (m) =%16.4f\nXS (m) =%16.4f\nYS (m) =%16.4f\n", P->a, P->e, P->lam0, P->phi0, P->k0, P->x0, P->y0, P->lamc, P->phic, P->c, P->n1, P->n2, P->XS +P->x0, P->YS + P->y0);*/
 ENDENTRY(P)
diff --git a/src/PJ_hammer.c b/src/PJ_hammer.c
index 31e7a127..8b15fe50 100644
--- a/src/PJ_hammer.c
+++ b/src/PJ_hammer.c
@@ -1,43 +1,43 @@
 #define PROJ_PARMS__ \
-	double w; \
-	double m, rm;
+    double w; \
+    double m, rm;
 #define PJ_LIB__
 #define EPS 1.0e-10
-# include	<projects.h>
+# include   <projects.h>
 PROJ_HEAD(hammer, "Hammer & Eckert-Greifendorff")
-	"\n\tMisc Sph, \n\tW= M=";
+    "\n\tMisc Sph, \n\tW= M=";
 FORWARD(s_forward); /* spheroid */
-	double cosphi, d;
+    double cosphi, d;
 
-	d = sqrt(2./(1. + (cosphi = cos(lp.phi)) * cos(lp.lam *= P->w)));
-	xy.x = P->m * d * cosphi * sin(lp.lam);
-	xy.y = P->rm * d * sin(lp.phi);
-	return (xy);
+    d = sqrt(2./(1. + (cosphi = cos(lp.phi)) * cos(lp.lam *= P->w)));
+    xy.x = P->m * d * cosphi * sin(lp.lam);
+    xy.y = P->rm * d * sin(lp.phi);
+    return (xy);
 }
 INVERSE(s_inverse); /* spheroid */
-	double z;
-	z = sqrt(1. - 0.25*P->w*P->w*xy.x*xy.x - 0.25*xy.y*xy.y);
-	if (fabs(2.*z*z-1.) < EPS) {
-		lp.lam = HUGE_VAL;
-		lp.phi = HUGE_VAL;
-		pj_errno = -14;
-	} else {
-		lp.lam = aatan2(P->w * xy.x * z,2. * z * z - 1)/P->w;
-		lp.phi = aasin(P->ctx,z * xy.y);
-	}
-	return (lp);
+    double z;
+    z = sqrt(1. - 0.25*P->w*P->w*xy.x*xy.x - 0.25*xy.y*xy.y);
+    if (fabs(2.*z*z-1.) < EPS) {
+        lp.lam = HUGE_VAL;
+        lp.phi = HUGE_VAL;
+        pj_errno = -14;
+    } else {
+        lp.lam = aatan2(P->w * xy.x * z,2. * z * z - 1)/P->w;
+        lp.phi = aasin(P->ctx,z * xy.y);
+    }
+    return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(hammer)
-	if (pj_param(P->ctx, P->params, "tW").i) {
-		if ((P->w = fabs(pj_param(P->ctx, P->params, "dW").f)) <= 0.) E_ERROR(-27);
-	} else
-		P->w = .5;
-	if (pj_param(P->ctx, P->params, "tM").i) {
-		if ((P->m = fabs(pj_param(P->ctx, P->params, "dM").f)) <= 0.) E_ERROR(-27);
-	} else
-		P->m = 1.;
-	P->rm = 1. / P->m;
-	P->m /= P->w;
-	P->es = 0.; P->fwd = s_forward; P->inv = s_inverse;
+    if (pj_param(P->ctx, P->params, "tW").i) {
+        if ((P->w = fabs(pj_param(P->ctx, P->params, "dW").f)) <= 0.) E_ERROR(-27);
+    } else
+        P->w = .5;
+    if (pj_param(P->ctx, P->params, "tM").i) {
+        if ((P->m = fabs(pj_param(P->ctx, P->params, "dM").f)) <= 0.) E_ERROR(-27);
+    } else
+        P->m = 1.;
+    P->rm = 1. / P->m;
+    P->m /= P->w;
+    P->es = 0.; P->fwd = s_forward; P->inv = s_inverse;
 ENDENTRY(P)
diff --git a/src/PJ_hatano.c b/src/PJ_hatano.c
index 7516ba6e..4ba9d6e8 100644
--- a/src/PJ_hatano.c
+++ b/src/PJ_hatano.c
@@ -1,51 +1,51 @@
 #define PJ_LIB__
-#include	<projects.h>
+#include    <projects.h>
 PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area") "\n\tPCyl, Sph.";
-#define NITER	20
-#define EPS	1e-7
+#define NITER   20
+#define EPS 1e-7
 #define ONETOL 1.000001
-#define CN	2.67595
-#define CS	2.43763
-#define RCN	0.37369906014686373063
-#define RCS	0.41023453108141924738
-#define FYCN	1.75859
-#define FYCS	1.93052
-#define RYCN	0.56863737426006061674
-#define RYCS	0.51799515156538134803
-#define FXC	0.85
-#define RXC	1.17647058823529411764
+#define CN  2.67595
+#define CS  2.43763
+#define RCN 0.37369906014686373063
+#define RCS 0.41023453108141924738
+#define FYCN    1.75859
+#define FYCS    1.93052
+#define RYCN    0.56863737426006061674
+#define RYCS    0.51799515156538134803
+#define FXC 0.85
+#define RXC 1.17647058823529411764
 FORWARD(s_forward); /* spheroid */
-	double th1, c;
-	int i;
-	(void) P;
+    double th1, c;
+    int i;
+    (void) P;
 
-	c = sin(lp.phi) * (lp.phi < 0. ? CS : CN);
-	for (i = NITER; i; --i) {
-		lp.phi -= th1 = (lp.phi + sin(lp.phi) - c) / (1. + cos(lp.phi));
-		if (fabs(th1) < EPS) break;
-	}
-	xy.x = FXC * lp.lam * cos(lp.phi *= .5);
-	xy.y = sin(lp.phi) * (lp.phi < 0. ? FYCS : FYCN);
-	return (xy);
+    c = sin(lp.phi) * (lp.phi < 0. ? CS : CN);
+    for (i = NITER; i; --i) {
+        lp.phi -= th1 = (lp.phi + sin(lp.phi) - c) / (1. + cos(lp.phi));
+        if (fabs(th1) < EPS) break;
+    }
+    xy.x = FXC * lp.lam * cos(lp.phi *= .5);
+    xy.y = sin(lp.phi) * (lp.phi < 0. ? FYCS : FYCN);
+    return (xy);
 }
 INVERSE(s_inverse); /* spheroid */
-	double th;
+    double th;
 
-	th = xy.y * ( xy.y < 0. ? RYCS : RYCN);
-	if (fabs(th) > 1.)
-		if (fabs(th) > ONETOL)	I_ERROR
-		else			th = th > 0. ? HALFPI : - HALFPI;
-	else
-		th = asin(th);
-	lp.lam = RXC * xy.x / cos(th);
-	th += th;
-	lp.phi = (th + sin(th)) * (xy.y < 0. ? RCS : RCN);
-	if (fabs(lp.phi) > 1.)
-		if (fabs(lp.phi) > ONETOL)	I_ERROR
-		else			lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;
-	else
-		lp.phi = asin(lp.phi);
-	return (lp);
+    th = xy.y * ( xy.y < 0. ? RYCS : RYCN);
+    if (fabs(th) > 1.)
+        if (fabs(th) > ONETOL)  I_ERROR
+        else            th = th > 0. ? HALFPI : - HALFPI;
+    else
+        th = asin(th);
+    lp.lam = RXC * xy.x / cos(th);
+    th += th;
+    lp.phi = (th + sin(th)) * (xy.y < 0. ? RCS : RCN);
+    if (fabs(lp.phi) > 1.)
+        if (fabs(lp.phi) > ONETOL)  I_ERROR
+        else            lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;
+    else
+        lp.phi = asin(lp.phi);
+    return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(hatano) P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)
diff --git a/src/PJ_healpix.c b/src/PJ_healpix.c
index 12d57ff2..a645ba34 100644
--- a/src/PJ_healpix.c
+++ b/src/PJ_healpix.c
@@ -2,10 +2,10 @@
  * Project: PROJ.4
  * Purpose: Implementation of the HEALPix and rHEALPix projections.
  *          For background see <http://code.scenzgrid.org/index.php/p/scenzgrid-py/source/tree/master/docs/rhealpix_dggs.pdf>.
- * Authors: Alex Raichev (raichev@cs.auckland.ac.nz) 
+ * Authors: Alex Raichev (raichev@cs.auckland.ac.nz)
  *          Michael Speth (spethm@landcareresearch.co.nz)
- * Notes:   Raichev implemented these projections in Python and 
- *          Speth translated them into C here.   
+ * Notes:   Raichev implemented these projections in Python and
+ *          Speth translated them into C here.
  ******************************************************************************
  * Copyright (c) 2001, Thomas Flemming, tf@ttqv.com
  *
@@ -33,17 +33,17 @@
     int south_square; \
     double qp; \
     double *apa;
-# define PJ_LIB__    
-# include	<projects.h>
+# define PJ_LIB__
+# 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>
+# include   <stdio.h>
 /* Matrix for counterclockwise rotation by pi/2: */
-# define R1 {{ 0,-1},{ 1, 0}}	
+# define R1 {{ 0,-1},{ 1, 0}}
 /* Matrix for counterclockwise rotation by pi: */
-# define R2 {{-1, 0},{ 0,-1}}	
+# define R2 {{-1, 0},{ 0,-1}}
 /* Matrix for counterclockwise rotation by 3*pi/2:  */
-# define R3 {{ 0, 1},{-1, 0}}	
+# define R3 {{ 0, 1},{-1, 0}}
 /* Identity matrix */
 # define IDENT {{1, 0},{0, 1}}
 /* IDENT, R1, R2, R3, R1 inverse, R2 inverse, R3 inverse:*/
@@ -74,25 +74,25 @@ double pj_sign (double v) {
  */
 static int get_rotate_index(int index) {
     switch(index) {
-	case 0:
-	    return 0;
-	case 1:
-	    return 1;
-	case 2: 
-	    return 2;
-	case 3: 
-	    return 3;
-	case -1:
-	    return 4;
-	case -2:
-	    return 5;
-	case -3:
-	    return 6;
+    case 0:
+        return 0;
+    case 1:
+        return 1;
+    case 2:
+        return 2;
+    case 3:
+        return 3;
+    case -1:
+        return 4;
+    case -2:
+        return 5;
+    case -3:
+        return 6;
     }
     return 0;
 }
 /**
- * Return 1 if point (testx, testy) lies in the interior of the polygon 
+ * Return 1 if point (testx, testy) lies in the interior of the polygon
  * determined by the vertices in vert, and return 0 otherwise.
  * See http://paulbourke.net/geometry/polygonmesh/ for more details.
  * @param nvert the number of vertices in the polygon.
@@ -137,69 +137,69 @@ static int pnpoly(int nvert, double vert[][2], double testx, double testy) {
 }
 /**
  * Return 1 if (x, y) lies in (the interior or boundary of) the image of the
- * HEALPix projection (in case proj=0) or in the image the rHEALPix projection 
+ * HEALPix projection (in case proj=0) or in the image the rHEALPix projection
  * (in case proj=1), and return 0 otherwise.
- * @param north_square the position of the north polar square (rHEALPix only)       
- * @param south_square the position of the south polar square (rHEALPix only)  
+ * @param north_square the position of the north polar square (rHEALPix only)
+ * @param south_square the position of the south polar square (rHEALPix only)
  **/
 int in_image(double x, double y, int proj, int north_square, int south_square) {
     if (proj == 0) {
-    	double healpixVertsJit[][2] = {
-    	    {-1.0*PI- EPS, PI/4.0},
-    	    {-3.0*PI/4.0, PI/2.0 + EPS},
-    	    {-1.0*PI/2.0, PI/4.0 + EPS},
-    	    {-1.0*PI/4.0, PI/2.0 + EPS},
-    	    {0.0, PI/4.0 + EPS},
-    	    {PI/4.0, PI/2.0 + EPS},
-    	    {PI/2.0, PI/4.0 + EPS},
-    	    {3.0*PI/4.0, PI/2.0 + EPS},
-    	    {PI+ EPS, PI/4.0},
-    	    {PI+ EPS, -1.0*PI/4.0},
-    	    {3.0*PI/4.0, -1.0*PI/2.0 - EPS},
-    	    {PI/2.0, -1.0*PI/4.0 - EPS},
-    	    {PI/4.0, -1.0*PI/2.0 - EPS},
-    	    {0.0, -1.0*PI/4.0 - EPS},
-    	    {-1.0*PI/4.0, -1.0*PI/2.0 - EPS},
-    	    {-1.0*PI/2.0, -1.0*PI/4.0 - EPS},
-    	    {-3.0*PI/4.0, -1.0*PI/2.0 - EPS},
-    	    {-1.0*PI - EPS, -1.0*PI/4.0}
-    	};
-    	return pnpoly((int)sizeof(healpixVertsJit)/
-    	              sizeof(healpixVertsJit[0]), healpixVertsJit, x, y);
+        double healpixVertsJit[][2] = {
+            {-1.0*PI- EPS, PI/4.0},
+            {-3.0*PI/4.0, PI/2.0 + EPS},
+            {-1.0*PI/2.0, PI/4.0 + EPS},
+            {-1.0*PI/4.0, PI/2.0 + EPS},
+            {0.0, PI/4.0 + EPS},
+            {PI/4.0, PI/2.0 + EPS},
+            {PI/2.0, PI/4.0 + EPS},
+            {3.0*PI/4.0, PI/2.0 + EPS},
+            {PI+ EPS, PI/4.0},
+            {PI+ EPS, -1.0*PI/4.0},
+            {3.0*PI/4.0, -1.0*PI/2.0 - EPS},
+            {PI/2.0, -1.0*PI/4.0 - EPS},
+            {PI/4.0, -1.0*PI/2.0 - EPS},
+            {0.0, -1.0*PI/4.0 - EPS},
+            {-1.0*PI/4.0, -1.0*PI/2.0 - EPS},
+            {-1.0*PI/2.0, -1.0*PI/4.0 - EPS},
+            {-3.0*PI/4.0, -1.0*PI/2.0 - EPS},
+            {-1.0*PI - EPS, -1.0*PI/4.0}
+        };
+        return pnpoly((int)sizeof(healpixVertsJit)/
+                      sizeof(healpixVertsJit[0]), healpixVertsJit, x, y);
     } else {
-    	double rhealpixVertsJit[][2] = {
-    	    {-1.0*PI - EPS, PI/4.0 + EPS},
-    	    {-1.0*PI + north_square*PI/2.0- EPS, PI/4.0 + EPS},
-    	    {-1.0*PI + north_square*PI/2.0- EPS, 3*PI/4.0 + EPS},
-    	    {-1.0*PI + (north_square + 1.0)*PI/2.0 + EPS, 3*PI/4.0 + EPS},
-    	    {-1.0*PI + (north_square + 1.0)*PI/2.0 + EPS, PI/4.0 + EPS},
-    	    {PI + EPS, PI/4.0 + EPS},
-    	    {PI + EPS, -1.0*PI/4.0 - EPS},
-    	    {-1.0*PI + (south_square + 1.0)*PI/2.0 + EPS, -1.0*PI/4.0 - EPS},
-    	    {-1.0*PI + (south_square + 1.0)*PI/2.0 + EPS, -3.0*PI/4.0 - EPS},
-    	    {-1.0*PI + south_square*PI/2.0 - EPS, -3.0*PI/4.0 - EPS},
-    	    {-1.0*PI + south_square*PI/2.0 - EPS, -1.0*PI/4.0 - EPS},
-    	    {-1.0*PI - EPS, -1.0*PI/4.0 - EPS}};
-    	return pnpoly((int)sizeof(rhealpixVertsJit)/
-    	              sizeof(rhealpixVertsJit[0]), rhealpixVertsJit, x, y);
+        double rhealpixVertsJit[][2] = {
+            {-1.0*PI - EPS, PI/4.0 + EPS},
+            {-1.0*PI + north_square*PI/2.0- EPS, PI/4.0 + EPS},
+            {-1.0*PI + north_square*PI/2.0- EPS, 3*PI/4.0 + EPS},
+            {-1.0*PI + (north_square + 1.0)*PI/2.0 + EPS, 3*PI/4.0 + EPS},
+            {-1.0*PI + (north_square + 1.0)*PI/2.0 + EPS, PI/4.0 + EPS},
+            {PI + EPS, PI/4.0 + EPS},
+            {PI + EPS, -1.0*PI/4.0 - EPS},
+            {-1.0*PI + (south_square + 1.0)*PI/2.0 + EPS, -1.0*PI/4.0 - EPS},
+            {-1.0*PI + (south_square + 1.0)*PI/2.0 + EPS, -3.0*PI/4.0 - EPS},
+            {-1.0*PI + south_square*PI/2.0 - EPS, -3.0*PI/4.0 - EPS},
+            {-1.0*PI + south_square*PI/2.0 - EPS, -1.0*PI/4.0 - EPS},
+            {-1.0*PI - EPS, -1.0*PI/4.0 - EPS}};
+        return pnpoly((int)sizeof(rhealpixVertsJit)/
+                      sizeof(rhealpixVertsJit[0]), rhealpixVertsJit, x, y);
     }
 }
 /**
  * Return the authalic latitude of latitude alpha (if inverse=0) or
- * return the approximate latitude of authalic latitude alpha (if inverse=1).    
- * P contains the relavent ellipsoid parameters. 
+ * return the approximate latitude of authalic latitude alpha (if inverse=1).
+ * P contains the relavent ellipsoid parameters.
  **/
 double auth_lat(PJ *P, double alpha, int inverse) {
     if (inverse == 0) {
         /* Authalic latitude. */
         double q = pj_qsfn(sin(alpha), P->e, 1.0 - P->es);
-        double qp = P->qp; 
-	    double ratio = q/qp;
-	    if (fabsl(ratio) > 1) {
-	        /* Rounding error. */
-	        ratio = pj_sign(ratio);
-	    }
-	    return asin(ratio);
+        double qp = P->qp;
+        double ratio = q/qp;
+        if (fabsl(ratio) > 1) {
+            /* Rounding error. */
+            ratio = pj_sign(ratio);
+        }
+        return asin(ratio);
     } else {
         /* Approximation to inverse authalic latitude. */
         return pj_authlat(alpha, P->apa);
@@ -216,46 +216,46 @@ XY healpix_sphere(LP lp) {
     XY xy;
     /* equatorial region */
     if ( fabsl(phi) <= phi0) {
-	    xy.x = lam;
-	    xy.y = 3.0*PI/8.0*sin(phi);
+        xy.x = lam;
+        xy.y = 3.0*PI/8.0*sin(phi);
     } else {
-	    double lamc;
-	    double sigma = sqrt(3.0*(1 - fabsl(sin(phi))));
-	    double cn = floor(2*lam / PI + 2);
-	    if (cn >= 4) {
-	        cn = 3;
-	    }
-	    lamc = -3*PI/4 + (PI/2)*cn;
-	    xy.x = lamc + (lam - lamc)*sigma;
-	    xy.y = pj_sign(phi)*PI/4*(2 - sigma);
+        double lamc;
+        double sigma = sqrt(3.0*(1 - fabsl(sin(phi))));
+        double cn = floor(2*lam / PI + 2);
+        if (cn >= 4) {
+            cn = 3;
+        }
+        lamc = -3*PI/4 + (PI/2)*cn;
+        xy.x = lamc + (lam - lamc)*sigma;
+        xy.y = pj_sign(phi)*PI/4*(2 - sigma);
     }
     return xy;
 }
 /**
- * Return the inverse of healpix_sphere(). 
+ * Return the inverse of healpix_sphere().
 **/
 LP healpix_sphere_inverse(XY xy) {
-    LP lp; 
+    LP lp;
     double x = xy.x;
     double y = xy.y;
     double y0 = PI/4.0;
     /* Equatorial region. */
     if (fabsl(y) <= y0) {
-	    lp.lam = x;
-	    lp.phi = asin(8.0*y/(3.0*PI));	
+        lp.lam = x;
+        lp.phi = asin(8.0*y/(3.0*PI));
     } else if (fabsl(y) < PI/2.0) {
-	    double cn = floor(2.0*x/PI + 2.0);
+        double cn = floor(2.0*x/PI + 2.0);
         double xc, tau;
-	    if (cn >= 4) {
-	        cn = 3;
-	    }
-	    xc = -3.0*PI/4.0 + (PI/2.0)*cn;
-	    tau = 2.0 - 4.0*fabsl(y)/PI;
-	    lp.lam = xc + (x - xc)/tau;	
-	    lp.phi = pj_sign(y)*asin(1.0 - pow(tau , 2.0)/3.0);
+        if (cn >= 4) {
+            cn = 3;
+        }
+        xc = -3.0*PI/4.0 + (PI/2.0)*cn;
+        tau = 2.0 - 4.0*fabsl(y)/PI;
+        lp.lam = xc + (x - xc)/tau;
+        lp.phi = pj_sign(y)*asin(1.0 - pow(tau , 2.0)/3.0);
     } else {
-	    lp.lam = -1.0*PI;
-	    lp.phi = pj_sign(y)*PI/2.0;
+        lp.lam = -1.0*PI;
+        lp.phi = pj_sign(y)*PI/2.0;
     }
     return (lp);
 }
@@ -266,7 +266,7 @@ LP healpix_sphere_inverse(XY xy) {
 static void vector_add(double a[2], double b[2], double *ret) {
     int i;
     for(i = 0; i < 2; i++) {
-	    ret[i] = a[i] + b[i];
+        ret[i] = a[i] + b[i];
     }
 }
 /**
@@ -275,12 +275,12 @@ static void vector_add(double a[2], double b[2], double *ret) {
  **/
 static void vector_sub(double a[2], double b[2], double*ret) {
     int i;
-    for(i = 0; i < 2; i++) { 
-	    ret[i] = a[i] - b[i];
+    for(i = 0; i < 2; i++) {
+        ret[i] = a[i] - b[i];
     }
 }
 /**
- * Return the 2 x 1 matrix product a*b, where a is a 2 x 2 matrix and 
+ * Return the 2 x 1 matrix product a*b, where a is a 2 x 2 matrix and
  * b is a 2 x 1 matrix.
  * @param ret holds a*b.
  **/
@@ -288,18 +288,18 @@ static void dot_product(double a[2][2], double b[2], double *ret) {
     int i, j;
     int length = 2;
     for(i = 0; i < length; i++) {
-	    ret[i] = 0;
-	    for(j = 0; j < length; j++) {
-	        ret[i] += a[i][j]*b[j];
-	    }
+        ret[i] = 0;
+        for(j = 0; j < length; j++) {
+            ret[i] += a[i][j]*b[j];
+        }
     }
 }
 /**
- * Return the number of the polar cap, the pole point coordinates, and 
+ * Return the number of the polar cap, the pole point coordinates, and
  * the region that (x, y) lies in.
- * If inverse=0, then assume (x,y) lies in the image of the HEALPix 
+ * If inverse=0, then assume (x,y) lies in the image of the HEALPix
  * projection of the unit sphere.
- * If inverse=1, then assume (x,y) lies in the image of the 
+ * If inverse=1, then assume (x,y) lies in the image of the
  * (north_square, south_square)-rHEALPix projection of the unit sphere.
  **/
 static CapMap get_cap(double x, double y, int north_square, int south_square,
@@ -309,84 +309,84 @@ static CapMap get_cap(double x, double y, int north_square, int south_square,
     capmap.x = x;
     capmap.y = y;
     if (inverse == 0) {
-    	if (y > PI/4.0) {
-    	    capmap.region = north;
-    	    c = PI/2.0; 
-    	} else if (y < -1*PI/4.0) {
-    	    capmap.region = south;
-    	    c = -1*PI/2.0;
-    	} else {
-    	    capmap.region = equatorial;
-    	    capmap.cn = 0;
-    	    return capmap;
-    	}
-    	/* polar region */
-    	if (x < -1*PI/2.0) {
-    	    capmap.cn = 0;
-    	    capmap.x = (-1*3.0*PI/4.0);
-    	    capmap.y = c;
-    	} else if (x >= -1*PI/2.0 && x < 0) {
-    	    capmap.cn = 1;
-    	    capmap.x = -1*PI/4.0;
-    	    capmap.y = c;
-    	} else if (x >= 0 && x < PI/2.0) {
-    	    capmap.cn = 2;
-    	    capmap.x = PI/4.0;
-    	    capmap.y = c;
-    	} else {
-    	    capmap.cn = 3;
-    	    capmap.x = 3.0*PI/4.0;
-    	    capmap.y = c;
-    	}
-    	return capmap;
+        if (y > PI/4.0) {
+            capmap.region = north;
+            c = PI/2.0;
+        } else if (y < -1*PI/4.0) {
+            capmap.region = south;
+            c = -1*PI/2.0;
+        } else {
+            capmap.region = equatorial;
+            capmap.cn = 0;
+            return capmap;
+        }
+        /* polar region */
+        if (x < -1*PI/2.0) {
+            capmap.cn = 0;
+            capmap.x = (-1*3.0*PI/4.0);
+            capmap.y = c;
+        } else if (x >= -1*PI/2.0 && x < 0) {
+            capmap.cn = 1;
+            capmap.x = -1*PI/4.0;
+            capmap.y = c;
+        } else if (x >= 0 && x < PI/2.0) {
+            capmap.cn = 2;
+            capmap.x = PI/4.0;
+            capmap.y = c;
+        } else {
+            capmap.cn = 3;
+            capmap.x = 3.0*PI/4.0;
+            capmap.y = c;
+        }
+        return capmap;
     } else {
-    	double eps;
-    	if (y > PI/4.0) {
-    	    capmap.region = north;
-    	    capmap.x = (-3.0*PI/4.0 + north_square*PI/2.0); 
-    	    capmap.y = PI/2.0;
-    	    x = x - north_square*PI/2.0;
-    	} else if (y < -1*PI/4.0) {
-    	    capmap.region = south;
-    	    capmap.x = (-3.0*PI/4.0 + south_square*PI/2); 
-    	    capmap.y = -1*PI/2.0;
-    	    x = x - south_square*PI/2.0;
-    	} else {
-    	    capmap.region = equatorial;
-    	    capmap.cn = 0;
-    	    return capmap;
-    	}
-    	/* Polar Region, find the HEALPix polar cap number that 
-    	   x, y moves to when rHEALPix polar square is disassembled. */
-    	eps = 1e-15; /* Kludge.  Fuzz to avoid some rounding errors. */
-    	if (capmap.region == north) {
-    	    if (y >= -1*x - PI/4.0 - eps && y < x + 5.0*PI/4.0 - eps) {
-    		    capmap.cn = (north_square + 1) % 4;
-    	    } else if (y > -1*x -1*PI/4.0 + eps && y >= x + 5.0*PI/4.0 - eps) {
-    		    capmap.cn = (north_square + 2) % 4;
-    	    } else if (y <= -1*x -1*PI/4.0 + eps && y > x + 5.0*PI/4.0 + eps) {
-    		    capmap.cn = (north_square + 3) % 4;
-    	    } else {
-    		    capmap.cn = north_square;
-    	    }
-    	} else if (capmap.region == south) {
-    	    if (y <= x + PI/4.0 + eps && y > -1*x - 5.0*PI/4 + eps) {
-    		    capmap.cn = (south_square + 1) % 4;
-    	    } else if (y < x + PI/4.0 - eps && y <= -1*x - 5.0*PI/4.0 + eps) {
-    		    capmap.cn = (south_square + 2) % 4;
-    	    } else if (y >= x + PI/4.0 - eps && y < -1*x - 5.0*PI/4.0 - eps) {
-    		    capmap.cn = (south_square + 3) % 4;
-    	    } else {
-    		    capmap.cn = south_square;
-    	    }
-    	}
-    	return capmap;
+        double eps;
+        if (y > PI/4.0) {
+            capmap.region = north;
+            capmap.x = (-3.0*PI/4.0 + north_square*PI/2.0);
+            capmap.y = PI/2.0;
+            x = x - north_square*PI/2.0;
+        } else if (y < -1*PI/4.0) {
+            capmap.region = south;
+            capmap.x = (-3.0*PI/4.0 + south_square*PI/2);
+            capmap.y = -1*PI/2.0;
+            x = x - south_square*PI/2.0;
+        } else {
+            capmap.region = equatorial;
+            capmap.cn = 0;
+            return capmap;
+        }
+        /* Polar Region, find the HEALPix polar cap number that
+           x, y moves to when rHEALPix polar square is disassembled. */
+        eps = 1e-15; /* Kludge.  Fuzz to avoid some rounding errors. */
+        if (capmap.region == north) {
+            if (y >= -1*x - PI/4.0 - eps && y < x + 5.0*PI/4.0 - eps) {
+                capmap.cn = (north_square + 1) % 4;
+            } else if (y > -1*x -1*PI/4.0 + eps && y >= x + 5.0*PI/4.0 - eps) {
+                capmap.cn = (north_square + 2) % 4;
+            } else if (y <= -1*x -1*PI/4.0 + eps && y > x + 5.0*PI/4.0 + eps) {
+                capmap.cn = (north_square + 3) % 4;
+            } else {
+                capmap.cn = north_square;
+            }
+        } else if (capmap.region == south) {
+            if (y <= x + PI/4.0 + eps && y > -1*x - 5.0*PI/4 + eps) {
+                capmap.cn = (south_square + 1) % 4;
+            } else if (y < x + PI/4.0 - eps && y <= -1*x - 5.0*PI/4.0 + eps) {
+                capmap.cn = (south_square + 2) % 4;
+            } else if (y >= x + PI/4.0 - eps && y < -1*x - 5.0*PI/4.0 - eps) {
+                capmap.cn = (south_square + 3) % 4;
+            } else {
+                capmap.cn = south_square;
+            }
+        }
+        return capmap;
     }
 }
 /**
- * Rearrange point (x, y) in the HEALPix projection by 
+ * Rearrange point (x, y) in the HEALPix projection by
  * combining the polar caps into two polar squares.
- * Put the north polar square in position north_square and 
+ * Put the north polar square in position north_square and
  * the south polar square in position south_square.
  * If inverse=1, then uncombine the polar caps.
  * @param north_square integer between 0 and 3.
@@ -402,65 +402,65 @@ static XY combine_caps(double x, double y, int north_square, int south_square,
     double ret_dot[2];
     CapMap capmap = get_cap(x, y, north_square, south_square, inverse);
     if (capmap.region == equatorial) {
-	    xy.x = capmap.x;
-	    xy.y = capmap.y;
-	    return xy;
+        xy.x = capmap.x;
+        xy.y = capmap.y;
+        return xy;
     }
     v[0] = x;
     v[1] = y;
     if (inverse == 0) {
-        /* Rotate (x, y) about its polar cap tip and then translate it to  
+        /* Rotate (x, y) about its polar cap tip and then translate it to
            north_square or south_square. */
-    	int pole = 0;
-    	double (*tmpRot)[2];
-    	double c[2] = {capmap.x, capmap.y};
-    	if (capmap.region == north) {
-    	    pole = north_square;
-    	    a[0] =  (-3.0*PI/4.0 + pole*PI/2);
-    	    a[1] =  (PI/2.0 + pole*0);
-    	    tmpRot = rot[get_rotate_index(capmap.cn - pole)];
-    	    vector_sub(v, c, v_min_c);	
-    	    dot_product(tmpRot, v_min_c, ret_dot);
-    	    vector_add(ret_dot, a, vector);
-    	} else {
-    	    pole = south_square;
-    	    a[0] =  (-3.0*PI/4.0 + pole*PI/2);
-    	    a[1] =  (PI/-2.0 + pole*0);
-    	    tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
-    	    vector_sub(v, c, v_min_c);	
-    	    dot_product(tmpRot, v_min_c, ret_dot);
-    	    vector_add(ret_dot, a, vector);
-    	}
-    	xy.x = vector[0];
-    	xy.y = vector[1];
-    	return xy;
+        int pole = 0;
+        double (*tmpRot)[2];
+        double c[2] = {capmap.x, capmap.y};
+        if (capmap.region == north) {
+            pole = north_square;
+            a[0] =  (-3.0*PI/4.0 + pole*PI/2);
+            a[1] =  (PI/2.0 + pole*0);
+            tmpRot = rot[get_rotate_index(capmap.cn - pole)];
+            vector_sub(v, c, v_min_c);
+            dot_product(tmpRot, v_min_c, ret_dot);
+            vector_add(ret_dot, a, vector);
+        } else {
+            pole = south_square;
+            a[0] =  (-3.0*PI/4.0 + pole*PI/2);
+            a[1] =  (PI/-2.0 + pole*0);
+            tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
+            vector_sub(v, c, v_min_c);
+            dot_product(tmpRot, v_min_c, ret_dot);
+            vector_add(ret_dot, a, vector);
+        }
+        xy.x = vector[0];
+        xy.y = vector[1];
+        return xy;
     } else {
         /* Inverse function.
          Unrotate (x, y) and then translate it back. */
-    	int pole = 0;
-    	double (*tmpRot)[2];
-    	double c[2] = {capmap.x, capmap.y};
-    	/* disassemble */
-    	if (capmap.region == north) {
-    	    pole = north_square;
-    	    a[0] =  (-3.0*PI/4.0 + capmap.cn*PI/2);
-    	    a[1] =  (PI/2.0 + capmap.cn*0);
-    	    tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
-    	    vector_sub(v, c, v_min_c);	
-    	    dot_product(tmpRot, v_min_c, ret_dot);
-    	    vector_add(ret_dot, a, vector);
-    	} else {
-    	    pole = south_square;
-    	    a[0] =  (-3.0*PI/4.0 + capmap.cn*PI/2);
-    	    a[1] =  (PI/-2.0 + capmap.cn*0);
-    	    tmpRot = rot[get_rotate_index(capmap.cn - pole)];
-    	    vector_sub(v, c, v_min_c);	
-    	    dot_product(tmpRot, v_min_c, ret_dot);
-    	    vector_add(ret_dot, a, vector);
-    	}
-    	xy.x = vector[0];
-    	xy.y = vector[1];
-    	return xy;
+        int pole = 0;
+        double (*tmpRot)[2];
+        double c[2] = {capmap.x, capmap.y};
+        /* disassemble */
+        if (capmap.region == north) {
+            pole = north_square;
+            a[0] =  (-3.0*PI/4.0 + capmap.cn*PI/2);
+            a[1] =  (PI/2.0 + capmap.cn*0);
+            tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
+            vector_sub(v, c, v_min_c);
+            dot_product(tmpRot, v_min_c, ret_dot);
+            vector_add(ret_dot, a, vector);
+        } else {
+            pole = south_square;
+            a[0] =  (-3.0*PI/4.0 + capmap.cn*PI/2);
+            a[1] =  (PI/-2.0 + capmap.cn*0);
+            tmpRot = rot[get_rotate_index(capmap.cn - pole)];
+            vector_sub(v, c, v_min_c);
+            dot_product(tmpRot, v_min_c, ret_dot);
+            vector_add(ret_dot, a, vector);
+        }
+        xy.x = vector[0];
+        xy.y = vector[1];
+        return xy;
     }
 }
 FORWARD(s_healpix_forward); /* sphere  */
@@ -476,20 +476,20 @@ FORWARD(e_healpix_forward); /* ellipsoid  */
 INVERSE(s_healpix_inverse); /* sphere */
     /* Check whether (x, y) lies in the HEALPix image */
     if (in_image(xy.x, xy.y, 0, 0, 0) == 0) {
-	    lp.lam = HUGE_VAL;
-	    lp.phi = HUGE_VAL;
-	    pj_ctx_set_errno(P->ctx, -15);
-	    return lp;
+        lp.lam = HUGE_VAL;
+        lp.phi = HUGE_VAL;
+        pj_ctx_set_errno(P->ctx, -15);
+        return lp;
     }
     return healpix_sphere_inverse(xy);
 }
 INVERSE(e_healpix_inverse); /* ellipsoid */
     /* Check whether (x, y) lies in the HEALPix image. */
     if (in_image(xy.x, xy.y, 0, 0, 0) == 0) {
-	    lp.lam = HUGE_VAL;
-	    lp.phi = HUGE_VAL;
-	    pj_ctx_set_errno(P->ctx, -15);
-	    return lp;
+        lp.lam = HUGE_VAL;
+        lp.phi = HUGE_VAL;
+        pj_ctx_set_errno(P->ctx, -15);
+        return lp;
     }
     lp = healpix_sphere_inverse(xy);
     lp.phi = auth_lat(P, lp.phi, 1);
@@ -507,10 +507,10 @@ FORWARD(e_rhealpix_forward); /* ellipsoid */
 INVERSE(s_rhealpix_inverse); /* sphere */
     /* Check whether (x, y) lies in the rHEALPix image. */
     if (in_image(xy.x, xy.y, 1, P->north_square, P->south_square) == 0) {
-	    lp.lam = HUGE_VAL;
-	    lp.phi = HUGE_VAL;
-	    pj_ctx_set_errno(P->ctx, -15);
-	    return lp;
+        lp.lam = HUGE_VAL;
+        lp.phi = HUGE_VAL;
+        pj_ctx_set_errno(P->ctx, -15);
+        return lp;
     }
     xy = combine_caps(xy.x, xy.y, P->north_square, P->south_square, 1);
     return healpix_sphere_inverse(xy);
@@ -518,10 +518,10 @@ INVERSE(s_rhealpix_inverse); /* sphere */
 INVERSE(e_rhealpix_inverse); /* ellipsoid */
     /* Check whether (x, y) lies in the rHEALPix image. */
     if (in_image(xy.x, xy.y, 1, P->north_square, P->south_square) == 0) {
-	    lp.lam = HUGE_VAL;
-	    lp.phi = HUGE_VAL;
-	    pj_ctx_set_errno(P->ctx, -15);
-	    return lp;
+        lp.lam = HUGE_VAL;
+        lp.phi = HUGE_VAL;
+        pj_ctx_set_errno(P->ctx, -15);
+        return lp;
     }
     xy = combine_caps(xy.x, xy.y, P->north_square, P->south_square, 1);
     lp = healpix_sphere_inverse(xy);
@@ -529,23 +529,23 @@ INVERSE(e_rhealpix_inverse); /* ellipsoid */
     return lp;
 }
 FREEUP;
-	if (P) {
-		if (P->apa)
-			pj_dalloc(P->apa);
-		pj_dalloc(P);
-	}
+    if (P) {
+        if (P->apa)
+            pj_dalloc(P->apa);
+        pj_dalloc(P);
+    }
 }
 ENTRY1(healpix, apa)
     if (P->es) {
         P->apa = pj_authset(P->es); /* For auth_lat(). */
         P->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */
-    	P->a = P->a*sqrt(0.5*P->qp); /* Set P->a to authalic radius. */
+        P->a = P->a*sqrt(0.5*P->qp); /* Set P->a to authalic radius. */
         P->ra = 1.0/P->a;
-    	P->fwd = e_healpix_forward;
-    	P->inv = e_healpix_inverse; 
+        P->fwd = e_healpix_forward;
+        P->inv = e_healpix_inverse;
     } else {
-    	P->fwd = s_healpix_forward;
-    	P->inv = s_healpix_inverse; 
+        P->fwd = s_healpix_forward;
+        P->inv = s_healpix_inverse;
     }
 ENDENTRY(P)
 ENTRY1(rhealpix, apa)
@@ -553,20 +553,20 @@ ENTRY1(rhealpix, apa)
     P->south_square = pj_param(P->ctx, P->params,"isouth_square").i;
     /* Check for valid north_square and south_square inputs. */
     if (P->north_square < 0 || P->north_square > 3) {
-	    E_ERROR(-47);
+        E_ERROR(-47);
     }
     if (P->south_square < 0 || P->south_square > 3) {
-	    E_ERROR(-47);
+        E_ERROR(-47);
     }
     if (P->es) {
         P->apa = pj_authset(P->es); /* For auth_lat(). */
         P->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */
-	    P->a = P->a*sqrt(0.5*P->qp); /* Set P->a to authalic radius. */
+        P->a = P->a*sqrt(0.5*P->qp); /* Set P->a to authalic radius. */
         P->ra = 1.0/P->a;
-	    P->fwd = e_rhealpix_forward;
-	    P->inv = e_rhealpix_inverse; 
+        P->fwd = e_rhealpix_forward;
+        P->inv = e_rhealpix_inverse;
     } else {
-	    P->fwd = s_rhealpix_forward;
-	    P->inv = s_rhealpix_inverse; 
+        P->fwd = s_rhealpix_forward;
+        P->inv = s_rhealpix_inverse;
     }
 ENDENTRY(P)
diff --git a/src/PJ_igh.c b/src/PJ_igh.c
index 4155c856..deca75b4 100644
--- a/src/PJ_igh.c
+++ b/src/PJ_igh.c
@@ -2,9 +2,9 @@
         struct PJconsts* pj[12]; \
         double dy0;
 #define PJ_LIB__
-#include	<projects.h>
+#include    <projects.h>
 PROJ_HEAD(igh, "Interrupted Goode Homolosine") "\n\tPCyl, Sph.";
-	C_NAMESPACE PJ
+    C_NAMESPACE PJ
 *pj_sinu(PJ *), *pj_moll(PJ *);
 
 static const double d4044118 = (40 + 44/60. + 11.8/3600.) * DEG_TO_RAD; // 40d 44' 11.8" [degrees]
@@ -116,9 +116,9 @@ FREEUP;
         if (P) {
                 int i;
                 for (i = 0; i < 12; ++i)
-		{
-			if (P->pj[i]) 
-				(*(P->pj[i]->pfree))(P->pj[i]); 
+        {
+            if (P->pj[i])
+                (*(P->pj[i]->pfree))(P->pj[i]);
                 }
                 pj_dalloc(P);
         }
@@ -149,7 +149,7 @@ ENTRY0(igh)
     if (!(P->pj[n-1] = pj_##proj(P->pj[n-1]))) E_ERROR_0; \
     P->pj[n-1]->x0 = x_0; \
     P->pj[n-1]->y0 = y_0; \
-    P->pj[n-1]->lam0 = lon_0; 
+    P->pj[n-1]->lam0 = lon_0;
 
         LP lp = { 0, d4044118 };
         XY xy1;
diff --git a/src/PJ_imw_p.c b/src/PJ_imw_p.c
index ae411116..1f209172 100644
--- a/src/PJ_imw_p.c
+++ b/src/PJ_imw_p.c
@@ -1,151 +1,151 @@
 #define PROJ_PARMS__ \
-	double	P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2; \
-	double	phi_1, phi_2, lam_1; \
-	double	*en; \
-	int	mode; /* = 0, phi_1 and phi_2 != 0, = 1, phi_1 = 0, = -1 phi_2 = 0 */
+    double  P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2; \
+    double  phi_1, phi_2, lam_1; \
+    double  *en; \
+    int mode; /* = 0, phi_1 and phi_2 != 0, = 1, phi_1 = 0, = -1 phi_2 = 0 */
 #define PJ_LIB__
-#include	<projects.h>
+#include    <projects.h>
 PROJ_HEAD(imw_p, "International Map of the World Polyconic")
-	"\n\tMod. Polyconic, Ell\n\tlat_1= and lat_2= [lon_1=]";
+    "\n\tMod. Polyconic, Ell\n\tlat_1= and lat_2= [lon_1=]";
 #define TOL 1e-10
 #define EPS 1e-10
-	static int
+    static int
 phi12(PJ *P, double *del, double *sig) {
-	int err = 0;
+    int err = 0;
 
-	if (!pj_param(P->ctx, P->params, "tlat_1").i ||
-		!pj_param(P->ctx, P->params, "tlat_2").i) {
-		err = -41;
-	} else {
-		P->phi_1 = pj_param(P->ctx, P->params, "rlat_1").f;
-		P->phi_2 = pj_param(P->ctx, P->params, "rlat_2").f;
-		*del = 0.5 * (P->phi_2 - P->phi_1);
-		*sig = 0.5 * (P->phi_2 + P->phi_1);
-		err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? -42 : 0;
-	}
-	return err;
+    if (!pj_param(P->ctx, P->params, "tlat_1").i ||
+        !pj_param(P->ctx, P->params, "tlat_2").i) {
+        err = -41;
+    } else {
+        P->phi_1 = pj_param(P->ctx, P->params, "rlat_1").f;
+        P->phi_2 = pj_param(P->ctx, P->params, "rlat_2").f;
+        *del = 0.5 * (P->phi_2 - P->phi_1);
+        *sig = 0.5 * (P->phi_2 + P->phi_1);
+        err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? -42 : 0;
+    }
+    return err;
 }
-	static XY
+    static XY
 loc_for(LP lp, PJ *P, double *yc) {
-	XY xy;
+    XY xy;
 
-	if (! lp.phi) {
-		xy.x = lp.lam;
-		xy.y = 0.;
-	} else {
-		double xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;
+    if (! lp.phi) {
+        xy.x = lp.lam;
+        xy.y = 0.;
+    } else {
+        double xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;
 
-		sp = sin(lp.phi);
-		m = pj_mlfn(lp.phi, sp, cos(lp.phi), P->en);
-		xa = P->Pp + P->Qp * m;
-		ya = P->P + P->Q * m;
-		R = 1. / (tan(lp.phi) * sqrt(1. - P->es * sp * sp));
-		C = sqrt(R * R - xa * xa);
-		if (lp.phi < 0.) C = - C;
-		C += ya - R;
-		if (P->mode < 0) {
-			xb = lp.lam;
-			yb = P->C2;
-		} else {
-			t = lp.lam * P->sphi_2;
-			xb = P->R_2 * sin(t);
-			yb = P->C2 + P->R_2 * (1. - cos(t));
-		}
-		if (P->mode > 0) {
-			xc = lp.lam;
-			*yc = 0.;
-		} else {
-			t = lp.lam * P->sphi_1;
-			xc = P->R_1 * sin(t);
-			*yc = P->R_1 * (1. - cos(t));
-		}
-		D = (xb - xc)/(yb - *yc);
-		B = xc + D * (C + R - *yc);
-		xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
-		if (lp.phi > 0)
-			xy.x = - xy.x;
-		xy.x = (B + xy.x) / (1. + D * D);
-		xy.y = sqrt(R * R - xy.x * xy.x);
-		if (lp.phi > 0)
-			xy.y = - xy.y;
-		xy.y += C + R;
-	}
-	return (xy);
+        sp = sin(lp.phi);
+        m = pj_mlfn(lp.phi, sp, cos(lp.phi), P->en);
+        xa = P->Pp + P->Qp * m;
+        ya = P->P + P->Q * m;
+        R = 1. / (tan(lp.phi) * sqrt(1. - P->es * sp * sp));
+        C = sqrt(R * R - xa * xa);
+        if (lp.phi < 0.) C = - C;
+        C += ya - R;
+        if (P->mode < 0) {
+            xb = lp.lam;
+            yb = P->C2;
+        } else {
+            t = lp.lam * P->sphi_2;
+            xb = P->R_2 * sin(t);
+            yb = P->C2 + P->R_2 * (1. - cos(t));
+        }
+        if (P->mode > 0) {
+            xc = lp.lam;
+            *yc = 0.;
+        } else {
+            t = lp.lam * P->sphi_1;
+            xc = P->R_1 * sin(t);
+            *yc = P->R_1 * (1. - cos(t));
+        }
+        D = (xb - xc)/(yb - *yc);
+        B = xc + D * (C + R - *yc);
+        xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
+        if (lp.phi > 0)
+            xy.x = - xy.x;
+        xy.x = (B + xy.x) / (1. + D * D);
+        xy.y = sqrt(R * R - xy.x * xy.x);
+        if (lp.phi > 0)
+            xy.y = - xy.y;
+        xy.y += C + R;
+    }
+    return (xy);
 }
 FORWARD(e_forward); /* ellipsoid */
-	double yc;
-	xy = loc_for(lp, P, &yc);
-	return (xy);
+    double yc;
+    xy = loc_for(lp, P, &yc);
+    return (xy);
 }
 INVERSE(e_inverse); /* ellipsoid */
-	XY t;
-	double yc;
+    XY t;
+    double yc;
 
-	lp.phi = P->phi_2;
-	lp.lam = xy.x / cos(lp.phi);
-	do {
-		t = loc_for(lp, P, &yc);
-		lp.phi = ((lp.phi - P->phi_1) * (xy.y - yc) / (t.y - yc)) + P->phi_1;
-		lp.lam = lp.lam * xy.x / t.x;
-	} while (fabs(t.x - xy.x) > TOL || fabs(t.y - xy.y) > TOL);
-	return (lp);
+    lp.phi = P->phi_2;
+    lp.lam = xy.x / cos(lp.phi);
+    do {
+        t = loc_for(lp, P, &yc);
+        lp.phi = ((lp.phi - P->phi_1) * (xy.y - yc) / (t.y - yc)) + P->phi_1;
+        lp.lam = lp.lam * xy.x / t.x;
+    } while (fabs(t.x - xy.x) > TOL || fabs(t.y - xy.y) > TOL);
+    return (lp);
 }
-	static void
+    static void
 xy(PJ *P, double phi, double *x, double *y, double *sp, double *R) {
-	double F;
+    double F;
 
-	*sp = sin(phi);
-	*R = 1./(tan(phi) * sqrt(1. - P->es * *sp * *sp ));
-	F = P->lam_1 * *sp;
-	*y = *R * (1 - cos(F));
-	*x = *R * sin(F);
+    *sp = sin(phi);
+    *R = 1./(tan(phi) * sqrt(1. - P->es * *sp * *sp ));
+    F = P->lam_1 * *sp;
+    *y = *R * (1 - cos(F));
+    *x = *R * sin(F);
 }
 FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
 ENTRY1(imw_p, en)
-	double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
-	int i;
+    double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
+    int i;
 
-	if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
-	if( (i = phi12(P, &del, &sig)) != 0)
-		E_ERROR(i);
-	if (P->phi_2 < P->phi_1) { /* make sure P->phi_1 most southerly */
-		del = P->phi_1;
-		P->phi_1 = P->phi_2;
-		P->phi_2 = del;
-	}
-	if (pj_param(P->ctx, P->params, "tlon_1").i)
-		P->lam_1 = pj_param(P->ctx, P->params, "rlon_1").f;
-	else { /* use predefined based upon latitude */
-		sig = fabs(sig * RAD_TO_DEG);
-		if (sig <= 60)		sig = 2.;
-		else if (sig <= 76) sig = 4.;
-		else				sig = 8.;
-		P->lam_1 = sig * DEG_TO_RAD;
-	}
-	P->mode = 0;
-	if (P->phi_1) xy(P, P->phi_1, &x1, &y1, &P->sphi_1, &P->R_1);
-	else {
-		P->mode = 1;
-		y1 = 0.;
-		x1 = P->lam_1;
-	}
-	if (P->phi_2) xy(P, P->phi_2, &x2, &T2, &P->sphi_2, &P->R_2);
-	else {
-		P->mode = -1;
-		T2 = 0.;
-		x2 = P->lam_1;
-	}
-	m1 = pj_mlfn(P->phi_1, P->sphi_1, cos(P->phi_1), P->en);
-	m2 = pj_mlfn(P->phi_2, P->sphi_2, cos(P->phi_2), P->en);
-	t = m2 - m1;
-	s = x2 - x1;
-	y2 = sqrt(t * t - s * s) + y1;
-	P->C2 = y2 - T2;
-	t = 1. / t;
-	P->P = (m2 * y1 - m1 * y2) * t;
-	P->Q = (y2 - y1) * t;
-	P->Pp = (m2 * x1 - m1 * x2) * t;
-	P->Qp = (x2 - x1) * t;
-	P->fwd = e_forward;
-	P->inv = e_inverse;
+    if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
+    if( (i = phi12(P, &del, &sig)) != 0)
+        E_ERROR(i);
+    if (P->phi_2 < P->phi_1) { /* make sure P->phi_1 most southerly */
+        del = P->phi_1;
+        P->phi_1 = P->phi_2;
+        P->phi_2 = del;
+    }
+    if (pj_param(P->ctx, P->params, "tlon_1").i)
+        P->lam_1 = pj_param(P->ctx, P->params, "rlon_1").f;
+    else { /* use predefined based upon latitude */
+        sig = fabs(sig * RAD_TO_DEG);
+        if (sig <= 60)      sig = 2.;
+        else if (sig <= 76) sig = 4.;
+        else                sig = 8.;
+        P->lam_1 = sig * DEG_TO_RAD;
+    }
+    P->mode = 0;
+    if (P->phi_1) xy(P, P->phi_1, &x1, &y1, &P->sphi_1, &P->R_1);
+    else {
+        P->mode = 1;
+        y1 = 0.;
+        x1 = P->lam_1;
+    }
+    if (P->phi_2) xy(P, P->phi_2, &x2, &T2, &P->sphi_2, &P->R_2);
+    else {
+        P->mode = -1;
+        T2 = 0.;
+        x2 = P->lam_1;
+    }
+    m1 = pj_mlfn(P->phi_1, P->sphi_1, cos(P->phi_1), P->en);
+    m2 = pj_mlfn(P->phi_2, P->sphi_2, cos(P->phi_2), P->en);
+    t = m2 - m1;
+    s = x2 - x1;
+    y2 = sqrt(t * t - s * s) + y1;
+    P->C2 = y2 - T2;
+    t = 1. / t;
+    P->P = (m2 * y1 - m1 * y2) * t;
+    P->Q = (y2 - y1) * t;
+    P->Pp = (m2 * x1 - m1 * x2) * t;
+    P->Qp = (x2 - x1) * t;
+    P->fwd = e_forward;
+    P->inv = e_inverse;
 ENDENTRY(P)
diff --git a/src/PJ_isea.c b/src/PJ_isea.c
index 7c09189c..0ba17807 100644
--- a/src/PJ_isea.c
+++ b/src/PJ_isea.c
@@ -30,78 +30,78 @@ struct hex {
 /* y *must* be positive down as the xy /iso conversion assumes this */
 ISEA_STATIC
 int hex_xy(struct hex *h) {
-	if (!h->iso) return 1;
-	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;
+    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;
+    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;
+    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
@@ -111,60 +111,60 @@ int hexbin2(double width, double x, double y,
 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
+    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;
+    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;
+    double x, y;
 };
 
 struct isea_geo {
-	double lon, lat;
+    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 */
+    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
+    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;
+    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},
+    {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
@@ -189,18 +189,18 @@ struct snyder_constants constants[] = {
 
 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}
+    {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 */
@@ -215,46 +215,46 @@ static int      tri_v1[] = {0, 0, 0, 0, 0, 0, 6, 7, 8, 9, 10, 2, 3, 4, 5, 1, 11,
 
 /* 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},
+    {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;
+    double          adj;
 
-	struct isea_geo v;
-	struct isea_geo c;
+    struct isea_geo v;
+    struct isea_geo c;
 
-	v = vertex[tri_v1[triangle]];
-	c = icostriangles[triangle];
+    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;
+    /* 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) */
@@ -269,49 +269,49 @@ 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;
+    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;
+    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;
+    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 */
@@ -319,173 +319,173 @@ 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));
+    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 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 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 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);
+        /* 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
-		 */
+        /*
+         * TODO
+         * translate coordinates to the origin for the particular
+         * hexagon on the flattened polyhedral map plot
+         */
 
-		out->x = x;
-		out->y = y;
-
-		return i;
-	}
+        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 */
+    /*
+     * 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 */
 }
 
 /*
@@ -509,73 +509,73 @@ 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;
+    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;
+    phi = pt->lat;
+    lambda = pt->lon;
+    alpha = np->lat;
+    beta = np->lon;
+    lambda0 = beta;
 
-	cos_p = cos(phi);
-	sin_a = sin(alpha);
+    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-7 */
+    sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0);
 
-	/* mpawm 5-8b */
+    /* 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)));
+    /* 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;
+    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;
+    /* 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);
+    phip = asin(sin_phip);
 
-	npt.lat = phip;
-	npt.lon = lambdap;
+    npt.lat = phip;
+    npt.lon = lambdap;
 
-	return npt;
+    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;
+    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 */
@@ -588,64 +588,64 @@ 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;
+    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;
+    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;
+    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_pt * out)
 {
-	struct isea_geo i, pole;
-	int             tri;
+    struct isea_geo i, pole;
+    int             tri;
 
-	pole.lat = g->o_lat;
-	pole.lon = g->o_lon;
+    pole.lat = g->o_lat;
+    pole.lon = g->o_lon;
 
-	i = isea_ctran(&pole, in, g->o_az);
+    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;
+    tri = isea_snyder_forward(&i, out);
+    out->x *= g->radius;
+    out->y *= g->radius;
+    g->triangle = tri;
 
-	return tri;
+    return tri;
 }
 
 #define DOWNTRI(tri) (((tri - 1) / 5) % 2 == 1)
@@ -654,249 +654,249 @@ ISEA_STATIC
 void
 isea_rotate(struct isea_pt * pt, double degrees)
 {
-	double          rad;
+    double          rad;
 
-	double          x, y;
+    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;
+    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);
+    x = pt->x * cos(rad) + pt->y * sin(rad);
+    y = -pt->x * sin(rad) + pt->y * cos(rad);
 
-	pt->x = x;
-	pt->y = y;
+    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;
+    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;
+    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;
+    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;
+    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;
+            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;
+    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
@@ -906,118 +906,118 @@ int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
 /* 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;
+        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;
+    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;
+    struct isea_dgg dgg;
 
 #define PJ_LIB__
 #include <projects.h>
@@ -1025,95 +1025,95 @@ isea_forward(struct isea_dgg *g, struct isea_geo *in)
 PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph";
 
 FORWARD(s_forward);
-	struct isea_pt out;
-	struct isea_geo in;
+    struct isea_pt out;
+    struct isea_geo in;
 
-	in.lon = lp.lam;
-	in.lat = lp.phi;
+    in.lon = lp.lam;
+    in.lat = lp.phi;
 
-	out = isea_forward(&P->dgg, &in);
-	
-	xy.x = out.x;
-	xy.y = out.y;
+    out = isea_forward(&P->dgg, &in);
 
-	return xy;
+    xy.x = out.x;
+    xy.y = out.y;
+
+    return xy;
 }
 FREEUP; if (P) pj_dalloc(P); }
 
 ENTRY0(isea)
-	char *opt;
+    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;
-	}
+/*      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/PJ_krovak.c b/src/PJ_krovak.c
index dd250134..09a8311e 100644
--- a/src/PJ_krovak.c
+++ b/src/PJ_krovak.c
@@ -30,7 +30,7 @@
  *****************************************************************************/
 
 #define PROJ_PARMS__ \
-	double	C_x;
+    double  C_x;
 #define PJ_LIB__
 
 #include <projects.h>
@@ -42,22 +42,22 @@ PROJ_HEAD(krovak, "Krovak") "\n\tPCyl., Ellps.";
 /**
    NOTES: According to EPSG the full Krovak projection method should have
           the following parameters.  Within PROJ.4 the azimuth, and pseudo
-          standard parallel are hardcoded in the algorithm and can't be 
+          standard parallel are hardcoded in the algorithm and can't be
           altered from outside.  The others all have defaults to match the
           common usage with Krovak projection.
 
   lat_0 = latitude of centre of the projection
-         
+
   lon_0 = longitude of centre of the projection
-  
+
   ** = azimuth (true) of the centre line passing through the centre of the projection
 
   ** = latitude of pseudo standard parallel
-   
+
   k  = scale factor on the pseudo standard parallel
-  
+
   x_0 = False Easting of the centre of the projection at the apex of the cone
-  
+
   y_0 = False Northing of the centre of the projection at the apex of the cone
 
  **/
@@ -68,127 +68,127 @@ FORWARD(e_forward); /* ellipsoid */
 /* calculate xy from lat/lon */
 
 /* Constants, identical to inverse transform function */
-	double s45, s90, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
-	double gfi, u, fi0, deltav, s, d, eps, ro;
+    double s45, s90, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
+    double gfi, u, fi0, deltav, s, d, eps, ro;
 
 
-	s45 = 0.785398163397448;    /* 45deg */
-	s90 = 2 * s45;
-	fi0 = P->phi0;    /* Latitude of projection centre 49deg 30' */
+    s45 = 0.785398163397448;    /* 45deg */
+    s90 = 2 * s45;
+    fi0 = P->phi0;    /* Latitude of projection centre 49deg 30' */
 
-   /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must 
+   /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
       be set to 1 here.
       Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
       e2=0.006674372230614;
    */
-	a =  1; /* 6377397.155; */
-	/* e2 = P->es;*/       /* 0.006674372230614; */
-	e2 = 0.006674372230614;
-	e = sqrt(e2);
+    a =  1; /* 6377397.155; */
+    /* e2 = P->es;*/       /* 0.006674372230614; */
+    e2 = 0.006674372230614;
+    e = sqrt(e2);
 
-	alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
+    alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
 
-	uq = 1.04216856380474;      /* DU(2, 59, 42, 42.69689) */
-	u0 = asin(sin(fi0) / alfa);
-	g = pow(   (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2.  );
+    uq = 1.04216856380474;      /* DU(2, 59, 42, 42.69689) */
+    u0 = asin(sin(fi0) / alfa);
+    g = pow(   (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2.  );
 
-	k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
+    k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
 
-	k1 = P->k0;
-	n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
-	s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78deg 30'00" N */
-	n = sin(s0);
-	ro0 = k1 * n0 / tan(s0);
-	ad = s90 - uq;
+    k1 = P->k0;
+    n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
+    s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78deg 30'00" N */
+    n = sin(s0);
+    ro0 = k1 * n0 / tan(s0);
+    ad = s90 - uq;
 
 /* Transformation */
 
-	gfi =pow ( ((1. + e * sin(lp.phi)) /
+    gfi =pow ( ((1. + e * sin(lp.phi)) /
                (1. - e * sin(lp.phi))) , (alfa * e / 2.));
 
-	u= 2. * (atan(k * pow( tan(lp.phi / 2. + s45), alfa) / gfi)-s45);
+    u= 2. * (atan(k * pow( tan(lp.phi / 2. + s45), alfa) / gfi)-s45);
 
-	deltav = - lp.lam * alfa;
+    deltav = - lp.lam * alfa;
 
-	s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav));
-	d = asin(cos(u) * sin(deltav) / cos(s));
-	eps = n * d;
-	ro = ro0 * pow(tan(s0 / 2. + s45) , n) / pow(tan(s / 2. + s45) , n)   ;
+    s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav));
+    d = asin(cos(u) * sin(deltav) / cos(s));
+    eps = n * d;
+    ro = ro0 * pow(tan(s0 / 2. + s45) , n) / pow(tan(s / 2. + s45) , n)   ;
 
    /* x and y are reverted! */
-	xy.y = ro * cos(eps) / a;
-	xy.x = ro * sin(eps) / a;
+    xy.y = ro * cos(eps) / a;
+    xy.x = ro * sin(eps) / a;
 
         if( !pj_param(P->ctx, P->params, "tczech").i )
-	  {
-	    xy.y *= -1.0;
-	    xy.x *= -1.0;
-	  }
+      {
+        xy.y *= -1.0;
+        xy.x *= -1.0;
+      }
 
-	return (xy);
+    return (xy);
 }
 
 
 
 INVERSE(e_inverse); /* ellipsoid */
-	/* calculate lat/lon from xy */
+    /* calculate lat/lon from xy */
 
 /* Constants, identisch wie in der Umkehrfunktion */
-	double s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
-	double u, deltav, s, d, eps, ro, fi1, xy0;
-	int ok;
+    double s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
+    double u, deltav, s, d, eps, ro, fi1, xy0;
+    int ok;
 
-	s45 = 0.785398163397448;    /* 45deg */
-	s90 = 2 * s45;
-	fi0 = P->phi0;    /* Latitude of projection centre 49deg 30' */
+    s45 = 0.785398163397448;    /* 45deg */
+    s90 = 2 * s45;
+    fi0 = P->phi0;    /* Latitude of projection centre 49deg 30' */
 
 
-   /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must 
+   /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
       be set to 1 here.
       Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
       e2=0.006674372230614;
    */
-	a = 1; /* 6377397.155; */
-	/* e2 = P->es; */      /* 0.006674372230614; */
-	e2 = 0.006674372230614;
-	e = sqrt(e2);
+    a = 1; /* 6377397.155; */
+    /* e2 = P->es; */      /* 0.006674372230614; */
+    e2 = 0.006674372230614;
+    e = sqrt(e2);
 
-	alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
-	uq = 1.04216856380474;      /* DU(2, 59, 42, 42.69689) */
-	u0 = asin(sin(fi0) / alfa);
-	g = pow(   (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2.  );
+    alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
+    uq = 1.04216856380474;      /* DU(2, 59, 42, 42.69689) */
+    u0 = asin(sin(fi0) / alfa);
+    g = pow(   (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2.  );
 
-	k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
+    k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
 
-	k1 = P->k0;
-	n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
-	s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78deg 30'00" N */
-	n = sin(s0);
-	ro0 = k1 * n0 / tan(s0);
-	ad = s90 - uq;
+    k1 = P->k0;
+    n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
+    s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78deg 30'00" N */
+    n = sin(s0);
+    ro0 = k1 * n0 / tan(s0);
+    ad = s90 - uq;
 
 
 /* Transformation */
    /* revert y, x*/
-	xy0=xy.x;
-	xy.x=xy.y;
-	xy.y=xy0;
+    xy0=xy.x;
+    xy.x=xy.y;
+    xy.y=xy0;
 
         if( !pj_param(P->ctx, P->params, "tczech").i )
-	  {
-	    xy.x *= -1.0;
-	    xy.y *= -1.0;
-	  }
+      {
+        xy.x *= -1.0;
+        xy.y *= -1.0;
+      }
 
-	ro = sqrt(xy.x * xy.x + xy.y * xy.y);
-	eps = atan2(xy.y, xy.x);
-	d = eps / sin(s0);
-	s = 2. * (atan(  pow(ro0 / ro, 1. / n) * tan(s0 / 2. + s45)) - s45);
+    ro = sqrt(xy.x * xy.x + xy.y * xy.y);
+    eps = atan2(xy.y, xy.x);
+    d = eps / sin(s0);
+    s = 2. * (atan(  pow(ro0 / ro, 1. / n) * tan(s0 / 2. + s45)) - s45);
 
-	u = asin(cos(ad) * sin(s) - sin(ad) * cos(s) * cos(d));
-	deltav = asin(cos(s) * sin(d) / cos(u));
+    u = asin(cos(ad) * sin(s) - sin(ad) * cos(s) * cos(d));
+    deltav = asin(cos(s) * sin(d) / cos(u));
 
-	lp.lam = P->lam0 - deltav / alfa;
+    lp.lam = P->lam0 - deltav / alfa;
 
 /* ITERATION FOR lp.phi */
    fi1 = u;
@@ -196,7 +196,7 @@ INVERSE(e_inverse); /* ellipsoid */
    ok = 0;
    do
    {
-   	lp.phi = 2. * ( atan( pow( k, -1. / alfa)  *
+    lp.phi = 2. * ( atan( pow( k, -1. / alfa)  *
                             pow( tan(u / 2. + s45) , 1. / alfa)  *
                             pow( (1. + e * sin(fi1)) / (1. - e * sin(fi1)) , e / 2.)
                            )  - s45);
@@ -215,34 +215,34 @@ INVERSE(e_inverse); /* ellipsoid */
 FREEUP; if (P) pj_dalloc(P); }
 
 ENTRY0(krovak)
-	double ts;
-	/* read some Parameters,
-	 * here Latitude Truescale */
+    double ts;
+    /* read some Parameters,
+     * here Latitude Truescale */
+
+    ts = pj_param(P->ctx, P->params, "rlat_ts").f;
+    P->C_x = ts;
 
-	ts = pj_param(P->ctx, P->params, "rlat_ts").f;
-	P->C_x = ts;
-	
-	/* we want Bessel as fixed ellipsoid */
-	P->a = 6377397.155;
-	P->e = sqrt(P->es = 0.006674372230614);
+    /* we want Bessel as fixed ellipsoid */
+    P->a = 6377397.155;
+    P->e = sqrt(P->es = 0.006674372230614);
 
         /* if latitude of projection center is not set, use 49d30'N */
-	if (!pj_param(P->ctx, P->params, "tlat_0").i)
-            P->phi0 = 0.863937979737193; 
+    if (!pj_param(P->ctx, P->params, "tlat_0").i)
+            P->phi0 = 0.863937979737193;
 
         /* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */
         /* that will correspond to using longitudes relative to greenwich    */
         /* as input and output, instead of lat/long relative to Ferro */
-	if (!pj_param(P->ctx, P->params, "tlon_0").i)
+    if (!pj_param(P->ctx, P->params, "tlon_0").i)
             P->lam0 = 0.7417649320975901 - 0.308341501185665;
 
         /* if scale not set default to 0.9999 */
-	if (!pj_param(P->ctx, P->params, "tk").i)
+    if (!pj_param(P->ctx, P->params, "tk").i)
             P->k0 = 0.9999;
 
-	/* always the same */
-        P->inv = e_inverse; 
-	P->fwd = e_forward;
+    /* always the same */
+        P->inv = e_inverse;
+    P->fwd = e_forward;
 
 ENDENTRY(P)