1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
|
#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)
|