/Proj4/pj_factors.c

http://github.com/route-me/route-me · C · 88 lines · 79 code · 1 blank · 8 comment · 18 complexity · 9a8bc6a72a6027262a74ee55f3e1fdda MD5 · raw file

  1. /* projection scale factors */
  2. #ifndef lint
  3. static const char SCCSID[]="@(#)pj_factors.c 4.9 94/03/17 GIE REL";
  4. #endif
  5. #define PJ_LIB__
  6. #include "projects.h"
  7. #include <errno.h>
  8. #ifndef DEFAULT_H
  9. #define DEFAULT_H 1e-5 /* radian default for numeric h */
  10. #endif
  11. #define EPS 1.0e-12
  12. int
  13. pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
  14. struct DERIVS der;
  15. double cosphi, t, n, r;
  16. /* check for forward and latitude or longitude overange */
  17. t = fabs(lp.phi)-HALFPI;
  18. if (t > EPS || fabs(lp.lam) > 10.) {
  19. pj_errno = -14;
  20. return 1;
  21. } else { /* proceed */
  22. errno = pj_errno = 0;
  23. if (h < EPS)
  24. h = DEFAULT_H;
  25. if (fabs(lp.phi) > (HALFPI - h))
  26. /* adjust to value around pi/2 where derived still exists*/
  27. lp.phi = lp.phi < 0. ? (-HALFPI+h) : (HALFPI-h);
  28. else if (P->geoc)
  29. lp.phi = atan(P->rone_es * tan(lp.phi));
  30. lp.lam -= P->lam0; /* compute del lp.lam */
  31. if (!P->over)
  32. lp.lam = adjlon(lp.lam); /* adjust del longitude */
  33. if (P->spc) /* get what projection analytic values */
  34. P->spc(lp, P, fac);
  35. if (((fac->code & (IS_ANAL_XL_YL+IS_ANAL_XP_YP)) !=
  36. (IS_ANAL_XL_YL+IS_ANAL_XP_YP)) &&
  37. pj_deriv(lp, h, P, &der))
  38. return 1;
  39. if (!(fac->code & IS_ANAL_XL_YL)) {
  40. fac->der.x_l = der.x_l;
  41. fac->der.y_l = der.y_l;
  42. }
  43. if (!(fac->code & IS_ANAL_XP_YP)) {
  44. fac->der.x_p = der.x_p;
  45. fac->der.y_p = der.y_p;
  46. }
  47. cosphi = cos(lp.phi);
  48. if (!(fac->code & IS_ANAL_HK)) {
  49. fac->h = hypot(fac->der.x_p, fac->der.y_p);
  50. fac->k = hypot(fac->der.x_l, fac->der.y_l) / cosphi;
  51. if (P->es) {
  52. t = sin(lp.phi);
  53. t = 1. - P->es * t * t;
  54. n = sqrt(t);
  55. fac->h *= t * n / P->one_es;
  56. fac->k *= n;
  57. r = t * t / P->one_es;
  58. } else
  59. r = 1.;
  60. } else if (P->es) {
  61. r = sin(lp.phi);
  62. r = 1. - P->es * r * r;
  63. r = r * r / P->one_es;
  64. } else
  65. r = 1.;
  66. /* convergence */
  67. if (!(fac->code & IS_ANAL_CONV)) {
  68. fac->conv = - atan2(fac->der.y_l, fac->der.x_l);
  69. if (fac->code & IS_ANAL_XL_YL)
  70. fac->code |= IS_ANAL_CONV;
  71. }
  72. /* areal scale factor */
  73. fac->s = (fac->der.y_p * fac->der.x_l - fac->der.x_p * fac->der.y_l) *
  74. r / cosphi;
  75. /* meridian-parallel angle theta prime */
  76. fac->thetap = aasin(fac->s / (fac->h * fac->k));
  77. /* Tissot ellips axis */
  78. t = fac->k * fac->k + fac->h * fac->h;
  79. fac->a = sqrt(t + 2. * fac->s);
  80. t = (t = t - 2. * fac->s) <= 0. ? 0. : sqrt(t);
  81. fac->b = 0.5 * (fac->a - t);
  82. fac->a = 0.5 * (fac->a + t);
  83. /* omega */
  84. fac->omega = 2. * aasin((fac->a - fac->b)/(fac->a + fac->b));
  85. }
  86. return 0;
  87. }