/Proj4/PJ_chamb.c

  1#ifndef lint
  2static const char SCCSID[]="@(#)PJ_chamb.c	4.1	94/02/15	GIE	REL";
  3#endif
  4typedef struct { double r, Az; } VECT;
  5#define PROJ_PARMS__ \
  6	struct { /* control point data */ \
  7		double phi, lam; \
  8		double cosphi, sinphi; \
  9		VECT v; \
 10		XY	p; \
 11		double Az; \
 12	} c[3]; \
 13	XY p; \
 14	double beta_0, beta_1, beta_2;
 15#define PJ_LIB__
 16#include	"projects.h"
 17PROJ_HEAD(chamb, "Chamberlin Trimetric") "\n\tMisc Sph, no inv."
 18"\n\tlat_1= lon_1= lat_2= lon_2= lat_3= lon_3=";
 19#include	<stdio.h>
 20#define THIRD 0.333333333333333333
 21#define TOL 1e-9
 22	static VECT /* distance and azimuth from point 1 to point 2 */
 23vect(double dphi, double c1, double s1, double c2, double s2, double dlam) {
 24	VECT v;
 25	double cdl, dp, dl;
 26
 27	cdl = cos(dlam);
 28	if (fabs(dphi) > 1. || fabs(dlam) > 1.)
 29		v.r = aacos(s1 * s2 + c1 * c2 * cdl);
 30	else { /* more accurate for smaller distances */
 31		dp = sin(.5 * dphi);
 32		dl = sin(.5 * dlam);
 33		v.r = 2. * aasin(sqrt(dp * dp + c1 * c2 * dl * dl));
 34	}
 35	if (fabs(v.r) > TOL)
 36		v.Az = atan2(c2 * sin(dlam), c1 * s2 - s1 * c2 * cdl);
 37	else
 38		v.r = v.Az = 0.;
 39	return v;
 40}
 41	static double /* law of cosines */
 42lc(double b,double c,double a) {
 43	return aacos(.5 * (b * b + c * c - a * a) / (b * c));
 44}
 45FORWARD(s_forward); /* spheroid */
 46	double sinphi, cosphi, a;
 47	VECT v[3];
 48	int i, j;
 49
 50	sinphi = sin(lp.phi);
 51	cosphi = cos(lp.phi);
 52	for (i = 0; i < 3; ++i) { /* dist/azimiths from control */
 53		v[i] = vect(lp.phi - P->c[i].phi, P->c[i].cosphi, P->c[i].sinphi,
 54			cosphi, sinphi, lp.lam - P->c[i].lam);
 55		if ( ! v[i].r)
 56			break;
 57		v[i].Az = adjlon(v[i].Az - P->c[i].v.Az);
 58	}
 59	if (i < 3) /* current point at control point */
 60		xy = P->c[i].p;
 61	else { /* point mean of intersepts */
 62		xy = P->p;
 63		for (i = 0; i < 3; ++i) {
 64			j = i == 2 ? 0 : i + 1;
 65			a = lc(P->c[i].v.r, v[i].r, v[j].r);
 66			if (v[i].Az < 0.)
 67				a = -a;
 68			if (! i) { /* coord comp unique to each arc */
 69				xy.x += v[i].r * cos(a);
 70				xy.y -= v[i].r * sin(a);
 71			} else if (i == 1) {
 72				a = P->beta_1 - a;
 73				xy.x -= v[i].r * cos(a);
 74				xy.y -= v[i].r * sin(a);
 75			} else {
 76				a = P->beta_2 - a;
 77				xy.x += v[i].r * cos(a);
 78				xy.y += v[i].r * sin(a);
 79			}
 80		}
 81		xy.x *= THIRD; /* mean of arc intercepts */
 82		xy.y *= THIRD;
 83	}
 84	return xy;
 85}
 86FREEUP; if (P) pj_dalloc(P); }
 87ENTRY0(chamb)
 88	int i, j;
 89	char line[10];
 90
 91	for (i = 0; i < 3; ++i) { /* get control point locations */
 92		(void)sprintf(line, "rlat_%d", i+1);
 93		P->c[i].phi = pj_param(P->params, line).f;
 94		(void)sprintf(line, "rlon_%d", i+1);
 95		P->c[i].lam = pj_param(P->params, line).f;
 96		P->c[i].lam = adjlon(P->c[i].lam - P->lam0);
 97		P->c[i].cosphi = cos(P->c[i].phi);
 98		P->c[i].sinphi = sin(P->c[i].phi);
 99	}
100	for (i = 0; i < 3; ++i) { /* inter ctl pt. distances and azimuths */
101		j = i == 2 ? 0 : i + 1;
102		P->c[i].v = vect(P->c[j].phi - P->c[i].phi, P->c[i].cosphi, P->c[i].sinphi,
103			P->c[j].cosphi, P->c[j].sinphi, P->c[j].lam - P->c[i].lam);
104		if (! P->c[i].v.r) E_ERROR(-25);
105		/* co-linearity problem ignored for now */
106	}
107	P->beta_0 = lc(P->c[0].v.r, P->c[2].v.r, P->c[1].v.r);
108	P->beta_1 = lc(P->c[0].v.r, P->c[1].v.r, P->c[2].v.r);
109	P->beta_2 = PI - P->beta_0;
110	P->p.y = 2. * (P->c[0].p.y = P->c[1].p.y = P->c[2].v.r * sin(P->beta_0));
111	P->c[2].p.y = 0.;
112	P->c[0].p.x = - (P->c[1].p.x = 0.5 * P->c[0].v.r);
113	P->p.x = P->c[2].p.x = P->c[0].p.x + P->c[2].v.r * cos(P->beta_0);
114	P->es = 0.; P->fwd = s_forward;
115ENDENTRY(P)