/Proj4/PJ_aea.c
C | 171 lines | 124 code | 9 blank | 38 comment | 22 complexity | f2b6cff19e17b55170f35cd81aac6ea4 MD5 | raw file
1/****************************************************************************** 2 * $Id: PJ_aea.c,v 1.4 2003/08/18 15:21:23 warmerda Exp $ 3 * 4 * Project: PROJ.4 5 * Purpose: Implementation of the aea (Albers Equal Area) projection. 6 * Author: Gerald Evenden 7 * 8 ****************************************************************************** 9 * Copyright (c) 1995, Gerald Evenden 10 * 11 * Permission is hereby granted, free of charge, to any person obtaining a 12 * copy of this software and associated documentation files (the "Software"), 13 * to deal in the Software without restriction, including without limitation 14 * the rights to use, copy, modify, merge, publish, distribute, sublicense, 15 * and/or sell copies of the Software, and to permit persons to whom the 16 * Software is furnished to do so, subject to the following conditions: 17 * 18 * The above copyright notice and this permission notice shall be included 19 * in all copies or substantial portions of the Software. 20 * 21 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 22 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 23 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 24 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 25 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 26 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 27 * DEALINGS IN THE SOFTWARE. 28 ****************************************************************************** 29 * 30 * $Log: PJ_aea.c,v $ 31 * Revision 1.4 2003/08/18 15:21:23 warmerda 32 * fixed initialization of en variable 33 * 34 * Revision 1.3 2002/12/14 19:27:06 warmerda 35 * updated header 36 * 37 */ 38 39#define PROJ_PARMS__ \ 40 double ec; \ 41 double n; \ 42 double c; \ 43 double dd; \ 44 double n2; \ 45 double rho0; \ 46 double rho; \ 47 double phi1; \ 48 double phi2; \ 49 double *en; \ 50 int ellips; 51 52#define PJ_LIB__ 53#include "projects.h" 54 55PJ_CVSID("$Id: PJ_aea.c,v 1.4 2003/08/18 15:21:23 warmerda Exp $"); 56 57# define EPS10 1.e-10 58# define TOL7 1.e-7 59 60PROJ_HEAD(aea, "Albers Equal Area") 61 "\n\tConic Sph&Ell\n\tlat_1= lat_2="; 62PROJ_HEAD(leac, "Lambert Equal Area Conic") 63 "\n\tConic, Sph&Ell\n\tlat_1= south"; 64/* determine latitude angle phi-1 */ 65# define N_ITER 15 66# define EPSILON 1.0e-7 67# define TOL 1.0e-10 68 static double 69phi1_(double qs, double Te, double Tone_es) { 70 int i; 71 double Phi, sinpi, cospi, con, com, dphi; 72 73 Phi = asin (.5 * qs); 74 if (Te < EPSILON) 75 return( Phi ); 76 i = N_ITER; 77 do { 78 sinpi = sin (Phi); 79 cospi = cos (Phi); 80 con = Te * sinpi; 81 com = 1. - con * con; 82 dphi = .5 * com * com / cospi * (qs / Tone_es - 83 sinpi / com + .5 / Te * log ((1. - con) / 84 (1. + con))); 85 Phi += dphi; 86 } while (fabs(dphi) > TOL && --i); 87 return( i ? Phi : HUGE_VAL ); 88} 89FORWARD(e_forward); /* ellipsoid & spheroid */ 90 if ((P->rho = P->c - (P->ellips ? P->n * pj_qsfn(sin(lp.phi), 91 P->e, P->one_es) : P->n2 * sin(lp.phi))) < 0.) F_ERROR 92 P->rho = P->dd * sqrt(P->rho); 93 xy.x = P->rho * sin( lp.lam *= P->n ); 94 xy.y = P->rho0 - P->rho * cos(lp.lam); 95 return (xy); 96} 97INVERSE(e_inverse) /* ellipsoid & spheroid */; 98 if( (P->rho = hypot(xy.x, xy.y = P->rho0 - xy.y)) != 0.0 ) { 99 if (P->n < 0.) { 100 P->rho = -P->rho; 101 xy.x = -xy.x; 102 xy.y = -xy.y; 103 } 104 lp.phi = P->rho / P->dd; 105 if (P->ellips) { 106 lp.phi = (P->c - lp.phi * lp.phi) / P->n; 107 if (fabs(P->ec - fabs(lp.phi)) > TOL7) { 108 if ((lp.phi = phi1_(lp.phi, P->e, P->one_es)) == HUGE_VAL) 109 I_ERROR 110 } else 111 lp.phi = lp.phi < 0. ? -HALFPI : HALFPI; 112 } else if (fabs(lp.phi = (P->c - lp.phi * lp.phi) / P->n2) <= 1.) 113 lp.phi = asin(lp.phi); 114 else 115 lp.phi = lp.phi < 0. ? -HALFPI : HALFPI; 116 lp.lam = atan2(xy.x, xy.y) / P->n; 117 } else { 118 lp.lam = 0.; 119 lp.phi = P->n > 0. ? HALFPI : - HALFPI; 120 } 121 return (lp); 122} 123FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } } 124 static PJ * 125setup(PJ *P) { 126 double cosphi, sinphi; 127 int secant; 128 129 if (fabs(P->phi1 + P->phi2) < EPS10) E_ERROR(-21); 130 P->n = sinphi = sin(P->phi1); 131 cosphi = cos(P->phi1); 132 secant = fabs(P->phi1 - P->phi2) >= EPS10; 133 if( (P->ellips = (P->es > 0.))) { 134 double ml1, m1; 135 136 if (!(P->en = pj_enfn(P->es))) E_ERROR_0; 137 m1 = pj_msfn(sinphi, cosphi, P->es); 138 ml1 = pj_qsfn(sinphi, P->e, P->one_es); 139 if (secant) { /* secant cone */ 140 double ml2, m2; 141 142 sinphi = sin(P->phi2); 143 cosphi = cos(P->phi2); 144 m2 = pj_msfn(sinphi, cosphi, P->es); 145 ml2 = pj_qsfn(sinphi, P->e, P->one_es); 146 P->n = (m1 * m1 - m2 * m2) / (ml2 - ml1); 147 } 148 P->ec = 1. - .5 * P->one_es * log((1. - P->e) / 149 (1. + P->e)) / P->e; 150 P->c = m1 * m1 + P->n * ml1; 151 P->dd = 1. / P->n; 152 P->rho0 = P->dd * sqrt(P->c - P->n * pj_qsfn(sin(P->phi0), 153 P->e, P->one_es)); 154 } else { 155 if (secant) P->n = .5 * (P->n + sin(P->phi2)); 156 P->n2 = P->n + P->n; 157 P->c = cosphi * cosphi + P->n2 * sinphi; 158 P->dd = 1. / P->n; 159 P->rho0 = P->dd * sqrt(P->c - P->n2 * sin(P->phi0)); 160 } 161 P->inv = e_inverse; P->fwd = e_forward; 162 return P; 163} 164ENTRY1(aea,en) 165 P->phi1 = pj_param(P->params, "rlat_1").f; 166 P->phi2 = pj_param(P->params, "rlat_2").f; 167ENDENTRY(setup(P)) 168ENTRY1(leac,en) 169 P->phi2 = pj_param(P->params, "rlat_1").f; 170 P->phi1 = pj_param(P->params, "bsouth").i ? - HALFPI: HALFPI; 171ENDENTRY(setup(P))