aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorHoward Butler <howard@hobu.co>2016-08-18 09:45:08 -0500
committerHoward Butler <howard@hobu.co>2016-08-18 09:45:08 -0500
commit77104aae973db34053e20163b73b511129652f73 (patch)
tree0b0a9d45a27e379dcca181556feea97a260572c3 /src
parent91e7ff28812d7c06e642a3e5e247610eb6e1988e (diff)
downloadPROJ-77104aae973db34053e20163b73b511129652f73.tar.gz
PROJ-77104aae973db34053e20163b73b511129652f73.zip
whitespace normalization for proj 4.9.3RC3
Diffstat (limited to 'src')
-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
12 files changed, 1927 insertions, 1929 deletions
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 20:31:10 +0000