/Proj4/PJ_imw_p.c

http://github.com/route-me/route-me · C · 154 lines · 148 code · 6 blank · 0 comment · 29 complexity · 2931594e2aeda4538c48b44525d88e5c MD5 · raw file

  1. #ifndef lint
  2. static const char SCCSID[]="@(#)PJ_imw_p.c 4.1 94/05/22 GIE REL";
  3. #endif
  4. #define PROJ_PARMS__ \
  5. double P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2; \
  6. double phi_1, phi_2, lam_1; \
  7. double *en; \
  8. int mode; /* = 0, phi_1 and phi_2 != 0, = 1, phi_1 = 0, = -1 phi_2 = 0 */
  9. #define PJ_LIB__
  10. #include "projects.h"
  11. PROJ_HEAD(imw_p, "International Map of the World Polyconic")
  12. "\n\tMod. Polyconic, Ell\n\tlat_1= and lat_2= [lon_1=]";
  13. #define TOL 1e-10
  14. #define EPS 1e-10
  15. static int
  16. phi12(PJ *P, double *del, double *sig) {
  17. int err = 0;
  18. if (!pj_param(P->params, "tlat_1").i ||
  19. !pj_param(P->params, "tlat_2").i) {
  20. err = -41;
  21. } else {
  22. P->phi_1 = pj_param(P->params, "rlat_1").f;
  23. P->phi_2 = pj_param(P->params, "rlat_2").f;
  24. *del = 0.5 * (P->phi_2 - P->phi_1);
  25. *sig = 0.5 * (P->phi_2 + P->phi_1);
  26. err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? -42 : 0;
  27. }
  28. return err;
  29. }
  30. static XY
  31. loc_for(LP lp, PJ *P, double *yc) {
  32. XY xy;
  33. if (! lp.phi) {
  34. xy.x = lp.lam;
  35. xy.y = 0.;
  36. } else {
  37. double xa, ya, xb, yb, xc, yc, D, B, m, sp, t, R, C;
  38. sp = sin(lp.phi);
  39. m = pj_mlfn(lp.phi, sp, cos(lp.phi), P->en);
  40. xa = P->Pp + P->Qp * m;
  41. ya = P->P + P->Q * m;
  42. R = 1. / (tan(lp.phi) * sqrt(1. - P->es * sp * sp));
  43. C = sqrt(R * R - xa * xa);
  44. if (lp.phi < 0.) C = - C;
  45. C += ya - R;
  46. if (P->mode < 0) {
  47. xb = lp.lam;
  48. yb = P->C2;
  49. } else {
  50. t = lp.lam * P->sphi_2;
  51. xb = P->R_2 * sin(t);
  52. yb = P->C2 + P->R_2 * (1. - cos(t));
  53. }
  54. if (P->mode > 0) {
  55. xc = lp.lam;
  56. yc = 0.;
  57. } else {
  58. t = lp.lam * P->sphi_1;
  59. xc = P->R_1 * sin(t);
  60. yc = P->R_1 * (1. - cos(t));
  61. }
  62. D = (xb - xc)/(yb - yc);
  63. B = xc + D * (C + R - yc);
  64. xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
  65. if (lp.phi > 0)
  66. xy.x = - xy.x;
  67. xy.x = (B + xy.x) / (1. + D * D);
  68. xy.y = sqrt(R * R - xy.x * xy.x);
  69. if (lp.phi > 0)
  70. xy.y = - xy.y;
  71. xy.y += C + R;
  72. }
  73. return (xy);
  74. }
  75. FORWARD(e_forward); /* ellipsoid */
  76. double yc;
  77. xy = loc_for(lp, P, &yc);
  78. return (xy);
  79. }
  80. INVERSE(e_inverse); /* ellipsoid */
  81. XY t;
  82. double yc;
  83. lp.phi = P->phi_2;
  84. lp.lam = xy.x / cos(lp.phi);
  85. do {
  86. t = loc_for(lp, P, &yc);
  87. lp.phi = ((lp.phi - P->phi_1) * (xy.y - yc) / (t.y - yc)) + P->phi_1;
  88. lp.lam = lp.lam * xy.x / t.x;
  89. } while (fabs(t.x - xy.x) > TOL || fabs(t.y - xy.y) > TOL);
  90. return (lp);
  91. }
  92. static void
  93. xy(PJ *P, double phi, double *x, double *y, double *sp, double *R) {
  94. double F;
  95. *sp = sin(phi);
  96. *R = 1./(tan(phi) * sqrt(1. - P->es * *sp * *sp ));
  97. F = P->lam_1 * *sp;
  98. *y = *R * (1 - cos(F));
  99. *x = *R * sin(F);
  100. }
  101. FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
  102. ENTRY1(imw_p, en)
  103. double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
  104. int i;
  105. if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
  106. if( (i = phi12(P, &del, &sig)) != 0)
  107. E_ERROR(i);
  108. if (P->phi_2 < P->phi_1) { /* make sure P->phi_1 most southerly */
  109. del = P->phi_1;
  110. P->phi_1 = P->phi_2;
  111. P->phi_2 = del;
  112. }
  113. if (pj_param(P->params, "tlon_1").i)
  114. P->lam_1 = pj_param(P->params, "rlon_1").f;
  115. else { /* use predefined based upon latitude */
  116. sig = fabs(sig * RAD_TO_DEG);
  117. if (sig <= 60) sig = 2.;
  118. else if (sig <= 76) sig = 4.;
  119. else sig = 8.;
  120. P->lam_1 = sig * DEG_TO_RAD;
  121. }
  122. P->mode = 0;
  123. if (P->phi_1) xy(P, P->phi_1, &x1, &y1, &P->sphi_1, &P->R_1);
  124. else {
  125. P->mode = 1;
  126. y1 = 0.;
  127. x1 = P->lam_1;
  128. }
  129. if (P->phi_2) xy(P, P->phi_2, &x2, &T2, &P->sphi_2, &P->R_2);
  130. else {
  131. P->mode = -1;
  132. T2 = 0.;
  133. x2 = P->lam_1;
  134. }
  135. m1 = pj_mlfn(P->phi_1, P->sphi_1, cos(P->phi_1), P->en);
  136. m2 = pj_mlfn(P->phi_2, P->sphi_2, cos(P->phi_2), P->en);
  137. t = m2 - m1;
  138. s = x2 - x1;
  139. y2 = sqrt(t * t - s * s) + y1;
  140. P->C2 = y2 - T2;
  141. t = 1. / t;
  142. P->P = (m2 * y1 - m1 * y2) * t;
  143. P->Q = (y2 - y1) * t;
  144. P->Pp = (m2 * x1 - m1 * x2) * t;
  145. P->Qp = (x2 - x1) * t;
  146. P->fwd = e_forward;
  147. P->inv = e_inverse;
  148. ENDENTRY(P)