/Proj4/PJ_eck4.c

http://github.com/route-me/route-me · C · 47 lines · 45 code · 2 blank · 0 comment · 5 complexity · a6d24e3fbfe0f994db7dff15595f2caf MD5 · raw file

  1. #ifndef lint
  2. static const char SCCSID[]="@(#)PJ_eck4.c 4.1 94/02/15 GIE REL";
  3. #endif
  4. #define PJ_LIB__
  5. #include "projects.h"
  6. PROJ_HEAD(eck4, "Eckert IV") "\n\tPCyl, Sph.";
  7. #define C_x .42223820031577120149
  8. #define C_y 1.32650042817700232218
  9. #define RC_y .75386330736002178205
  10. #define C_p 3.57079632679489661922
  11. #define RC_p .28004957675577868795
  12. #define EPS 1e-7
  13. #define NITER 6
  14. FORWARD(s_forward); /* spheroid */
  15. double p, V, s, c;
  16. int i;
  17. p = C_p * sin(lp.phi);
  18. V = lp.phi * lp.phi;
  19. lp.phi *= 0.895168 + V * ( 0.0218849 + V * 0.00826809 );
  20. for (i = NITER; i ; --i) {
  21. c = cos(lp.phi);
  22. s = sin(lp.phi);
  23. lp.phi -= V = (lp.phi + s * (c + 2.) - p) /
  24. (1. + c * (c + 2.) - s * s);
  25. if (fabs(V) < EPS)
  26. break;
  27. }
  28. if (!i) {
  29. xy.x = C_x * lp.lam;
  30. xy.y = lp.phi < 0. ? -C_y : C_y;
  31. } else {
  32. xy.x = C_x * lp.lam * (1. + cos(lp.phi));
  33. xy.y = C_y * sin(lp.phi);
  34. }
  35. return (xy);
  36. }
  37. INVERSE(s_inverse); /* spheroid */
  38. double c;
  39. lp.phi = aasin(xy.y / C_y);
  40. lp.lam = xy.x / (C_x * (1. + (c = cos(lp.phi))));
  41. lp.phi = aasin((lp.phi + sin(lp.phi) * (c + 2.)) / C_p);
  42. return (lp);
  43. }
  44. FREEUP; if (P) pj_dalloc(P); }
  45. ENTRY0(eck4); P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)