/Proj4/PJ_mbtfpq.c

 1#ifndef lint
 2static const char SCCSID[]="@(#)PJ_mbtfpq.c	4.1 94/02/15     GIE     REL";
 3#endif
 4#define PJ_LIB__
 5#include	"projects.h"
 6PROJ_HEAD(mbtfpq, "McBryde-Thomas Flat-Polar Quartic") "\n\tCyl., Sph.";
 7#define NITER	20
 8#define EPS	1e-7
 9#define ONETOL 1.000001
10#define C	1.70710678118654752440
11#define RC	0.58578643762690495119
12#define FYC	1.87475828462269495505
13#define RYC	0.53340209679417701685
14#define FXC	0.31245971410378249250
15#define RXC	3.20041258076506210122
16FORWARD(s_forward); /* spheroid */
17	double th1, c;
18	int i;
19
20	c = C * sin(lp.phi);
21	for (i = NITER; i; --i) {
22		lp.phi -= th1 = (sin(.5*lp.phi) + sin(lp.phi) - c) /
23			(.5*cos(.5*lp.phi)  + cos(lp.phi));
24		if (fabs(th1) < EPS) break;
25	}
26	xy.x = FXC * lp.lam * (1.0 + 2. * cos(lp.phi)/cos(0.5 * lp.phi));
27	xy.y = FYC * sin(0.5 * lp.phi);
28	return (xy);
29}
30INVERSE(s_inverse); /* spheroid */
31	double t;
32
33	lp.phi = RYC * xy.y;
34	if (fabs(lp.phi) > 1.) {
35		if (fabs(lp.phi) > ONETOL)	I_ERROR
36		else if (lp.phi < 0.) { t = -1.; lp.phi = -PI; }
37		else { t = 1.; lp.phi = PI; }
38	} else
39		lp.phi = 2. * asin(t = lp.phi);
40	lp.lam = RXC * xy.x / (1. + 2. * cos(lp.phi)/cos(0.5 * lp.phi));
41	lp.phi = RC * (t + sin(lp.phi));
42	if (fabs(lp.phi) > 1.)
43		if (fabs(lp.phi) > ONETOL)	I_ERROR
44		else			lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
45	else
46		lp.phi = asin(lp.phi);
47	return (lp);
48}
49FREEUP; if (P) pj_dalloc(P); }
50ENTRY0(mbtfpq) P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)