New

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@776 4e78687f-474d-0410-85f9-8d5e500ac6b2

1 files changed, 60 insertions, 0 deletions

diff --git a/src/pj_mlfn.c b/src/pj_mlfn.c
new file mode 100644
index 00000000..0b980c06
--- /dev/null
+++ b/src/pj_mlfn.c
@@ -0,0 +1,60 @@
+#ifndef lint
+static const char SCCSID[]="@(#)pj_mlfn.c	4.5	95/07/06	GIE	REL";
+#endif
+#include <projects.h>
+/* meridinal distance for ellipsoid and inverse
+**	8th degree - accurate to < 1e-5 meters when used in conjuction
+**		with typical major axis values.
+**	Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
+*/
+#define C00 1.
+#define C02 .25
+#define C04 .046875
+#define C06 .01953125
+#define C08 .01068115234375
+#define C22 .75
+#define C44 .46875
+#define C46 .01302083333333333333
+#define C48 .00712076822916666666
+#define C66 .36458333333333333333
+#define C68 .00569661458333333333
+#define C88 .3076171875
+#define EPS 1e-11
+#define MAX_ITER 10
+#define EN_SIZE 5
+	double *
+pj_enfn(double es) {
+	double t, *en;
+
+	if (en = (double *)pj_malloc(EN_SIZE * sizeof(double))) {
+		en[0] = C00 - es * (C02 + es * (C04 + es * (C06 + es * C08)));
+		en[1] = es * (C22 - es * (C04 + es * (C06 + es * C08)));
+		en[2] = (t = es * es) * (C44 - es * (C46 + es * C48));
+		en[3] = (t *= es) * (C66 - es * C68);
+		en[4] = t * es * C88;
+	} /* else return NULL if unable to allocate memory */
+	return en;
+}
+	double
+pj_mlfn(double phi, double sphi, double cphi, double *en) {
+	cphi *= sphi;
+	sphi *= sphi;
+	return(en[0] * phi - cphi * (en[1] + sphi*(en[2]
+		+ sphi*(en[3] + sphi*en[4]))));
+}
+	double
+pj_inv_mlfn(double arg, double es, double *en) {
+	double s, t, phi, k = 1./(1.-es);
+	int i;
+
+	phi = arg;
+	for (i = MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */
+		s = sin(phi);
+		t = 1. - es * s * s;
+		phi -= t = (pj_mlfn(phi, s, cos(phi), en) - arg) * (t * sqrt(t)) * k;
+		if (fabs(t) < EPS)
+			return phi;
+	}
+	pj_errno = -17;
+	return phi;
+}