/Proj4/PJ_putp6.c

http://github.com/route-me/route-me · C · 62 lines · 60 code · 2 blank · 0 comment · 4 complexity · 97fa2093a30a8aaf02d04453ecf707cf MD5 · raw file

  1. #ifndef lint
  2. static const char SCCSID[]="@(#)PJ_putp6.c 4.1 94/02/15 GIE REL";
  3. #endif
  4. #define PROJ_PARMS__ \
  5. double C_x, C_y, A, B, D;
  6. #define PJ_LIB__
  7. #include "projects.h"
  8. PROJ_HEAD(putp6, "Putnins P6") "\n\tPCyl., Sph.";
  9. PROJ_HEAD(putp6p, "Putnins P6'") "\n\tPCyl., Sph.";
  10. #define EPS 1e-10
  11. #define NITER 10
  12. #define CON_POLE 1.732050807568877
  13. FORWARD(s_forward); /* spheroid */
  14. double p, r, V;
  15. int i;
  16. p = P->B * sin(lp.phi);
  17. lp.phi *= 1.10265779;
  18. for (i = NITER; i ; --i) {
  19. r = sqrt(1. + lp.phi * lp.phi);
  20. lp.phi -= V = ( (P->A - r) * lp.phi - log(lp.phi + r) - p ) /
  21. (P->A - 2. * r);
  22. if (fabs(V) < EPS)
  23. break;
  24. }
  25. if (!i)
  26. lp.phi = p < 0. ? -CON_POLE : CON_POLE;
  27. xy.x = P->C_x * lp.lam * (P->D - sqrt(1. + lp.phi * lp.phi));
  28. xy.y = P->C_y * lp.phi;
  29. return (xy);
  30. }
  31. INVERSE(s_inverse); /* spheroid */
  32. double r;
  33. lp.phi = xy.y / P->C_y;
  34. r = sqrt(1. + lp.phi * lp.phi);
  35. lp.lam = xy.x / (P->C_x * (P->D - r));
  36. lp.phi = aasin( ( (P->A - r) * lp.phi - log(lp.phi + r) ) / P->B);
  37. return (lp);
  38. }
  39. FREEUP; if (P) pj_dalloc(P); }
  40. static PJ *
  41. setup(PJ *P) {
  42. P->es = 0.;
  43. P->inv = s_inverse;
  44. P->fwd = s_forward;
  45. return P;
  46. }
  47. ENTRY0(putp6)
  48. P->C_x = 1.01346;
  49. P->C_y = 0.91910;
  50. P->A = 4.;
  51. P->B = 2.1471437182129378784;
  52. P->D = 2.;
  53. ENDENTRY(setup(P))
  54. ENTRY0(putp6p)
  55. P->C_x = 0.44329;
  56. P->C_y = 0.80404;
  57. P->A = 6.;
  58. P->B = 5.61125;
  59. P->D = 3.;
  60. ENDENTRY(setup(P))