PageRenderTime 43ms CodeModel.GetById 25ms RepoModel.GetById 0ms app.codeStats 0ms

/Proj4/PJ_bipc.c

http://github.com/route-me/route-me
C | 135 lines | 133 code | 2 blank | 0 comment | 28 complexity | 6130e12f99d4cc039c55d1d2f4e6bdc3 MD5 | raw file
  1. #ifndef lint
  2. static 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"
  8. PROJ_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
  29. FORWARD(s_forward); /* spheroid */
  30. double cphi, sphi, tphi, t, al, Az, z, Av, cdlam, sdlam, r;
  31. int tag;
  32. cphi = cos(lp.phi);
  33. sphi = sin(lp.phi);
  34. cdlam = cos(sdlam = lamB - lp.lam);
  35. sdlam = sin(sdlam);
  36. if (fabs(fabs(lp.phi) - HALFPI) < EPS10) {
  37. Az = lp.phi < 0. ? PI : 0.;
  38. tphi = HUGE_VAL;
  39. } else {
  40. tphi = sphi / cphi;
  41. Az = atan2(sdlam , C45 * (tphi - cdlam));
  42. }
  43. if( (tag = (Az > Azba)) ) {
  44. cdlam = cos(sdlam = lp.lam + R110);
  45. sdlam = sin(sdlam);
  46. z = S20 * sphi + C20 * cphi * cdlam;
  47. if (fabs(z) > 1.) {
  48. if (fabs(z) > ONEEPS) F_ERROR
  49. else z = z < 0. ? -1. : 1.;
  50. } else
  51. z = acos(z);
  52. if (tphi != HUGE_VAL)
  53. Az = atan2(sdlam, (C20 * tphi - S20 * cdlam));
  54. Av = Azab;
  55. xy.y = rhoc;
  56. } else {
  57. z = S45 * (sphi + cphi * cdlam);
  58. if (fabs(z) > 1.) {
  59. if (fabs(z) > ONEEPS) F_ERROR
  60. else z = z < 0. ? -1. : 1.;
  61. } else
  62. z = acos(z);
  63. Av = Azba;
  64. xy.y = -rhoc;
  65. }
  66. if (z < 0.) F_ERROR;
  67. r = F * (t = pow(tan(.5 * z), n));
  68. if ((al = .5 * (R104 - z)) < 0.) F_ERROR;
  69. al = (t + pow(al, n)) / T;
  70. if (fabs(al) > 1.) {
  71. if (fabs(al) > ONEEPS) F_ERROR
  72. else al = al < 0. ? -1. : 1.;
  73. } else
  74. al = acos(al);
  75. if (fabs(t = n * (Av - Az)) < al)
  76. r /= cos(al + (tag ? t : -t));
  77. xy.x = r * sin(t);
  78. xy.y += (tag ? -r : r) * cos(t);
  79. if (P->noskew) {
  80. t = xy.x;
  81. xy.x = -xy.x * cAzc - xy.y * sAzc;
  82. xy.y = -xy.y * cAzc + t * sAzc;
  83. }
  84. return (xy);
  85. }
  86. INVERSE(s_inverse); /* spheroid */
  87. double t, r, rp, rl, al, z, fAz, Az, s, c, Av;
  88. int neg, i;
  89. if (P->noskew) {
  90. t = xy.x;
  91. xy.x = -xy.x * cAzc + xy.y * sAzc;
  92. xy.y = -xy.y * cAzc - t * sAzc;
  93. }
  94. if( (neg = (xy.x < 0.)) ) {
  95. xy.y = rhoc - xy.y;
  96. s = S20;
  97. c = C20;
  98. Av = Azab;
  99. } else {
  100. xy.y += rhoc;
  101. s = S45;
  102. c = C45;
  103. Av = Azba;
  104. }
  105. rl = rp = r = hypot(xy.x, xy.y);
  106. fAz = fabs(Az = atan2(xy.x, xy.y));
  107. for (i = NITER; i ; --i) {
  108. z = 2. * atan(pow(r / F,1 / n));
  109. al = acos((pow(tan(.5 * z), n) +
  110. pow(tan(.5 * (R104 - z)), n)) / T);
  111. if (fAz < al)
  112. r = rp * cos(al + (neg ? Az : -Az));
  113. if (fabs(rl - r) < EPS)
  114. break;
  115. rl = r;
  116. }
  117. if (! i) I_ERROR;
  118. Az = Av - Az / n;
  119. lp.phi = asin(s * cos(z) + c * sin(z) * cos(Az));
  120. lp.lam = atan2(sin(Az), c / tan(z) - s * cos(Az));
  121. if (neg)
  122. lp.lam -= R110;
  123. else
  124. lp.lam = lamB - lp.lam;
  125. return (lp);
  126. }
  127. FREEUP; if (P) pj_dalloc(P); }
  128. ENTRY0(bipc)
  129. P->noskew = pj_param(P->params, "bns").i;
  130. P->inv = s_inverse;
  131. P->fwd = s_forward;
  132. P->es = 0.;
  133. ENDENTRY(P)