/Proj4/PJ_putp6.c

 1#ifndef lint
 2static const char SCCSID[]="@(#)PJ_putp6.c	4.1	94/02/15	GIE	REL";
 3#endif
 4#define PROJ_PARMS__ \
 5	double C_x, C_y, A, B, D;
 6#define PJ_LIB__
 7#include	"projects.h"
 8PROJ_HEAD(putp6, "Putnins P6") "\n\tPCyl., Sph.";
 9PROJ_HEAD(putp6p, "Putnins P6'") "\n\tPCyl., Sph.";
10#define EPS	1e-10
11#define NITER	10
12#define CON_POLE 1.732050807568877
13FORWARD(s_forward); /* spheroid */
14	double p, r, V;
15	int i;
16
17	p = P->B * sin(lp.phi);
18	lp.phi *=  1.10265779;
19	for (i = NITER; i ; --i) {
20		r = sqrt(1. + lp.phi * lp.phi);
21		lp.phi -= V = ( (P->A - r) * lp.phi - log(lp.phi + r) - p ) /
22			(P->A - 2. * r);
23		if (fabs(V) < EPS)
24			break;
25	}
26	if (!i)
27		lp.phi = p < 0. ? -CON_POLE : CON_POLE;
28	xy.x = P->C_x * lp.lam * (P->D - sqrt(1. + lp.phi * lp.phi));
29	xy.y = P->C_y * lp.phi;
30	return (xy);
31}
32INVERSE(s_inverse); /* spheroid */
33	double r;
34
35	lp.phi = xy.y / P->C_y;
36	r = sqrt(1. + lp.phi * lp.phi);
37	lp.lam = xy.x / (P->C_x * (P->D - r));
38	lp.phi = aasin( ( (P->A - r) * lp.phi - log(lp.phi + r) ) / P->B);
39	return (lp);
40}
41FREEUP; if (P) pj_dalloc(P); }
42	static PJ *
43setup(PJ *P) {
44	P->es = 0.;
45	P->inv = s_inverse;
46	P->fwd = s_forward;
47	return P;
48}
49ENTRY0(putp6)
50	P->C_x = 1.01346;
51	P->C_y = 0.91910;
52	P->A   = 4.;
53	P->B   = 2.1471437182129378784;
54	P->D   = 2.;
55ENDENTRY(setup(P))
56ENTRY0(putp6p)
57	P->C_x = 0.44329;
58	P->C_y = 0.80404;
59	P->A   = 6.;
60	P->B   = 5.61125;
61	P->D   = 3.;
62ENDENTRY(setup(P))