/Proj4/PJ_aea.c

http://github.com/route-me/route-me · C · 171 lines · 124 code · 9 blank · 38 comment · 22 complexity · f2b6cff19e17b55170f35cd81aac6ea4 MD5 · raw file

  1. /******************************************************************************
  2. * $Id: PJ_aea.c,v 1.4 2003/08/18 15:21:23 warmerda Exp $
  3. *
  4. * Project: PROJ.4
  5. * Purpose: Implementation of the aea (Albers Equal Area) 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_aea.c,v $
  31. * Revision 1.4 2003/08/18 15:21:23 warmerda
  32. * fixed initialization of en variable
  33. *
  34. * Revision 1.3 2002/12/14 19:27:06 warmerda
  35. * updated header
  36. *
  37. */
  38. #define PROJ_PARMS__ \
  39. double ec; \
  40. double n; \
  41. double c; \
  42. double dd; \
  43. double n2; \
  44. double rho0; \
  45. double rho; \
  46. double phi1; \
  47. double phi2; \
  48. double *en; \
  49. int ellips;
  50. #define PJ_LIB__
  51. #include "projects.h"
  52. PJ_CVSID("$Id: PJ_aea.c,v 1.4 2003/08/18 15:21:23 warmerda Exp $");
  53. # define EPS10 1.e-10
  54. # define TOL7 1.e-7
  55. PROJ_HEAD(aea, "Albers Equal Area")
  56. "\n\tConic Sph&Ell\n\tlat_1= lat_2=";
  57. PROJ_HEAD(leac, "Lambert Equal Area Conic")
  58. "\n\tConic, Sph&Ell\n\tlat_1= south";
  59. /* determine latitude angle phi-1 */
  60. # define N_ITER 15
  61. # define EPSILON 1.0e-7
  62. # define TOL 1.0e-10
  63. static double
  64. phi1_(double qs, double Te, double Tone_es) {
  65. int i;
  66. double Phi, sinpi, cospi, con, com, dphi;
  67. Phi = asin (.5 * qs);
  68. if (Te < EPSILON)
  69. return( Phi );
  70. i = N_ITER;
  71. do {
  72. sinpi = sin (Phi);
  73. cospi = cos (Phi);
  74. con = Te * sinpi;
  75. com = 1. - con * con;
  76. dphi = .5 * com * com / cospi * (qs / Tone_es -
  77. sinpi / com + .5 / Te * log ((1. - con) /
  78. (1. + con)));
  79. Phi += dphi;
  80. } while (fabs(dphi) > TOL && --i);
  81. return( i ? Phi : HUGE_VAL );
  82. }
  83. FORWARD(e_forward); /* ellipsoid & spheroid */
  84. if ((P->rho = P->c - (P->ellips ? P->n * pj_qsfn(sin(lp.phi),
  85. P->e, P->one_es) : P->n2 * sin(lp.phi))) < 0.) F_ERROR
  86. P->rho = P->dd * sqrt(P->rho);
  87. xy.x = P->rho * sin( lp.lam *= P->n );
  88. xy.y = P->rho0 - P->rho * cos(lp.lam);
  89. return (xy);
  90. }
  91. INVERSE(e_inverse) /* ellipsoid & spheroid */;
  92. if( (P->rho = hypot(xy.x, xy.y = P->rho0 - xy.y)) != 0.0 ) {
  93. if (P->n < 0.) {
  94. P->rho = -P->rho;
  95. xy.x = -xy.x;
  96. xy.y = -xy.y;
  97. }
  98. lp.phi = P->rho / P->dd;
  99. if (P->ellips) {
  100. lp.phi = (P->c - lp.phi * lp.phi) / P->n;
  101. if (fabs(P->ec - fabs(lp.phi)) > TOL7) {
  102. if ((lp.phi = phi1_(lp.phi, P->e, P->one_es)) == HUGE_VAL)
  103. I_ERROR
  104. } else
  105. lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
  106. } else if (fabs(lp.phi = (P->c - lp.phi * lp.phi) / P->n2) <= 1.)
  107. lp.phi = asin(lp.phi);
  108. else
  109. lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
  110. lp.lam = atan2(xy.x, xy.y) / P->n;
  111. } else {
  112. lp.lam = 0.;
  113. lp.phi = P->n > 0. ? HALFPI : - HALFPI;
  114. }
  115. return (lp);
  116. }
  117. FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
  118. static PJ *
  119. setup(PJ *P) {
  120. double cosphi, sinphi;
  121. int secant;
  122. if (fabs(P->phi1 + P->phi2) < EPS10) E_ERROR(-21);
  123. P->n = sinphi = sin(P->phi1);
  124. cosphi = cos(P->phi1);
  125. secant = fabs(P->phi1 - P->phi2) >= EPS10;
  126. if( (P->ellips = (P->es > 0.))) {
  127. double ml1, m1;
  128. if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
  129. m1 = pj_msfn(sinphi, cosphi, P->es);
  130. ml1 = pj_qsfn(sinphi, P->e, P->one_es);
  131. if (secant) { /* secant cone */
  132. double ml2, m2;
  133. sinphi = sin(P->phi2);
  134. cosphi = cos(P->phi2);
  135. m2 = pj_msfn(sinphi, cosphi, P->es);
  136. ml2 = pj_qsfn(sinphi, P->e, P->one_es);
  137. P->n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
  138. }
  139. P->ec = 1. - .5 * P->one_es * log((1. - P->e) /
  140. (1. + P->e)) / P->e;
  141. P->c = m1 * m1 + P->n * ml1;
  142. P->dd = 1. / P->n;
  143. P->rho0 = P->dd * sqrt(P->c - P->n * pj_qsfn(sin(P->phi0),
  144. P->e, P->one_es));
  145. } else {
  146. if (secant) P->n = .5 * (P->n + sin(P->phi2));
  147. P->n2 = P->n + P->n;
  148. P->c = cosphi * cosphi + P->n2 * sinphi;
  149. P->dd = 1. / P->n;
  150. P->rho0 = P->dd * sqrt(P->c - P->n2 * sin(P->phi0));
  151. }
  152. P->inv = e_inverse; P->fwd = e_forward;
  153. return P;
  154. }
  155. ENTRY1(aea,en)
  156. P->phi1 = pj_param(P->params, "rlat_1").f;
  157. P->phi2 = pj_param(P->params, "rlat_2").f;
  158. ENDENTRY(setup(P))
  159. ENTRY1(leac,en)
  160. P->phi2 = pj_param(P->params, "rlat_1").f;
  161. P->phi1 = pj_param(P->params, "bsouth").i ? - HALFPI: HALFPI;
  162. ENDENTRY(setup(P))