diff options
| author | Even Rouault <even.rouault@spatialys.com> | 2018-12-19 13:00:37 +0100 |
|---|---|---|
| committer | Even Rouault <even.rouault@spatialys.com> | 2018-12-26 10:08:55 +0100 |
| commit | df574ae332d57f556fd56314883b3354cab1d0ff (patch) | |
| tree | 63a68d40d7ed7932d6329d9c7baa340bb6294f7f /src/projections/mod_ster.cpp | |
| parent | e6de172371ea203f6393d745641d66c82b5b13e2 (diff) | |
| download | PROJ-df574ae332d57f556fd56314883b3354cab1d0ff.tar.gz PROJ-df574ae332d57f556fd56314883b3354cab1d0ff.zip | |
cpp conversion: remove useless pj_, PJ_ and proj_ filename prefixes
Diffstat (limited to 'src/projections/mod_ster.cpp')
| -rw-r--r-- | src/projections/mod_ster.cpp | 282 |
1 files changed, 282 insertions, 0 deletions
diff --git a/src/projections/mod_ster.cpp b/src/projections/mod_ster.cpp new file mode 100644 index 00000000..7c4f363b --- /dev/null +++ b/src/projections/mod_ster.cpp @@ -0,0 +1,282 @@ +/* based upon Snyder and Linck, USGS-NMD */ +#define PJ_LIB__ +#include <errno.h> +#include "projects.h" +#include "proj_math.h" + +PROJ_HEAD(mil_os, "Miller Oblated Stereographic") "\n\tAzi(mod)"; +PROJ_HEAD(lee_os, "Lee Oblated Stereographic") "\n\tAzi(mod)"; +PROJ_HEAD(gs48, "Mod. Stereographic of 48 U.S.") "\n\tAzi(mod)"; +PROJ_HEAD(alsk, "Mod. Stereographic of Alaska") "\n\tAzi(mod)"; +PROJ_HEAD(gs50, "Mod. Stereographic of 50 U.S.") "\n\tAzi(mod)"; + +#define EPSLN 1e-12 + +namespace { // anonymous namespace +struct pj_opaque { + const COMPLEX *zcoeff; \ + double cchio, schio; \ + int n; +}; +} // anonymous namespace + + +static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ + XY xy = {0.0,0.0}; + struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); + double sinlon, coslon, esphi, chi, schi, cchi, s; + COMPLEX p; + + sinlon = sin(lp.lam); + coslon = cos(lp.lam); + esphi = P->e * sin(lp.phi); + chi = 2. * atan(tan((M_HALFPI + lp.phi) * .5) * + pow((1. - esphi) / (1. + esphi), P->e * .5)) - M_HALFPI; + schi = sin(chi); + cchi = cos(chi); + s = 2. / (1. + Q->schio * schi + Q->cchio * cchi * coslon); + p.r = s * cchi * sinlon; + p.i = s * (Q->cchio * schi - Q->schio * cchi * coslon); + p = pj_zpoly1(p, Q->zcoeff, Q->n); + xy.x = p.r; + xy.y = p.i; + + return xy; +} + + +static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ + LP lp = {0.0,0.0}; + struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); + int nn; + COMPLEX p, fxy, fpxy, dp; + double den, rh = 0.0, z, sinz = 0.0, cosz = 0.0, chi, phi = 0.0, esphi; + + p.r = xy.x; + p.i = xy.y; + for (nn = 20; nn ;--nn) { + fxy = pj_zpolyd1(p, Q->zcoeff, Q->n, &fpxy); + fxy.r -= xy.x; + fxy.i -= xy.y; + den = fpxy.r * fpxy.r + fpxy.i * fpxy.i; + dp.r = -(fxy.r * fpxy.r + fxy.i * fpxy.i) / den; + dp.i = -(fxy.i * fpxy.r - fxy.r * fpxy.i) / den; + p.r += dp.r; + p.i += dp.i; + if ((fabs(dp.r) + fabs(dp.i)) <= EPSLN) + break; + } + if (nn) { + rh = hypot(p.r, p.i); + z = 2. * atan(.5 * rh); + sinz = sin(z); + cosz = cos(z); + lp.lam = P->lam0; + if (fabs(rh) <= EPSLN) { + /* if we end up here input coordinates were (0,0). + * pj_inv() adds P->lam0 to lp.lam, this way we are + * sure to get the correct offset */ + lp.lam = 0.0; + lp.phi = P->phi0; + return lp; + } + chi = aasin(P->ctx, cosz * Q->schio + p.i * sinz * Q->cchio / rh); + phi = chi; + for (nn = 20; nn ;--nn) { + double dphi; + esphi = P->e * sin(phi); + dphi = 2. * atan(tan((M_HALFPI + chi) * .5) * + pow((1. + esphi) / (1. - esphi), P->e * .5)) - M_HALFPI - phi; + phi += dphi; + if (fabs(dphi) <= EPSLN) + break; + } + } + if (nn) { + lp.phi = phi; + lp.lam = atan2(p.r * sinz, rh * Q->cchio * cosz - p.i * + Q->schio * sinz); + } else + lp.lam = lp.phi = HUGE_VAL; + return lp; +} + + +static PJ *setup(PJ *P) { /* general initialization */ + struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); + double esphi, chio; + + if (P->es != 0.0) { + esphi = P->e * sin(P->phi0); + chio = 2. * atan(tan((M_HALFPI + P->phi0) * .5) * + pow((1. - esphi) / (1. + esphi), P->e * .5)) - M_HALFPI; + } else + chio = P->phi0; + Q->schio = sin(chio); + Q->cchio = cos(chio); + P->inv = e_inverse; + P->fwd = e_forward; + + return P; +} + + +/* Miller Oblated Stereographic */ +PJ *PROJECTION(mil_os) { + static const COMPLEX AB[] = { + {0.924500, 0.}, + {0., 0.}, + {0.019430, 0.} + }; + + struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + + Q->n = 2; + P->lam0 = DEG_TO_RAD * 20.; + P->phi0 = DEG_TO_RAD * 18.; + Q->zcoeff = AB; + P->es = 0.; + + return setup(P); +} + + +/* Lee Oblated Stereographic */ +PJ *PROJECTION(lee_os) { + static const COMPLEX AB[] = { + {0.721316, 0.}, + {0., 0.}, + {-0.0088162, -0.00617325} + }; + + struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + + Q->n = 2; + P->lam0 = DEG_TO_RAD * -165.; + P->phi0 = DEG_TO_RAD * -10.; + Q->zcoeff = AB; + P->es = 0.; + + return setup(P); +} + + +PJ *PROJECTION(gs48) { + static const COMPLEX /* 48 United States */ + AB[] = { + {0.98879, 0.}, + {0., 0.}, + {-0.050909, 0.}, + {0., 0.}, + {0.075528, 0.} + }; + + struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + + Q->n = 4; + P->lam0 = DEG_TO_RAD * -96.; + P->phi0 = DEG_TO_RAD * 39.; + Q->zcoeff = AB; + P->es = 0.; + P->a = 6370997.; + + return setup(P); +} + + +PJ *PROJECTION(alsk) { + static const COMPLEX ABe[] = { /* Alaska ellipsoid */ + { .9945303, 0.}, + { .0052083, -.0027404}, + { .0072721, .0048181}, + {-.0151089, -.1932526}, + { .0642675, -.1381226}, + { .3582802, -.2884586}, + }; + + static const COMPLEX ABs[] = { /* Alaska sphere */ + { .9972523, 0.}, + { .0052513, -.0041175}, + { .0074606, .0048125}, + {-.0153783, -.1968253}, + { .0636871, -.1408027}, + { .3660976, -.2937382} + }; + + struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + + Q->n = 5; + P->lam0 = DEG_TO_RAD * -152.; + P->phi0 = DEG_TO_RAD * 64.; + if (P->es != 0.0) { /* fixed ellipsoid/sphere */ + Q->zcoeff = ABe; + P->a = 6378206.4; + P->e = sqrt(P->es = 0.00676866); + } else { + Q->zcoeff = ABs; + P->a = 6370997.; + } + + return setup(P); +} + + +PJ *PROJECTION(gs50) { + static const COMPLEX ABe[] = { /* GS50 ellipsoid */ + { .9827497, 0.}, + { .0210669, .0053804}, + {-.1031415, -.0571664}, + {-.0323337, -.0322847}, + { .0502303, .1211983}, + { .0251805, .0895678}, + {-.0012315, -.1416121}, + { .0072202, -.1317091}, + {-.0194029, .0759677}, + {-.0210072, .0834037} + }; + + static const COMPLEX ABs[] = { /* GS50 sphere */ + { .9842990, 0.}, + { .0211642, .0037608}, + {-.1036018, -.0575102}, + {-.0329095, -.0320119}, + { .0499471, .1223335}, + { .0260460, .0899805}, + { .0007388, -.1435792}, + { .0075848, -.1334108}, + {-.0216473, .0776645}, + {-.0225161, .0853673} + }; + + struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + + Q->n = 9; + P->lam0 = DEG_TO_RAD * -120.; + P->phi0 = DEG_TO_RAD * 45.; + if (P->es != 0.0) { /* fixed ellipsoid/sphere */ + Q->zcoeff = ABe; + P->a = 6378206.4; + P->e = sqrt(P->es = 0.00676866); + } else { + Q->zcoeff = ABs; + P->a = 6370997.; + } + + return setup(P); +} + |
