/Proj4/PJ_bipc.c
C | 135 lines | 133 code | 2 blank | 0 comment | 28 complexity | 6130e12f99d4cc039c55d1d2f4e6bdc3 MD5 | raw file
- #ifndef lint
- static const char SCCSID[]="@(#)PJ_bipc.c 4.1 94/02/15 GIE REL";
- #endif
- #define PROJ_PARMS__ \
- int noskew;
- #define PJ_LIB__
- # include "projects.h"
- PROJ_HEAD(bipc, "Bipolar conic of western hemisphere")
- "\n\tConic Sph.";
- # define EPS 1e-10
- # define EPS10 1e-10
- # define ONEEPS 1.000000001
- # define NITER 10
- # define lamB -.34894976726250681539
- # define n .63055844881274687180
- # define F 1.89724742567461030582
- # define Azab .81650043674686363166
- # define Azba 1.82261843856185925133
- # define T 1.27246578267089012270
- # define rhoc 1.20709121521568721927
- # define cAzc .69691523038678375519
- # define sAzc .71715351331143607555
- # define C45 .70710678118654752469
- # define S45 .70710678118654752410
- # define C20 .93969262078590838411
- # define S20 -.34202014332566873287
- # define R110 1.91986217719376253360
- # define R104 1.81514242207410275904
- FORWARD(s_forward); /* spheroid */
- double cphi, sphi, tphi, t, al, Az, z, Av, cdlam, sdlam, r;
- int tag;
- cphi = cos(lp.phi);
- sphi = sin(lp.phi);
- cdlam = cos(sdlam = lamB - lp.lam);
- sdlam = sin(sdlam);
- if (fabs(fabs(lp.phi) - HALFPI) < EPS10) {
- Az = lp.phi < 0. ? PI : 0.;
- tphi = HUGE_VAL;
- } else {
- tphi = sphi / cphi;
- Az = atan2(sdlam , C45 * (tphi - cdlam));
- }
- if( (tag = (Az > Azba)) ) {
- cdlam = cos(sdlam = lp.lam + R110);
- sdlam = sin(sdlam);
- z = S20 * sphi + C20 * cphi * cdlam;
- if (fabs(z) > 1.) {
- if (fabs(z) > ONEEPS) F_ERROR
- else z = z < 0. ? -1. : 1.;
- } else
- z = acos(z);
- if (tphi != HUGE_VAL)
- Az = atan2(sdlam, (C20 * tphi - S20 * cdlam));
- Av = Azab;
- xy.y = rhoc;
- } else {
- z = S45 * (sphi + cphi * cdlam);
- if (fabs(z) > 1.) {
- if (fabs(z) > ONEEPS) F_ERROR
- else z = z < 0. ? -1. : 1.;
- } else
- z = acos(z);
- Av = Azba;
- xy.y = -rhoc;
- }
- if (z < 0.) F_ERROR;
- r = F * (t = pow(tan(.5 * z), n));
- if ((al = .5 * (R104 - z)) < 0.) F_ERROR;
- al = (t + pow(al, n)) / T;
- if (fabs(al) > 1.) {
- if (fabs(al) > ONEEPS) F_ERROR
- else al = al < 0. ? -1. : 1.;
- } else
- al = acos(al);
- if (fabs(t = n * (Av - Az)) < al)
- r /= cos(al + (tag ? t : -t));
- xy.x = r * sin(t);
- xy.y += (tag ? -r : r) * cos(t);
- if (P->noskew) {
- t = xy.x;
- xy.x = -xy.x * cAzc - xy.y * sAzc;
- xy.y = -xy.y * cAzc + t * sAzc;
- }
- return (xy);
- }
- INVERSE(s_inverse); /* spheroid */
- double t, r, rp, rl, al, z, fAz, Az, s, c, Av;
- int neg, i;
- if (P->noskew) {
- t = xy.x;
- xy.x = -xy.x * cAzc + xy.y * sAzc;
- xy.y = -xy.y * cAzc - t * sAzc;
- }
- if( (neg = (xy.x < 0.)) ) {
- xy.y = rhoc - xy.y;
- s = S20;
- c = C20;
- Av = Azab;
- } else {
- xy.y += rhoc;
- s = S45;
- c = C45;
- Av = Azba;
- }
- rl = rp = r = hypot(xy.x, xy.y);
- fAz = fabs(Az = atan2(xy.x, xy.y));
- for (i = NITER; i ; --i) {
- z = 2. * atan(pow(r / F,1 / n));
- al = acos((pow(tan(.5 * z), n) +
- pow(tan(.5 * (R104 - z)), n)) / T);
- if (fAz < al)
- r = rp * cos(al + (neg ? Az : -Az));
- if (fabs(rl - r) < EPS)
- break;
- rl = r;
- }
- if (! i) I_ERROR;
- Az = Av - Az / n;
- lp.phi = asin(s * cos(z) + c * sin(z) * cos(Az));
- lp.lam = atan2(sin(Az), c / tan(z) - s * cos(Az));
- if (neg)
- lp.lam -= R110;
- else
- lp.lam = lamB - lp.lam;
- return (lp);
- }
- FREEUP; if (P) pj_dalloc(P); }
- ENTRY0(bipc)
- P->noskew = pj_param(P->params, "bns").i;
- P->inv = s_inverse;
- P->fwd = s_forward;
- P->es = 0.;
- ENDENTRY(P)