/Proj4/PJ_imw_p.c

  1#ifndef lint
  2static const char SCCSID[]="@(#)PJ_imw_p.c	4.1	94/05/22	GIE	REL";
  3#endif
  4#define PROJ_PARMS__ \
  5	double	P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2; \
  6	double	phi_1, phi_2, lam_1; \
  7	double	*en; \
  8	int	mode; /* = 0, phi_1 and phi_2 != 0, = 1, phi_1 = 0, = -1 phi_2 = 0 */
  9#define PJ_LIB__
 10#include	"projects.h"
 11PROJ_HEAD(imw_p, "International Map of the World Polyconic")
 12	"\n\tMod. Polyconic, Ell\n\tlat_1= and lat_2= [lon_1=]";
 13#define TOL 1e-10
 14#define EPS 1e-10
 15	static int
 16phi12(PJ *P, double *del, double *sig) {
 17	int err = 0;
 18
 19	if (!pj_param(P->params, "tlat_1").i ||
 20		!pj_param(P->params, "tlat_2").i) {
 21		err = -41;
 22	} else {
 23		P->phi_1 = pj_param(P->params, "rlat_1").f;
 24		P->phi_2 = pj_param(P->params, "rlat_2").f;
 25		*del = 0.5 * (P->phi_2 - P->phi_1);
 26		*sig = 0.5 * (P->phi_2 + P->phi_1);
 27		err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? -42 : 0;
 28	}
 29	return err;
 30}
 31	static XY
 32loc_for(LP lp, PJ *P, double *yc) {
 33	XY xy;
 34
 35	if (! lp.phi) {
 36		xy.x = lp.lam;
 37		xy.y = 0.;
 38	} else {
 39		double xa, ya, xb, yb, xc, yc, D, B, m, sp, t, R, C;
 40
 41		sp = sin(lp.phi);
 42		m = pj_mlfn(lp.phi, sp, cos(lp.phi), P->en);
 43		xa = P->Pp + P->Qp * m;
 44		ya = P->P + P->Q * m;
 45		R = 1. / (tan(lp.phi) * sqrt(1. - P->es * sp * sp));
 46		C = sqrt(R * R - xa * xa);
 47		if (lp.phi < 0.) C = - C;
 48		C += ya - R;
 49		if (P->mode < 0) {
 50			xb = lp.lam;
 51			yb = P->C2;
 52		} else {
 53			t = lp.lam * P->sphi_2;
 54			xb = P->R_2 * sin(t);
 55			yb = P->C2 + P->R_2 * (1. - cos(t));
 56		}
 57		if (P->mode > 0) {
 58			xc = lp.lam;
 59			yc = 0.;
 60		} else {
 61			t = lp.lam * P->sphi_1;
 62			xc = P->R_1 * sin(t);
 63			yc = P->R_1 * (1. - cos(t));
 64		}
 65		D = (xb - xc)/(yb - yc);
 66		B = xc + D * (C + R - yc);
 67		xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
 68		if (lp.phi > 0)
 69			xy.x = - xy.x;
 70		xy.x = (B + xy.x) / (1. + D * D);
 71		xy.y = sqrt(R * R - xy.x * xy.x);
 72		if (lp.phi > 0)
 73			xy.y = - xy.y;
 74		xy.y += C + R;
 75	}
 76	return (xy);
 77}
 78FORWARD(e_forward); /* ellipsoid */
 79	double yc;
 80	xy = loc_for(lp, P, &yc);
 81	return (xy);
 82}
 83INVERSE(e_inverse); /* ellipsoid */
 84	XY t;
 85	double yc;
 86
 87	lp.phi = P->phi_2;
 88	lp.lam = xy.x / cos(lp.phi);
 89	do {
 90		t = loc_for(lp, P, &yc);
 91		lp.phi = ((lp.phi - P->phi_1) * (xy.y - yc) / (t.y - yc)) + P->phi_1;
 92		lp.lam = lp.lam * xy.x / t.x;
 93	} while (fabs(t.x - xy.x) > TOL || fabs(t.y - xy.y) > TOL);
 94	return (lp);
 95}
 96	static void
 97xy(PJ *P, double phi, double *x, double *y, double *sp, double *R) {
 98	double F;
 99
100	*sp = sin(phi);
101	*R = 1./(tan(phi) * sqrt(1. - P->es * *sp * *sp ));
102	F = P->lam_1 * *sp;
103	*y = *R * (1 - cos(F));
104	*x = *R * sin(F);
105}
106FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
107ENTRY1(imw_p, en)
108	double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
109	int i;
110
111	if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
112	if( (i = phi12(P, &del, &sig)) != 0)
113		E_ERROR(i);
114	if (P->phi_2 < P->phi_1) { /* make sure P->phi_1 most southerly */
115		del = P->phi_1;
116		P->phi_1 = P->phi_2;
117		P->phi_2 = del;
118	}
119	if (pj_param(P->params, "tlon_1").i)
120		P->lam_1 = pj_param(P->params, "rlon_1").f;
121	else { /* use predefined based upon latitude */
122		sig = fabs(sig * RAD_TO_DEG);
123		if (sig <= 60)		sig = 2.;
124		else if (sig <= 76) sig = 4.;
125		else				sig = 8.;
126		P->lam_1 = sig * DEG_TO_RAD;
127	}
128	P->mode = 0;
129	if (P->phi_1) xy(P, P->phi_1, &x1, &y1, &P->sphi_1, &P->R_1);
130	else {
131		P->mode = 1;
132		y1 = 0.;
133		x1 = P->lam_1;
134	}
135	if (P->phi_2) xy(P, P->phi_2, &x2, &T2, &P->sphi_2, &P->R_2);
136	else {
137		P->mode = -1;
138		T2 = 0.;
139		x2 = P->lam_1;
140	}
141	m1 = pj_mlfn(P->phi_1, P->sphi_1, cos(P->phi_1), P->en);
142	m2 = pj_mlfn(P->phi_2, P->sphi_2, cos(P->phi_2), P->en);
143	t = m2 - m1;
144	s = x2 - x1;
145	y2 = sqrt(t * t - s * s) + y1;
146	P->C2 = y2 - T2;
147	t = 1. / t;
148	P->P = (m2 * y1 - m1 * y2) * t;
149	P->Q = (y2 - y1) * t;
150	P->Pp = (m2 * x1 - m1 * x2) * t;
151	P->Qp = (x2 - x1) * t;
152	P->fwd = e_forward;
153	P->inv = e_inverse;
154ENDENTRY(P)