PageRenderTime 17ms CodeModel.GetById 0ms RepoModel.GetById 0ms app.codeStats 0ms

/Proj4/PJ_cass.c

http://github.com/route-me/route-me
C | 82 lines | 81 code | 1 blank | 0 comment | 5 complexity | ae3dae96ab70eb78789fe671ff78d8dd MD5 | raw file
  1. #ifndef lint
  2. static const char SCCSID[]="@(#)PJ_cass.c 4.1 94/02/15 GIE REL";
  3. #endif
  4. #define PROJ_PARMS__ \
  5. double m0; \
  6. double n; \
  7. double t; \
  8. double a1; \
  9. double c; \
  10. double r; \
  11. double dd; \
  12. double d2; \
  13. double a2; \
  14. double tn; \
  15. double *en;
  16. #define PJ_LIB__
  17. # include "projects.h"
  18. PROJ_HEAD(cass, "Cassini") "\n\tCyl, Sph&Ell";
  19. # define EPS10 1e-10
  20. # define C1 .16666666666666666666
  21. # define C2 .00833333333333333333
  22. # define C3 .04166666666666666666
  23. # define C4 .33333333333333333333
  24. # define C5 .06666666666666666666
  25. FORWARD(e_forward); /* ellipsoid */
  26. xy.y = pj_mlfn(lp.phi, P->n = sin(lp.phi), P->c = cos(lp.phi), P->en);
  27. P->n = 1./sqrt(1. - P->es * P->n * P->n);
  28. P->tn = tan(lp.phi); P->t = P->tn * P->tn;
  29. P->a1 = lp.lam * P->c;
  30. P->c *= P->es * P->c / (1 - P->es);
  31. P->a2 = P->a1 * P->a1;
  32. xy.x = P->n * P->a1 * (1. - P->a2 * P->t *
  33. (C1 - (8. - P->t + 8. * P->c) * P->a2 * C2));
  34. xy.y -= P->m0 - P->n * P->tn * P->a2 *
  35. (.5 + (5. - P->t + 6. * P->c) * P->a2 * C3);
  36. return (xy);
  37. }
  38. FORWARD(s_forward); /* spheroid */
  39. xy.x = asin(cos(lp.phi) * sin(lp.lam));
  40. xy.y = atan2(tan(lp.phi) , cos(lp.lam)) - P->phi0;
  41. return (xy);
  42. }
  43. INVERSE(e_inverse); /* ellipsoid */
  44. double ph1;
  45. ph1 = pj_inv_mlfn(P->m0 + xy.y, P->es, P->en);
  46. P->tn = tan(ph1); P->t = P->tn * P->tn;
  47. P->n = sin(ph1);
  48. P->r = 1. / (1. - P->es * P->n * P->n);
  49. P->n = sqrt(P->r);
  50. P->r *= (1. - P->es) * P->n;
  51. P->dd = xy.x / P->n;
  52. P->d2 = P->dd * P->dd;
  53. lp.phi = ph1 - (P->n * P->tn / P->r) * P->d2 *
  54. (.5 - (1. + 3. * P->t) * P->d2 * C3);
  55. lp.lam = P->dd * (1. + P->t * P->d2 *
  56. (-C4 + (1. + 3. * P->t) * P->d2 * C5)) / cos(ph1);
  57. return (lp);
  58. }
  59. INVERSE(s_inverse); /* spheroid */
  60. lp.phi = asin(sin(P->dd = xy.y + P->phi0) * cos(xy.x));
  61. lp.lam = atan2(tan(xy.x), cos(P->dd));
  62. return (lp);
  63. }
  64. FREEUP;
  65. if (P) {
  66. if (P->en)
  67. pj_dalloc(P->en);
  68. pj_dalloc(P);
  69. }
  70. }
  71. ENTRY1(cass, en)
  72. if (P->es) {
  73. if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
  74. P->m0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), P->en);
  75. P->inv = e_inverse;
  76. P->fwd = e_forward;
  77. } else {
  78. P->inv = s_inverse;
  79. P->fwd = s_forward;
  80. }
  81. ENDENTRY(P)