PageRenderTime 208ms CodeModel.GetById 18ms RepoModel.GetById 1ms app.codeStats 0ms

/Proj4/geod_for.c

http://github.com/route-me/route-me
C | 107 lines | 106 code | 1 blank | 0 comment | 23 complexity | 92f137bd6e14fa48b50ed8ed04bc33f6 MD5 | raw file
  1. #ifndef lint
  2. static const char SCCSID[]="@(#)geod_for.c 4.6 95/09/23 GIE REL";
  3. #endif
  4. # include "projects.h"
  5. # include "geodesic.h"
  6. # define MERI_TOL 1e-9
  7. static double
  8. th1,costh1,sinth1,sina12,cosa12,M,N,c1,c2,D,P,s1;
  9. static int
  10. merid, signS;
  11. void
  12. geod_pre(void) {
  13. al12 = adjlon(al12); /* reduce to +- 0-PI */
  14. signS = fabs(al12) > HALFPI ? 1 : 0;
  15. th1 = ellipse ? atan(onef * tan(phi1)) : phi1;
  16. costh1 = cos(th1);
  17. sinth1 = sin(th1);
  18. if ((merid = fabs(sina12 = sin(al12)) < MERI_TOL)) {
  19. sina12 = 0.;
  20. cosa12 = fabs(al12) < HALFPI ? 1. : -1.;
  21. M = 0.;
  22. } else {
  23. cosa12 = cos(al12);
  24. M = costh1 * sina12;
  25. }
  26. N = costh1 * cosa12;
  27. if (ellipse) {
  28. if (merid) {
  29. c1 = 0.;
  30. c2 = f4;
  31. D = 1. - c2;
  32. D *= D;
  33. P = c2 / D;
  34. } else {
  35. c1 = geod_f * M;
  36. c2 = f4 * (1. - M * M);
  37. D = (1. - c2)*(1. - c2 - c1 * M);
  38. P = (1. + .5 * c1 * M) * c2 / D;
  39. }
  40. }
  41. if (merid) s1 = HALFPI - th1;
  42. else {
  43. s1 = (fabs(M) >= 1.) ? 0. : acos(M);
  44. s1 = sinth1 / sin(s1);
  45. s1 = (fabs(s1) >= 1.) ? 0. : acos(s1);
  46. }
  47. }
  48. void
  49. geod_for(void) {
  50. double d,sind,u,V,X,ds,cosds,sinds,ss,de;
  51. ss = 0.;
  52. if (ellipse) {
  53. d = geod_S / (D * geod_a);
  54. if (signS) d = -d;
  55. u = 2. * (s1 - d);
  56. V = cos(u + d);
  57. X = c2 * c2 * (sind = sin(d)) * cos(d) * (2. * V * V - 1.);
  58. ds = d + X - 2. * P * V * (1. - 2. * P * cos(u)) * sind;
  59. ss = s1 + s1 - ds;
  60. } else {
  61. ds = geod_S / geod_a;
  62. if (signS) ds = - ds;
  63. }
  64. cosds = cos(ds);
  65. sinds = sin(ds);
  66. if (signS) sinds = - sinds;
  67. al21 = N * cosds - sinth1 * sinds;
  68. if (merid) {
  69. phi2 = atan( tan(HALFPI + s1 - ds) / onef);
  70. if (al21 > 0.) {
  71. al21 = PI;
  72. if (signS)
  73. de = PI;
  74. else {
  75. phi2 = - phi2;
  76. de = 0.;
  77. }
  78. } else {
  79. al21 = 0.;
  80. if (signS) {
  81. phi2 = - phi2;
  82. de = 0;
  83. } else
  84. de = PI;
  85. }
  86. } else {
  87. al21 = atan(M / al21);
  88. if (al21 > 0)
  89. al21 += PI;
  90. if (al12 < 0.)
  91. al21 -= PI;
  92. al21 = adjlon(al21);
  93. phi2 = atan(-(sinth1 * cosds + N * sinds) * sin(al21) /
  94. (ellipse ? onef * M : M));
  95. de = atan2(sinds * sina12 ,
  96. (costh1 * cosds - sinth1 * sinds * cosa12));
  97. if (ellipse)
  98. if (signS)
  99. de += c1 * ((1. - c2) * ds +
  100. c2 * sinds * cos(ss));
  101. else
  102. de -= c1 * ((1. - c2) * ds -
  103. c2 * sinds * cos(ss));
  104. }
  105. lam2 = adjlon( lam1 + de );
  106. }