/Proj4/PJ_labrd.c

  1#ifndef lint
  2static const char SCCSID[]="@(#)PJ_labrd.c	4.1	94/02/15	GIE	REL";
  3#endif
  4#define PROJ_PARMS__ \
  5	double	Az, kRg, p0s, A, C, Ca, Cb, Cc, Cd; \
  6	int		rot;
  7#define PJ_LIB__
  8#include	"projects.h"
  9PROJ_HEAD(labrd, "Laborde") "\n\tCyl, Sph\n\tSpecial for Madagascar";
 10#define EPS	1.e-10
 11FORWARD(e_forward);
 12	double V1, V2, ps, sinps, cosps, sinps2, cosps2, I1, I2, I3, I4, I5, I6,
 13		x2, y2, t;
 14
 15	V1 = P->A * log( tan(FORTPI + .5 * lp.phi) );
 16	t = P->e * sin(lp.phi);
 17	V2 = .5 * P->e * P->A * log ((1. + t)/(1. - t));
 18	ps = 2. * (atan(exp(V1 - V2 + P->C)) - FORTPI);
 19	I1 = ps - P->p0s;
 20	cosps = cos(ps);	cosps2 = cosps * cosps;
 21	sinps = sin(ps);	sinps2 = sinps * sinps;
 22	I4 = P->A * cosps;
 23	I2 = .5 * P->A * I4 * sinps;
 24	I3 = I2 * P->A * P->A * (5. * cosps2 - sinps2) / 12.;
 25	I6 = I4 * P->A * P->A;
 26	I5 = I6 * (cosps2 - sinps2) / 6.;
 27	I6 *= P->A * P->A *
 28		(5. * cosps2 * cosps2 + sinps2 * (sinps2 - 18. * cosps2)) / 120.;
 29	t = lp.lam * lp.lam;
 30	xy.x = P->kRg * lp.lam * (I4 + t * (I5 + t * I6));
 31	xy.y = P->kRg * (I1 + t * (I2 + t * I3));
 32	x2 = xy.x * xy.x;
 33	y2 = xy.y * xy.y;
 34	V1 = 3. * xy.x * y2 - xy.x * x2;
 35	V2 = xy.y * y2 - 3. * x2 * xy.y;
 36	xy.x += P->Ca * V1 + P->Cb * V2;
 37	xy.y += P->Ca * V2 - P->Cb * V1;
 38	return (xy);
 39}
 40INVERSE(e_inverse); /* ellipsoid & spheroid */
 41	double x2, y2, V1, V2, V3, V4, t, t2, ps, pe, tpe, s,
 42		I7, I8, I9, I10, I11, d, Re;
 43	int i;
 44
 45	x2 = xy.x * xy.x;
 46	y2 = xy.y * xy.y;
 47	V1 = 3. * xy.x * y2 - xy.x * x2;
 48	V2 = xy.y * y2 - 3. * x2 * xy.y;
 49	V3 = xy.x * (5. * y2 * y2 + x2 * (-10. * y2 + x2 ));
 50	V4 = xy.y * (5. * x2 * x2 + y2 * (-10. * x2 + y2 ));
 51	xy.x += - P->Ca * V1 - P->Cb * V2 + P->Cc * V3 + P->Cd * V4;
 52	xy.y +=   P->Cb * V1 - P->Ca * V2 - P->Cd * V3 + P->Cc * V4;
 53	ps = P->p0s + xy.y / P->kRg;
 54	pe = ps + P->phi0 - P->p0s;
 55	for ( i = 20; i; --i) {
 56		V1 = P->A * log(tan(FORTPI + .5 * pe));
 57		tpe = P->e * sin(pe);
 58		V2 = .5 * P->e * P->A * log((1. + tpe)/(1. - tpe));
 59		t = ps - 2. * (atan(exp(V1 - V2 + P->C)) - FORTPI);
 60		pe += t;
 61		if (fabs(t) < EPS)
 62			break;
 63	}
 64/*
 65	if (!i) {
 66	} else {
 67	}
 68*/
 69	t = P->e * sin(pe);
 70	t = 1. - t * t;
 71	Re = P->one_es / ( t * sqrt(t) );
 72	t = tan(ps);
 73	t2 = t * t;
 74	s = P->kRg * P->kRg;
 75	d = Re * P->k0 * P->kRg;
 76	I7 = t / (2. * d);
 77	I8 = t * (5. + 3. * t2) / (24. * d * s);
 78	d = cos(ps) * P->kRg * P->A;
 79	I9 = 1. / d;
 80	d *= s;
 81	I10 = (1. + 2. * t2) / (6. * d);
 82	I11 = (5. + t2 * (28. + 24. * t2)) / (120. * d * s);
 83	x2 = xy.x * xy.x;
 84	lp.phi = pe + x2 * (-I7 + I8 * x2);
 85	lp.lam = xy.x * (I9 + x2 * (-I10 + x2 * I11));
 86	return (lp);
 87}
 88FREEUP; if (P) pj_dalloc(P); }
 89ENTRY0(labrd)
 90	double Az, sinp, R, N, t;
 91
 92	P->rot	= pj_param(P->params, "bno_rot").i == 0;
 93	Az = pj_param(P->params, "razi").f;
 94	sinp = sin(P->phi0);
 95	t = 1. - P->es * sinp * sinp;
 96	N = 1. / sqrt(t);
 97	R = P->one_es * N / t;
 98	P->kRg = P->k0 * sqrt( N * R );
 99	P->p0s = atan( sqrt(R / N) * tan(P->phi0) );
100	P->A = sinp / sin(P->p0s);
101	t = P->e * sinp;
102	P->C = .5 * P->e * P->A * log((1. + t)/(1. - t)) +
103		- P->A * log( tan(FORTPI + .5 * P->phi0))
104		+ log( tan(FORTPI + .5 * P->p0s));
105	t = Az + Az;
106	P->Ca = (1. - cos(t)) * ( P->Cb = 1. / (12. * P->kRg * P->kRg) );
107	P->Cb *= sin(t);
108	P->Cc = 3. * (P->Ca * P->Ca - P->Cb * P->Cb);
109	P->Cd = 6. * P->Ca * P->Cb;
110	P->inv = e_inverse;
111	P->fwd = e_forward;
112ENDENTRY(P)