/Proj4/PJ_bipc.c

  1#ifndef lint
  2static const char SCCSID[]="@(#)PJ_bipc.c	4.1	94/02/15	GIE	REL";
  3#endif
  4#define PROJ_PARMS__ \
  5	int	noskew;
  6#define PJ_LIB__
  7# include	"projects.h"
  8PROJ_HEAD(bipc, "Bipolar conic of western hemisphere")
  9	"\n\tConic Sph.";
 10# define EPS	1e-10
 11# define EPS10	1e-10
 12# define ONEEPS 1.000000001
 13# define NITER	10
 14# define lamB	-.34894976726250681539
 15# define n	.63055844881274687180
 16# define F	1.89724742567461030582
 17# define Azab	.81650043674686363166
 18# define Azba	1.82261843856185925133
 19# define T	1.27246578267089012270
 20# define rhoc	1.20709121521568721927
 21# define cAzc	.69691523038678375519
 22# define sAzc	.71715351331143607555
 23# define C45	.70710678118654752469
 24# define S45	.70710678118654752410
 25# define C20	.93969262078590838411
 26# define S20	-.34202014332566873287
 27# define R110	1.91986217719376253360
 28# define R104	1.81514242207410275904
 29FORWARD(s_forward); /* spheroid */
 30	double cphi, sphi, tphi, t, al, Az, z, Av, cdlam, sdlam, r;
 31	int tag;
 32
 33	cphi = cos(lp.phi);
 34	sphi = sin(lp.phi);
 35	cdlam = cos(sdlam = lamB - lp.lam);
 36	sdlam = sin(sdlam);
 37	if (fabs(fabs(lp.phi) - HALFPI) < EPS10) {
 38		Az = lp.phi < 0. ? PI : 0.;
 39		tphi = HUGE_VAL;
 40	} else {
 41		tphi = sphi / cphi;
 42		Az = atan2(sdlam , C45 * (tphi - cdlam));
 43	}
 44	if( (tag = (Az > Azba)) ) {
 45		cdlam = cos(sdlam = lp.lam + R110);
 46		sdlam = sin(sdlam);
 47		z = S20 * sphi + C20 * cphi * cdlam;
 48		if (fabs(z) > 1.) {
 49			if (fabs(z) > ONEEPS) F_ERROR
 50			else z = z < 0. ? -1. : 1.;
 51		} else
 52			z = acos(z);
 53		if (tphi != HUGE_VAL)
 54			Az = atan2(sdlam, (C20 * tphi - S20 * cdlam));
 55		Av = Azab;
 56		xy.y = rhoc;
 57	} else {
 58		z = S45 * (sphi + cphi * cdlam);
 59		if (fabs(z) > 1.) {
 60			if (fabs(z) > ONEEPS) F_ERROR
 61			else z = z < 0. ? -1. : 1.;
 62		} else
 63			z = acos(z);
 64		Av = Azba;
 65		xy.y = -rhoc;
 66	}
 67	if (z < 0.) F_ERROR;
 68	r = F * (t = pow(tan(.5 * z), n));
 69	if ((al = .5 * (R104 - z)) < 0.) F_ERROR;
 70	al = (t + pow(al, n)) / T;
 71	if (fabs(al) > 1.) {
 72		if (fabs(al) > ONEEPS) F_ERROR
 73		else al = al < 0. ? -1. : 1.;
 74	} else
 75		al = acos(al);
 76	if (fabs(t = n * (Av - Az)) < al)
 77		r /= cos(al + (tag ? t : -t));
 78	xy.x = r * sin(t);
 79	xy.y += (tag ? -r : r) * cos(t);
 80	if (P->noskew) {
 81		t = xy.x;
 82		xy.x = -xy.x * cAzc - xy.y * sAzc; 
 83		xy.y = -xy.y * cAzc + t * sAzc; 
 84	}
 85	return (xy);
 86}
 87INVERSE(s_inverse); /* spheroid */
 88	double t, r, rp, rl, al, z, fAz, Az, s, c, Av;
 89	int neg, i;
 90
 91	if (P->noskew) {
 92		t = xy.x;
 93		xy.x = -xy.x * cAzc + xy.y * sAzc; 
 94		xy.y = -xy.y * cAzc - t * sAzc; 
 95	}
 96	if( (neg = (xy.x < 0.)) ) {
 97		xy.y = rhoc - xy.y;
 98		s = S20;
 99		c = C20;
100		Av = Azab;
101	} else {
102		xy.y += rhoc;
103		s = S45;
104		c = C45;
105		Av = Azba;
106	}
107	rl = rp = r = hypot(xy.x, xy.y);
108	fAz = fabs(Az = atan2(xy.x, xy.y));
109	for (i = NITER; i ; --i) {
110		z = 2. * atan(pow(r / F,1 / n));
111		al = acos((pow(tan(.5 * z), n) +
112		   pow(tan(.5 * (R104 - z)), n)) / T);
113		if (fAz < al)
114			r = rp * cos(al + (neg ? Az : -Az));
115		if (fabs(rl - r) < EPS)
116			break;
117		rl = r;
118	}
119	if (! i) I_ERROR;
120	Az = Av - Az / n;
121	lp.phi = asin(s * cos(z) + c * sin(z) * cos(Az));
122	lp.lam = atan2(sin(Az), c / tan(z) - s * cos(Az));
123	if (neg)
124		lp.lam -= R110;
125	else
126		lp.lam = lamB - lp.lam;
127	return (lp);
128}
129FREEUP; if (P) pj_dalloc(P); }
130ENTRY0(bipc)
131	P->noskew = pj_param(P->params, "bns").i;
132	P->inv = s_inverse;
133	P->fwd = s_forward;
134	P->es = 0.;
135ENDENTRY(P)