/*
The Patterson projection was designed by Tom Patterson, US National Park
Service, in 2014. The polynomial equation was developed by Bojan Savric 
and Bernhard Jenny, College of Earth, Ocean, and Atmospheric Sciences, 
Oregon State University.
Port to PROJ.4 by Bernhard Jenny, 22 February 2015
*/

#define PJ_LIB__
#include	<projects.h>
PROJ_HEAD(patterson, "Patterson") "\n\tCyl., Sph.";

#define K1 1.0148
#define K2 0.23185
#define K3 -0.14499
#define K4 0.02406
#define C1 K1
#define C2 (5 * K2)
#define C3 (7 * K3)
#define C4 (9 * K4)
#define EPS 1e-11
#define MAX_Y 1.790857183

FORWARD(s_forward); /* spheroid */
	double lat_sq;

	lat_sq = lp.phi * lp.phi;
	xy.x = lp.lam;
	xy.y = lp.phi * (K1 + lat_sq * lat_sq * (K2 + lat_sq * (K3 + K4 * lat_sq)));
    return (xy);
}
INVERSE(s_inverse); /* spheroid */
	double yc, tol, y2, y4, f, fder;

    /* make sure y is inside valid range */
	if (xy.y > MAX_Y) {
		xy.y = MAX_Y;
	} else if (xy.y < -MAX_Y) {
		xy.y = -MAX_Y;
	}

    /* latitude */
	yc = xy.y;
    for (;;) { /* Newton-Raphson */
		y2 = yc * yc;
		y4 = y2 * y2;
		f = (yc * (K1 + y2 * y2 * (K2 + y2 * (K3 + K4 * y2)))) - xy.y;
		fder = C1 + y2 * y2 * (C2 + y2 * (C3 + C4 * y2));
		yc -= tol = f / fder;
		if (fabs(tol) < EPS) {
			break;
		}
	}
	lp.phi = yc;

    /* longitude */
	lp.lam = xy.x;

	return (lp);
}
FREEUP; if (P) pj_dalloc(P); }
ENTRY0(patterson) P->es = 0; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)