/Proj4/PJ_laea.c

http://github.com/route-me/route-me · C · 236 lines · 227 code · 6 blank · 3 comment · 43 complexity · 36d75befc602f80ec7a2bb1e389ce102 MD5 · raw file

  1. #ifndef lint
  2. static const char SCCSID[]="@(#)PJ_laea.c 4.1 94/02/15 GIE REL";
  3. #endif
  4. #define PROJ_PARMS__ \
  5. double sinb1; \
  6. double cosb1; \
  7. double xmf; \
  8. double ymf; \
  9. double mmf; \
  10. double qp; \
  11. double dd; \
  12. double rq; \
  13. double *apa; \
  14. int mode;
  15. #define PJ_LIB__
  16. #include "projects.h"
  17. PROJ_HEAD(laea, "Lambert Azimuthal Equal Area") "\n\tAzi, Sph&Ell";
  18. #define sinph0 P->sinb1
  19. #define cosph0 P->cosb1
  20. #define EPS10 1.e-10
  21. #define NITER 20
  22. #define CONV 1.e-10
  23. #define N_POLE 0
  24. #define S_POLE 1
  25. #define EQUIT 2
  26. #define OBLIQ 3
  27. FORWARD(e_forward); /* ellipsoid */
  28. double coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;
  29. coslam = cos(lp.lam);
  30. sinlam = sin(lp.lam);
  31. sinphi = sin(lp.phi);
  32. q = pj_qsfn(sinphi, P->e, P->one_es);
  33. if (P->mode == OBLIQ || P->mode == EQUIT) {
  34. sinb = q / P->qp;
  35. cosb = sqrt(1. - sinb * sinb);
  36. }
  37. switch (P->mode) {
  38. case OBLIQ:
  39. b = 1. + P->sinb1 * sinb + P->cosb1 * cosb * coslam;
  40. break;
  41. case EQUIT:
  42. b = 1. + cosb * coslam;
  43. break;
  44. case N_POLE:
  45. b = HALFPI + lp.phi;
  46. q = P->qp - q;
  47. break;
  48. case S_POLE:
  49. b = lp.phi - HALFPI;
  50. q = P->qp + q;
  51. break;
  52. }
  53. if (fabs(b) < EPS10) F_ERROR;
  54. switch (P->mode) {
  55. case OBLIQ:
  56. xy.y = P->ymf * ( b = sqrt(2. / b) )
  57. * (P->cosb1 * sinb - P->sinb1 * cosb * coslam);
  58. goto eqcon;
  59. break;
  60. case EQUIT:
  61. xy.y = (b = sqrt(2. / (1. + cosb * coslam))) * sinb * P->ymf;
  62. eqcon:
  63. xy.x = P->xmf * b * cosb * sinlam;
  64. break;
  65. case N_POLE:
  66. case S_POLE:
  67. if (q >= 0.) {
  68. xy.x = (b = sqrt(q)) * sinlam;
  69. xy.y = coslam * (P->mode == S_POLE ? b : -b);
  70. } else
  71. xy.x = xy.y = 0.;
  72. break;
  73. }
  74. return (xy);
  75. }
  76. FORWARD(s_forward); /* spheroid */
  77. double coslam, cosphi, sinphi;
  78. sinphi = sin(lp.phi);
  79. cosphi = cos(lp.phi);
  80. coslam = cos(lp.lam);
  81. switch (P->mode) {
  82. case EQUIT:
  83. xy.y = 1. + cosphi * coslam;
  84. goto oblcon;
  85. case OBLIQ:
  86. xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;
  87. oblcon:
  88. if (xy.y <= EPS10) F_ERROR;
  89. xy.x = (xy.y = sqrt(2. / xy.y)) * cosphi * sin(lp.lam);
  90. xy.y *= P->mode == EQUIT ? sinphi :
  91. cosph0 * sinphi - sinph0 * cosphi * coslam;
  92. break;
  93. case N_POLE:
  94. coslam = -coslam;
  95. case S_POLE:
  96. if (fabs(lp.phi + P->phi0) < EPS10) F_ERROR;
  97. xy.y = FORTPI - lp.phi * .5;
  98. xy.y = 2. * (P->mode == S_POLE ? cos(xy.y) : sin(xy.y));
  99. xy.x = xy.y * sin(lp.lam);
  100. xy.y *= coslam;
  101. break;
  102. }
  103. return (xy);
  104. }
  105. INVERSE(e_inverse); /* ellipsoid */
  106. double cCe, sCe, q, rho, ab=0.0;
  107. switch (P->mode) {
  108. case EQUIT:
  109. case OBLIQ:
  110. if ((rho = hypot(xy.x /= P->dd, xy.y *= P->dd)) < EPS10) {
  111. lp.lam = 0.;
  112. lp.phi = P->phi0;
  113. return (lp);
  114. }
  115. cCe = cos(sCe = 2. * asin(.5 * rho / P->rq));
  116. xy.x *= (sCe = sin(sCe));
  117. if (P->mode == OBLIQ) {
  118. q = P->qp * (ab = cCe * P->sinb1 + xy.y * sCe * P->cosb1 / rho);
  119. xy.y = rho * P->cosb1 * cCe - xy.y * P->sinb1 * sCe;
  120. } else {
  121. q = P->qp * (ab = xy.y * sCe / rho);
  122. xy.y = rho * cCe;
  123. }
  124. break;
  125. case N_POLE:
  126. xy.y = -xy.y;
  127. case S_POLE:
  128. if (!(q = (xy.x * xy.x + xy.y * xy.y)) ) {
  129. lp.lam = 0.;
  130. lp.phi = P->phi0;
  131. return (lp);
  132. }
  133. /*
  134. q = P->qp - q;
  135. */
  136. ab = 1. - q / P->qp;
  137. if (P->mode == S_POLE)
  138. ab = - ab;
  139. break;
  140. }
  141. lp.lam = atan2(xy.x, xy.y);
  142. lp.phi = pj_authlat(asin(ab), P->apa);
  143. return (lp);
  144. }
  145. INVERSE(s_inverse); /* spheroid */
  146. double cosz=0.0, rh, sinz=0.0;
  147. rh = hypot(xy.x, xy.y);
  148. if ((lp.phi = rh * .5 ) > 1.) I_ERROR;
  149. lp.phi = 2. * asin(lp.phi);
  150. if (P->mode == OBLIQ || P->mode == EQUIT) {
  151. sinz = sin(lp.phi);
  152. cosz = cos(lp.phi);
  153. }
  154. switch (P->mode) {
  155. case EQUIT:
  156. lp.phi = fabs(rh) <= EPS10 ? 0. : asin(xy.y * sinz / rh);
  157. xy.x *= sinz;
  158. xy.y = cosz * rh;
  159. break;
  160. case OBLIQ:
  161. lp.phi = fabs(rh) <= EPS10 ? P->phi0 :
  162. asin(cosz * sinph0 + xy.y * sinz * cosph0 / rh);
  163. xy.x *= sinz * cosph0;
  164. xy.y = (cosz - sin(lp.phi) * sinph0) * rh;
  165. break;
  166. case N_POLE:
  167. xy.y = -xy.y;
  168. lp.phi = HALFPI - lp.phi;
  169. break;
  170. case S_POLE:
  171. lp.phi -= HALFPI;
  172. break;
  173. }
  174. lp.lam = (xy.y == 0. && (P->mode == EQUIT || P->mode == OBLIQ)) ?
  175. 0. : atan2(xy.x, xy.y);
  176. return (lp);
  177. }
  178. FREEUP;
  179. if (P) {
  180. if (P->apa)
  181. pj_dalloc(P->apa);
  182. pj_dalloc(P);
  183. }
  184. }
  185. ENTRY1(laea,apa)
  186. double t;
  187. if (fabs((t = fabs(P->phi0)) - HALFPI) < EPS10)
  188. P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
  189. else if (fabs(t) < EPS10)
  190. P->mode = EQUIT;
  191. else
  192. P->mode = OBLIQ;
  193. if (P->es) {
  194. double sinphi;
  195. P->e = sqrt(P->es);
  196. P->qp = pj_qsfn(1., P->e, P->one_es);
  197. P->mmf = .5 / (1. - P->es);
  198. P->apa = pj_authset(P->es);
  199. switch (P->mode) {
  200. case N_POLE:
  201. case S_POLE:
  202. P->dd = 1.;
  203. break;
  204. case EQUIT:
  205. P->dd = 1. / (P->rq = sqrt(.5 * P->qp));
  206. P->xmf = 1.;
  207. P->ymf = .5 * P->qp;
  208. break;
  209. case OBLIQ:
  210. P->rq = sqrt(.5 * P->qp);
  211. sinphi = sin(P->phi0);
  212. P->sinb1 = pj_qsfn(sinphi, P->e, P->one_es) / P->qp;
  213. P->cosb1 = sqrt(1. - P->sinb1 * P->sinb1);
  214. P->dd = cos(P->phi0) / (sqrt(1. - P->es * sinphi * sinphi) *
  215. P->rq * P->cosb1);
  216. P->ymf = (P->xmf = P->rq) / P->dd;
  217. P->xmf *= P->dd;
  218. break;
  219. }
  220. P->inv = e_inverse;
  221. P->fwd = e_forward;
  222. } else {
  223. if (P->mode == OBLIQ) {
  224. sinph0 = sin(P->phi0);
  225. cosph0 = cos(P->phi0);
  226. }
  227. P->inv = s_inverse;
  228. P->fwd = s_forward;
  229. }
  230. ENDENTRY(P)