diff options
| author | Frank Warmerdam <warmerdam@pobox.com> | 2000-04-04 17:29:53 +0000 |
|---|---|---|
| committer | Frank Warmerdam <warmerdam@pobox.com> | 2000-04-04 17:29:53 +0000 |
| commit | 067352cad73074fe83ea5a779de03e68d36e5783 (patch) | |
| tree | 25bded117cfd8f2e1d3420e524e02c0c0ab7f6c2 /src | |
| parent | 060d898a3a1338b06b82c9c1d948253708085ca1 (diff) | |
| download | PROJ-067352cad73074fe83ea5a779de03e68d36e5783.tar.gz PROJ-067352cad73074fe83ea5a779de03e68d36e5783.zip | |
Applied patch to INVERSE() from Criag Bruce (cbruce@cubewerx.com)
in order to make it work proprely with coordinates close to zero.
git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@832 4e78687f-474d-0410-85f9-8d5e500ac6b2
Diffstat (limited to 'src')
| -rw-r--r-- | src/PJ_ortho.c | 76 |
1 files changed, 41 insertions, 35 deletions
diff --git a/src/PJ_ortho.c b/src/PJ_ortho.c index 98abf731..185e6ce6 100644 --- a/src/PJ_ortho.c +++ b/src/PJ_ortho.c @@ -38,44 +38,50 @@ FORWARD(s_forward); /* spheroid */ xy.x = cosphi * sin(lp.lam); return (xy); } + INVERSE(s_inverse); /* spheroid */ - double rh, cosc, sinc; + double rh, cosc, sinc; - if ((sinc = (rh = hypot(xy.x, xy.y))) > 1.) { - if ((sinc - 1.) > EPS10) I_ERROR; - sinc = 1.; - } - cosc = sqrt(1. - sinc * sinc); /* in this range OK */ - if (fabs(rh) <= EPS10) - lp.phi = P->phi0; - else switch (P->mode) { - case N_POLE: - xy.y = -xy.y; - lp.phi = acos(sinc); - break; - case S_POLE: - lp.phi = - acos(sinc); - break; - case EQUIT: - lp.phi = xy.y * sinc / rh; - xy.x *= sinc; - xy.y = cosc * rh; - goto sinchk; - case OBLIQ: - lp.phi = cosc * P->sinph0 + xy.y * sinc * P->cosph0 / rh; - xy.y = (cosc - P->sinph0 * lp.phi) * rh; - xy.x *= sinc * P->cosph0; -sinchk: - if (fabs(lp.phi) >= 1.) - lp.phi = lp.phi < 0. ? -HALFPI : HALFPI; - else - lp.phi = asin(lp.phi); - break; - } - lp.lam = (xy.y == 0. && (P->mode == OBLIQ || P->mode == EQUIT)) ? - (xy.x == 0. ? 0. : xy.x < 0. ? -HALFPI : HALFPI) : atan2(xy.x, xy.y); - return (lp); + if ((sinc = (rh = hypot(xy.x, xy.y))) > 1.) { + if ((sinc - 1.) > EPS10) I_ERROR; + sinc = 1.; + } + cosc = sqrt(1. - sinc * sinc); /* in this range OK */ + if (fabs(rh) <= EPS10) { + lp.phi = P->phi0; + lp.lam = 0.0; + } else { + switch (P->mode) { + case N_POLE: + xy.y = -xy.y; + lp.phi = acos(sinc); + break; + case S_POLE: + lp.phi = - acos(sinc); + break; + case EQUIT: + lp.phi = xy.y * sinc / rh; + xy.x *= sinc; + xy.y = cosc * rh; + goto sinchk; + case OBLIQ: + lp.phi = cosc * P->sinph0 + xy.y * sinc * P->cosph0 /rh; + xy.y = (cosc - P->sinph0 * lp.phi) * rh; + xy.x *= sinc * P->cosph0; + sinchk: + if (fabs(lp.phi) >= 1.) + lp.phi = lp.phi < 0. ? -HALFPI : HALFPI; + else + lp.phi = asin(lp.phi); + break; + } + lp.lam = (xy.y == 0. && (P->mode == OBLIQ || P->mode == EQUIT)) + ? (xy.x == 0. ? 0. : xy.x < 0. ? -HALFPI : HALFPI) + : atan2(xy.x, xy.y); + } + return (lp); } + FREEUP; if (P) pj_dalloc(P); } ENTRY0(ortho) if (fabs(fabs(P->phi0) - HALFPI) <= EPS10) |