/Proj4/PJ_bipc.c
C | 135 lines | 133 code | 2 blank | 0 comment | 28 complexity | 6130e12f99d4cc039c55d1d2f4e6bdc3 MD5 | raw file
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)