/Proj4/PJ_tmerc.c
C | 148 lines | 143 code | 5 blank | 0 comment | 22 complexity | e34753a3c0172d7c2f4850c9462f3595 MD5 | raw file
1#ifndef lint 2static const char SCCSID[]="@(#)PJ_tmerc.c 4.2 94/06/02 GIE REL"; 3#endif 4#define PROJ_PARMS__ \ 5 double esp; \ 6 double ml0; \ 7 double *en; 8#define PJ_LIB__ 9#include "projects.h" 10PROJ_HEAD(tmerc, "Transverse Mercator") "\n\tCyl, Sph&Ell"; 11PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)") 12 "\n\tCyl, Sph\n\tzone= south"; 13#define EPS10 1.e-10 14#define aks0 P->esp 15#define aks5 P->ml0 16#define FC1 1. 17#define FC2 .5 18#define FC3 .16666666666666666666 19#define FC4 .08333333333333333333 20#define FC5 .05 21#define FC6 .03333333333333333333 22#define FC7 .02380952380952380952 23#define FC8 .01785714285714285714 24FORWARD(e_forward); /* ellipse */ 25 double al, als, n, cosphi, sinphi, t; 26 27 sinphi = sin(lp.phi); cosphi = cos(lp.phi); 28 t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.; 29 t *= t; 30 al = cosphi * lp.lam; 31 als = al * al; 32 al /= sqrt(1. - P->es * sinphi * sinphi); 33 n = P->esp * cosphi * cosphi; 34 xy.x = P->k0 * al * (FC1 + 35 FC3 * als * (1. - t + n + 36 FC5 * als * (5. + t * (t - 18.) + n * (14. - 58. * t) 37 + FC7 * als * (61. + t * ( t * (179. - t) - 479. ) ) 38 ))); 39 xy.y = P->k0 * (pj_mlfn(lp.phi, sinphi, cosphi, P->en) - P->ml0 + 40 sinphi * al * lp.lam * FC2 * ( 1. + 41 FC4 * als * (5. - t + n * (9. + 4. * n) + 42 FC6 * als * (61. + t * (t - 58.) + n * (270. - 330 * t) 43 + FC8 * als * (1385. + t * ( t * (543. - t) - 3111.) ) 44 )))); 45 return (xy); 46} 47FORWARD(s_forward); /* sphere */ 48 double b, cosphi; 49 50 b = (cosphi = cos(lp.phi)) * sin(lp.lam); 51 if (fabs(fabs(b) - 1.) <= EPS10) F_ERROR; 52 xy.x = aks5 * log((1. + b) / (1. - b)); 53 if ((b = fabs( xy.y = cosphi * cos(lp.lam) / sqrt(1. - b * b) )) >= 1.) { 54 if ((b - 1.) > EPS10) F_ERROR 55 else xy.y = 0.; 56 } else 57 xy.y = acos(xy.y); 58 if (lp.phi < 0.) xy.y = -xy.y; 59 xy.y = aks0 * (xy.y - P->phi0); 60 return (xy); 61} 62INVERSE(e_inverse); /* ellipsoid */ 63 double n, con, cosphi, d, ds, sinphi, t; 64 65 lp.phi = pj_inv_mlfn(P->ml0 + xy.y / P->k0, P->es, P->en); 66 if (fabs(lp.phi) >= HALFPI) { 67 lp.phi = xy.y < 0. ? -HALFPI : HALFPI; 68 lp.lam = 0.; 69 } else { 70 sinphi = sin(lp.phi); 71 cosphi = cos(lp.phi); 72 t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.; 73 n = P->esp * cosphi * cosphi; 74 d = xy.x * sqrt(con = 1. - P->es * sinphi * sinphi) / P->k0; 75 con *= t; 76 t *= t; 77 ds = d * d; 78 lp.phi -= (con * ds / (1.-P->es)) * FC2 * (1. - 79 ds * FC4 * (5. + t * (3. - 9. * n) + n * (1. - 4 * n) - 80 ds * FC6 * (61. + t * (90. - 252. * n + 81 45. * t) + 46. * n 82 - ds * FC8 * (1385. + t * (3633. + t * (4095. + 1574. * t)) ) 83 ))); 84 lp.lam = d*(FC1 - 85 ds*FC3*( 1. + 2.*t + n - 86 ds*FC5*(5. + t*(28. + 24.*t + 8.*n) + 6.*n 87 - ds * FC7 * (61. + t * (662. + t * (1320. + 720. * t)) ) 88 ))) / cosphi; 89 } 90 return (lp); 91} 92INVERSE(s_inverse); /* sphere */ 93 double h, g; 94 95 h = exp(xy.x / aks0); 96 g = .5 * (h - 1. / h); 97 h = cos(P->phi0 + xy.y / aks0); 98 lp.phi = asin(sqrt((1. - h * h) / (1. + g * g))); 99 if (xy.y < 0.) lp.phi = -lp.phi; 100 lp.lam = (g || h) ? atan2(g, h) : 0.; 101 return (lp); 102} 103FREEUP; 104 if (P) { 105 if (P->en) 106 pj_dalloc(P->en); 107 pj_dalloc(P); 108 } 109} 110 static PJ * 111setup(PJ *P) { /* general initialization */ 112 if (P->es) { 113 if (!(P->en = pj_enfn(P->es))) 114 E_ERROR_0; 115 P->ml0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), P->en); 116 P->esp = P->es / (1. - P->es); 117 P->inv = e_inverse; 118 P->fwd = e_forward; 119 } else { 120 aks0 = P->k0; 121 aks5 = .5 * aks0; 122 P->inv = s_inverse; 123 P->fwd = s_forward; 124 } 125 return P; 126} 127ENTRY1(tmerc, en) 128ENDENTRY(setup(P)) 129ENTRY1(utm, en) 130 int zone; 131 132 if (!P->es) E_ERROR(-34); 133 P->y0 = pj_param(P->params, "bsouth").i ? 10000000. : 0.; 134 P->x0 = 500000.; 135 if (pj_param(P->params, "tzone").i) /* zone input ? */ 136 if ((zone = pj_param(P->params, "izone").i) > 0 && zone <= 60) 137 --zone; 138 else 139 E_ERROR(-35) 140 else /* nearest central meridian input */ 141 if ((zone = floor((adjlon(P->lam0) + PI) * 30. / PI)) < 0) 142 zone = 0; 143 else if (zone >= 60) 144 zone = 59; 145 P->lam0 = (zone + .5) * PI / 30. - PI; 146 P->k0 = 0.9996; 147 P->phi0 = 0.; 148ENDENTRY(setup(P))