PageRenderTime 28ms CodeModel.GetById 22ms app.highlight 4ms RepoModel.GetById 1ms app.codeStats 0ms

/Proj4/PJ_fouc_s.c

http://github.com/route-me/route-me
C | 48 lines | 46 code | 2 blank | 0 comment | 7 complexity | c30ebab0f1da09140f96f3685b57ddec MD5 | raw file
 1#ifndef lint
 2static const char SCCSID[]="@(#)PJ_fouc_s.c	4.1	94/02/15	GIE	REL";
 3#endif
 4#define PROJ_PARMS__ \
 5	double n, n1;
 6#define PJ_LIB__
 7#include	"projects.h"
 8PROJ_HEAD(fouc_s, "Foucaut Sinusoidal") "\n\tPCyl., Sph.";
 9#define MAX_ITER    10
10#define LOOP_TOL    1e-7
11FORWARD(s_forward); /* spheroid */
12	double t;
13
14	t = cos(lp.phi);
15	xy.x = lp.lam * t / (P->n + P->n1 * t);
16	xy.y = P->n * lp.phi + P->n1 * sin(lp.phi);
17	return (xy);
18}
19INVERSE(s_inverse); /* spheroid */
20	double V;
21	int i;
22
23	if (P->n) {
24		lp.phi = xy.y;
25		for (i = MAX_ITER; i ; --i) {
26			lp.phi -= V = (P->n * lp.phi + P->n1 * sin(lp.phi) - xy.y ) /
27				(P->n + P->n1 * cos(lp.phi));
28			if (fabs(V) < LOOP_TOL)
29				break;
30		}
31		if (!i)
32			lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
33	} else
34		lp.phi = aasin(xy.y);
35	V = cos(lp.phi);
36	lp.lam = xy.x * (P->n + P->n1 * V) / V;
37	return (lp);
38}
39FREEUP; if (P) pj_dalloc(P); }
40ENTRY0(fouc_s)
41	P->n = pj_param(P->params, "dn").f;
42	if (P->n < 0. || P->n > 1.)
43		E_ERROR(-99)
44	P->n1 = 1. - P->n;
45	P->es = 0;
46	P->inv = s_inverse;
47	P->fwd = s_forward;
48ENDENTRY(P)