/Proj4/PJ_laea.c

  1#ifndef lint
  2static const char SCCSID[]="@(#)PJ_laea.c	4.1	94/02/15	GIE	REL";
  3#endif
  4#define PROJ_PARMS__ \
  5	double	sinb1; \
  6	double	cosb1; \
  7	double	xmf; \
  8	double	ymf; \
  9	double	mmf; \
 10	double	qp; \
 11	double	dd; \
 12	double	rq; \
 13	double	*apa; \
 14	int		mode;
 15#define PJ_LIB__
 16#include	"projects.h"
 17PROJ_HEAD(laea, "Lambert Azimuthal Equal Area") "\n\tAzi, Sph&Ell";
 18#define sinph0	P->sinb1
 19#define cosph0	P->cosb1
 20#define EPS10	1.e-10
 21#define NITER	20
 22#define CONV	1.e-10
 23#define N_POLE	0
 24#define S_POLE	1
 25#define EQUIT	2
 26#define OBLIQ	3
 27FORWARD(e_forward); /* ellipsoid */
 28	double coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;
 29
 30	coslam = cos(lp.lam);
 31	sinlam = sin(lp.lam);
 32	sinphi = sin(lp.phi);
 33	q = pj_qsfn(sinphi, P->e, P->one_es);
 34	if (P->mode == OBLIQ || P->mode == EQUIT) {
 35		sinb = q / P->qp;
 36		cosb = sqrt(1. - sinb * sinb);
 37	}
 38	switch (P->mode) {
 39	case OBLIQ:
 40		b = 1. + P->sinb1 * sinb + P->cosb1 * cosb * coslam;
 41		break;
 42	case EQUIT:
 43		b = 1. + cosb * coslam;
 44		break;
 45	case N_POLE:
 46		b = HALFPI + lp.phi;
 47		q = P->qp - q;
 48		break;
 49	case S_POLE:
 50		b = lp.phi - HALFPI;
 51		q = P->qp + q;
 52		break;
 53	}
 54	if (fabs(b) < EPS10) F_ERROR;
 55	switch (P->mode) {
 56	case OBLIQ:
 57		xy.y = P->ymf * ( b = sqrt(2. / b) )
 58		   * (P->cosb1 * sinb - P->sinb1 * cosb * coslam);
 59		goto eqcon;
 60		break;
 61	case EQUIT:
 62		xy.y = (b = sqrt(2. / (1. + cosb * coslam))) * sinb * P->ymf; 
 63eqcon:
 64		xy.x = P->xmf * b * cosb * sinlam;
 65		break;
 66	case N_POLE:
 67	case S_POLE:
 68		if (q >= 0.) {
 69			xy.x = (b = sqrt(q)) * sinlam;
 70			xy.y = coslam * (P->mode == S_POLE ? b : -b);
 71		} else
 72			xy.x = xy.y = 0.;
 73		break;
 74	}
 75	return (xy);
 76}
 77FORWARD(s_forward); /* spheroid */
 78	double  coslam, cosphi, sinphi;
 79
 80	sinphi = sin(lp.phi);
 81	cosphi = cos(lp.phi);
 82	coslam = cos(lp.lam);
 83	switch (P->mode) {
 84	case EQUIT:
 85		xy.y = 1. + cosphi * coslam;
 86		goto oblcon;
 87	case OBLIQ:
 88		xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;
 89oblcon:
 90		if (xy.y <= EPS10) F_ERROR;
 91		xy.x = (xy.y = sqrt(2. / xy.y)) * cosphi * sin(lp.lam);
 92		xy.y *= P->mode == EQUIT ? sinphi :
 93		   cosph0 * sinphi - sinph0 * cosphi * coslam;
 94		break;
 95	case N_POLE:
 96		coslam = -coslam;
 97	case S_POLE:
 98		if (fabs(lp.phi + P->phi0) < EPS10) F_ERROR;
 99		xy.y = FORTPI - lp.phi * .5;
100		xy.y = 2. * (P->mode == S_POLE ? cos(xy.y) : sin(xy.y));
101		xy.x = xy.y * sin(lp.lam);
102		xy.y *= coslam;
103		break;
104	}
105	return (xy);
106}
107INVERSE(e_inverse); /* ellipsoid */
108	double cCe, sCe, q, rho, ab=0.0;
109
110	switch (P->mode) {
111	case EQUIT:
112	case OBLIQ:
113		if ((rho = hypot(xy.x /= P->dd, xy.y *=  P->dd)) < EPS10) {
114			lp.lam = 0.;
115			lp.phi = P->phi0;
116			return (lp);
117		}
118		cCe = cos(sCe = 2. * asin(.5 * rho / P->rq));
119		xy.x *= (sCe = sin(sCe));
120		if (P->mode == OBLIQ) {
121			q = P->qp * (ab = cCe * P->sinb1 + xy.y * sCe * P->cosb1 / rho);
122			xy.y = rho * P->cosb1 * cCe - xy.y * P->sinb1 * sCe;
123		} else {
124			q = P->qp * (ab = xy.y * sCe / rho);
125			xy.y = rho * cCe;
126		}
127		break;
128	case N_POLE:
129		xy.y = -xy.y;
130	case S_POLE:
131		if (!(q = (xy.x * xy.x + xy.y * xy.y)) ) {
132			lp.lam = 0.;
133			lp.phi = P->phi0;
134			return (lp);
135		}
136		/*
137		q = P->qp - q;
138		*/
139		ab = 1. - q / P->qp;
140		if (P->mode == S_POLE)
141			ab = - ab;
142		break;
143	}
144	lp.lam = atan2(xy.x, xy.y);
145	lp.phi = pj_authlat(asin(ab), P->apa);
146	return (lp);
147}
148INVERSE(s_inverse); /* spheroid */
149	double  cosz=0.0, rh, sinz=0.0;
150
151	rh = hypot(xy.x, xy.y);
152	if ((lp.phi = rh * .5 ) > 1.) I_ERROR;
153	lp.phi = 2. * asin(lp.phi);
154	if (P->mode == OBLIQ || P->mode == EQUIT) {
155		sinz = sin(lp.phi);
156		cosz = cos(lp.phi);
157	}
158	switch (P->mode) {
159	case EQUIT:
160		lp.phi = fabs(rh) <= EPS10 ? 0. : asin(xy.y * sinz / rh);
161		xy.x *= sinz;
162		xy.y = cosz * rh;
163		break;
164	case OBLIQ:
165		lp.phi = fabs(rh) <= EPS10 ? P->phi0 :
166		   asin(cosz * sinph0 + xy.y * sinz * cosph0 / rh);
167		xy.x *= sinz * cosph0;
168		xy.y = (cosz - sin(lp.phi) * sinph0) * rh;
169		break;
170	case N_POLE:
171		xy.y = -xy.y;
172		lp.phi = HALFPI - lp.phi;
173		break;
174	case S_POLE:
175		lp.phi -= HALFPI;
176		break;
177	}
178	lp.lam = (xy.y == 0. && (P->mode == EQUIT || P->mode == OBLIQ)) ?
179		0. : atan2(xy.x, xy.y);
180	return (lp);
181}
182FREEUP;
183    if (P) {
184		if (P->apa)
185			pj_dalloc(P->apa);
186		pj_dalloc(P);
187	}
188}
189ENTRY1(laea,apa)
190	double t;
191
192	if (fabs((t = fabs(P->phi0)) - HALFPI) < EPS10)
193		P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
194	else if (fabs(t) < EPS10)
195		P->mode = EQUIT;
196	else
197		P->mode = OBLIQ;
198	if (P->es) {
199		double sinphi;
200
201		P->e = sqrt(P->es);
202		P->qp = pj_qsfn(1., P->e, P->one_es);
203		P->mmf = .5 / (1. - P->es);
204		P->apa = pj_authset(P->es);
205		switch (P->mode) {
206		case N_POLE:
207		case S_POLE:
208			P->dd = 1.;
209			break;
210		case EQUIT:
211			P->dd = 1. / (P->rq = sqrt(.5 * P->qp));
212			P->xmf = 1.;
213			P->ymf = .5 * P->qp;
214			break;
215		case OBLIQ:
216			P->rq = sqrt(.5 * P->qp);
217			sinphi = sin(P->phi0);
218			P->sinb1 = pj_qsfn(sinphi, P->e, P->one_es) / P->qp;
219			P->cosb1 = sqrt(1. - P->sinb1 * P->sinb1);
220			P->dd = cos(P->phi0) / (sqrt(1. - P->es * sinphi * sinphi) *
221			   P->rq * P->cosb1);
222			P->ymf = (P->xmf = P->rq) / P->dd;
223			P->xmf *= P->dd;
224			break;
225		}
226		P->inv = e_inverse;
227		P->fwd = e_forward;
228	} else {
229		if (P->mode == OBLIQ) {
230			sinph0 = sin(P->phi0);
231			cosph0 = cos(P->phi0);
232		}
233		P->inv = s_inverse;
234		P->fwd = s_forward;
235	}
236ENDENTRY(P)