aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--man/man1/cs2cs.116
-rw-r--r--src/PJ_eck1.c184
-rw-r--r--src/PJ_eck4.c234
-rw-r--r--src/PJ_eck5.c186
-rw-r--r--src/PJ_fahey.c190
-rw-r--r--src/PJ_goode.c244
-rw-r--r--src/PJ_hatano.c268
-rw-r--r--src/PJ_isea.c2350
-rw-r--r--src/PJ_vandg4.c194
-rw-r--r--src/pj_factors.c2
-rw-r--r--src/pj_gridinfo.c2
-rw-r--r--src/pj_list.c1
-rw-r--r--src/pj_list.h1
13 files changed, 1935 insertions, 1937 deletions
diff --git a/man/man1/cs2cs.1 b/man/man1/cs2cs.1
index 61221389..bfebe8ea 100644
--- a/man/man1/cs2cs.1
+++ b/man/man1/cs2cs.1
@@ -2,7 +2,7 @@
.nr LL 5.5i
.ad b
.hy 1
-.TH CS2CS 1 "2000/03/21 Rel. 4.4"
+.TH CS2CS 1 "2000/03/21 Rel. 4.4"
.SH NAME
cs2cs \- cartographic coordinate system filter
.SH SYNOPSIS
@@ -21,7 +21,7 @@ file[s]
.I Cs2cs
performs transformation between the source and destination cartographic
coordinate system on a set of input points. The coordinate system
-transformation can include translation between projected and geographic
+transformation can include translation between projected and geographic
coordinates as well as the application of datum shifts.
.PP
The following control parameters can appear in any order:
@@ -75,7 +75,7 @@ that can be selected with
.B +units
or
.B \-ld
-list of datums that can be selected with
+list of datums that can be selected with
.B +datum.
.TP
.BI \-r
@@ -121,13 +121,13 @@ and supplementary documentation for Release 4.
.PP
The \fIcs2cs\fR program requires two coordinate system definitions. The
first (or primary is defined based on all projection parameters not
-appearing after the \fB+to\fR argument. All projection parameters
+appearing after the \fB+to\fR argument. All projection parameters
appearing after the \fB+to\fR argument are considered the definition
of the second coordinate system. If there is no second coordinate system
defined, a geographic coordinate system based on the datum and ellipsoid of
the source coordinate system is assumed. Note that the source and destination
coordinate system can both be projections, both be geographic, or one of
-each and may have the same or different datums.
+each and may have the same or different datums.
.PP
Additional projection control parameters may be contained in two
auxiliary control files:
@@ -160,15 +160,15 @@ Input geographic data
(longitude and latitude) must be in DMS or decimal degrees format and input
cartesian data must be in units consistent with the ellipsoid
major axis or sphere radius units.
-Output geographic coordinates will normally be in DMS format (use
-.B \-f %.12f
+Output geographic coordinates will normally be in DMS format (use
+.B \-f %.12f
for decimal degrees with 12 decimal places), while
projected (cartesian) coordinates will be in linear (meter, feet) units.
.SH EXAMPLE
The following script
.RS 5
- \f(CWcs2cs +proj=latlong +datum=NAD83
+ \f(CWcs2cs +proj=latlong +datum=NAD83
+to +proj=utm +zone=10 +datum=NAD27 \-r <<EOF
45d15'33.1" 111.5W
45d15.551666667N \-111d30
diff --git a/src/PJ_eck1.c b/src/PJ_eck1.c
index da2c8685..f285e00b 100644
--- a/src/PJ_eck1.c
+++ b/src/PJ_eck1.c
@@ -1,92 +1,92 @@
-#define PJ_LIB__
-#include <projects.h>
-
-PROJ_HEAD(eck1, "Eckert I") "\n\tPCyl., Sph.";
-#define FC 0.92131773192356127802
-#define RP 0.31830988618379067154
-
-
-static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
- XY xy = {0.0,0.0};
- (void) P;
-
- xy.x = FC * lp.lam * (1. - RP * fabs(lp.phi));
- xy.y = FC * lp.phi;
-
- return xy;
-}
-
-
-static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
- LP lp = {0.0,0.0};
- (void) P;
-
- lp.phi = xy.y / FC;
- lp.lam = xy.x / (FC * (1. - RP * fabs(lp.phi)));
-
- return (lp);
-}
-
-
-static void *freeup_new (PJ *P) { /* Destructor */
- return pj_dealloc(P);
-}
-
-
-static void freeup (PJ *P) {
- freeup_new (P);
- return;
-}
-
-
-PJ *PROJECTION(eck1) {
- P->es = 0.0;
- P->inv = s_inverse;
- P->fwd = s_forward;
-
- return P ;
-}
-
-
-#ifdef PJ_OMIT_SELFTEST
-int pj_eck1_selftest (void) {return 0;}
-#else
-
-int pj_eck1_selftest (void) {
- double tolerance_lp = 1e-10;
- double tolerance_xy = 1e-7;
-
- char s_args[] = {"+proj=eck1 +a=6400000 +lat_1=0.5 +lat_2=2"};
-
- LP fwd_in[] = {
- { 2, 1},
- { 2,-1},
- {-2, 1},
- {-2,-1}
- };
- XY s_fwd_expect[] = {
- { 204680.88820295094, 102912.17842606473},
- { 204680.88820295094, -102912.17842606473},
- {-204680.88820295094, 102912.17842606473},
- {-204680.88820295094, -102912.17842606473},
- };
-
- XY inv_in[] = {
- { 200, 100},
- { 200,-100},
- {-200, 100},
- {-200,-100}
- };
-
- LP s_inv_expect[] = {
- { 0.0019434150820034624, 0.00097170229538813102},
- { 0.0019434150820034624, -0.00097170229538813102},
- {-0.0019434150820034624, 0.00097170229538813102},
- {-0.0019434150820034624, -0.00097170229538813102},
- };
-
- return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
-}
-
-
-#endif
+#define PJ_LIB__
+#include <projects.h>
+
+PROJ_HEAD(eck1, "Eckert I") "\n\tPCyl., Sph.";
+#define FC 0.92131773192356127802
+#define RP 0.31830988618379067154
+
+
+static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
+ XY xy = {0.0,0.0};
+ (void) P;
+
+ xy.x = FC * lp.lam * (1. - RP * fabs(lp.phi));
+ xy.y = FC * lp.phi;
+
+ return xy;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
+ LP lp = {0.0,0.0};
+ (void) P;
+
+ lp.phi = xy.y / FC;
+ lp.lam = xy.x / (FC * (1. - RP * fabs(lp.phi)));
+
+ return (lp);
+}
+
+
+static void *freeup_new (PJ *P) { /* Destructor */
+ return pj_dealloc(P);
+}
+
+
+static void freeup (PJ *P) {
+ freeup_new (P);
+ return;
+}
+
+
+PJ *PROJECTION(eck1) {
+ P->es = 0.0;
+ P->inv = s_inverse;
+ P->fwd = s_forward;
+
+ return P ;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_eck1_selftest (void) {return 0;}
+#else
+
+int pj_eck1_selftest (void) {
+ double tolerance_lp = 1e-10;
+ double tolerance_xy = 1e-7;
+
+ char s_args[] = {"+proj=eck1 +a=6400000 +lat_1=0.5 +lat_2=2"};
+
+ LP fwd_in[] = {
+ { 2, 1},
+ { 2,-1},
+ {-2, 1},
+ {-2,-1}
+ };
+ XY s_fwd_expect[] = {
+ { 204680.88820295094, 102912.17842606473},
+ { 204680.88820295094, -102912.17842606473},
+ {-204680.88820295094, 102912.17842606473},
+ {-204680.88820295094, -102912.17842606473},
+ };
+
+ XY inv_in[] = {
+ { 200, 100},
+ { 200,-100},
+ {-200, 100},
+ {-200,-100}
+ };
+
+ LP s_inv_expect[] = {
+ { 0.0019434150820034624, 0.00097170229538813102},
+ { 0.0019434150820034624, -0.00097170229538813102},
+ {-0.0019434150820034624, 0.00097170229538813102},
+ {-0.0019434150820034624, -0.00097170229538813102},
+ };
+
+ return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
+}
+
+
+#endif
diff --git a/src/PJ_eck4.c b/src/PJ_eck4.c
index 8a56f019..66eaa9d0 100644
--- a/src/PJ_eck4.c
+++ b/src/PJ_eck4.c
@@ -1,117 +1,117 @@
-#define PJ_LIB__
-#include <projects.h>
-
-PROJ_HEAD(eck4, "Eckert IV") "\n\tPCyl, Sph.";
-
-#define C_x .42223820031577120149
-#define C_y 1.32650042817700232218
-#define RC_y .75386330736002178205
-#define C_p 3.57079632679489661922
-#define RC_p .28004957675577868795
-#define EPS 1e-7
-#define NITER 6
-
-
-static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
- XY xy = {0.0,0.0};
- double p, V, s, c;
- int i;
- (void) P;
-
- p = C_p * sin(lp.phi);
- V = lp.phi * lp.phi;
- lp.phi *= 0.895168 + V * ( 0.0218849 + V * 0.00826809 );
- for (i = NITER; i ; --i) {
- c = cos(lp.phi);
- s = sin(lp.phi);
- lp.phi -= V = (lp.phi + s * (c + 2.) - p) /
- (1. + c * (c + 2.) - s * s);
- if (fabs(V) < EPS)
- break;
- }
- if (!i) {
- xy.x = C_x * lp.lam;
- xy.y = lp.phi < 0. ? -C_y : C_y;
- } else {
- xy.x = C_x * lp.lam * (1. + cos(lp.phi));
- xy.y = C_y * sin(lp.phi);
- }
- return xy;
-}
-
-
-static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
- LP lp = {0.0,0.0};
- double c;
-
- lp.phi = aasin(P->ctx,xy.y / C_y);
- lp.lam = xy.x / (C_x * (1. + (c = cos(lp.phi))));
- lp.phi = aasin(P->ctx,(lp.phi + sin(lp.phi) * (c + 2.)) / C_p);
- return lp;
-}
-
-
-static void *freeup_new (PJ *P) { /* Destructor */
- if (0==P)
- return 0;
- return pj_dealloc(P);
-}
-
-static void freeup (PJ *P) {
- freeup_new (P);
- return;
-}
-
-
-PJ *PROJECTION(eck4) {
- P->es = 0.0;
- P->inv = s_inverse;
- P->fwd = s_forward;
-
- return P;
-}
-
-
-#ifdef PJ_OMIT_SELFTEST
-int pj_eck4_selftest (void) {return 0;}
-#else
-
-int pj_eck4_selftest (void) {
- double tolerance_lp = 1e-10;
- double tolerance_xy = 1e-7;
-
- char s_args[] = {"+proj=eck4 +a=6400000 +lat_1=0.5 +lat_2=2"};
-
- LP fwd_in[] = {
- { 2, 1},
- { 2,-1},
- {-2, 1},
- {-2,-1}
- };
-
- XY s_fwd_expect[] = {
- { 188646.38935641639, 132268.54017406539},
- { 188646.38935641639, -132268.54017406539},
- {-188646.38935641639, 132268.54017406539},
- {-188646.38935641639, -132268.54017406539},
- };
-
- XY inv_in[] = {
- { 200, 100},
- { 200,-100},
- {-200, 100},
- {-200,-100}
- };
-
- LP s_inv_expect[] = {
- { 0.0021202405520236059, 0.00075601458836610643},
- { 0.0021202405520236059, -0.00075601458836610643},
- {-0.0021202405520236059, 0.00075601458836610643},
- {-0.0021202405520236059, -0.00075601458836610643},
- };
-
- return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
-}
-
-
-#endif
+#define PJ_LIB__
+#include <projects.h>
+
+PROJ_HEAD(eck4, "Eckert IV") "\n\tPCyl, Sph.";
+
+#define C_x .42223820031577120149
+#define C_y 1.32650042817700232218
+#define RC_y .75386330736002178205
+#define C_p 3.57079632679489661922
+#define RC_p .28004957675577868795
+#define EPS 1e-7
+#define NITER 6
+
+
+static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
+ XY xy = {0.0,0.0};
+ double p, V, s, c;
+ int i;
+ (void) P;
+
+ p = C_p * sin(lp.phi);
+ V = lp.phi * lp.phi;
+ lp.phi *= 0.895168 + V * ( 0.0218849 + V * 0.00826809 );
+ for (i = NITER; i ; --i) {
+ c = cos(lp.phi);
+ s = sin(lp.phi);
+ lp.phi -= V = (lp.phi + s * (c + 2.) - p) /
+ (1. + c * (c + 2.) - s * s);
+ if (fabs(V) < EPS)
+ break;
+ }
+ if (!i) {
+ xy.x = C_x * lp.lam;
+ xy.y = lp.phi < 0. ? -C_y : C_y;
+ } else {
+ xy.x = C_x * lp.lam * (1. + cos(lp.phi));
+ xy.y = C_y * sin(lp.phi);
+ }
+ return xy;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
+ LP lp = {0.0,0.0};
+ double c;
+
+ lp.phi = aasin(P->ctx,xy.y / C_y);
+ lp.lam = xy.x / (C_x * (1. + (c = cos(lp.phi))));
+ lp.phi = aasin(P->ctx,(lp.phi + sin(lp.phi) * (c + 2.)) / C_p);
+ return lp;
+}
+
+
+static void *freeup_new (PJ *P) { /* Destructor */
+ if (0==P)
+ return 0;
+ return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+ freeup_new (P);
+ return;
+}
+
+
+PJ *PROJECTION(eck4) {
+ P->es = 0.0;
+ P->inv = s_inverse;
+ P->fwd = s_forward;
+
+ return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_eck4_selftest (void) {return 0;}
+#else
+
+int pj_eck4_selftest (void) {
+ double tolerance_lp = 1e-10;
+ double tolerance_xy = 1e-7;
+
+ char s_args[] = {"+proj=eck4 +a=6400000 +lat_1=0.5 +lat_2=2"};
+
+ LP fwd_in[] = {
+ { 2, 1},
+ { 2,-1},
+ {-2, 1},
+ {-2,-1}
+ };
+
+ XY s_fwd_expect[] = {
+ { 188646.38935641639, 132268.54017406539},
+ { 188646.38935641639, -132268.54017406539},
+ {-188646.38935641639, 132268.54017406539},
+ {-188646.38935641639, -132268.54017406539},
+ };
+
+ XY inv_in[] = {
+ { 200, 100},
+ { 200,-100},
+ {-200, 100},
+ {-200,-100}
+ };
+
+ LP s_inv_expect[] = {
+ { 0.0021202405520236059, 0.00075601458836610643},
+ { 0.0021202405520236059, -0.00075601458836610643},
+ {-0.0021202405520236059, 0.00075601458836610643},
+ {-0.0021202405520236059, -0.00075601458836610643},
+ };
+
+ return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
+}
+
+
+#endif
diff --git a/src/PJ_eck5.c b/src/PJ_eck5.c
index d7626939..4d1eeb17 100644
--- a/src/PJ_eck5.c
+++ b/src/PJ_eck5.c
@@ -1,93 +1,93 @@
-#define PJ_LIB__
-#include <projects.h>
-
-PROJ_HEAD(eck5, "Eckert V") "\n\tPCyl, Sph.";
-
-#define XF 0.44101277172455148219
-#define RXF 2.26750802723822639137
-#define YF 0.88202554344910296438
-#define RYF 1.13375401361911319568
-
-
-static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
- XY xy = {0.0,0.0};
- (void) P;
- xy.x = XF * (1. + cos(lp.phi)) * lp.lam;
- xy.y = YF * lp.phi;
-
- return xy;
-}
-
-
-static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
- LP lp = {0.0,0.0};
- (void) P;
- lp.lam = RXF * xy.x / (1. + cos( lp.phi = RYF * xy.y));
-
- return lp;
-}
-
-
-static void *freeup_new (PJ *P) { /* Destructor */
- if (0==P)
- return 0;
- return pj_dealloc(P);
-}
-
-static void freeup (PJ *P) {
- freeup_new (P);
- return;
-}
-
-
-PJ *PROJECTION(eck5) {
- P->es = 0.0;
- P->inv = s_inverse;
- P->fwd = s_forward;
-
- return P;
-}
-
-#ifdef PJ_OMIT_SELFTEST
-int pj_eck5_selftest (void) {return 0;}
-#else
-
-int pj_eck5_selftest (void) {
- double tolerance_lp = 1e-10;
- double tolerance_xy = 1e-7;
-
- char s_args[] = {"+proj=eck5 +a=6400000 +lat_1=0.5 +lat_2=2"};
-
- LP fwd_in[] = {
- { 2, 1},
- { 2,-1},
- {-2, 1},
- {-2,-1}
- };
-
- XY s_fwd_expect[] = {
- { 197031.39213406085, 98523.198847226551},
- { 197031.39213406085, -98523.198847226551},
- {-197031.39213406085, 98523.198847226551},
- {-197031.39213406085, -98523.198847226551},
- };
-
- XY inv_in[] = {
- { 200, 100},
- { 200,-100},
- {-200, 100},
- {-200,-100}
- };
-
- LP s_inv_expect[] = {
- {0.002029978749734037, 0.001014989374787388},
- {0.002029978749734037, -0.001014989374787388},
- {-0.002029978749734037, 0.001014989374787388},
- {-0.002029978749734037, -0.001014989374787388},
- };
-
- return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
-}
-
-
-#endif
+#define PJ_LIB__
+#include <projects.h>
+
+PROJ_HEAD(eck5, "Eckert V") "\n\tPCyl, Sph.";
+
+#define XF 0.44101277172455148219
+#define RXF 2.26750802723822639137
+#define YF 0.88202554344910296438
+#define RYF 1.13375401361911319568
+
+
+static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
+ XY xy = {0.0,0.0};
+ (void) P;
+ xy.x = XF * (1. + cos(lp.phi)) * lp.lam;
+ xy.y = YF * lp.phi;
+
+ return xy;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
+ LP lp = {0.0,0.0};
+ (void) P;
+ lp.lam = RXF * xy.x / (1. + cos( lp.phi = RYF * xy.y));
+
+ return lp;
+}
+
+
+static void *freeup_new (PJ *P) { /* Destructor */
+ if (0==P)
+ return 0;
+ return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+ freeup_new (P);
+ return;
+}
+
+
+PJ *PROJECTION(eck5) {
+ P->es = 0.0;
+ P->inv = s_inverse;
+ P->fwd = s_forward;
+
+ return P;
+}
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_eck5_selftest (void) {return 0;}
+#else
+
+int pj_eck5_selftest (void) {
+ double tolerance_lp = 1e-10;
+ double tolerance_xy = 1e-7;
+
+ char s_args[] = {"+proj=eck5 +a=6400000 +lat_1=0.5 +lat_2=2"};
+
+ LP fwd_in[] = {
+ { 2, 1},
+ { 2,-1},
+ {-2, 1},
+ {-2,-1}
+ };
+
+ XY s_fwd_expect[] = {
+ { 197031.39213406085, 98523.198847226551},
+ { 197031.39213406085, -98523.198847226551},
+ {-197031.39213406085, 98523.198847226551},
+ {-197031.39213406085, -98523.198847226551},
+ };
+
+ XY inv_in[] = {
+ { 200, 100},
+ { 200,-100},
+ {-200, 100},
+ {-200,-100}
+ };
+
+ LP s_inv_expect[] = {
+ {0.002029978749734037, 0.001014989374787388},
+ {0.002029978749734037, -0.001014989374787388},
+ {-0.002029978749734037, 0.001014989374787388},
+ {-0.002029978749734037, -0.001014989374787388},
+ };
+
+ return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
+}
+
+
+#endif
diff --git a/src/PJ_fahey.c b/src/PJ_fahey.c
index b8841c96..69fd40b7 100644
--- a/src/PJ_fahey.c
+++ b/src/PJ_fahey.c
@@ -1,95 +1,95 @@
-#define PJ_LIB__
-#include <projects.h>
-
-PROJ_HEAD(fahey, "Fahey") "\n\tPcyl, Sph.";
-
-#define TOL 1e-6
-
-
-static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
- XY xy = {0.0,0.0};
- (void) P;
-
- xy.x = tan(0.5 * lp.phi);
- xy.y = 1.819152 * xy.x;
- xy.x = 0.819152 * lp.lam * asqrt(1 - xy.x * xy.x);
- return xy;
-}
-
-
-static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
- LP lp = {0.0,0.0};
- (void) P;
-
- xy.y /= 1.819152;
- lp.phi = 2. * atan(xy.y);
- xy.y = 1. - xy.y * xy.y;
- lp.lam = fabs(xy.y) < TOL ? 0. : xy.x / (0.819152 * sqrt(xy.y));
- return lp;
-}
-
-
-static void *freeup_new (PJ *P) { /* Destructor */
- if (0==P)
- return 0;
- return pj_dealloc(P);
-}
-
-
-static void freeup (PJ *P) {
- freeup_new (P);
- return;
-}
-
-
-PJ *PROJECTION(fahey) {
- P->es = 0.;
- P->inv = s_inverse;
- P->fwd = s_forward;
-
- return P;
-}
-
-#ifdef PJ_OMIT_SELFTEST
-int pj_fahey_selftest (void) {return 0;}
-#else
-
-int pj_fahey_selftest (void) {
- double tolerance_lp = 1e-10;
- double tolerance_xy = 1e-7;
-
- char s_args[] = {"+proj=fahey +a=6400000 +lat_1=0.5 +lat_2=2"};
-
- LP fwd_in[] = {
- { 2, 1},
- { 2,-1},
- {-2, 1},
- {-2,-1}
- };
-
- XY s_fwd_expect[] = {
- { 182993.34464912376, 101603.19356988439},
- { 182993.34464912376, -101603.19356988439},
- {-182993.34464912376, 101603.19356988439},
- {-182993.34464912376, -101603.19356988439},
- };
-
- XY inv_in[] = {
- { 200, 100},
- { 200,-100},
- {-200, 100},
- {-200,-100}
- };
-
- LP s_inv_expect[] = {
- {0.0021857886080359551, 0.00098424601668238403},
- {0.0021857886080359551, -0.00098424601668238403},
- {-0.0021857886080359551, 0.00098424601668238403},
- {-0.0021857886080359551, -0.00098424601668238403},
- };
-
- return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
-}
-
-
-#endif
+#define PJ_LIB__
+#include <projects.h>
+
+PROJ_HEAD(fahey, "Fahey") "\n\tPcyl, Sph.";
+
+#define TOL 1e-6
+
+
+static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
+ XY xy = {0.0,0.0};
+ (void) P;
+
+ xy.x = tan(0.5 * lp.phi);
+ xy.y = 1.819152 * xy.x;
+ xy.x = 0.819152 * lp.lam * asqrt(1 - xy.x * xy.x);
+ return xy;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
+ LP lp = {0.0,0.0};
+ (void) P;
+
+ xy.y /= 1.819152;
+ lp.phi = 2. * atan(xy.y);
+ xy.y = 1. - xy.y * xy.y;
+ lp.lam = fabs(xy.y) < TOL ? 0. : xy.x / (0.819152 * sqrt(xy.y));
+ return lp;
+}
+
+
+static void *freeup_new (PJ *P) { /* Destructor */
+ if (0==P)
+ return 0;
+ return pj_dealloc(P);
+}
+
+
+static void freeup (PJ *P) {
+ freeup_new (P);
+ return;
+}
+
+
+PJ *PROJECTION(fahey) {
+ P->es = 0.;
+ P->inv = s_inverse;
+ P->fwd = s_forward;
+
+ return P;
+}
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_fahey_selftest (void) {return 0;}
+#else
+
+int pj_fahey_selftest (void) {
+ double tolerance_lp = 1e-10;
+ double tolerance_xy = 1e-7;
+
+ char s_args[] = {"+proj=fahey +a=6400000 +lat_1=0.5 +lat_2=2"};
+
+ LP fwd_in[] = {
+ { 2, 1},
+ { 2,-1},
+ {-2, 1},
+ {-2,-1}
+ };
+
+ XY s_fwd_expect[] = {
+ { 182993.34464912376, 101603.19356988439},
+ { 182993.34464912376, -101603.19356988439},
+ {-182993.34464912376, 101603.19356988439},
+ {-182993.34464912376, -101603.19356988439},
+ };
+
+ XY inv_in[] = {
+ { 200, 100},
+ { 200,-100},
+ {-200, 100},
+ {-200,-100}
+ };
+
+ LP s_inv_expect[] = {
+ {0.0021857886080359551, 0.00098424601668238403},
+ {0.0021857886080359551, -0.00098424601668238403},
+ {-0.0021857886080359551, 0.00098424601668238403},
+ {-0.0021857886080359551, -0.00098424601668238403},
+ };
+
+ return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
+}
+
+
+#endif
diff --git a/src/PJ_goode.c b/src/PJ_goode.c
index 48fc9ad5..d44bd2a1 100644
--- a/src/PJ_goode.c
+++ b/src/PJ_goode.c
@@ -1,122 +1,122 @@
-#define PJ_LIB__
-#include <projects.h>
-
-PROJ_HEAD(goode, "Goode Homolosine") "\n\tPCyl, Sph.";
-
-#define Y_COR 0.05280
-#define PHI_LIM 0.71093078197902358062
-
-C_NAMESPACE PJ *pj_sinu(PJ *), *pj_moll(PJ *);
-
-struct pj_opaque {
- PJ *sinu;
- PJ *moll;
-};
-
-
-static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
- XY xy = {0.0,0.0};
- struct pj_opaque *Q = P->opaque;
-
- if (fabs(lp.phi) <= PHI_LIM)
- xy = Q->sinu->fwd(lp, Q->sinu);
- else {
- xy = Q->moll->fwd(lp, Q->moll);
- xy.y -= lp.phi >= 0.0 ? Y_COR : -Y_COR;
- }
- return xy;
-}
-
-
-static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
- LP lp = {0.0,0.0};
- struct pj_opaque *Q = P->opaque;
-
- if (fabs(xy.y) <= PHI_LIM)
- lp = Q->sinu->inv(xy, Q->sinu);
- else {
- xy.y += xy.y >= 0.0 ? Y_COR : -Y_COR;
- lp = Q->moll->inv(xy, Q->moll);
- }
- return lp;
-}
-
-
-static void *freeup_new (PJ *P) { /* Destructor */
- if (0==P)
- return 0;
- if (0==P->opaque)
- return pj_dealloc(P);
- if (P->opaque->sinu)
- pj_dealloc(P->opaque->sinu);
- if (P->opaque->moll)
- pj_dealloc(P->opaque->moll);
- pj_dealloc (P->opaque);
- return pj_dealloc(P);
-
-}
-
-
-static void freeup (PJ *P) {
- freeup_new (P);
- return;
-}
-
-
-PJ *PROJECTION(goode) {
- struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
- if (0==Q)
- return freeup_new (P);
- P->opaque = Q;
-
- P->es = 0.;
- if (!(Q->sinu = pj_sinu(0)) || !(Q->moll = pj_moll(0)))
- E_ERROR_0;
- Q->sinu->es = 0.;
- Q->sinu->ctx = P->ctx;
- Q->moll->ctx = P->ctx;
- if (!(Q->sinu = pj_sinu(Q->sinu)) || !(Q->moll = pj_moll(Q->moll)))
- E_ERROR_0;
-
- P->fwd = s_forward;
- P->inv = s_inverse;
-
- return P;
-}
-
-
-#ifdef PJ_OMIT_SELFTEST
-int pj_goode_selftest (void) {return 0;}
-#else
-
-int pj_goode_selftest (void) {
- double tolerance_lp = 1e-10;
- double tolerance_xy = 1e-7;
-
- char s_args[] = {"+proj=goode +a=6400000 +lat_1=0.5 +lat_2=2"};
-
- LP fwd_in[] = {
- { 2, 1},
- { 2,-1},
- {-2, 1},
- {-2,-1}
- };
-
- XY s_fwd_expect[] = {
- { 223368.11902663155, 111701.07212763709}, { 223368.11902663155, -111701.07212763709}, {-223368.11902663155, 111701.07212763709}, {-223368.11902663155, -111701.07212763709}, };
-
- XY inv_in[] = {
- { 200, 100},
- { 200,-100},
- {-200, 100},
- {-200,-100}
- };
-
- LP s_inv_expect[] = {
- { 0.0017904931100023887, 0.00089524655489191132}, { 0.0017904931100023887, -0.00089524655489191132}, {-0.0017904931100023887, 0.00089524655489191132}, {-0.0017904931100023887, -0.00089524655489191132}, };
-
- return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
-}
-
-
-#endif
+#define PJ_LIB__
+#include <projects.h>
+
+PROJ_HEAD(goode, "Goode Homolosine") "\n\tPCyl, Sph.";
+
+#define Y_COR 0.05280
+#define PHI_LIM 0.71093078197902358062
+
+C_NAMESPACE PJ *pj_sinu(PJ *), *pj_moll(PJ *);
+
+struct pj_opaque {
+ PJ *sinu;
+ PJ *moll;
+};
+
+
+static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
+ XY xy = {0.0,0.0};
+ struct pj_opaque *Q = P->opaque;
+
+ if (fabs(lp.phi) <= PHI_LIM)
+ xy = Q->sinu->fwd(lp, Q->sinu);
+ else {
+ xy = Q->moll->fwd(lp, Q->moll);
+ xy.y -= lp.phi >= 0.0 ? Y_COR : -Y_COR;
+ }
+ return xy;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
+ LP lp = {0.0,0.0};
+ struct pj_opaque *Q = P->opaque;
+
+ if (fabs(xy.y) <= PHI_LIM)
+ lp = Q->sinu->inv(xy, Q->sinu);
+ else {
+ xy.y += xy.y >= 0.0 ? Y_COR : -Y_COR;
+ lp = Q->moll->inv(xy, Q->moll);
+ }
+ return lp;
+}
+
+
+static void *freeup_new (PJ *P) { /* Destructor */
+ if (0==P)
+ return 0;
+ if (0==P->opaque)
+ return pj_dealloc(P);
+ if (P->opaque->sinu)
+ pj_dealloc(P->opaque->sinu);
+ if (P->opaque->moll)
+ pj_dealloc(P->opaque->moll);
+ pj_dealloc (P->opaque);
+ return pj_dealloc(P);
+
+}
+
+
+static void freeup (PJ *P) {
+ freeup_new (P);
+ return;
+}
+
+
+PJ *PROJECTION(goode) {
+ struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+ if (0==Q)
+ return freeup_new (P);
+ P->opaque = Q;
+
+ P->es = 0.;
+ if (!(Q->sinu = pj_sinu(0)) || !(Q->moll = pj_moll(0)))
+ E_ERROR_0;
+ Q->sinu->es = 0.;
+ Q->sinu->ctx = P->ctx;
+ Q->moll->ctx = P->ctx;
+ if (!(Q->sinu = pj_sinu(Q->sinu)) || !(Q->moll = pj_moll(Q->moll)))
+ E_ERROR_0;
+
+ P->fwd = s_forward;
+ P->inv = s_inverse;
+
+ return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_goode_selftest (void) {return 0;}
+#else
+
+int pj_goode_selftest (void) {
+ double tolerance_lp = 1e-10;
+ double tolerance_xy = 1e-7;
+
+ char s_args[] = {"+proj=goode +a=6400000 +lat_1=0.5 +lat_2=2"};
+
+ LP fwd_in[] = {
+ { 2, 1},
+ { 2,-1},
+ {-2, 1},
+ {-2,-1}
+ };
+
+ XY s_fwd_expect[] = {
+ { 223368.11902663155, 111701.07212763709}, { 223368.11902663155, -111701.07212763709}, {-223368.11902663155, 111701.07212763709}, {-223368.11902663155, -111701.07212763709}, };
+
+ XY inv_in[] = {
+ { 200, 100},
+ { 200,-100},
+ {-200, 100},
+ {-200,-100}
+ };
+
+ LP s_inv_expect[] = {
+ { 0.0017904931100023887, 0.00089524655489191132}, { 0.0017904931100023887, -0.00089524655489191132}, {-0.0017904931100023887, 0.00089524655489191132}, {-0.0017904931100023887, -0.00089524655489191132}, };
+
+ return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
+}
+
+
+#endif
diff --git a/src/PJ_hatano.c b/src/PJ_hatano.c
index 34c70c0f..67973b17 100644
--- a/src/PJ_hatano.c
+++ b/src/PJ_hatano.c
@@ -1,134 +1,134 @@
-#define PJ_LIB__
-#include <projects.h>
-
-PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area") "\n\tPCyl, Sph.";
-
-#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
-
-
-static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
- XY xy = {0.0,0.0};
- 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;
-}
-
-
-static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
- LP lp = {0.0,0.0};
- double th;
-
- th = xy.y * ( xy.y < 0. ? RYCS : RYCN);
- if (fabs(th) > 1.) {
- if (fabs(th) > ONETOL) {
- I_ERROR;
- } else {
- th = th > 0. ? M_HALFPI : - M_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. ? M_HALFPI : - M_HALFPI;
- }
- } else {
- lp.phi = asin(lp.phi);
- }
-
- return (lp);
-}
-
-
-static void *freeup_new (PJ *P) { /* Destructor */
- if (0==P)
- return 0;
-
- return pj_dealloc(P);
-}
-
-static void freeup (PJ *P) {
- freeup_new (P);
- return;
-}
-
-
-PJ *PROJECTION(hatano) {
- P->es = 0.;
- P->inv = s_inverse;
- P->fwd = s_forward;
-
- return P;
-}
-
-#ifdef PJ_OMIT_SELFTEST
-int pj_hatano_selftest (void) {return 0;}
-#else
-
-int pj_hatano_selftest (void) {
- double tolerance_lp = 1e-10;
- double tolerance_xy = 1e-7;
-
- char s_args[] = {"+proj=hatano +a=6400000 +lat_1=0.5 +lat_2=2"};
-
- LP fwd_in[] = {
- { 2, 1},
- { 2,-1},
- {-2, 1},
- {-2,-1}
- };
-
- XY s_fwd_expect[] = {
- { 189878.87894652804, 131409.8024406255 },
- { 189881.08195244463, -131409.14227607418 },
- {-189878.87894652804, 131409.8024406255 },
- {-189881.08195244463, -131409.14227607418 },
- };
-
- XY inv_in[] = {
- { 200, 100},
- { 200,-100},
- {-200, 100},
- {-200,-100}
- };
-
- LP s_inv_expect[] = {
- { 0.0021064624821817597, 0.00076095689425791926 },
- { 0.0021064624821676096, -0.00076095777439265377 },
- {-0.0021064624821817597, 0.00076095689425791926 },
- {-0.0021064624821676096, -0.00076095777439265377 },
- };
-
- return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
-}
-
-
-#endif
+#define PJ_LIB__
+#include <projects.h>
+
+PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area") "\n\tPCyl, Sph.";
+
+#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
+
+
+static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
+ XY xy = {0.0,0.0};
+ 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;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
+ LP lp = {0.0,0.0};
+ double th;
+
+ th = xy.y * ( xy.y < 0. ? RYCS : RYCN);
+ if (fabs(th) > 1.) {
+ if (fabs(th) > ONETOL) {
+ I_ERROR;
+ } else {
+ th = th > 0. ? M_HALFPI : - M_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. ? M_HALFPI : - M_HALFPI;
+ }
+ } else {
+ lp.phi = asin(lp.phi);
+ }
+
+ return (lp);
+}
+
+
+static void *freeup_new (PJ *P) { /* Destructor */
+ if (0==P)
+ return 0;
+
+ return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+ freeup_new (P);
+ return;
+}
+
+
+PJ *PROJECTION(hatano) {
+ P->es = 0.;
+ P->inv = s_inverse;
+ P->fwd = s_forward;
+
+ return P;
+}
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_hatano_selftest (void) {return 0;}
+#else
+
+int pj_hatano_selftest (void) {
+ double tolerance_lp = 1e-10;
+ double tolerance_xy = 1e-7;
+
+ char s_args[] = {"+proj=hatano +a=6400000 +lat_1=0.5 +lat_2=2"};
+
+ LP fwd_in[] = {
+ { 2, 1},
+ { 2,-1},
+ {-2, 1},
+ {-2,-1}
+ };
+
+ XY s_fwd_expect[] = {
+ { 189878.87894652804, 131409.8024406255 },
+ { 189881.08195244463, -131409.14227607418 },
+ {-189878.87894652804, 131409.8024406255 },
+ {-189881.08195244463, -131409.14227607418 },
+ };
+
+ XY inv_in[] = {
+ { 200, 100},
+ { 200,-100},
+ {-200, 100},
+ {-200,-100}
+ };
+
+ LP s_inv_expect[] = {
+ { 0.0021064624821817597, 0.00076095689425791926 },
+ { 0.0021064624821676096, -0.00076095777439265377 },
+ {-0.0021064624821817597, 0.00076095689425791926 },
+ {-0.0021064624821676096, -0.00076095777439265377 },
+ };
+
+ return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
+}
+
+
+#endif
diff --git a/src/PJ_isea.c b/src/PJ_isea.c
index ba9d3c62..1eaec1f9 100644
--- a/src/PJ_isea.c
+++ b/src/PJ_isea.c
@@ -1,1176 +1,1176 @@
-/*
- * This code was entirely written by Nathan Wagner
- * and is in the public domain.
- */
-
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <float.h>
-
-#ifndef M_PI
-# define M_PI 3.14159265358979323846
-#endif
-
-/*
- * Proj 4 provides its own entry points into
- * the code, so none of the library functions
- * need to be global
- */
-#define ISEA_STATIC static
-#ifndef ISEA_STATIC
-#define ISEA_STATIC
-#endif
-
-struct hex {
- int iso;
- int x, y, z;
+/*
+ * This code was entirely written by Nathan Wagner
+ * and is in the public domain.
+ */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <float.h>
+
+#ifndef M_PI
+# define M_PI 3.14159265358979323846
+#endif
+
+/*
+ * Proj 4 provides its own entry points into
+ * the code, so none of the library functions
+ * need to be global
+ */
+#define ISEA_STATIC static
+#ifndef ISEA_STATIC
+#define ISEA_STATIC
+#endif
+
+struct hex {
+ int iso;
+ int x, y, z;
+};
+
+/* y *must* be positive down as the xy /iso conversion assumes this */
+ISEA_STATIC
+int hex_xy(struct hex *h) {
+ if (!h->iso) return 1;
+ if (h->x >= 0) {
+ h->y = -h->y - (h->x+1)/2;
+ } else {
+ /* need to round toward -inf, not toward zero, so x-1 */
+ h->y = -h->y - h->x/2;
+ }
+ h->iso = 0;
+
+ return 1;
+}
+
+ISEA_STATIC
+int hex_iso(struct hex *h) {
+ if (h->iso) return 1;
+
+ if (h->x >= 0) {
+ h->y = (-h->y - (h->x+1)/2);
+ } else {
+ /* need to round toward -inf, not toward zero, so x-1 */
+ h->y = (-h->y - (h->x)/2);
+ }
+
+ h->z = -h->x - h->y;
+ h->iso = 1;
+ return 1;
+}
+
+ISEA_STATIC
+int hexbin2(double width, double x, double y,
+ int *i, int *j) {
+ double z, rx, ry, rz;
+ double abs_dx, abs_dy, abs_dz;
+ int ix, iy, iz, s;
+ struct hex h;
+
+ x = x / cos(30 * M_PI / 180.0); /* rotated X coord */
+ y = y - x / 2.0; /* adjustment for rotated X */
+
+ /* adjust for actual hexwidth */
+ x /= width;
+ y /= width;
+
+ z = -x - y;
+
+ ix = rx = floor(x + 0.5);
+ iy = ry = floor(y + 0.5);
+ iz = rz = floor(z + 0.5);
+
+ s = ix + iy + iz;
+
+ if (s) {
+ abs_dx = fabs(rx - x);
+ abs_dy = fabs(ry - y);
+ abs_dz = fabs(rz - z);
+
+ if (abs_dx >= abs_dy && abs_dx >= abs_dz) {
+ ix -= s;
+ } else if (abs_dy >= abs_dx && abs_dy >= abs_dz) {
+ iy -= s;
+ } else {
+ iz -= s;
+ }
+ }
+ h.x = ix;
+ h.y = iy;
+ h.z = iz;
+ h.iso = 1;
+
+ hex_xy(&h);
+ *i = h.x;
+ *j = h.y;
+ return ix * 100 + iy;
+}
+#ifndef ISEA_STATIC
+#define ISEA_STATIC
+#endif
+
+enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = 20 };
+enum isea_topology { ISEA_HEXAGON=6, ISEA_TRIANGLE=3, ISEA_DIAMOND=4 };
+enum isea_address_form { ISEA_GEO, ISEA_Q2DI, ISEA_SEQNUM, ISEA_INTERLEAVE,
+ ISEA_PLANE, ISEA_Q2DD, ISEA_PROJTRI, ISEA_VERTEX2DD, ISEA_HEX
+};
+
+struct isea_dgg {
+ int polyhedron; /* ignored, icosahedron */
+ double o_lat, o_lon, o_az; /* orientation, radians */
+ int pole; /* true if standard snyder */
+ int topology; /* ignored, hexagon */
+ int aperture; /* valid values depend on partitioning method */
+ int resolution;
+ double radius; /* radius of the earth in meters, ignored 1.0 */
+ int output; /* an isea_address_form */
+ int triangle; /* triangle of last transformed point */
+ int quad; /* quad of last transformed point */
+ unsigned long serial;
};
-
-/* y *must* be positive down as the xy /iso conversion assumes this */
-ISEA_STATIC
-int hex_xy(struct hex *h) {
- if (!h->iso) return 1;
- if (h->x >= 0) {
- h->y = -h->y - (h->x+1)/2;
- } else {
- /* need to round toward -inf, not toward zero, so x-1 */
- h->y = -h->y - h->x/2;
- }
- h->iso = 0;
-
- return 1;
-}
-
-ISEA_STATIC
-int hex_iso(struct hex *h) {
- if (h->iso) return 1;
-
- if (h->x >= 0) {
- h->y = (-h->y - (h->x+1)/2);
- } else {
- /* need to round toward -inf, not toward zero, so x-1 */
- h->y = (-h->y - (h->x)/2);
- }
-
- h->z = -h->x - h->y;
- h->iso = 1;
- return 1;
-}
-
-ISEA_STATIC
-int hexbin2(double width, double x, double y,
- int *i, int *j) {
- double z, rx, ry, rz;
- double abs_dx, abs_dy, abs_dz;
- int ix, iy, iz, s;
- struct hex h;
-
- x = x / cos(30 * M_PI / 180.0); /* rotated X coord */
- y = y - x / 2.0; /* adjustment for rotated X */
-
- /* adjust for actual hexwidth */
- x /= width;
- y /= width;
-
- z = -x - y;
-
- ix = rx = floor(x + 0.5);
- iy = ry = floor(y + 0.5);
- iz = rz = floor(z + 0.5);
-
- s = ix + iy + iz;
-
- if (s) {
- abs_dx = fabs(rx - x);
- abs_dy = fabs(ry - y);
- abs_dz = fabs(rz - z);
-
- if (abs_dx >= abs_dy && abs_dx >= abs_dz) {
- ix -= s;
- } else if (abs_dy >= abs_dx && abs_dy >= abs_dz) {
- iy -= s;
- } else {
- iz -= s;
- }
- }
- h.x = ix;
- h.y = iy;
- h.z = iz;
- h.iso = 1;
-
- hex_xy(&h);
- *i = h.x;
- *j = h.y;
- return ix * 100 + iy;
-}
-#ifndef ISEA_STATIC
-#define ISEA_STATIC
-#endif
-
-enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = 20 };
-enum isea_topology { ISEA_HEXAGON=6, ISEA_TRIANGLE=3, ISEA_DIAMOND=4 };
-enum isea_address_form { ISEA_GEO, ISEA_Q2DI, ISEA_SEQNUM, ISEA_INTERLEAVE,
- ISEA_PLANE, ISEA_Q2DD, ISEA_PROJTRI, ISEA_VERTEX2DD, ISEA_HEX
-};
-
-struct isea_dgg {
- int polyhedron; /* ignored, icosahedron */
- double o_lat, o_lon, o_az; /* orientation, radians */
- int pole; /* true if standard snyder */
- int topology; /* ignored, hexagon */
- int aperture; /* valid values depend on partitioning method */
- int resolution;
- double radius; /* radius of the earth in meters, ignored 1.0 */
- int output; /* an isea_address_form */
- int triangle; /* triangle of last transformed point */
- int quad; /* quad of last transformed point */
- unsigned long serial;
-};
-
-struct isea_pt {
- double x, y;
-};
-
-struct isea_geo {
- double lon, lat;
-};
-
-struct isea_address {
- int type; /* enum isea_address_form */
- int number;
- double x,y; /* or i,j or lon,lat depending on type */
-};
-
-/* ENDINC */
-
-enum snyder_polyhedron {
- SNYDER_POLY_HEXAGON, SNYDER_POLY_PENTAGON,
- SNYDER_POLY_TETRAHEDRON, SNYDER_POLY_CUBE,
- SNYDER_POLY_OCTAHEDRON, SNYDER_POLY_DODECAHEDRON,
- SNYDER_POLY_ICOSAHEDRON
-};
-
-struct snyder_constants {
- double g, G, theta, ea_w, ea_a, ea_b, g_w, g_a, g_b;
-};
-
-/* TODO put these in radians to avoid a later conversion */
-ISEA_STATIC
-struct snyder_constants constants[] = {
- {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0},
- {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027},
- {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
- {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
- {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
- {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
- {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0},
-};
-
-#define E 52.62263186
-#define F 10.81231696
-
-#define DEG60 1.04719755119659774614
-#define DEG120 2.09439510239319549229
-#define DEG72 1.25663706143591729537
-#define DEG90 1.57079632679489661922
-#define DEG144 2.51327412287183459075
-#define DEG36 0.62831853071795864768
-#define DEG108 1.88495559215387594306
-#define DEG180 M_PI
-/* sqrt(5)/M_PI */
-#define ISEA_SCALE 0.8301572857837594396028083
-
-/* 26.565051177 degrees */
-#define V_LAT 0.46364760899944494524
-
-#define RAD2DEG (180.0/M_PI)
-#define DEG2RAD (M_PI/180.0)
-
-ISEA_STATIC
-struct isea_geo vertex[] = {
- {0.0, DEG90},
- {DEG180, V_LAT},
- {-DEG108, V_LAT},
- {-DEG36, V_LAT},
- {DEG36, V_LAT},
- {DEG108, V_LAT},
- {-DEG144, -V_LAT},
- {-DEG72, -V_LAT},
- {0.0, -V_LAT},
- {DEG72, -V_LAT},
- {DEG144, -V_LAT},
- {0.0, -DEG90}
-};
-
-/* TODO make an isea_pt array of the vertices as well */
-
-static int tri_v1[] = {0, 0, 0, 0, 0, 0, 6, 7, 8, 9, 10, 2, 3, 4, 5, 1, 11, 11, 11, 11, 11};
-
-/* 52.62263186 */
-#define E_RAD 0.91843818702186776133
-
-/* 10.81231696 */
-#define F_RAD 0.18871053072122403508
-
-/* triangle Centers */
-struct isea_geo icostriangles[] = {
- {0.0, 0.0},
- {-DEG144, E_RAD},
- {-DEG72, E_RAD},
- {0.0, E_RAD},
- {DEG72, E_RAD},
- {DEG144, E_RAD},
- {-DEG144, F_RAD},
- {-DEG72, F_RAD},
- {0.0, F_RAD},
- {DEG72, F_RAD},
- {DEG144, F_RAD},
- {-DEG108, -F_RAD},
- {-DEG36, -F_RAD},
- {DEG36, -F_RAD},
- {DEG108, -F_RAD},
- {DEG180, -F_RAD},
- {-DEG108, -E_RAD},
- {-DEG36, -E_RAD},
- {DEG36, -E_RAD},
- {DEG108, -E_RAD},
- {DEG180, -E_RAD},
-};
-
-static double
-az_adjustment(int triangle)
-{
- double adj;
-
- struct isea_geo v;
- struct isea_geo c;
-
- v = vertex[tri_v1[triangle]];
- c = icostriangles[triangle];
-
- /* TODO looks like the adjustment is always either 0 or 180 */
- /* at least if you pick your vertex carefully */
- adj = atan2(cos(v.lat) * sin(v.lon - c.lon),
- cos(c.lat) * sin(v.lat)
- - sin(c.lat) * cos(v.lat) * cos(v.lon - c.lon));
- return adj;
-}
-
-/* R tan(g) sin(60) */
-#define TABLE_G 0.6615845383
-
-/* H = 0.25 R tan g = */
-#define TABLE_H 0.1909830056
-
-#define RPRIME 0.91038328153090290025
-
-ISEA_STATIC
-struct isea_pt
-isea_triangle_xy(int triangle)
-{
- struct isea_pt c;
- double Rprime = 0.91038328153090290025;
-
- triangle = (triangle - 1) % 20;
-
- c.x = TABLE_G * ((triangle % 5) - 2) * 2.0;
- if (triangle > 9) {
- c.x += TABLE_G;
- }
- switch (triangle / 5) {
- case 0:
- c.y = 5.0 * TABLE_H;
- break;
- case 1:
- c.y = TABLE_H;
- break;
- case 2:
- c.y = -TABLE_H;
- break;
- case 3:
- c.y = -5.0 * TABLE_H;
- break;
- default:
- /* should be impossible */
- exit(EXIT_FAILURE);
- };
- c.x *= Rprime;
- c.y *= Rprime;
-
- return c;
-}
-
-/* snyder eq 14 */
-static double
-sph_azimuth(double f_lon, double f_lat, double t_lon, double t_lat)
-{
- double az;
-
- az = atan2(cos(t_lat) * sin(t_lon - f_lon),
- cos(f_lat) * sin(t_lat)
- - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon)
- );
- return az;
-}
-
-/* coord needs to be in radians */
-ISEA_STATIC
-int
-isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)
-{
- int i;
-
- /*
- * spherical distance from center of polygon face to any of its
- * vertexes on the globe
- */
- double g;
-
- /*
- * spherical angle between radius vector to center and adjacent edge
- * of spherical polygon on the globe
- */
- double G;
-
- /*
- * plane angle between radius vector to center and adjacent edge of
- * plane polygon
- */
- double theta;
-
- /* additional variables from snyder */
- double q, Rprime, H, Ag, Azprime, Az, dprime, f, rho,
- x, y;
-
- /* variables used to store intermediate results */
- double cot_theta, tan_g, az_offset;
-
- /* how many multiples of 60 degrees we adjust the azimuth */
- int Az_adjust_multiples;
-
- struct snyder_constants c;
-
- /*
- * TODO by locality of reference, start by trying the same triangle
- * as last time
- */
-
- /* TODO put these constants in as radians to begin with */
- c = constants[SNYDER_POLY_ICOSAHEDRON];
- theta = c.theta * DEG2RAD;
- g = c.g * DEG2RAD;
- G = c.G * DEG2RAD;
-
- for (i = 1; i <= 20; i++) {
- double z;
- struct isea_geo center;
-
- center = icostriangles[i];
-
- /* step 1 */
- z = acos(sin(center.lat) * sin(ll->lat)
- + cos(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon));
- /* not on this triangle */
- if (z > g + 0.000005) { /* TODO DBL_EPSILON */
- continue;
- }
-
- Az = sph_azimuth(center.lon, center.lat, ll->lon, ll->lat);
-
- /* step 2 */
-
- /* This calculates "some" vertex coordinate */
- az_offset = az_adjustment(i);
-
- Az -= az_offset;
-
- /* TODO I don't know why we do this. It's not in snyder */
- /* maybe because we should have picked a better vertex */
- if (Az < 0.0) {
- Az += 2.0 * M_PI;
- }
- /*
- * adjust Az for the point to fall within the range of 0 to
- * 2(90 - theta) or 60 degrees for the hexagon, by
- * and therefore 120 degrees for the triangle
- * of the icosahedron
- * subtracting or adding multiples of 60 degrees to Az and
- * recording the amount of adjustment
- */
-
- Az_adjust_multiples = 0;
- while (Az < 0.0) {
- Az += DEG120;
- Az_adjust_multiples--;
- }
- while (Az > DEG120 + DBL_EPSILON) {
- Az -= DEG120;
- Az_adjust_multiples++;
- }
-
- /* step 3 */
- cot_theta = 1.0 / tan(theta);
- tan_g = tan(g); /* TODO this is a constant */
-
- /* Calculate q from eq 9. */
- /* TODO cot_theta is cot(30) */
- q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta);
-
- /* not in this triangle */
- if (z > q + 0.000005) {
- continue;
- }
- /* step 4 */
-
- /* Apply equations 5-8 and 10-12 in order */
-
- /* eq 5 */
- /* Rprime = 0.9449322893 * R; */
- /* R' in the paper is for the truncated */
- Rprime = 0.91038328153090290025;
-
- /* eq 6 */
- H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G));
-
- /* eq 7 */
- /* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */
- Ag = Az + G + H - DEG180;
-
- /* eq 8 */
- Azprime = atan2(2.0 * Ag, Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta);
-
- /* eq 10 */
- /* cot(theta) = 1.73205080756887729355 */
- dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta);
-
- /* eq 11 */
- f = dprime / (2.0 * Rprime * sin(q / 2.0));
-
- /* eq 12 */
- rho = 2.0 * Rprime * f * sin(z / 2.0);
-
- /*
- * add back the same 60 degree multiple adjustment from step
- * 2 to Azprime
- */
-
- Azprime += DEG120 * Az_adjust_multiples;
-
- /* calculate rectangular coordinates */
-
- x = rho * sin(Azprime);
- y = rho * cos(Azprime);
-
- /*
- * TODO
- * translate coordinates to the origin for the particular
- * hexagon on the flattened polyhedral map plot
- */
-
- out->x = x;
- out->y = y;
-
- return i;
- }
-
- /*
- * should be impossible, this implies that the coordinate is not on
- * any triangle
- */
-
- fprintf(stderr, "impossible transform: %f %f is not on any triangle\n",
- ll->lon * RAD2DEG, ll->lat * RAD2DEG);
-
- exit(EXIT_FAILURE);
-
- /* not reached */
- return 0; /* supresses a warning */
-}
-
-/*
- * return the new coordinates of any point in orginal coordinate system.
- * Define a point (newNPold) in orginal coordinate system as the North Pole in
- * new coordinate system, and the great circle connect the original and new
- * North Pole as the lon0 longitude in new coordinate system, given any point
- * in orginal coordinate system, this function return the new coordinates.
- */
-
-#define PRECISION 0.0000000000005
-
-/* formula from Snyder, Map Projections: A working manual, p31 */
-/*
- * old north pole at np in new coordinates
- * could be simplified a bit with fewer intermediates
- *
- * TODO take a result pointer
- */
-ISEA_STATIC
-struct isea_geo
-snyder_ctran(struct isea_geo * np, struct isea_geo * pt)
-{
- struct isea_geo npt;
- double alpha, phi, lambda, lambda0, beta, lambdap, phip;
- double sin_phip;
- double lp_b; /* lambda prime minus beta */
- double cos_p, sin_a;
-
- phi = pt->lat;
- lambda = pt->lon;
- alpha = np->lat;
- beta = np->lon;
- lambda0 = beta;
-
- cos_p = cos(phi);
- sin_a = sin(alpha);
-
- /* mpawm 5-7 */
- sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0);
-
- /* mpawm 5-8b */
-
- /* use the two argument form so we end up in the right quadrant */
- lp_b = atan2(cos_p * sin(lambda - lambda0),
- (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));
-
- lambdap = lp_b + beta;
-
- /* normalize longitude */
- /* TODO can we just do a modulus ? */
- lambdap = fmod(lambdap, 2 * M_PI);
- while (lambdap > M_PI)
- lambdap -= 2 * M_PI;
- while (lambdap < -M_PI)
- lambdap += 2 * M_PI;
-
- phip = asin(sin_phip);
-
- npt.lat = phip;
- npt.lon = lambdap;
-
- return npt;
-}
-
-ISEA_STATIC
-struct isea_geo
-isea_ctran(struct isea_geo * np, struct isea_geo * pt, double lon0)
-{
- struct isea_geo npt;
-
- np->lon += M_PI;
- npt = snyder_ctran(np, pt);
- np->lon -= M_PI;
-
- npt.lon -= (M_PI - lon0 + np->lon);
-
- /*
- * snyder is down tri 3, isea is along side of tri1 from vertex 0 to
- * vertex 1 these are 180 degrees apart
- */
- npt.lon += M_PI;
- /* normalize longitude */
- npt.lon = fmod(npt.lon, 2 * M_PI);
- while (npt.lon > M_PI)
- npt.lon -= 2 * M_PI;
- while (npt.lon < -M_PI)
- npt.lon += 2 * M_PI;
-
- return npt;
-}
-
-/* in radians */
-#define ISEA_STD_LAT 1.01722196792335072101
-#define ISEA_STD_LON .19634954084936207740
-
-/* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */
-
-ISEA_STATIC
-int
-isea_grid_init(struct isea_dgg * g)
-{
- if (!g)
- return 0;
-
- g->polyhedron = 20;
- g->o_lat = ISEA_STD_LAT;
- g->o_lon = ISEA_STD_LON;
- g->o_az = 0.0;
- g->aperture = 4;
- g->resolution = 6;
- g->radius = 1.0;
- g->topology = 6;
-
- return 1;
-}
-
-ISEA_STATIC
-int
-isea_orient_isea(struct isea_dgg * g)
-{
- if (!g)
- return 0;
- g->o_lat = ISEA_STD_LAT;
- g->o_lon = ISEA_STD_LON;
- g->o_az = 0.0;
- return 1;
-}
-
-ISEA_STATIC
-int
-isea_orient_pole(struct isea_dgg * g)
-{
- if (!g)
- return 0;
- g->o_lat = M_PI / 2.0;
- g->o_lon = 0.0;
- g->o_az = 0;
- return 1;
-}
-
-ISEA_STATIC
-int
-isea_transform(struct isea_dgg * g, struct isea_geo * in,
- struct isea_pt * out)
-{
- struct isea_geo i, pole;
- int tri;
-
- pole.lat = g->o_lat;
- pole.lon = g->o_lon;
-
- i = isea_ctran(&pole, in, g->o_az);
-
- tri = isea_snyder_forward(&i, out);
- out->x *= g->radius;
- out->y *= g->radius;
- g->triangle = tri;
-
- return tri;
-}
-
-#define DOWNTRI(tri) (((tri - 1) / 5) % 2 == 1)
-
-ISEA_STATIC
-void
-isea_rotate(struct isea_pt * pt, double degrees)
-{
- double rad;
-
- double x, y;
-
- rad = -degrees * M_PI / 180.0;
- while (rad >= 2.0 * M_PI) rad -= 2.0 * M_PI;
- while (rad <= -2.0 * M_PI) rad += 2.0 * M_PI;
-
- x = pt->x * cos(rad) + pt->y * sin(rad);
- y = -pt->x * sin(rad) + pt->y * cos(rad);
-
- pt->x = x;
- pt->y = y;
-}
-
-ISEA_STATIC
-int isea_tri_plane(int tri, struct isea_pt *pt, double radius) {
- struct isea_pt tc; /* center of triangle */
-
- if (DOWNTRI(tri)) {
- isea_rotate(pt, 180.0);
- }
- tc = isea_triangle_xy(tri);
- tc.x *= radius;
- tc.y *= radius;
- pt->x += tc.x;
- pt->y += tc.y;
-
- return tri;
-}
-
-/* convert projected triangle coords to quad xy coords, return quad number */
-ISEA_STATIC
-int
-isea_ptdd(int tri, struct isea_pt *pt) {
- int downtri, quad;
-
- downtri = (((tri - 1) / 5) % 2 == 1);
- quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
-
- isea_rotate(pt, downtri ? 240.0 : 60.0);
- if (downtri) {
- pt->x += 0.5;
- /* pt->y += cos(30.0 * M_PI / 180.0); */
- pt->y += .86602540378443864672;
- }
- return quad;
-}
-
-ISEA_STATIC
-int
-isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di)
-{
- struct isea_pt v;
- double hexwidth;
- double sidelength; /* in hexes */
- int d, i;
- int maxcoord;
- struct hex h;
-
- /* This is the number of hexes from apex to base of a triangle */
- sidelength = (pow(2.0, g->resolution) + 1.0) / 2.0;
-
- /* apex to base is cos(30deg) */
- hexwidth = cos(M_PI / 6.0) / sidelength;
-
- /* TODO I think sidelength is always x.5, so
- * (int)sidelength * 2 + 1 might be just as good
- */
- maxcoord = (int) (sidelength * 2.0 + 0.5);
-
- v = *pt;
- hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);
- h.iso = 0;
- hex_iso(&h);
-
- d = h.x - h.z;
- i = h.x + h.y + h.y;
-
- /*
- * you want to test for max coords for the next quad in the same
- * "row" first to get the case where both are max
- */
- if (quad <= 5) {
- if (d == 0 && i == maxcoord) {
- /* north pole */
- quad = 0;
- d = 0;
- i = 0;
- } else if (i == maxcoord) {
- /* upper right in next quad */
- quad += 1;
- if (quad == 6)
- quad = 1;
- i = maxcoord - d;
- d = 0;
- } else if (d == maxcoord) {
- /* lower right in quad to lower right */
- quad += 5;
- d = 0;
- }
- } else if (quad >= 6) {
- if (i == 0 && d == maxcoord) {
- /* south pole */
- quad = 11;
- d = 0;
- i = 0;
- } else if (d == maxcoord) {
- /* lower right in next quad */
- quad += 1;
- if (quad == 11)
- quad = 6;
- d = maxcoord - i;
- i = 0;
- } else if (i == maxcoord) {
- /* upper right in quad to upper right */
- quad = (quad - 4) % 5;
- i = 0;
- }
- }
-
- di->x = d;
- di->y = i;
-
- g->quad = quad;
- return quad;
-}
-
-ISEA_STATIC
-int
-isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di) {
- struct isea_pt v;
- double hexwidth;
- int sidelength; /* in hexes */
- struct hex h;
-
- if (g->aperture == 3 && g->resolution % 2 != 0) {
- return isea_dddi_ap3odd(g, quad, pt, di);
- }
- /* todo might want to do this as an iterated loop */
- if (g->aperture >0) {
- sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
- } else {
- sidelength = g->resolution;
- }
-
- hexwidth = 1.0 / sidelength;
-
- v = *pt;
- isea_rotate(&v, -30.0);
- hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);
- h.iso = 0;
- hex_iso(&h);
-
- /* we may actually be on another quad */
- if (quad <= 5) {
- if (h.x == 0 && h.z == -sidelength) {
- /* north pole */
- quad = 0;
- h.z = 0;
- h.y = 0;
- h.x = 0;
- } else if (h.z == -sidelength) {
- quad = quad + 1;
- if (quad == 6)
- quad = 1;
- h.y = sidelength - h.x;
- h.z = h.x - sidelength;
- h.x = 0;
- } else if (h.x == sidelength) {
- quad += 5;
- h.y = -h.z;
- h.x = 0;
- }
- } else if (quad >= 6) {
- if (h.z == 0 && h.x == sidelength) {
- /* south pole */
- quad = 11;
- h.x = 0;
- h.y = 0;
- h.z = 0;
- } else if (h.x == sidelength) {
- quad = quad + 1;
- if (quad == 11)
- quad = 6;
- h.x = h.y + sidelength;
- h.y = 0;
- h.z = -h.x;
- } else if (h.y == -sidelength) {
- quad -= 4;
- h.y = 0;
- h.z = -h.x;
- }
- }
- di->x = h.x;
- di->y = -h.z;
-
- g->quad = quad;
- return quad;
-}
-
-ISEA_STATIC
-int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,
- struct isea_pt *di) {
- struct isea_pt v;
- int quad;
-
- v = *pt;
- quad = isea_ptdd(tri, &v);
- quad = isea_dddi(g, quad, &v, di);
- return quad;
-}
-
-/* q2di to seqnum */
-ISEA_STATIC
-int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
- int sidelength;
- int sn, height;
- int hexes;
-
- if (quad == 0) {
- g->serial = 1;
- return g->serial;
- }
- /* hexes in a quad */
- hexes = (int) (pow(g->aperture, g->resolution) + 0.5);
- if (quad == 11) {
- g->serial = 1 + 10 * hexes + 1;
- return g->serial;
- }
- if (g->aperture == 3 && g->resolution % 2 == 1) {
- height = (int) (pow(g->aperture, (g->resolution - 1) / 2.0));
- sn = ((int) di->x) * height;
- sn += ((int) di->y) / height;
- sn += (quad - 1) * hexes;
- sn += 2;
- } else {
- sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
- sn = (quad - 1) * hexes + sidelength * di->x + di->y + 2;
- }
-
- g->serial = sn;
- return sn;
-}
-
-/* TODO just encode the quad in the d or i coordinate
- * quad is 0-11, which can be four bits.
- * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf
- */
-/* convert a q2di to global hex coord */
-ISEA_STATIC
-int isea_hex(struct isea_dgg *g, int tri,
- struct isea_pt *pt, struct isea_pt *hex) {
- struct isea_pt v;
- int sidelength;
- int d, i, x, y, quad;
-
- quad = isea_ptdi(g, tri, pt, &v);
-
- hex->x = ((int)v.x << 4) + quad;
- hex->y = v.y;
-
- return 1;
-
- d = v.x;
- i = v.y;
-
- /* Aperture 3 odd resolutions */
- if (g->aperture == 3 && g->resolution % 2 != 0) {
- int offset = (int)(pow(3.0, g->resolution - 1) + 0.5);
-
- d += offset * ((g->quad-1) % 5);
- i += offset * ((g->quad-1) % 5);
-
- if (quad == 0) {
- d = 0;
- i = offset;
- } else if (quad == 11) {
- d = 2 * offset;
- i = 0;
- } else if (quad > 5) {
- d += offset;
- }
-
- x = (2*d - i) /3;
- y = (2*i - d) /3;
-
- hex->x = x + offset / 3;
- hex->y = y + 2 * offset / 3;
- return 1;
- }
-
- /* aperture 3 even resolutions and aperture 4 */
- sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
- if (g->quad == 0) {
- hex->x = 0;
- hex->y = sidelength;
- } else if (g->quad == 11) {
- hex->x = sidelength * 2;
- hex->y = 0;
- } else {
- hex->x = d + sidelength * ((g->quad-1) % 5);
- if (g->quad > 5) hex->x += sidelength;
- hex->y = i + sidelength * ((g->quad-1) % 5);
- }
-
- return 1;
-}
-
-ISEA_STATIC
-struct isea_pt
-isea_forward(struct isea_dgg *g, struct isea_geo *in)
-{
- int tri;
- struct isea_pt out, coord;
-
- tri = isea_transform(g, in, &out);
-
- if (g->output == ISEA_PLANE) {
- isea_tri_plane(tri, &out, g->radius);
- return out;
- }
-
- /* convert to isea standard triangle size */
- out.x = out.x / g->radius * ISEA_SCALE;
- out.y = out.y / g->radius * ISEA_SCALE;
- out.x += 0.5;
- out.y += 2.0 * .14433756729740644112;
-
- switch (g->output) {
- case ISEA_PROJTRI:
- /* nothing to do, already in projected triangle */
- break;
- case ISEA_VERTEX2DD:
- g->quad = isea_ptdd(tri, &out);
- break;
- case ISEA_Q2DD:
- /* Same as above, we just don't print as much */
- g->quad = isea_ptdd(tri, &out);
- break;
- case ISEA_Q2DI:
- g->quad = isea_ptdi(g, tri, &out, &coord);
- return coord;
- break;
- case ISEA_SEQNUM:
- isea_ptdi(g, tri, &out, &coord);
- /* disn will set g->serial */
- isea_disn(g, g->quad, &coord);
- return coord;
- break;
- case ISEA_HEX:
- isea_hex(g, tri, &out, &coord);
- return coord;
- break;
- }
-
- return out;
-}
-/*
- * Proj 4 integration code follows
- */
-
-#define PJ_LIB__
-#include <projects.h>
-
-PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph";
-
-struct pj_opaque {
- struct isea_dgg dgg;
-};
-
-
-static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
- XY xy = {0.0,0.0};
- struct pj_opaque *Q = P->opaque;
- struct isea_pt out;
- struct isea_geo in;
-
- in.lon = lp.lam;
- in.lat = lp.phi;
-
- out = isea_forward(&Q->dgg, &in);
-
- xy.x = out.x;
- xy.y = out.y;
-
- return xy;
-}
-
-
-static void *freeup_new (PJ *P) { /* Destructor */
- if (0==P)
- return 0;
- if (0==P->opaque)
- return pj_dealloc (P);
-
- pj_dealloc (P->opaque);
- return pj_dealloc(P);
-}
-
-static void freeup (PJ *P) {
- freeup_new (P);
- return;
-}
-
-
-PJ *PROJECTION(isea) {
- char *opt;
- struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
- if (0==Q)
- return freeup_new (P);
- P->opaque = Q;
-
-
- P->fwd = s_forward;
- isea_grid_init(&Q->dgg);
-
- Q->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(&Q->dgg);
- } else if (!strcmp(opt, "pole")) {
- isea_orient_pole(&Q->dgg);
- } else {
- E_ERROR(-34);
- }
- }
-
- if (pj_param(P->ctx,P->params, "tazi").i) {
- Q->dgg.o_az = pj_param(P->ctx,P->params, "razi").f;
- }
-
- if (pj_param(P->ctx,P->params, "tlon_0").i) {
- Q->dgg.o_lon = pj_param(P->ctx,P->params, "rlon_0").f;
- }
-
- if (pj_param(P->ctx,P->params, "tlat_0").i) {
- Q->dgg.o_lat = pj_param(P->ctx,P->params, "rlat_0").f;
- }
-
- if (pj_param(P->ctx,P->params, "taperture").i) {
- Q->dgg.aperture = pj_param(P->ctx,P->params, "iaperture").i;
- }
-
- if (pj_param(P->ctx,P->params, "tresolution").i) {
- Q->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")) {
- Q->dgg.output = ISEA_PLANE;
- } else if (!strcmp(opt, "di")) {
- Q->dgg.output = ISEA_Q2DI;
- }
- else if (!strcmp(opt, "dd")) {
- Q->dgg.output = ISEA_Q2DD;
- }
- else if (!strcmp(opt, "hex")) {
- Q->dgg.output = ISEA_HEX;
- }
- else {
- /* TODO verify error code. Possibly eliminate magic */
- E_ERROR(-34);
- }
- }
-
- if (pj_param(P->ctx,P->params, "trescale").i) {
- Q->dgg.radius = ISEA_SCALE;
- }
-
- if (pj_param(P->ctx,P->params, "tresolution").i) {
- Q->dgg.resolution = pj_param(P->ctx,P->params, "iresolution").i;
- } else {
- Q->dgg.resolution = 4;
- }
-
- if (pj_param(P->ctx,P->params, "taperture").i) {
- Q->dgg.aperture = pj_param(P->ctx,P->params, "iaperture").i;
- } else {
- Q->dgg.aperture = 3;
- }
-
- return P;
-}
-
-
-#ifdef PJ_OMIT_SELFTEST
-int pj_isea_selftest (void) {return 0;}
-#else
-
-int pj_isea_selftest (void) {
- double tolerance_lp = 1e-10;
- double tolerance_xy = 1e-7;
-
- char s_args[] = {"+proj=isea +a=6400000 +lat_1=0.5 +lat_2=2"};
-
- LP fwd_in[] = {
- { 2, 1},
- { 2,-1},
- {-2, 1},
- {-2,-1}
- };
-
- XY s_fwd_expect[] = {
- {-1097074.9480224741, 3442909.3090371834},
- {-1097074.9482647954, 3233611.7285857084},
- {-1575486.3536415542, 3442168.3420281881},
- {-1575486.353880283, 3234352.6955947056},
- };
-
- return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, 0, 0, 0);
-}
-
-
-#endif
+
+struct isea_pt {
+ double x, y;
+};
+
+struct isea_geo {
+ double lon, lat;
+};
+
+struct isea_address {
+ int type; /* enum isea_address_form */
+ int number;
+ double x,y; /* or i,j or lon,lat depending on type */
+};
+
+/* ENDINC */
+
+enum snyder_polyhedron {
+ SNYDER_POLY_HEXAGON, SNYDER_POLY_PENTAGON,
+ SNYDER_POLY_TETRAHEDRON, SNYDER_POLY_CUBE,
+ SNYDER_POLY_OCTAHEDRON, SNYDER_POLY_DODECAHEDRON,
+ SNYDER_POLY_ICOSAHEDRON
+};
+
+struct snyder_constants {
+ double g, G, theta, ea_w, ea_a, ea_b, g_w, g_a, g_b;
+};
+
+/* TODO put these in radians to avoid a later conversion */
+ISEA_STATIC
+struct snyder_constants constants[] = {
+ {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0},
+ {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027},
+ {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+ {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0},
+};
+
+#define E 52.62263186
+#define F 10.81231696
+
+#define DEG60 1.04719755119659774614
+#define DEG120 2.09439510239319549229
+#define DEG72 1.25663706143591729537
+#define DEG90 1.57079632679489661922
+#define DEG144 2.51327412287183459075
+#define DEG36 0.62831853071795864768
+#define DEG108 1.88495559215387594306
+#define DEG180 M_PI
+/* sqrt(5)/M_PI */
+#define ISEA_SCALE 0.8301572857837594396028083
+
+/* 26.565051177 degrees */
+#define V_LAT 0.46364760899944494524
+
+#define RAD2DEG (180.0/M_PI)
+#define DEG2RAD (M_PI/180.0)
+
+ISEA_STATIC
+struct isea_geo vertex[] = {
+ {0.0, DEG90},
+ {DEG180, V_LAT},
+ {-DEG108, V_LAT},
+ {-DEG36, V_LAT},
+ {DEG36, V_LAT},
+ {DEG108, V_LAT},
+ {-DEG144, -V_LAT},
+ {-DEG72, -V_LAT},
+ {0.0, -V_LAT},
+ {DEG72, -V_LAT},
+ {DEG144, -V_LAT},
+ {0.0, -DEG90}
+};
+
+/* TODO make an isea_pt array of the vertices as well */
+
+static int tri_v1[] = {0, 0, 0, 0, 0, 0, 6, 7, 8, 9, 10, 2, 3, 4, 5, 1, 11, 11, 11, 11, 11};
+
+/* 52.62263186 */
+#define E_RAD 0.91843818702186776133
+
+/* 10.81231696 */
+#define F_RAD 0.18871053072122403508
+
+/* triangle Centers */
+struct isea_geo icostriangles[] = {
+ {0.0, 0.0},
+ {-DEG144, E_RAD},
+ {-DEG72, E_RAD},
+ {0.0, E_RAD},
+ {DEG72, E_RAD},
+ {DEG144, E_RAD},
+ {-DEG144, F_RAD},
+ {-DEG72, F_RAD},
+ {0.0, F_RAD},
+ {DEG72, F_RAD},
+ {DEG144, F_RAD},
+ {-DEG108, -F_RAD},
+ {-DEG36, -F_RAD},
+ {DEG36, -F_RAD},
+ {DEG108, -F_RAD},
+ {DEG180, -F_RAD},
+ {-DEG108, -E_RAD},
+ {-DEG36, -E_RAD},
+ {DEG36, -E_RAD},
+ {DEG108, -E_RAD},
+ {DEG180, -E_RAD},
+};
+
+static double
+az_adjustment(int triangle)
+{
+ double adj;
+
+ struct isea_geo v;
+ struct isea_geo c;
+
+ v = vertex[tri_v1[triangle]];
+ c = icostriangles[triangle];
+
+ /* TODO looks like the adjustment is always either 0 or 180 */
+ /* at least if you pick your vertex carefully */
+ adj = atan2(cos(v.lat) * sin(v.lon - c.lon),
+ cos(c.lat) * sin(v.lat)
+ - sin(c.lat) * cos(v.lat) * cos(v.lon - c.lon));
+ return adj;
+}
+
+/* R tan(g) sin(60) */
+#define TABLE_G 0.6615845383
+
+/* H = 0.25 R tan g = */
+#define TABLE_H 0.1909830056
+
+#define RPRIME 0.91038328153090290025
+
+ISEA_STATIC
+struct isea_pt
+isea_triangle_xy(int triangle)
+{
+ struct isea_pt c;
+ double Rprime = 0.91038328153090290025;
+
+ triangle = (triangle - 1) % 20;
+
+ c.x = TABLE_G * ((triangle % 5) - 2) * 2.0;
+ if (triangle > 9) {
+ c.x += TABLE_G;
+ }
+ switch (triangle / 5) {
+ case 0:
+ c.y = 5.0 * TABLE_H;
+ break;
+ case 1:
+ c.y = TABLE_H;
+ break;
+ case 2:
+ c.y = -TABLE_H;
+ break;
+ case 3:
+ c.y = -5.0 * TABLE_H;
+ break;
+ default:
+ /* should be impossible */
+ exit(EXIT_FAILURE);
+ };
+ c.x *= Rprime;
+ c.y *= Rprime;
+
+ return c;
+}
+
+/* snyder eq 14 */
+static double
+sph_azimuth(double f_lon, double f_lat, double t_lon, double t_lat)
+{
+ double az;
+
+ az = atan2(cos(t_lat) * sin(t_lon - f_lon),
+ cos(f_lat) * sin(t_lat)
+ - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon)
+ );
+ return az;
+}
+
+/* coord needs to be in radians */
+ISEA_STATIC
+int
+isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)
+{
+ int i;
+
+ /*
+ * spherical distance from center of polygon face to any of its
+ * vertexes on the globe
+ */
+ double g;
+
+ /*
+ * spherical angle between radius vector to center and adjacent edge
+ * of spherical polygon on the globe
+ */
+ double G;
+
+ /*
+ * plane angle between radius vector to center and adjacent edge of
+ * plane polygon
+ */
+ double theta;
+
+ /* additional variables from snyder */
+ double q, Rprime, H, Ag, Azprime, Az, dprime, f, rho,
+ x, y;
+
+ /* variables used to store intermediate results */
+ double cot_theta, tan_g, az_offset;
+
+ /* how many multiples of 60 degrees we adjust the azimuth */
+ int Az_adjust_multiples;
+
+ struct snyder_constants c;
+
+ /*
+ * TODO by locality of reference, start by trying the same triangle
+ * as last time
+ */
+
+ /* TODO put these constants in as radians to begin with */
+ c = constants[SNYDER_POLY_ICOSAHEDRON];
+ theta = c.theta * DEG2RAD;
+ g = c.g * DEG2RAD;
+ G = c.G * DEG2RAD;
+
+ for (i = 1; i <= 20; i++) {
+ double z;
+ struct isea_geo center;
+
+ center = icostriangles[i];
+
+ /* step 1 */
+ z = acos(sin(center.lat) * sin(ll->lat)
+ + cos(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon));
+ /* not on this triangle */
+ if (z > g + 0.000005) { /* TODO DBL_EPSILON */
+ continue;
+ }
+
+ Az = sph_azimuth(center.lon, center.lat, ll->lon, ll->lat);
+
+ /* step 2 */
+
+ /* This calculates "some" vertex coordinate */
+ az_offset = az_adjustment(i);
+
+ Az -= az_offset;
+
+ /* TODO I don't know why we do this. It's not in snyder */
+ /* maybe because we should have picked a better vertex */
+ if (Az < 0.0) {
+ Az += 2.0 * M_PI;
+ }
+ /*
+ * adjust Az for the point to fall within the range of 0 to
+ * 2(90 - theta) or 60 degrees for the hexagon, by
+ * and therefore 120 degrees for the triangle
+ * of the icosahedron
+ * subtracting or adding multiples of 60 degrees to Az and
+ * recording the amount of adjustment
+ */
+
+ Az_adjust_multiples = 0;
+ while (Az < 0.0) {
+ Az += DEG120;
+ Az_adjust_multiples--;
+ }
+ while (Az > DEG120 + DBL_EPSILON) {
+ Az -= DEG120;
+ Az_adjust_multiples++;
+ }
+
+ /* step 3 */
+ cot_theta = 1.0 / tan(theta);
+ tan_g = tan(g); /* TODO this is a constant */
+
+ /* Calculate q from eq 9. */
+ /* TODO cot_theta is cot(30) */
+ q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta);
+
+ /* not in this triangle */
+ if (z > q + 0.000005) {
+ continue;
+ }
+ /* step 4 */
+
+ /* Apply equations 5-8 and 10-12 in order */
+
+ /* eq 5 */
+ /* Rprime = 0.9449322893 * R; */
+ /* R' in the paper is for the truncated */
+ Rprime = 0.91038328153090290025;
+
+ /* eq 6 */
+ H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G));
+
+ /* eq 7 */
+ /* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */
+ Ag = Az + G + H - DEG180;
+
+ /* eq 8 */
+ Azprime = atan2(2.0 * Ag, Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta);
+
+ /* eq 10 */
+ /* cot(theta) = 1.73205080756887729355 */
+ dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta);
+
+ /* eq 11 */
+ f = dprime / (2.0 * Rprime * sin(q / 2.0));
+
+ /* eq 12 */
+ rho = 2.0 * Rprime * f * sin(z / 2.0);
+
+ /*
+ * add back the same 60 degree multiple adjustment from step
+ * 2 to Azprime
+ */
+
+ Azprime += DEG120 * Az_adjust_multiples;
+
+ /* calculate rectangular coordinates */
+
+ x = rho * sin(Azprime);
+ y = rho * cos(Azprime);
+
+ /*
+ * TODO
+ * translate coordinates to the origin for the particular
+ * hexagon on the flattened polyhedral map plot
+ */
+
+ out->x = x;
+ out->y = y;
+
+ return i;
+ }
+
+ /*
+ * should be impossible, this implies that the coordinate is not on
+ * any triangle
+ */
+
+ fprintf(stderr, "impossible transform: %f %f is not on any triangle\n",
+ ll->lon * RAD2DEG, ll->lat * RAD2DEG);
+
+ exit(EXIT_FAILURE);
+
+ /* not reached */
+ return 0; /* supresses a warning */
+}
+
+/*
+ * return the new coordinates of any point in orginal coordinate system.
+ * Define a point (newNPold) in orginal coordinate system as the North Pole in
+ * new coordinate system, and the great circle connect the original and new
+ * North Pole as the lon0 longitude in new coordinate system, given any point
+ * in orginal coordinate system, this function return the new coordinates.
+ */
+
+#define PRECISION 0.0000000000005
+
+/* formula from Snyder, Map Projections: A working manual, p31 */
+/*
+ * old north pole at np in new coordinates
+ * could be simplified a bit with fewer intermediates
+ *
+ * TODO take a result pointer
+ */
+ISEA_STATIC
+struct isea_geo
+snyder_ctran(struct isea_geo * np, struct isea_geo * pt)
+{
+ struct isea_geo npt;
+ double alpha, phi, lambda, lambda0, beta, lambdap, phip;
+ double sin_phip;
+ double lp_b; /* lambda prime minus beta */
+ double cos_p, sin_a;
+
+ phi = pt->lat;
+ lambda = pt->lon;
+ alpha = np->lat;
+ beta = np->lon;
+ lambda0 = beta;
+
+ cos_p = cos(phi);
+ sin_a = sin(alpha);
+
+ /* mpawm 5-7 */
+ sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0);
+
+ /* mpawm 5-8b */
+
+ /* use the two argument form so we end up in the right quadrant */
+ lp_b = atan2(cos_p * sin(lambda - lambda0),
+ (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));
+
+ lambdap = lp_b + beta;
+
+ /* normalize longitude */
+ /* TODO can we just do a modulus ? */
+ lambdap = fmod(lambdap, 2 * M_PI);
+ while (lambdap > M_PI)
+ lambdap -= 2 * M_PI;
+ while (lambdap < -M_PI)
+ lambdap += 2 * M_PI;
+
+ phip = asin(sin_phip);
+
+ npt.lat = phip;
+ npt.lon = lambdap;
+
+ return npt;
+}
+
+ISEA_STATIC
+struct isea_geo
+isea_ctran(struct isea_geo * np, struct isea_geo * pt, double lon0)
+{
+ struct isea_geo npt;
+
+ np->lon += M_PI;
+ npt = snyder_ctran(np, pt);
+ np->lon -= M_PI;
+
+ npt.lon -= (M_PI - lon0 + np->lon);
+
+ /*
+ * snyder is down tri 3, isea is along side of tri1 from vertex 0 to
+ * vertex 1 these are 180 degrees apart
+ */
+ npt.lon += M_PI;
+ /* normalize longitude */
+ npt.lon = fmod(npt.lon, 2 * M_PI);
+ while (npt.lon > M_PI)
+ npt.lon -= 2 * M_PI;
+ while (npt.lon < -M_PI)
+ npt.lon += 2 * M_PI;
+
+ return npt;
+}
+
+/* in radians */
+#define ISEA_STD_LAT 1.01722196792335072101
+#define ISEA_STD_LON .19634954084936207740
+
+/* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */
+
+ISEA_STATIC
+int
+isea_grid_init(struct isea_dgg * g)
+{
+ if (!g)
+ return 0;
+
+ g->polyhedron = 20;
+ g->o_lat = ISEA_STD_LAT;
+ g->o_lon = ISEA_STD_LON;
+ g->o_az = 0.0;
+ g->aperture = 4;
+ g->resolution = 6;
+ g->radius = 1.0;
+ g->topology = 6;
+
+ return 1;
+}
+
+ISEA_STATIC
+int
+isea_orient_isea(struct isea_dgg * g)
+{
+ if (!g)
+ return 0;
+ g->o_lat = ISEA_STD_LAT;
+ g->o_lon = ISEA_STD_LON;
+ g->o_az = 0.0;
+ return 1;
+}
+
+ISEA_STATIC
+int
+isea_orient_pole(struct isea_dgg * g)
+{
+ if (!g)
+ return 0;
+ g->o_lat = M_PI / 2.0;
+ g->o_lon = 0.0;
+ g->o_az = 0;
+ return 1;
+}
+
+ISEA_STATIC
+int
+isea_transform(struct isea_dgg * g, struct isea_geo * in,
+ struct isea_pt * out)
+{
+ struct isea_geo i, pole;
+ int tri;
+
+ pole.lat = g->o_lat;
+ pole.lon = g->o_lon;
+
+ i = isea_ctran(&pole, in, g->o_az);
+
+ tri = isea_snyder_forward(&i, out);
+ out->x *= g->radius;
+ out->y *= g->radius;
+ g->triangle = tri;
+
+ return tri;
+}
+
+#define DOWNTRI(tri) (((tri - 1) / 5) % 2 == 1)
+
+ISEA_STATIC
+void
+isea_rotate(struct isea_pt * pt, double degrees)
+{
+ double rad;
+
+ double x, y;
+
+ rad = -degrees * M_PI / 180.0;
+ while (rad >= 2.0 * M_PI) rad -= 2.0 * M_PI;
+ while (rad <= -2.0 * M_PI) rad += 2.0 * M_PI;
+
+ x = pt->x * cos(rad) + pt->y * sin(rad);
+ y = -pt->x * sin(rad) + pt->y * cos(rad);
+
+ pt->x = x;
+ pt->y = y;
+}
+
+ISEA_STATIC
+int isea_tri_plane(int tri, struct isea_pt *pt, double radius) {
+ struct isea_pt tc; /* center of triangle */
+
+ if (DOWNTRI(tri)) {
+ isea_rotate(pt, 180.0);
+ }
+ tc = isea_triangle_xy(tri);
+ tc.x *= radius;
+ tc.y *= radius;
+ pt->x += tc.x;
+ pt->y += tc.y;
+
+ return tri;
+}
+
+/* convert projected triangle coords to quad xy coords, return quad number */
+ISEA_STATIC
+int
+isea_ptdd(int tri, struct isea_pt *pt) {
+ int downtri, quad;
+
+ downtri = (((tri - 1) / 5) % 2 == 1);
+ quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
+
+ isea_rotate(pt, downtri ? 240.0 : 60.0);
+ if (downtri) {
+ pt->x += 0.5;
+ /* pt->y += cos(30.0 * M_PI / 180.0); */
+ pt->y += .86602540378443864672;
+ }
+ return quad;
+}
+
+ISEA_STATIC
+int
+isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di)
+{
+ struct isea_pt v;
+ double hexwidth;
+ double sidelength; /* in hexes */
+ int d, i;
+ int maxcoord;
+ struct hex h;
+
+ /* This is the number of hexes from apex to base of a triangle */
+ sidelength = (pow(2.0, g->resolution) + 1.0) / 2.0;
+
+ /* apex to base is cos(30deg) */
+ hexwidth = cos(M_PI / 6.0) / sidelength;
+
+ /* TODO I think sidelength is always x.5, so
+ * (int)sidelength * 2 + 1 might be just as good
+ */
+ maxcoord = (int) (sidelength * 2.0 + 0.5);
+
+ v = *pt;
+ hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);
+ h.iso = 0;
+ hex_iso(&h);
+
+ d = h.x - h.z;
+ i = h.x + h.y + h.y;
+
+ /*
+ * you want to test for max coords for the next quad in the same
+ * "row" first to get the case where both are max
+ */
+ if (quad <= 5) {
+ if (d == 0 && i == maxcoord) {
+ /* north pole */
+ quad = 0;
+ d = 0;
+ i = 0;
+ } else if (i == maxcoord) {
+ /* upper right in next quad */
+ quad += 1;
+ if (quad == 6)
+ quad = 1;
+ i = maxcoord - d;
+ d = 0;
+ } else if (d == maxcoord) {
+ /* lower right in quad to lower right */
+ quad += 5;
+ d = 0;
+ }
+ } else if (quad >= 6) {
+ if (i == 0 && d == maxcoord) {
+ /* south pole */
+ quad = 11;
+ d = 0;
+ i = 0;
+ } else if (d == maxcoord) {
+ /* lower right in next quad */
+ quad += 1;
+ if (quad == 11)
+ quad = 6;
+ d = maxcoord - i;
+ i = 0;
+ } else if (i == maxcoord) {
+ /* upper right in quad to upper right */
+ quad = (quad - 4) % 5;
+ i = 0;
+ }
+ }
+
+ di->x = d;
+ di->y = i;
+
+ g->quad = quad;
+ return quad;
+}
+
+ISEA_STATIC
+int
+isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di) {
+ struct isea_pt v;
+ double hexwidth;
+ int sidelength; /* in hexes */
+ struct hex h;
+
+ if (g->aperture == 3 && g->resolution % 2 != 0) {
+ return isea_dddi_ap3odd(g, quad, pt, di);
+ }
+ /* todo might want to do this as an iterated loop */
+ if (g->aperture >0) {
+ sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
+ } else {
+ sidelength = g->resolution;
+ }
+
+ hexwidth = 1.0 / sidelength;
+
+ v = *pt;
+ isea_rotate(&v, -30.0);
+ hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);
+ h.iso = 0;
+ hex_iso(&h);
+
+ /* we may actually be on another quad */
+ if (quad <= 5) {
+ if (h.x == 0 && h.z == -sidelength) {
+ /* north pole */
+ quad = 0;
+ h.z = 0;
+ h.y = 0;
+ h.x = 0;
+ } else if (h.z == -sidelength) {
+ quad = quad + 1;
+ if (quad == 6)
+ quad = 1;
+ h.y = sidelength - h.x;
+ h.z = h.x - sidelength;
+ h.x = 0;
+ } else if (h.x == sidelength) {
+ quad += 5;
+ h.y = -h.z;
+ h.x = 0;
+ }
+ } else if (quad >= 6) {
+ if (h.z == 0 && h.x == sidelength) {
+ /* south pole */
+ quad = 11;
+ h.x = 0;
+ h.y = 0;
+ h.z = 0;
+ } else if (h.x == sidelength) {
+ quad = quad + 1;
+ if (quad == 11)
+ quad = 6;
+ h.x = h.y + sidelength;
+ h.y = 0;
+ h.z = -h.x;
+ } else if (h.y == -sidelength) {
+ quad -= 4;
+ h.y = 0;
+ h.z = -h.x;
+ }
+ }
+ di->x = h.x;
+ di->y = -h.z;
+
+ g->quad = quad;
+ return quad;
+}
+
+ISEA_STATIC
+int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,
+ struct isea_pt *di) {
+ struct isea_pt v;
+ int quad;
+
+ v = *pt;
+ quad = isea_ptdd(tri, &v);
+ quad = isea_dddi(g, quad, &v, di);
+ return quad;
+}
+
+/* q2di to seqnum */
+ISEA_STATIC
+int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
+ int sidelength;
+ int sn, height;
+ int hexes;
+
+ if (quad == 0) {
+ g->serial = 1;
+ return g->serial;
+ }
+ /* hexes in a quad */
+ hexes = (int) (pow(g->aperture, g->resolution) + 0.5);
+ if (quad == 11) {
+ g->serial = 1 + 10 * hexes + 1;
+ return g->serial;
+ }
+ if (g->aperture == 3 && g->resolution % 2 == 1) {
+ height = (int) (pow(g->aperture, (g->resolution - 1) / 2.0));
+ sn = ((int) di->x) * height;
+ sn += ((int) di->y) / height;
+ sn += (quad - 1) * hexes;
+ sn += 2;
+ } else {
+ sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
+ sn = (quad - 1) * hexes + sidelength * di->x + di->y + 2;
+ }
+
+ g->serial = sn;
+ return sn;
+}
+
+/* TODO just encode the quad in the d or i coordinate
+ * quad is 0-11, which can be four bits.
+ * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf
+ */
+/* convert a q2di to global hex coord */
+ISEA_STATIC
+int isea_hex(struct isea_dgg *g, int tri,
+ struct isea_pt *pt, struct isea_pt *hex) {
+ struct isea_pt v;
+ int sidelength;
+ int d, i, x, y, quad;
+
+ quad = isea_ptdi(g, tri, pt, &v);
+
+ hex->x = ((int)v.x << 4) + quad;
+ hex->y = v.y;
+
+ return 1;
+
+ d = v.x;
+ i = v.y;
+
+ /* Aperture 3 odd resolutions */
+ if (g->aperture == 3 && g->resolution % 2 != 0) {
+ int offset = (int)(pow(3.0, g->resolution - 1) + 0.5);
+
+ d += offset * ((g->quad-1) % 5);
+ i += offset * ((g->quad-1) % 5);
+
+ if (quad == 0) {
+ d = 0;
+ i = offset;
+ } else if (quad == 11) {
+ d = 2 * offset;
+ i = 0;
+ } else if (quad > 5) {
+ d += offset;
+ }
+
+ x = (2*d - i) /3;
+ y = (2*i - d) /3;
+
+ hex->x = x + offset / 3;
+ hex->y = y + 2 * offset / 3;
+ return 1;
+ }
+
+ /* aperture 3 even resolutions and aperture 4 */
+ sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
+ if (g->quad == 0) {
+ hex->x = 0;
+ hex->y = sidelength;
+ } else if (g->quad == 11) {
+ hex->x = sidelength * 2;
+ hex->y = 0;
+ } else {
+ hex->x = d + sidelength * ((g->quad-1) % 5);
+ if (g->quad > 5) hex->x += sidelength;
+ hex->y = i + sidelength * ((g->quad-1) % 5);
+ }
+
+ return 1;
+}
+
+ISEA_STATIC
+struct isea_pt
+isea_forward(struct isea_dgg *g, struct isea_geo *in)
+{
+ int tri;
+ struct isea_pt out, coord;
+
+ tri = isea_transform(g, in, &out);
+
+ if (g->output == ISEA_PLANE) {
+ isea_tri_plane(tri, &out, g->radius);
+ return out;
+ }
+
+ /* convert to isea standard triangle size */
+ out.x = out.x / g->radius * ISEA_SCALE;
+ out.y = out.y / g->radius * ISEA_SCALE;
+ out.x += 0.5;
+ out.y += 2.0 * .14433756729740644112;
+
+ switch (g->output) {
+ case ISEA_PROJTRI:
+ /* nothing to do, already in projected triangle */
+ break;
+ case ISEA_VERTEX2DD:
+ g->quad = isea_ptdd(tri, &out);
+ break;
+ case ISEA_Q2DD:
+ /* Same as above, we just don't print as much */
+ g->quad = isea_ptdd(tri, &out);
+ break;
+ case ISEA_Q2DI:
+ g->quad = isea_ptdi(g, tri, &out, &coord);
+ return coord;
+ break;
+ case ISEA_SEQNUM:
+ isea_ptdi(g, tri, &out, &coord);
+ /* disn will set g->serial */
+ isea_disn(g, g->quad, &coord);
+ return coord;
+ break;
+ case ISEA_HEX:
+ isea_hex(g, tri, &out, &coord);
+ return coord;
+ break;
+ }
+
+ return out;
+}
+/*
+ * Proj 4 integration code follows
+ */
+
+#define PJ_LIB__
+#include <projects.h>
+
+PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph";
+
+struct pj_opaque {
+ struct isea_dgg dgg;
+};
+
+
+static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
+ XY xy = {0.0,0.0};
+ struct pj_opaque *Q = P->opaque;
+ struct isea_pt out;
+ struct isea_geo in;
+
+ in.lon = lp.lam;
+ in.lat = lp.phi;
+
+ out = isea_forward(&Q->dgg, &in);
+
+ xy.x = out.x;
+ xy.y = out.y;
+
+ return xy;
+}
+
+
+static void *freeup_new (PJ *P) { /* Destructor */
+ if (0==P)
+ return 0;
+ if (0==P->opaque)
+ return pj_dealloc (P);
+
+ pj_dealloc (P->opaque);
+ return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+ freeup_new (P);
+ return;
+}
+
+
+PJ *PROJECTION(isea) {
+ char *opt;
+ struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+ if (0==Q)
+ return freeup_new (P);
+ P->opaque = Q;
+
+
+ P->fwd = s_forward;
+ isea_grid_init(&Q->dgg);
+
+ Q->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(&Q->dgg);
+ } else if (!strcmp(opt, "pole")) {
+ isea_orient_pole(&Q->dgg);
+ } else {
+ E_ERROR(-34);
+ }
+ }
+
+ if (pj_param(P->ctx,P->params, "tazi").i) {
+ Q->dgg.o_az = pj_param(P->ctx,P->params, "razi").f;
+ }
+
+ if (pj_param(P->ctx,P->params, "tlon_0").i) {
+ Q->dgg.o_lon = pj_param(P->ctx,P->params, "rlon_0").f;
+ }
+
+ if (pj_param(P->ctx,P->params, "tlat_0").i) {
+ Q->dgg.o_lat = pj_param(P->ctx,P->params, "rlat_0").f;
+ }
+
+ if (pj_param(P->ctx,P->params, "taperture").i) {
+ Q->dgg.aperture = pj_param(P->ctx,P->params, "iaperture").i;
+ }
+
+ if (pj_param(P->ctx,P->params, "tresolution").i) {
+ Q->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")) {
+ Q->dgg.output = ISEA_PLANE;
+ } else if (!strcmp(opt, "di")) {
+ Q->dgg.output = ISEA_Q2DI;
+ }
+ else if (!strcmp(opt, "dd")) {
+ Q->dgg.output = ISEA_Q2DD;
+ }
+ else if (!strcmp(opt, "hex")) {
+ Q->dgg.output = ISEA_HEX;
+ }
+ else {
+ /* TODO verify error code. Possibly eliminate magic */
+ E_ERROR(-34);
+ }
+ }
+
+ if (pj_param(P->ctx,P->params, "trescale").i) {
+ Q->dgg.radius = ISEA_SCALE;
+ }
+
+ if (pj_param(P->ctx,P->params, "tresolution").i) {
+ Q->dgg.resolution = pj_param(P->ctx,P->params, "iresolution").i;
+ } else {
+ Q->dgg.resolution = 4;
+ }
+
+ if (pj_param(P->ctx,P->params, "taperture").i) {
+ Q->dgg.aperture = pj_param(P->ctx,P->params, "iaperture").i;
+ } else {
+ Q->dgg.aperture = 3;
+ }
+
+ return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_isea_selftest (void) {return 0;}
+#else
+
+int pj_isea_selftest (void) {
+ double tolerance_lp = 1e-10;
+ double tolerance_xy = 1e-7;
+
+ char s_args[] = {"+proj=isea +a=6400000 +lat_1=0.5 +lat_2=2"};
+
+ LP fwd_in[] = {
+ { 2, 1},
+ { 2,-1},
+ {-2, 1},
+ {-2,-1}
+ };
+
+ XY s_fwd_expect[] = {
+ {-1097074.9480224741, 3442909.3090371834},
+ {-1097074.9482647954, 3233611.7285857084},
+ {-1575486.3536415542, 3442168.3420281881},
+ {-1575486.353880283, 3234352.6955947056},
+ };
+
+ return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, 0, 0, 0);
+}
+
+
+#endif
diff --git a/src/PJ_vandg4.c b/src/PJ_vandg4.c
index 4d4032b3..ca079000 100644
--- a/src/PJ_vandg4.c
+++ b/src/PJ_vandg4.c
@@ -1,97 +1,97 @@
-#define PJ_LIB__
-#include <projects.h>
-
-PROJ_HEAD(vandg4, "van der Grinten IV") "\n\tMisc Sph, no inv.";
-
-#define TOL 1e-10
-
-
-static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
- XY xy = {0.0,0.0};
- double x1, t, bt, ct, ft, bt2, ct2, dt, dt2;
- (void) P;
-
- if (fabs(lp.phi) < TOL) {
- xy.x = lp.lam;
- xy.y = 0.;
- } else if (fabs(lp.lam) < TOL || fabs(fabs(lp.phi) - M_HALFPI) < TOL) {
- xy.x = 0.;
- xy.y = lp.phi;
- } else {
- bt = fabs(M_TWO_D_PI * lp.phi);
- bt2 = bt * bt;
- ct = 0.5 * (bt * (8. - bt * (2. + bt2)) - 5.)
- / (bt2 * (bt - 1.));
- ct2 = ct * ct;
- dt = M_TWO_D_PI * lp.lam;
- dt = dt + 1. / dt;
- dt = sqrt(dt * dt - 4.);
- if ((fabs(lp.lam) - M_HALFPI) < 0.) dt = -dt;
- dt2 = dt * dt;
- x1 = bt + ct; x1 *= x1;
- t = bt + 3.*ct;
- ft = x1 * (bt2 + ct2 * dt2 - 1.) + (1.-bt2) * (
- bt2 * (t * t + 4. * ct2) +
- ct2 * (12. * bt * ct + 4. * ct2) );
- x1 = (dt*(x1 + ct2 - 1.) + 2.*sqrt(ft)) /
- (4.* x1 + dt2);
- xy.x = M_HALFPI * x1;
- xy.y = M_HALFPI * sqrt(1. + dt * fabs(x1) - x1 * x1);
- if (lp.lam < 0.) xy.x = -xy.x;
- if (lp.phi < 0.) xy.y = -xy.y;
- }
- return xy;
-}
-
-
-static void *freeup_new (PJ *P) { /* Destructor */
- if (0==P)
- return 0;
-
- return pj_dealloc(P);
-}
-
-
-static void freeup (PJ *P) {
- freeup_new (P);
- return;
-}
-
-
-PJ *PROJECTION(vandg4) {
- P->es = 0.;
- P->fwd = s_forward;
-
- return P;
-}
-
-
-#ifdef PJ_OMIT_SELFTEST
-int pj_vandg4_selftest (void) {return 0;}
-#else
-
-int pj_vandg4_selftest (void) {
- double tolerance_lp = 1e-10;
- double tolerance_xy = 1e-7;
-
- char s_args[] = {"+proj=vandg4 +a=6400000 +lat_1=0.5 +lat_2=2"};
-
- LP fwd_in[] = {
- { 2, 1},
- { 2,-1},
- {-2, 1},
- {-2,-1}
- };
-
- XY s_fwd_expect[] = {
- { 223374.57729435508, 111701.19548415358 },
- { 223374.57729435508, -111701.19548415358 },
- {-223374.57729435508, 111701.19548415358 },
- {-223374.57729435508, -111701.19548415358 },
- };
-
- return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, 0, 0, 0);
-}
-
-
-#endif
+#define PJ_LIB__
+#include <projects.h>
+
+PROJ_HEAD(vandg4, "van der Grinten IV") "\n\tMisc Sph, no inv.";
+
+#define TOL 1e-10
+
+
+static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
+ XY xy = {0.0,0.0};
+ double x1, t, bt, ct, ft, bt2, ct2, dt, dt2;
+ (void) P;
+
+ if (fabs(lp.phi) < TOL) {
+ xy.x = lp.lam;
+ xy.y = 0.;
+ } else if (fabs(lp.lam) < TOL || fabs(fabs(lp.phi) - M_HALFPI) < TOL) {
+ xy.x = 0.;
+ xy.y = lp.phi;
+ } else {
+ bt = fabs(M_TWO_D_PI * lp.phi);
+ bt2 = bt * bt;
+ ct = 0.5 * (bt * (8. - bt * (2. + bt2)) - 5.)
+ / (bt2 * (bt - 1.));
+ ct2 = ct * ct;
+ dt = M_TWO_D_PI * lp.lam;
+ dt = dt + 1. / dt;
+ dt = sqrt(dt * dt - 4.);
+ if ((fabs(lp.lam) - M_HALFPI) < 0.) dt = -dt;
+ dt2 = dt * dt;
+ x1 = bt + ct; x1 *= x1;
+ t = bt + 3.*ct;
+ ft = x1 * (bt2 + ct2 * dt2 - 1.) + (1.-bt2) * (
+ bt2 * (t * t + 4. * ct2) +
+ ct2 * (12. * bt * ct + 4. * ct2) );
+ x1 = (dt*(x1 + ct2 - 1.) + 2.*sqrt(ft)) /
+ (4.* x1 + dt2);
+ xy.x = M_HALFPI * x1;
+ xy.y = M_HALFPI * sqrt(1. + dt * fabs(x1) - x1 * x1);
+ if (lp.lam < 0.) xy.x = -xy.x;
+ if (lp.phi < 0.) xy.y = -xy.y;
+ }
+ return xy;
+}
+
+
+static void *freeup_new (PJ *P) { /* Destructor */
+ if (0==P)
+ return 0;
+
+ return pj_dealloc(P);
+}
+
+
+static void freeup (PJ *P) {
+ freeup_new (P);
+ return;
+}
+
+
+PJ *PROJECTION(vandg4) {
+ P->es = 0.;
+ P->fwd = s_forward;
+
+ return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_vandg4_selftest (void) {return 0;}
+#else
+
+int pj_vandg4_selftest (void) {
+ double tolerance_lp = 1e-10;
+ double tolerance_xy = 1e-7;
+
+ char s_args[] = {"+proj=vandg4 +a=6400000 +lat_1=0.5 +lat_2=2"};
+
+ LP fwd_in[] = {
+ { 2, 1},
+ { 2,-1},
+ {-2, 1},
+ {-2,-1}
+ };
+
+ XY s_fwd_expect[] = {
+ { 223374.57729435508, 111701.19548415358 },
+ { 223374.57729435508, -111701.19548415358 },
+ {-223374.57729435508, 111701.19548415358 },
+ {-223374.57729435508, -111701.19548415358 },
+ };
+
+ return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, 0, 0, 0);
+}
+
+
+#endif
diff --git a/src/pj_factors.c b/src/pj_factors.c
index f7cb75e6..ac23ed71 100644
--- a/src/pj_factors.c
+++ b/src/pj_factors.c
@@ -21,7 +21,7 @@ pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
if (h < EPS)
h = DEFAULT_H;
- if (fabs(lp.phi) > (M_HALFPI - h))
+ if (fabs(lp.phi) > (M_HALFPI - h))
/* adjust to value around pi/2 where derived still exists*/
lp.phi = lp.phi < 0. ? (-M_HALFPI+h) : (M_HALFPI-h);
else if (P->geoc)
diff --git a/src/pj_gridinfo.c b/src/pj_gridinfo.c
index 48ef7b87..5f528de1 100644
--- a/src/pj_gridinfo.c
+++ b/src/pj_gridinfo.c
@@ -443,7 +443,7 @@ static int pj_gridinfo_init_ntv2( projCtx ctx, PAFile fid, PJ_GRIDINFO *gilist )
pj_ctx_set_errno( ctx, -38 );
return 0;
}
-
+
if( header[8] == 11 )
must_swap = !IS_LSB;
else
diff --git a/src/pj_list.c b/src/pj_list.c
index a8f171ab..9fde000f 100644
--- a/src/pj_list.c
+++ b/src/pj_list.c
@@ -53,4 +53,3 @@ struct PJ_LIST *pj_get_list_ref (void) {
struct PJ_SELFTEST_LIST *pj_get_selftest_list_ref (void) {
return pj_selftest_list;
}
-
diff --git a/src/pj_list.h b/src/pj_list.h
index baacacc3..bc4c448a 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -147,4 +147,3 @@ PROJ_HEAD(weren, "Werenskiold I")
PROJ_HEAD(wink1, "Winkel I")
PROJ_HEAD(wink2, "Winkel II")
PROJ_HEAD(wintri, "Winkel Tripel")
-
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 11:16:05 +0000