diff options
Diffstat (limited to 'src/PJ_mod_ster.c')
| -rw-r--r-- | src/PJ_mod_ster.c | 214 |
1 files changed, 214 insertions, 0 deletions
diff --git a/src/PJ_mod_ster.c b/src/PJ_mod_ster.c new file mode 100644 index 00000000..04fecb68 --- /dev/null +++ b/src/PJ_mod_ster.c @@ -0,0 +1,214 @@ +#ifndef lint +static const char SCCSID[]="@(#)PJ_mod_ster.c 4.1 94/02/15 GIE REL"; +#endif +/* based upon Snyder and Linck, USGS-NMD */ +#define PROJ_PARMS__ \ + COMPLEX *zcoeff; \ + double cchio, schio; \ + int n; +#define PJ_LIB__ +#include <projects.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. Stererographics of 48 U.S.") "\n\tAzi(mod)"; +PROJ_HEAD(alsk, "Mod. Stererographics of Alaska") "\n\tAzi(mod)"; +PROJ_HEAD(gs50, "Mod. Stererographics of 50 U.S.") "\n\tAzi(mod)"; +#define EPSLN 1e-10 + +FORWARD(e_forward); /* ellipsoid */ + 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((HALFPI + lp.phi) * .5) * + pow((1. - esphi) / (1. + esphi), P->e * .5)) - HALFPI; + schi = sin(chi); + cchi = cos(chi); + s = 2. / (1. + P->schio * schi + P->cchio * cchi * coslon); + p.r = s * cchi * sinlon; + p.i = s * (P->cchio * schi - P->schio * cchi * coslon); + p = pj_zpoly1(p, P->zcoeff, P->n); + xy.x = p.r; + xy.y = p.i; + return xy; +} +INVERSE(e_inverse); /* ellipsoid */ + int nn; + COMPLEX p, fxy, fpxy, dp; + double den, rh, z, sinz, cosz, chi, phi, dphi, esphi; + + p.r = xy.x; + p.i = xy.y; + for (nn = 20; nn ;--nn) { + fxy = pj_zpolyd1(p, P->zcoeff, P->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) { + lp.phi = P->phi0; + return lp; + } + chi = aasin(cosz * P->schio + p.i * sinz * P->cchio / rh); + phi = chi; + for (nn = 20; nn ;--nn) { + esphi = P->e * sin(phi); + dphi = 2. * atan(tan((HALFPI + chi) * .5) * + pow((1. + esphi) / (1. - esphi), P->e * .5)) - HALFPI - phi; + phi += dphi; + if (fabs(dphi) <= EPSLN) + break; + } + } + if (nn) { + lp.phi = phi; + lp.lam = atan2(p.r * sinz, rh * P->cchio * cosz - p.i * + P->schio * sinz); + } else + lp.lam = lp.phi = HUGE_VAL; + return lp; +} +FREEUP; if (P) pj_dalloc(P); } + static PJ * +setup(PJ *P) { /* general initialization */ + double esphi, chio; + + if (P->es) { + esphi = P->e * sin(P->phi0); + chio = 2. * atan(tan((HALFPI + P->phi0) * .5) * + pow((1. - esphi) / (1. + esphi), P->e * .5)) - HALFPI; + } else + chio = P->phi0; + P->schio = sin(chio); + P->cchio = cos(chio); + P->inv = e_inverse; P->fwd = e_forward; + return P; +} +ENTRY0(mil_os) + static COMPLEX /* Miller Oblated Stereographic */ +AB[] = { + 0.924500, 0., + 0., 0., + 0.019430, 0. +}; + + P->n = 2; + P->lam0 = DEG_TO_RAD * 20.; + P->phi0 = DEG_TO_RAD * 18.; + P->zcoeff = AB; + P->es = 0.; +ENDENTRY(setup(P)) +ENTRY0(lee_os) + static COMPLEX /* Lee Oblated Stereographic */ +AB[] = { + 0.721316, 0., + 0., 0., + -0.0088162, -0.00617325 +}; + + P->n = 2; + P->lam0 = DEG_TO_RAD * -165.; + P->phi0 = DEG_TO_RAD * -10.; + P->zcoeff = AB; + P->es = 0.; +ENDENTRY(setup(P)) +ENTRY0(gs48) + static COMPLEX /* 48 United States */ +AB[] = { + 0.98879, 0., + 0., 0., + -0.050909, 0., + 0., 0., + 0.075528, 0. +}; + + P->n = 4; + P->lam0 = DEG_TO_RAD * -96.; + P->phi0 = DEG_TO_RAD * -39.; + P->zcoeff = AB; + P->es = 0.; + P->a = 6370997.; +ENDENTRY(setup(P)) +ENTRY0(alsk) + static COMPLEX +ABe[] = { /* Alaska ellipsoid */ + .9945303, 0., + .0052083, -.0027404, + .0072721, .0048181, + -.0151089, -.1932526, + .0642675, -.1381226, + .3582802, -.2884586}, +ABs[] = { /* Alaska sphere */ + .9972523, 0., + .0052513, -.0041175, + .0074606, .0048125, + -.0153783, -.1968253, + .0636871, -.1408027, + .3660976, -.2937382 +}; + + P->n = 5; + P->lam0 = DEG_TO_RAD * -152.; + P->phi0 = DEG_TO_RAD * 64.; + if (P->es) { /* fixed ellipsoid/sphere */ + P->zcoeff = ABe; + P->a = 6378206.4; + P->e = sqrt(P->es = 0.00676866); + } else { + P->zcoeff = ABs; + P->a = 6370997.; + } +ENDENTRY(setup(P)) +ENTRY0(gs50) + static COMPLEX +ABe[] = { /* GS50 ellipsoid */ + .9827497, 0., + .0210669, .0053804, + -.1031415, -.0571664, + -.0323337, -.0322847, + .0502303, .1211983, + .0251805, .0895678, + -.0012315, -.1416121, + .0072202, -.1317091, + -.0194029, .0759677, + -.0210072, .0834037 +}, +ABs[] = { /* GS50 sphere */ + .9842990, 0., + .0211642, .0037608, + -.1036018, -.0575102, + -.0329095, -.0320119, + .0499471, .1223335, + .0260460, .0899805, + .0007388, -.1435792, + .0075848, -.1334108, + -.0216473, .0776645, + -.0225161, .0853673 +}; + + P->n = 9; + P->lam0 = DEG_TO_RAD * -120.; + P->phi0 = DEG_TO_RAD * 45.; + if (P->es) { /* fixed ellipsoid/sphere */ + P->zcoeff = ABe; + P->a = 6378206.4; + P->e = sqrt(P->es = 0.00676866); + } else { + P->zcoeff = ABs; + P->a = 6370997.; + } +ENDENTRY(setup(P)) |