From 77104aae973db34053e20163b73b511129652f73 Mon Sep 17 00:00:00 2001 From: Howard Butler Date: Thu, 18 Aug 2016 09:45:08 -0500 Subject: whitespace normalization for proj 4.9.3RC3 --- src/PJ_eck1.c | 184 ++--- src/PJ_eck4.c | 234 +++--- src/PJ_eck5.c | 186 ++--- src/PJ_fahey.c | 190 ++--- src/PJ_goode.c | 244 +++--- src/PJ_hatano.c | 268 +++--- src/PJ_isea.c | 2350 ++++++++++++++++++++++++++--------------------------- src/PJ_vandg4.c | 194 ++--- src/pj_factors.c | 2 +- src/pj_gridinfo.c | 2 +- src/pj_list.c | 1 - src/pj_list.h | 1 - 12 files changed, 1927 insertions(+), 1929 deletions(-) (limited to 'src') 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 - -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 + +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 - -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 + +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 - -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 + +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 - -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 + +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 - -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 + +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 - -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 + +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 -#include -#include -#include - -#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 +#include +#include +#include + +#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 - -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 + +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 - -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 + +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") - -- cgit v1.2.3