/Proj4/PJ_cea.c

 1#ifndef lint
 2static const char SCCSID[]="@(#)PJ_cea.c	4.1	94/02/15	GIE	REL";
 3#endif
 4#define PROJ_PARMS__ \
 5	double qp; \
 6	double *apa;
 7#define PJ_LIB__
 8# include	"projects.h"
 9PROJ_HEAD(cea, "Equal Area Cylindrical") "\n\tCyl, Sph&Ell\n\tlat_ts=";
10# define EPS	1e-10
11FORWARD(e_forward); /* spheroid */
12	xy.x = P->k0 * lp.lam;
13	xy.y = .5 * pj_qsfn(sin(lp.phi), P->e, P->one_es) / P->k0;
14	return (xy);
15}
16FORWARD(s_forward); /* spheroid */
17	xy.x = P->k0 * lp.lam;
18	xy.y = sin(lp.phi) / P->k0;
19	return (xy);
20}
21INVERSE(e_inverse); /* spheroid */
22	lp.phi = pj_authlat(asin( 2. * xy.y * P->k0 / P->qp), P->apa);
23	lp.lam = xy.x / P->k0;
24	return (lp);
25}
26INVERSE(s_inverse); /* spheroid */
27	double t;
28
29	if ((t = fabs(xy.y *= P->k0)) - EPS <= 1.) {
30		if (t >= 1.)
31			lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
32		else
33			lp.phi = asin(xy.y);
34		lp.lam = xy.x / P->k0;
35	} else I_ERROR;
36	return (lp);
37}
38FREEUP;
39	if (P) {
40		if (P->apa)
41			pj_dalloc(P->apa);
42		pj_dalloc(P);
43	}
44}
45ENTRY1(cea, apa)
46	double t;
47
48	if (pj_param(P->params, "tlat_ts").i &&
49		(P->k0 = cos(t = pj_param(P->params, "rlat_ts").f)) < 0.) E_ERROR(-24)
50	else
51		t = 0.;
52	if (P->es) {
53		t = sin(t);
54		P->k0 /= sqrt(1. - P->es * t * t);
55		P->e = sqrt(P->es);
56		if (!(P->apa = pj_authset(P->es))) E_ERROR_0;
57		P->qp = pj_qsfn(1., P->e, P->one_es);
58		P->inv = e_inverse;
59		P->fwd = e_forward;
60	} else {
61		P->inv = s_inverse;
62		P->fwd = s_forward;
63	}
64ENDENTRY(P)