diff options
Diffstat (limited to 'src/PJ_igh.cpp')
| -rw-r--r-- | src/PJ_igh.cpp | 225 |
1 files changed, 225 insertions, 0 deletions
diff --git a/src/PJ_igh.cpp b/src/PJ_igh.cpp new file mode 100644 index 00000000..476d1c6b --- /dev/null +++ b/src/PJ_igh.cpp @@ -0,0 +1,225 @@ +#define PJ_LIB__ + +#include <errno.h> +#include <math.h> + +#include "projects.h" + +PROJ_HEAD(igh, "Interrupted Goode Homolosine") "\n\tPCyl, Sph"; + +C_NAMESPACE PJ *pj_sinu(PJ *), *pj_moll(PJ *); + +/* 40d 44' 11.8" [degrees] */ +/* +static const double d4044118 = (40 + 44/60. + 11.8/3600.) * DEG_TO_RAD; +has been replaced by this define, to eliminate portability issue: +Initializer element not computable at load time +*/ +#define d4044118 ((40 + 44/60. + 11.8/3600.) * DEG_TO_RAD) + +static const double d10 = 10 * DEG_TO_RAD; +static const double d20 = 20 * DEG_TO_RAD; +static const double d30 = 30 * DEG_TO_RAD; +static const double d40 = 40 * DEG_TO_RAD; +static const double d50 = 50 * DEG_TO_RAD; +static const double d60 = 60 * DEG_TO_RAD; +static const double d80 = 80 * DEG_TO_RAD; +static const double d90 = 90 * DEG_TO_RAD; +static const double d100 = 100 * DEG_TO_RAD; +static const double d140 = 140 * DEG_TO_RAD; +static const double d160 = 160 * DEG_TO_RAD; +static const double d180 = 180 * DEG_TO_RAD; + +static const double EPSLN = 1.e-10; /* allow a little 'slack' on zone edge positions */ + +struct pj_opaque { + struct PJconsts* pj[12]; \ + double dy0; +}; + + +static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ + XY xy; + struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); + int z; + + if (lp.phi >= d4044118) { /* 1|2 */ + z = (lp.lam <= -d40 ? 1: 2); + } + else if (lp.phi >= 0) { /* 3|4 */ + z = (lp.lam <= -d40 ? 3: 4); + } + else if (lp.phi >= -d4044118) { /* 5|6|7|8 */ + if (lp.lam <= -d100) z = 5; /* 5 */ + else if (lp.lam <= -d20) z = 6; /* 6 */ + else if (lp.lam <= d80) z = 7; /* 7 */ + else z = 8; /* 8 */ + } + else { /* 9|10|11|12 */ + if (lp.lam <= -d100) z = 9; /* 9 */ + else if (lp.lam <= -d20) z = 10; /* 10 */ + else if (lp.lam <= d80) z = 11; /* 11 */ + else z = 12; /* 12 */ + } + + lp.lam -= Q->pj[z-1]->lam0; + xy = Q->pj[z-1]->fwd(lp, Q->pj[z-1]); + xy.x += Q->pj[z-1]->x0; + xy.y += Q->pj[z-1]->y0; + + return xy; +} + + +static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ + LP lp = {0.0,0.0}; + struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); + const double y90 = Q->dy0 + sqrt(2); /* lt=90 corresponds to y=y0+sqrt(2) */ + + int z = 0; + if (xy.y > y90+EPSLN || xy.y < -y90+EPSLN) /* 0 */ + z = 0; + else if (xy.y >= d4044118) /* 1|2 */ + z = (xy.x <= -d40? 1: 2); + else if (xy.y >= 0) /* 3|4 */ + z = (xy.x <= -d40? 3: 4); + else if (xy.y >= -d4044118) { /* 5|6|7|8 */ + if (xy.x <= -d100) z = 5; /* 5 */ + else if (xy.x <= -d20) z = 6; /* 6 */ + else if (xy.x <= d80) z = 7; /* 7 */ + else z = 8; /* 8 */ + } + else { /* 9|10|11|12 */ + if (xy.x <= -d100) z = 9; /* 9 */ + else if (xy.x <= -d20) z = 10; /* 10 */ + else if (xy.x <= d80) z = 11; /* 11 */ + else z = 12; /* 12 */ + } + + if (z) { + int ok = 0; + + xy.x -= Q->pj[z-1]->x0; + xy.y -= Q->pj[z-1]->y0; + lp = Q->pj[z-1]->inv(xy, Q->pj[z-1]); + lp.lam += Q->pj[z-1]->lam0; + + switch (z) { + case 1: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d40+EPSLN) || + ((lp.lam >= -d40-EPSLN && lp.lam <= -d10+EPSLN) && + (lp.phi >= d60-EPSLN && lp.phi <= d90+EPSLN)); break; + case 2: ok = (lp.lam >= -d40-EPSLN && lp.lam <= d180+EPSLN) || + ((lp.lam >= -d180-EPSLN && lp.lam <= -d160+EPSLN) && + (lp.phi >= d50-EPSLN && lp.phi <= d90+EPSLN)) || + ((lp.lam >= -d50-EPSLN && lp.lam <= -d40+EPSLN) && + (lp.phi >= d60-EPSLN && lp.phi <= d90+EPSLN)); break; + case 3: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d40+EPSLN); break; + case 4: ok = (lp.lam >= -d40-EPSLN && lp.lam <= d180+EPSLN); break; + case 5: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d100+EPSLN); break; + case 6: ok = (lp.lam >= -d100-EPSLN && lp.lam <= -d20+EPSLN); break; + case 7: ok = (lp.lam >= -d20-EPSLN && lp.lam <= d80+EPSLN); break; + case 8: ok = (lp.lam >= d80-EPSLN && lp.lam <= d180+EPSLN); break; + case 9: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d100+EPSLN); break; + case 10: ok = (lp.lam >= -d100-EPSLN && lp.lam <= -d20+EPSLN); break; + case 11: ok = (lp.lam >= -d20-EPSLN && lp.lam <= d80+EPSLN); break; + case 12: ok = (lp.lam >= d80-EPSLN && lp.lam <= d180+EPSLN); break; + } + z = (!ok? 0: z); /* projectable? */ + } + + if (!z) lp.lam = HUGE_VAL; + if (!z) lp.phi = HUGE_VAL; + + return lp; +} + + +static PJ *destructor (PJ *P, int errlev) { + int i; + if (0==P) + return 0; + + if (0==P->opaque) + return pj_default_destructor (P, errlev); + + for (i = 0; i < 12; ++i) { + if (static_cast<struct pj_opaque*>(P->opaque)->pj[i]) + static_cast<struct pj_opaque*>(P->opaque)->pj[i]->destructor(static_cast<struct pj_opaque*>(P->opaque)->pj[i], errlev); + } + + return pj_default_destructor(P, errlev); +} + + + +/* + Zones: + + -180 -40 180 + +--------------+-------------------------+ Zones 1,2,9,10,11 & 12: + |1 |2 | Mollweide projection + | | | + +--------------+-------------------------+ Zones 3,4,5,6,7 & 8: + |3 |4 | Sinusoidal projection + | | | + 0 +-------+------+-+-----------+-----------+ + |5 |6 |7 |8 | + | | | | | + +-------+--------+-----------+-----------+ + |9 |10 |11 |12 | + | | | | | + +-------+--------+-----------+-----------+ + -180 -100 -20 80 180 +*/ + +#define SETUP(n, proj, x_0, y_0, lon_0) \ + if (!(Q->pj[n-1] = pj_##proj(0))) return destructor(P, ENOMEM); \ + if (!(Q->pj[n-1] = pj_##proj(Q->pj[n-1]))) return destructor(P, ENOMEM); \ + Q->pj[n-1]->ctx = P->ctx; \ + Q->pj[n-1]->x0 = x_0; \ + Q->pj[n-1]->y0 = y_0; \ + Q->pj[n-1]->lam0 = lon_0; + + +PJ *PROJECTION(igh) { + XY xy1, xy3; + LP lp = { 0, d4044118 }; + struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque))); + if (0==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + + + /* sinusoidal zones */ + SETUP(3, sinu, -d100, 0, -d100); + SETUP(4, sinu, d30, 0, d30); + SETUP(5, sinu, -d160, 0, -d160); + SETUP(6, sinu, -d60, 0, -d60); + SETUP(7, sinu, d20, 0, d20); + SETUP(8, sinu, d140, 0, d140); + + /* mollweide zones */ + SETUP(1, moll, -d100, 0, -d100); + + /* y0 ? */ + xy1 = Q->pj[0]->fwd(lp, Q->pj[0]); /* zone 1 */ + xy3 = Q->pj[2]->fwd(lp, Q->pj[2]); /* zone 3 */ + /* y0 + xy1.y = xy3.y for lt = 40d44'11.8" */ + Q->dy0 = xy3.y - xy1.y; + + Q->pj[0]->y0 = Q->dy0; + + /* mollweide zones (cont'd) */ + SETUP( 2, moll, d30, Q->dy0, d30); + SETUP( 9, moll, -d160, -Q->dy0, -d160); + SETUP(10, moll, -d60, -Q->dy0, -d60); + SETUP(11, moll, d20, -Q->dy0, d20); + SETUP(12, moll, d140, -Q->dy0, d140); + + P->inv = s_inverse; + P->fwd = s_forward; + P->destructor = destructor; + P->es = 0.; + + return P; +} |