/Proj4/PJ_labrd.c

http://github.com/route-me/route-me · C · 112 lines · 104 code · 3 blank · 5 comment · 4 complexity · 9020484bd657bb2090dba4c2afc96a70 MD5 · raw file

  1. #ifndef lint
  2. static const char SCCSID[]="@(#)PJ_labrd.c 4.1 94/02/15 GIE REL";
  3. #endif
  4. #define PROJ_PARMS__ \
  5. double Az, kRg, p0s, A, C, Ca, Cb, Cc, Cd; \
  6. int rot;
  7. #define PJ_LIB__
  8. #include "projects.h"
  9. PROJ_HEAD(labrd, "Laborde") "\n\tCyl, Sph\n\tSpecial for Madagascar";
  10. #define EPS 1.e-10
  11. FORWARD(e_forward);
  12. double V1, V2, ps, sinps, cosps, sinps2, cosps2, I1, I2, I3, I4, I5, I6,
  13. x2, y2, t;
  14. V1 = P->A * log( tan(FORTPI + .5 * lp.phi) );
  15. t = P->e * sin(lp.phi);
  16. V2 = .5 * P->e * P->A * log ((1. + t)/(1. - t));
  17. ps = 2. * (atan(exp(V1 - V2 + P->C)) - FORTPI);
  18. I1 = ps - P->p0s;
  19. cosps = cos(ps); cosps2 = cosps * cosps;
  20. sinps = sin(ps); sinps2 = sinps * sinps;
  21. I4 = P->A * cosps;
  22. I2 = .5 * P->A * I4 * sinps;
  23. I3 = I2 * P->A * P->A * (5. * cosps2 - sinps2) / 12.;
  24. I6 = I4 * P->A * P->A;
  25. I5 = I6 * (cosps2 - sinps2) / 6.;
  26. I6 *= P->A * P->A *
  27. (5. * cosps2 * cosps2 + sinps2 * (sinps2 - 18. * cosps2)) / 120.;
  28. t = lp.lam * lp.lam;
  29. xy.x = P->kRg * lp.lam * (I4 + t * (I5 + t * I6));
  30. xy.y = P->kRg * (I1 + t * (I2 + t * I3));
  31. x2 = xy.x * xy.x;
  32. y2 = xy.y * xy.y;
  33. V1 = 3. * xy.x * y2 - xy.x * x2;
  34. V2 = xy.y * y2 - 3. * x2 * xy.y;
  35. xy.x += P->Ca * V1 + P->Cb * V2;
  36. xy.y += P->Ca * V2 - P->Cb * V1;
  37. return (xy);
  38. }
  39. INVERSE(e_inverse); /* ellipsoid & spheroid */
  40. double x2, y2, V1, V2, V3, V4, t, t2, ps, pe, tpe, s,
  41. I7, I8, I9, I10, I11, d, Re;
  42. int i;
  43. x2 = xy.x * xy.x;
  44. y2 = xy.y * xy.y;
  45. V1 = 3. * xy.x * y2 - xy.x * x2;
  46. V2 = xy.y * y2 - 3. * x2 * xy.y;
  47. V3 = xy.x * (5. * y2 * y2 + x2 * (-10. * y2 + x2 ));
  48. V4 = xy.y * (5. * x2 * x2 + y2 * (-10. * x2 + y2 ));
  49. xy.x += - P->Ca * V1 - P->Cb * V2 + P->Cc * V3 + P->Cd * V4;
  50. xy.y += P->Cb * V1 - P->Ca * V2 - P->Cd * V3 + P->Cc * V4;
  51. ps = P->p0s + xy.y / P->kRg;
  52. pe = ps + P->phi0 - P->p0s;
  53. for ( i = 20; i; --i) {
  54. V1 = P->A * log(tan(FORTPI + .5 * pe));
  55. tpe = P->e * sin(pe);
  56. V2 = .5 * P->e * P->A * log((1. + tpe)/(1. - tpe));
  57. t = ps - 2. * (atan(exp(V1 - V2 + P->C)) - FORTPI);
  58. pe += t;
  59. if (fabs(t) < EPS)
  60. break;
  61. }
  62. /*
  63. if (!i) {
  64. } else {
  65. }
  66. */
  67. t = P->e * sin(pe);
  68. t = 1. - t * t;
  69. Re = P->one_es / ( t * sqrt(t) );
  70. t = tan(ps);
  71. t2 = t * t;
  72. s = P->kRg * P->kRg;
  73. d = Re * P->k0 * P->kRg;
  74. I7 = t / (2. * d);
  75. I8 = t * (5. + 3. * t2) / (24. * d * s);
  76. d = cos(ps) * P->kRg * P->A;
  77. I9 = 1. / d;
  78. d *= s;
  79. I10 = (1. + 2. * t2) / (6. * d);
  80. I11 = (5. + t2 * (28. + 24. * t2)) / (120. * d * s);
  81. x2 = xy.x * xy.x;
  82. lp.phi = pe + x2 * (-I7 + I8 * x2);
  83. lp.lam = xy.x * (I9 + x2 * (-I10 + x2 * I11));
  84. return (lp);
  85. }
  86. FREEUP; if (P) pj_dalloc(P); }
  87. ENTRY0(labrd)
  88. double Az, sinp, R, N, t;
  89. P->rot = pj_param(P->params, "bno_rot").i == 0;
  90. Az = pj_param(P->params, "razi").f;
  91. sinp = sin(P->phi0);
  92. t = 1. - P->es * sinp * sinp;
  93. N = 1. / sqrt(t);
  94. R = P->one_es * N / t;
  95. P->kRg = P->k0 * sqrt( N * R );
  96. P->p0s = atan( sqrt(R / N) * tan(P->phi0) );
  97. P->A = sinp / sin(P->p0s);
  98. t = P->e * sinp;
  99. P->C = .5 * P->e * P->A * log((1. + t)/(1. - t)) +
  100. - P->A * log( tan(FORTPI + .5 * P->phi0))
  101. + log( tan(FORTPI + .5 * P->p0s));
  102. t = Az + Az;
  103. P->Ca = (1. - cos(t)) * ( P->Cb = 1. / (12. * P->kRg * P->kRg) );
  104. P->Cb *= sin(t);
  105. P->Cc = 3. * (P->Ca * P->Ca - P->Cb * P->Cb);
  106. P->Cd = 6. * P->Ca * P->Cb;
  107. P->inv = e_inverse;
  108. P->fwd = e_forward;
  109. ENDENTRY(P)