/Proj4/pj_mlfn.c

 1#ifndef lint
 2static const char SCCSID[]="@(#)pj_mlfn.c	4.5	95/07/06	GIE	REL";
 3#endif
 4#include "projects.h"
 5/* meridinal distance for ellipsoid and inverse
 6**	8th degree - accurate to < 1e-5 meters when used in conjuction
 7**		with typical major axis values.
 8**	Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
 9*/
10#define C00 1.
11#define C02 .25
12#define C04 .046875
13#define C06 .01953125
14#define C08 .01068115234375
15#define C22 .75
16#define C44 .46875
17#define C46 .01302083333333333333
18#define C48 .00712076822916666666
19#define C66 .36458333333333333333
20#define C68 .00569661458333333333
21#define C88 .3076171875
22#define EPS 1e-11
23#define MAX_ITER 10
24#define EN_SIZE 5
25	double *
26pj_enfn(double es) {
27	double t, *en;
28
29	en = (double *)pj_malloc(EN_SIZE * sizeof(double));
30	if (en) {
31		en[0] = C00 - es * (C02 + es * (C04 + es * (C06 + es * C08)));
32		en[1] = es * (C22 - es * (C04 + es * (C06 + es * C08)));
33		en[2] = (t = es * es) * (C44 - es * (C46 + es * C48));
34		en[3] = (t *= es) * (C66 - es * C68);
35		en[4] = t * es * C88;
36	} /* else return NULL if unable to allocate memory */
37	return en;
38}
39	double
40pj_mlfn(double phi, double sphi, double cphi, double *en) {
41	cphi *= sphi;
42	sphi *= sphi;
43	return(en[0] * phi - cphi * (en[1] + sphi*(en[2]
44		+ sphi*(en[3] + sphi*en[4]))));
45}
46	double
47pj_inv_mlfn(double arg, double es, double *en) {
48	double s, t, phi, k = 1./(1.-es);
49	int i;
50
51	phi = arg;
52	for (i = MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */
53		s = sin(phi);
54		t = 1. - es * s * s;
55		phi -= t = (pj_mlfn(phi, s, cos(phi), en) - arg) * (t * sqrt(t)) * k;
56		if (fabs(t) < EPS)
57			return phi;
58	}
59	pj_errno = -17;
60	return phi;
61}