PageRenderTime 32ms CodeModel.GetById 19ms RepoModel.GetById 0ms app.codeStats 0ms

/Proj4/PJ_mbt_fps.c

http://github.com/route-me/route-me
C | 41 lines | 39 code | 2 blank | 0 comment | 3 complexity | c85925fedd16e409b2cb925dbee37078 MD5 | raw file
  1. #ifndef lint
  2. static const char SCCSID[]="@(#)PJ_mbt_fps.c 4.1 94/02/15 GIE REL";
  3. #endif
  4. #define PJ_LIB__
  5. #include "projects.h"
  6. PROJ_HEAD(mbt_fps, "McBryde-Thomas Flat-Pole Sine (No. 2)") "\n\tCyl., Sph.";
  7. #define MAX_ITER 10
  8. #define LOOP_TOL 1e-7
  9. #define C1 0.45503
  10. #define C2 1.36509
  11. #define C3 1.41546
  12. #define C_x 0.22248
  13. #define C_y 1.44492
  14. #define C1_2 0.33333333333333333333333333
  15. FORWARD(s_forward); /* spheroid */
  16. double k, V, t;
  17. int i;
  18. k = C3 * sin(lp.phi);
  19. for (i = MAX_ITER; i ; --i) {
  20. t = lp.phi / C2;
  21. lp.phi -= V = (C1 * sin(t) + sin(lp.phi) - k) /
  22. (C1_2 * cos(t) + cos(lp.phi));
  23. if (fabs(V) < LOOP_TOL)
  24. break;
  25. }
  26. t = lp.phi / C2;
  27. xy.x = C_x * lp.lam * (1. + 3. * cos(lp.phi)/cos(t) );
  28. xy.y = C_y * sin(t);
  29. return (xy);
  30. }
  31. INVERSE(s_inverse); /* spheroid */
  32. double t, s;
  33. lp.phi = C2 * (t = aasin(xy.y / C_y));
  34. lp.lam = xy.x / (C_x * (1. + 3. * cos(lp.phi)/cos(t)));
  35. lp.phi = aasin((C1 * sin(t) + sin(lp.phi)) / C3);
  36. return (lp);
  37. }
  38. FREEUP; if (P) pj_dalloc(P); }
  39. ENTRY0(mbt_fps) P->es = 0; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)