/Proj4/PJ_chamb.c

http://github.com/route-me/route-me · C · 115 lines · 111 code · 3 blank · 1 comment · 21 complexity · 6e8164807f3919f4c3becd986ade064e MD5 · raw file

  1. #ifndef lint
  2. static const char SCCSID[]="@(#)PJ_chamb.c 4.1 94/02/15 GIE REL";
  3. #endif
  4. typedef struct { double r, Az; } VECT;
  5. #define PROJ_PARMS__ \
  6. struct { /* control point data */ \
  7. double phi, lam; \
  8. double cosphi, sinphi; \
  9. VECT v; \
  10. XY p; \
  11. double Az; \
  12. } c[3]; \
  13. XY p; \
  14. double beta_0, beta_1, beta_2;
  15. #define PJ_LIB__
  16. #include "projects.h"
  17. PROJ_HEAD(chamb, "Chamberlin Trimetric") "\n\tMisc Sph, no inv."
  18. "\n\tlat_1= lon_1= lat_2= lon_2= lat_3= lon_3=";
  19. #include <stdio.h>
  20. #define THIRD 0.333333333333333333
  21. #define TOL 1e-9
  22. static VECT /* distance and azimuth from point 1 to point 2 */
  23. vect(double dphi, double c1, double s1, double c2, double s2, double dlam) {
  24. VECT v;
  25. double cdl, dp, dl;
  26. cdl = cos(dlam);
  27. if (fabs(dphi) > 1. || fabs(dlam) > 1.)
  28. v.r = aacos(s1 * s2 + c1 * c2 * cdl);
  29. else { /* more accurate for smaller distances */
  30. dp = sin(.5 * dphi);
  31. dl = sin(.5 * dlam);
  32. v.r = 2. * aasin(sqrt(dp * dp + c1 * c2 * dl * dl));
  33. }
  34. if (fabs(v.r) > TOL)
  35. v.Az = atan2(c2 * sin(dlam), c1 * s2 - s1 * c2 * cdl);
  36. else
  37. v.r = v.Az = 0.;
  38. return v;
  39. }
  40. static double /* law of cosines */
  41. lc(double b,double c,double a) {
  42. return aacos(.5 * (b * b + c * c - a * a) / (b * c));
  43. }
  44. FORWARD(s_forward); /* spheroid */
  45. double sinphi, cosphi, a;
  46. VECT v[3];
  47. int i, j;
  48. sinphi = sin(lp.phi);
  49. cosphi = cos(lp.phi);
  50. for (i = 0; i < 3; ++i) { /* dist/azimiths from control */
  51. v[i] = vect(lp.phi - P->c[i].phi, P->c[i].cosphi, P->c[i].sinphi,
  52. cosphi, sinphi, lp.lam - P->c[i].lam);
  53. if ( ! v[i].r)
  54. break;
  55. v[i].Az = adjlon(v[i].Az - P->c[i].v.Az);
  56. }
  57. if (i < 3) /* current point at control point */
  58. xy = P->c[i].p;
  59. else { /* point mean of intersepts */
  60. xy = P->p;
  61. for (i = 0; i < 3; ++i) {
  62. j = i == 2 ? 0 : i + 1;
  63. a = lc(P->c[i].v.r, v[i].r, v[j].r);
  64. if (v[i].Az < 0.)
  65. a = -a;
  66. if (! i) { /* coord comp unique to each arc */
  67. xy.x += v[i].r * cos(a);
  68. xy.y -= v[i].r * sin(a);
  69. } else if (i == 1) {
  70. a = P->beta_1 - a;
  71. xy.x -= v[i].r * cos(a);
  72. xy.y -= v[i].r * sin(a);
  73. } else {
  74. a = P->beta_2 - a;
  75. xy.x += v[i].r * cos(a);
  76. xy.y += v[i].r * sin(a);
  77. }
  78. }
  79. xy.x *= THIRD; /* mean of arc intercepts */
  80. xy.y *= THIRD;
  81. }
  82. return xy;
  83. }
  84. FREEUP; if (P) pj_dalloc(P); }
  85. ENTRY0(chamb)
  86. int i, j;
  87. char line[10];
  88. for (i = 0; i < 3; ++i) { /* get control point locations */
  89. (void)sprintf(line, "rlat_%d", i+1);
  90. P->c[i].phi = pj_param(P->params, line).f;
  91. (void)sprintf(line, "rlon_%d", i+1);
  92. P->c[i].lam = pj_param(P->params, line).f;
  93. P->c[i].lam = adjlon(P->c[i].lam - P->lam0);
  94. P->c[i].cosphi = cos(P->c[i].phi);
  95. P->c[i].sinphi = sin(P->c[i].phi);
  96. }
  97. for (i = 0; i < 3; ++i) { /* inter ctl pt. distances and azimuths */
  98. j = i == 2 ? 0 : i + 1;
  99. P->c[i].v = vect(P->c[j].phi - P->c[i].phi, P->c[i].cosphi, P->c[i].sinphi,
  100. P->c[j].cosphi, P->c[j].sinphi, P->c[j].lam - P->c[i].lam);
  101. if (! P->c[i].v.r) E_ERROR(-25);
  102. /* co-linearity problem ignored for now */
  103. }
  104. P->beta_0 = lc(P->c[0].v.r, P->c[2].v.r, P->c[1].v.r);
  105. P->beta_1 = lc(P->c[0].v.r, P->c[1].v.r, P->c[2].v.r);
  106. P->beta_2 = PI - P->beta_0;
  107. P->p.y = 2. * (P->c[0].p.y = P->c[1].p.y = P->c[2].v.r * sin(P->beta_0));
  108. P->c[2].p.y = 0.;
  109. P->c[0].p.x = - (P->c[1].p.x = 0.5 * P->c[0].v.r);
  110. P->p.x = P->c[2].p.x = P->c[0].p.x + P->c[2].v.r * cos(P->beta_0);
  111. P->es = 0.; P->fwd = s_forward;
  112. ENDENTRY(P)