/Proj4/PJ_aeqd.c

http://github.com/route-me/route-me · C · 277 lines · 232 code · 11 blank · 34 comment · 49 complexity · 734bf86c1287f4407c696e528eacc50b MD5 · raw file

  1. /******************************************************************************
  2. * $Id: PJ_aeqd.c,v 1.3 2002/12/14 19:27:06 warmerda Exp $
  3. *
  4. * Project: PROJ.4
  5. * Purpose: Implementation of the aeqd (Azimuthal Equidistant) projection.
  6. * Author: Gerald Evenden
  7. *
  8. ******************************************************************************
  9. * Copyright (c) 1995, Gerald Evenden
  10. *
  11. * Permission is hereby granted, free of charge, to any person obtaining a
  12. * copy of this software and associated documentation files (the "Software"),
  13. * to deal in the Software without restriction, including without limitation
  14. * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  15. * and/or sell copies of the Software, and to permit persons to whom the
  16. * Software is furnished to do so, subject to the following conditions:
  17. *
  18. * The above copyright notice and this permission notice shall be included
  19. * in all copies or substantial portions of the Software.
  20. *
  21. * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  22. * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  23. * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  24. * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  25. * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  26. * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  27. * DEALINGS IN THE SOFTWARE.
  28. ******************************************************************************
  29. *
  30. * $Log: PJ_aeqd.c,v $
  31. * Revision 1.3 2002/12/14 19:27:06 warmerda
  32. * updated header
  33. *
  34. */
  35. #define PROJ_PARMS__ \
  36. double sinph0; \
  37. double cosph0; \
  38. double *en; \
  39. double M1; \
  40. double N1; \
  41. double Mp; \
  42. double He; \
  43. double G; \
  44. int mode;
  45. #define PJ_LIB__
  46. #include "projects.h"
  47. PJ_CVSID("$Id: PJ_aeqd.c,v 1.3 2002/12/14 19:27:06 warmerda Exp $");
  48. PROJ_HEAD(aeqd, "Azimuthal Equidistant") "\n\tAzi, Sph&Ell\n\tlat_0 guam";
  49. #define EPS10 1.e-10
  50. #define TOL 1.e-14
  51. #define N_POLE 0
  52. #define S_POLE 1
  53. #define EQUIT 2
  54. #define OBLIQ 3
  55. FORWARD(e_guam_fwd); /* Guam elliptical */
  56. double cosphi, sinphi, t;
  57. cosphi = cos(lp.phi);
  58. sinphi = sin(lp.phi);
  59. t = 1. / sqrt(1. - P->es * sinphi * sinphi);
  60. xy.x = lp.lam * cosphi * t;
  61. xy.y = pj_mlfn(lp.phi, sinphi, cosphi, P->en) - P->M1 +
  62. .5 * lp.lam * lp.lam * cosphi * sinphi * t;
  63. return (xy);
  64. }
  65. FORWARD(e_forward); /* elliptical */
  66. double coslam, cosphi, sinphi, rho, s, H, H2, c, Az, t, ct, st, cA, sA;
  67. coslam = cos(lp.lam);
  68. cosphi = cos(lp.phi);
  69. sinphi = sin(lp.phi);
  70. switch (P->mode) {
  71. case N_POLE:
  72. coslam = - coslam;
  73. case S_POLE:
  74. xy.x = (rho = fabs(P->Mp - pj_mlfn(lp.phi, sinphi, cosphi, P->en))) *
  75. sin(lp.lam);
  76. xy.y = rho * coslam;
  77. break;
  78. case EQUIT:
  79. case OBLIQ:
  80. if (fabs(lp.lam) < EPS10 && fabs(lp.phi - P->phi0) < EPS10) {
  81. xy.x = xy.y = 0.;
  82. break;
  83. }
  84. t = atan2(P->one_es * sinphi + P->es * P->N1 * P->sinph0 *
  85. sqrt(1. - P->es * sinphi * sinphi), cosphi);
  86. ct = cos(t); st = sin(t);
  87. Az = atan2(sin(lp.lam) * ct, P->cosph0 * st - P->sinph0 * coslam * ct);
  88. cA = cos(Az); sA = sin(Az);
  89. s = aasin( fabs(sA) < TOL ?
  90. (P->cosph0 * st - P->sinph0 * coslam * ct) / cA :
  91. sin(lp.lam) * ct / sA );
  92. H = P->He * cA;
  93. H2 = H * H;
  94. c = P->N1 * s * (1. + s * s * (- H2 * (1. - H2)/6. +
  95. s * ( P->G * H * (1. - 2. * H2 * H2) / 8. +
  96. s * ((H2 * (4. - 7. * H2) - 3. * P->G * P->G * (1. - 7. * H2)) /
  97. 120. - s * P->G * H / 48.))));
  98. xy.x = c * sA;
  99. xy.y = c * cA;
  100. break;
  101. }
  102. return (xy);
  103. }
  104. FORWARD(s_forward); /* spherical */
  105. double coslam, cosphi, sinphi;
  106. sinphi = sin(lp.phi);
  107. cosphi = cos(lp.phi);
  108. coslam = cos(lp.lam);
  109. switch (P->mode) {
  110. case EQUIT:
  111. xy.y = cosphi * coslam;
  112. goto oblcon;
  113. case OBLIQ:
  114. xy.y = P->sinph0 * sinphi + P->cosph0 * cosphi * coslam;
  115. oblcon:
  116. if (fabs(fabs(xy.y) - 1.) < TOL)
  117. if (xy.y < 0.)
  118. F_ERROR
  119. else
  120. xy.x = xy.y = 0.;
  121. else {
  122. xy.y = acos(xy.y);
  123. xy.y /= sin(xy.y);
  124. xy.x = xy.y * cosphi * sin(lp.lam);
  125. xy.y *= (P->mode == EQUIT) ? sinphi :
  126. P->cosph0 * sinphi - P->sinph0 * cosphi * coslam;
  127. }
  128. break;
  129. case N_POLE:
  130. lp.phi = -lp.phi;
  131. coslam = -coslam;
  132. case S_POLE:
  133. if (fabs(lp.phi - HALFPI) < EPS10) F_ERROR;
  134. xy.x = (xy.y = (HALFPI + lp.phi)) * sin(lp.lam);
  135. xy.y *= coslam;
  136. break;
  137. }
  138. return (xy);
  139. }
  140. INVERSE(e_guam_inv); /* Guam elliptical */
  141. double x2, t;
  142. int i;
  143. x2 = 0.5 * xy.x * xy.x;
  144. lp.phi = P->phi0;
  145. for (i = 0; i < 3; ++i) {
  146. t = P->e * sin(lp.phi);
  147. lp.phi = pj_inv_mlfn(P->M1 + xy.y -
  148. x2 * tan(lp.phi) * (t = sqrt(1. - t * t)), P->es, P->en);
  149. }
  150. lp.lam = xy.x * t / cos(lp.phi);
  151. return (lp);
  152. }
  153. INVERSE(e_inverse); /* elliptical */
  154. double c, Az, cosAz, A, B, D, E, F, psi, t;
  155. if ((c = hypot(xy.x, xy.y)) < EPS10) {
  156. lp.phi = P->phi0;
  157. lp.lam = 0.;
  158. return (lp);
  159. }
  160. if (P->mode == OBLIQ || P->mode == EQUIT) {
  161. cosAz = cos(Az = atan2(xy.x, xy.y));
  162. t = P->cosph0 * cosAz;
  163. B = P->es * t / P->one_es;
  164. A = - B * t;
  165. B *= 3. * (1. - A) * P->sinph0;
  166. D = c / P->N1;
  167. E = D * (1. - D * D * (A * (1. + A) / 6. + B * (1. + 3.*A) * D / 24.));
  168. F = 1. - E * E * (A / 2. + B * E / 6.);
  169. psi = aasin(P->sinph0 * cos(E) + t * sin(E));
  170. lp.lam = aasin(sin(Az) * sin(E) / cos(psi));
  171. if ((t = fabs(psi)) < EPS10)
  172. lp.phi = 0.;
  173. else if (fabs(t - HALFPI) < 0.)
  174. lp.phi = HALFPI;
  175. else
  176. lp.phi = atan((1. - P->es * F * P->sinph0 / sin(psi)) * tan(psi) /
  177. P->one_es);
  178. } else { /* Polar */
  179. lp.phi = pj_inv_mlfn(P->mode == N_POLE ? P->Mp - c : P->Mp + c,
  180. P->es, P->en);
  181. lp.lam = atan2(xy.x, P->mode == N_POLE ? -xy.y : xy.y);
  182. }
  183. return (lp);
  184. }
  185. INVERSE(s_inverse); /* spherical */
  186. double cosc, c_rh, sinc;
  187. if ((c_rh = hypot(xy.x, xy.y)) > PI) {
  188. if (c_rh - EPS10 > PI) I_ERROR;
  189. c_rh = PI;
  190. } else if (c_rh < EPS10) {
  191. lp.phi = P->phi0;
  192. lp.lam = 0.;
  193. return (lp);
  194. }
  195. if (P->mode == OBLIQ || P->mode == EQUIT) {
  196. sinc = sin(c_rh);
  197. cosc = cos(c_rh);
  198. if (P->mode == EQUIT) {
  199. lp.phi = aasin(xy.y * sinc / c_rh);
  200. xy.x *= sinc;
  201. xy.y = cosc * c_rh;
  202. } else {
  203. lp.phi = aasin(cosc * P->sinph0 + xy.y * sinc * P->cosph0 /
  204. c_rh);
  205. xy.y = (cosc - P->sinph0 * sin(lp.phi)) * c_rh;
  206. xy.x *= sinc * P->cosph0;
  207. }
  208. lp.lam = xy.y == 0. ? 0. : atan2(xy.x, xy.y);
  209. } else if (P->mode == N_POLE) {
  210. lp.phi = HALFPI - c_rh;
  211. lp.lam = atan2(xy.x, -xy.y);
  212. } else {
  213. lp.phi = c_rh - HALFPI;
  214. lp.lam = atan2(xy.x, xy.y);
  215. }
  216. return (lp);
  217. }
  218. FREEUP;
  219. if (P) {
  220. if (P->en)
  221. pj_dalloc(P->en);
  222. pj_dalloc(P);
  223. }
  224. }
  225. ENTRY1(aeqd, en)
  226. P->phi0 = pj_param(P->params, "rlat_0").f;
  227. if (fabs(fabs(P->phi0) - HALFPI) < EPS10) {
  228. P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
  229. P->sinph0 = P->phi0 < 0. ? -1. : 1.;
  230. P->cosph0 = 0.;
  231. } else if (fabs(P->phi0) < EPS10) {
  232. P->mode = EQUIT;
  233. P->sinph0 = 0.;
  234. P->cosph0 = 1.;
  235. } else {
  236. P->mode = OBLIQ;
  237. P->sinph0 = sin(P->phi0);
  238. P->cosph0 = cos(P->phi0);
  239. }
  240. if (! P->es) {
  241. P->inv = s_inverse; P->fwd = s_forward;
  242. } else {
  243. if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
  244. if (pj_param(P->params, "bguam").i) {
  245. P->M1 = pj_mlfn(P->phi0, P->sinph0, P->cosph0, P->en);
  246. P->inv = e_guam_inv; P->fwd = e_guam_fwd;
  247. } else {
  248. switch (P->mode) {
  249. case N_POLE:
  250. P->Mp = pj_mlfn(HALFPI, 1., 0., P->en);
  251. break;
  252. case S_POLE:
  253. P->Mp = pj_mlfn(-HALFPI, -1., 0., P->en);
  254. break;
  255. case EQUIT:
  256. case OBLIQ:
  257. P->inv = e_inverse; P->fwd = e_forward;
  258. P->N1 = 1. / sqrt(1. - P->es * P->sinph0 * P->sinph0);
  259. P->G = P->sinph0 * (P->He = P->e / sqrt(P->one_es));
  260. P->He *= P->cosph0;
  261. break;
  262. }
  263. P->inv = e_inverse; P->fwd = e_forward;
  264. }
  265. }
  266. ENDENTRY(P)