diff options
Diffstat (limited to 'src/3rdparty/proj/PJ_tmerc.c')
| -rw-r--r-- | src/3rdparty/proj/PJ_tmerc.c | 175 |
1 files changed, 0 insertions, 175 deletions
diff --git a/src/3rdparty/proj/PJ_tmerc.c b/src/3rdparty/proj/PJ_tmerc.c deleted file mode 100644 index 14cf2e8e..00000000 --- a/src/3rdparty/proj/PJ_tmerc.c +++ /dev/null @@ -1,175 +0,0 @@ -#define PROJ_PARMS__ \ - double esp; \ - double ml0; \ - double *en; -#define PJ_LIB__ -#include <projects.h> -PROJ_HEAD(tmerc, "Transverse Mercator") "\n\tCyl, Sph&Ell"; -PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)") - "\n\tCyl, Sph\n\tzone= south"; -#define EPS10 1.e-10 -#define aks0 P->esp -#define aks5 P->ml0 -#define FC1 1. -#define FC2 .5 -#define FC3 .16666666666666666666 -#define FC4 .08333333333333333333 -#define FC5 .05 -#define FC6 .03333333333333333333 -#define FC7 .02380952380952380952 -#define FC8 .01785714285714285714 -FORWARD(e_forward); /* ellipse */ - double al, als, n, cosphi, sinphi, t; - - /* - * Fail if our longitude is more than 90 degrees from the - * central meridian since the results are essentially garbage. - * Is error -20 really an appropriate return value? - * - * http://trac.osgeo.org/proj/ticket/5 - */ - if( lp.lam < -HALFPI || lp.lam > HALFPI ) - { - xy.x = HUGE_VAL; - xy.y = HUGE_VAL; - pj_errno = -14; - return xy; - } - - sinphi = sin(lp.phi); cosphi = cos(lp.phi); - t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.; - t *= t; - al = cosphi * lp.lam; - als = al * al; - al /= sqrt(1. - P->es * sinphi * sinphi); - n = P->esp * cosphi * cosphi; - xy.x = P->k0 * al * (FC1 + - FC3 * als * (1. - t + n + - FC5 * als * (5. + t * (t - 18.) + n * (14. - 58. * t) - + FC7 * als * (61. + t * ( t * (179. - t) - 479. ) ) - ))); - xy.y = P->k0 * (pj_mlfn(lp.phi, sinphi, cosphi, P->en) - P->ml0 + - sinphi * al * lp.lam * FC2 * ( 1. + - FC4 * als * (5. - t + n * (9. + 4. * n) + - FC6 * als * (61. + t * (t - 58.) + n * (270. - 330 * t) - + FC8 * als * (1385. + t * ( t * (543. - t) - 3111.) ) - )))); - return (xy); -} -FORWARD(s_forward); /* sphere */ - double b, cosphi; - - /* - * Fail if our longitude is more than 90 degrees from the - * central meridian since the results are essentially garbage. - * Is error -20 really an appropriate return value? - * - * http://trac.osgeo.org/proj/ticket/5 - */ - if( lp.lam < -HALFPI || lp.lam > HALFPI ) - { - xy.x = HUGE_VAL; - xy.y = HUGE_VAL; - pj_errno = -14; - return xy; - } - - b = (cosphi = cos(lp.phi)) * sin(lp.lam); - if (fabs(fabs(b) - 1.) <= EPS10) F_ERROR; - xy.x = aks5 * log((1. + b) / (1. - b)); - if ((b = fabs( xy.y = cosphi * cos(lp.lam) / sqrt(1. - b * b) )) >= 1.) { - if ((b - 1.) > EPS10) F_ERROR - else xy.y = 0.; - } else - xy.y = acos(xy.y); - if (lp.phi < 0.) xy.y = -xy.y; - xy.y = aks0 * (xy.y - P->phi0); - return (xy); -} -INVERSE(e_inverse); /* ellipsoid */ - double n, con, cosphi, d, ds, sinphi, t; - - lp.phi = pj_inv_mlfn(P->ml0 + xy.y / P->k0, P->es, P->en); - if (fabs(lp.phi) >= HALFPI) { - lp.phi = xy.y < 0. ? -HALFPI : HALFPI; - lp.lam = 0.; - } else { - sinphi = sin(lp.phi); - cosphi = cos(lp.phi); - t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.; - n = P->esp * cosphi * cosphi; - d = xy.x * sqrt(con = 1. - P->es * sinphi * sinphi) / P->k0; - con *= t; - t *= t; - ds = d * d; - lp.phi -= (con * ds / (1.-P->es)) * FC2 * (1. - - ds * FC4 * (5. + t * (3. - 9. * n) + n * (1. - 4 * n) - - ds * FC6 * (61. + t * (90. - 252. * n + - 45. * t) + 46. * n - - ds * FC8 * (1385. + t * (3633. + t * (4095. + 1574. * t)) ) - ))); - lp.lam = d*(FC1 - - ds*FC3*( 1. + 2.*t + n - - ds*FC5*(5. + t*(28. + 24.*t + 8.*n) + 6.*n - - ds * FC7 * (61. + t * (662. + t * (1320. + 720. * t)) ) - ))) / cosphi; - } - return (lp); -} -INVERSE(s_inverse); /* sphere */ - double h, g; - - h = exp(xy.x / aks0); - g = .5 * (h - 1. / h); - h = cos(P->phi0 + xy.y / aks0); - lp.phi = asin(sqrt((1. - h * h) / (1. + g * g))); - if (xy.y < 0.) lp.phi = -lp.phi; - lp.lam = (g || h) ? atan2(g, h) : 0.; - return (lp); -} -FREEUP; - if (P) { - if (P->en) - pj_dalloc(P->en); - pj_dalloc(P); - } -} - static PJ * -setup(PJ *P) { /* general initialization */ - if (P->es) { - if (!(P->en = pj_enfn(P->es))) - E_ERROR_0; - P->ml0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), P->en); - P->esp = P->es / (1. - P->es); - P->inv = e_inverse; - P->fwd = e_forward; - } else { - aks0 = P->k0; - aks5 = .5 * aks0; - P->inv = s_inverse; - P->fwd = s_forward; - } - return P; -} -ENTRY1(tmerc, en) -ENDENTRY(setup(P)) -ENTRY1(utm, en) - int zone; - - if (!P->es) E_ERROR(-34); - P->y0 = pj_param(P->params, "bsouth").i ? 10000000. : 0.; - P->x0 = 500000.; - if (pj_param(P->params, "tzone").i) /* zone input ? */ - if ((zone = pj_param(P->params, "izone").i) > 0 && zone <= 60) - --zone; - else - E_ERROR(-35) - else /* nearest central meridian input */ - if ((zone = floor((adjlon(P->lam0) + PI) * 30. / PI)) < 0) - zone = 0; - else if (zone >= 60) - zone = 59; - P->lam0 = (zone + .5) * PI / 30. - PI; - P->k0 = 0.9996; - P->phi0 = 0.; -ENDENTRY(setup(P)) |