/Proj4/PJ_mbtfpq.c

http://github.com/route-me/route-me · C · 50 lines · 48 code · 2 blank · 0 comment · 10 complexity · 07b0ecb3ae08242322f67b5cb32108b4 MD5 · raw file

  1. #ifndef lint
  2. static const char SCCSID[]="@(#)PJ_mbtfpq.c 4.1 94/02/15 GIE REL";
  3. #endif
  4. #define PJ_LIB__
  5. #include "projects.h"
  6. PROJ_HEAD(mbtfpq, "McBryde-Thomas Flat-Polar Quartic") "\n\tCyl., Sph.";
  7. #define NITER 20
  8. #define EPS 1e-7
  9. #define ONETOL 1.000001
  10. #define C 1.70710678118654752440
  11. #define RC 0.58578643762690495119
  12. #define FYC 1.87475828462269495505
  13. #define RYC 0.53340209679417701685
  14. #define FXC 0.31245971410378249250
  15. #define RXC 3.20041258076506210122
  16. FORWARD(s_forward); /* spheroid */
  17. double th1, c;
  18. int i;
  19. c = C * sin(lp.phi);
  20. for (i = NITER; i; --i) {
  21. lp.phi -= th1 = (sin(.5*lp.phi) + sin(lp.phi) - c) /
  22. (.5*cos(.5*lp.phi) + cos(lp.phi));
  23. if (fabs(th1) < EPS) break;
  24. }
  25. xy.x = FXC * lp.lam * (1.0 + 2. * cos(lp.phi)/cos(0.5 * lp.phi));
  26. xy.y = FYC * sin(0.5 * lp.phi);
  27. return (xy);
  28. }
  29. INVERSE(s_inverse); /* spheroid */
  30. double t;
  31. lp.phi = RYC * xy.y;
  32. if (fabs(lp.phi) > 1.) {
  33. if (fabs(lp.phi) > ONETOL) I_ERROR
  34. else if (lp.phi < 0.) { t = -1.; lp.phi = -PI; }
  35. else { t = 1.; lp.phi = PI; }
  36. } else
  37. lp.phi = 2. * asin(t = lp.phi);
  38. lp.lam = RXC * xy.x / (1. + 2. * cos(lp.phi)/cos(0.5 * lp.phi));
  39. lp.phi = RC * (t + sin(lp.phi));
  40. if (fabs(lp.phi) > 1.)
  41. if (fabs(lp.phi) > ONETOL) I_ERROR
  42. else lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
  43. else
  44. lp.phi = asin(lp.phi);
  45. return (lp);
  46. }
  47. FREEUP; if (P) pj_dalloc(P); }
  48. ENTRY0(mbtfpq) P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)