/Proj4/PJ_aea.c

  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))