1 /* 2 ** libproj -- library of cartographic projections 3 ** 4 ** Copyright (c) 2003, 2006 Gerald I. Evenden 5 */ 6 static const char 7 LIBPROJ_ID[] = "Id"; 8 /* 9 ** Permission is hereby granted, free of charge, to any person obtaining 10 ** a copy of this software and associated documentation files (the 11 ** "Software"), to deal in the Software without restriction, including 12 ** without limitation the rights to use, copy, modify, merge, publish, 13 ** distribute, sublicense, and/or sell copies of the Software, and to 14 ** permit persons to whom the Software is furnished to do so, subject to 15 ** the following conditions: 16 ** 17 ** The above copyright notice and this permission notice shall be 18 ** included in all copies or substantial portions of the Software. 19 ** 20 ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 21 ** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 22 ** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 23 ** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY 24 ** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, 25 ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE 26 ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 27 */ 28 #define TOL 1e-10 29 #define PROJ_PARMS__ \ 30 double M0, MP; \ 31 double qp; \ 32 void *en; \ 33 void *apa; 34 #define PROJ_LIB__ 35 # include <lib_proj.h> 36 PROJ_HEAD(tcea, "Transverse Cylindrical Equal-Area") 37 "\n\tCyl, Sph&Ell\n\tk_0=(1)"; 38 # define EPS 1e-10 39 FORWARD(e_forward); /* ellipsoid */ 40 double beta, betac, phic, sc, M; 41 42 beta = proj_auth_lat(lp.phi, P->apa); 43 if (fabs(fabs(lp.lam) - HALFPI) > TOL) { 44 betac = atan2(tan(beta),cos(lp.lam)); 45 phic = proj_auth_inv(betac, P->apa); 46 sc = sin(phic); 47 xy.x = cos(beta)*cos(phic)*sin(lp.lam)/ 48 (P->k0*cos(betac)*sqrt(1.-P->es*sc*sc)); 49 } else { 50 betac = phic = lp.phi >=0. ? HALFPI : -HALFPI; 51 sc = lp.phi >= 0. ? 1. : -1.; 52 xy.x = cos(beta)*sin(lp.lam)/ (P->k0*sqrt(1.-P->es)); 53 } 54 M = proj_mdist(phic, sc, cos(phic), P->en); 55 xy.y = P->k0 * (M - P->M0); 56 return (xy); 57 } 58 FORWARD(s_forward); /* spheroid */ 59 xy.x = cos(lp.phi)*sin(lp.lam)/P->k0; 60 xy.y = P->k0 * (atan2(tan(lp.phi), cos(lp.lam)) - P->phi0); 61 return (xy); 62 } 63 INVERSE(e_inverse); /* ellipsoid */ 64 double phic, betac, betap, beta, sc, t; 65 66 t = P->M0 + xy.y/P->k0; 67 phic = proj_inv_mdist(t, P->en); 68 sc = sin(phic); 69 betac = proj_auth_lat(phic, P->apa); 70 betap = -asin(P->k0 * xy.x * cos(betac) * sqrt(1.-P->es*sc*sc)/cos(phic)); 71 beta = asin(cos(betap) * sin(betac)); 72 lp.lam = - atan2(tan(betap),cos(betac)); 73 if (fabs(t) > P->MP) 74 lp.lam += lp.lam < 0. ? PI : - PI; 75 lp.phi = proj_auth_inv(beta, P->apa); 76 return (lp); 77 } 78 INVERSE(s_inverse); /* spheroid */ 79 double t, D, r; 80 81 t = xy.x * P->k0; 82 r = sqrt(1. - t * t); 83 D = xy.y / P->k0 + P->phi0; 84 lp.phi = asin(r * sin(D)); 85 lp.lam = atan2(t,(r * cos(D))); 86 return (lp); 87 } 88 FREEUP; 89 if (P) { 90 if (P->apa) 91 free(P->apa); 92 if (P->en) 93 free(P->en); 94 free(P); 95 } 96 } 97 ENTRY2(tcea, apa, en) 98 double t; 99 100 if (P->es) { 101 if (!(P->apa = proj_auth_ini(P->es, &t))) E_ERROR_0; 102 if (!(P->en = proj_mdist_ini(P->es))) E_ERROR_0; 103 P->qp = proj_qsfn(HALFPI, P->apa); 104 P->M0 = proj_mdist(P->phi0, sin(P->phi0), cos(P->phi0), P->en); 105 P->MP = proj_mdist(HALFPI, 1., 0., P->en); 106 P->inv = e_inverse; 107 P->fwd = e_forward; 108 } else { 109 P->inv = s_inverse; 110 P->fwd = s_forward; 111 } 112 ENDENTRY(P) 113 /* 114 ** Log: proj_tcea.c 115 ** Revision 3.1 2006/01/11 01:38:18 gie 116 ** Initial 117 ** 118 */ 119