/Proj4/PJ_airy.c

http://github.com/route-me/route-me · C · 130 lines · 90 code · 6 blank · 34 comment · 22 complexity · b9aaab73accf5d2d28cc4e4fa911fcd2 MD5 · raw file

  1. /******************************************************************************
  2. * $Id: PJ_airy.c,v 1.2 2002/12/14 19:30:40 warmerda Exp $
  3. *
  4. * Project: PROJ.4
  5. * Purpose: Implementation of the airy (Airy) 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_airy.c,v $
  31. * Revision 1.2 2002/12/14 19:30:40 warmerda
  32. * updated header
  33. *
  34. */
  35. #define PROJ_PARMS__ \
  36. double p_halfpi; \
  37. double sinph0; \
  38. double cosph0; \
  39. double Cb; \
  40. int mode; \
  41. int no_cut; /* do not cut at hemisphere limit */
  42. #define PJ_LIB__
  43. #include "projects.h"
  44. PJ_CVSID("$Id: PJ_airy.c,v 1.2 2002/12/14 19:30:40 warmerda Exp $");
  45. PROJ_HEAD(airy, "Airy") "\n\tMisc Sph, no inv.\n\tno_cut lat_b=";
  46. # define EPS 1.e-10
  47. # define N_POLE 0
  48. # define S_POLE 1
  49. # define EQUIT 2
  50. # define OBLIQ 3
  51. FORWARD(s_forward); /* spheroid */
  52. double sinlam, coslam, cosphi, sinphi, t, s, Krho, cosz;
  53. sinlam = sin(lp.lam);
  54. coslam = cos(lp.lam);
  55. switch (P->mode) {
  56. case EQUIT:
  57. case OBLIQ:
  58. sinphi = sin(lp.phi);
  59. cosphi = cos(lp.phi);
  60. cosz = cosphi * coslam;
  61. if (P->mode == OBLIQ)
  62. cosz = P->sinph0 * sinphi + P->cosph0 * cosz;
  63. if (!P->no_cut && cosz < -EPS)
  64. F_ERROR;
  65. if (fabs(s = 1. - cosz) > EPS) {
  66. t = 0.5 * (1. + cosz);
  67. Krho = -log(t)/s - P->Cb / t;
  68. } else
  69. Krho = 0.5 - P->Cb;
  70. xy.x = Krho * cosphi * sinlam;
  71. if (P->mode == OBLIQ)
  72. xy.y = Krho * (P->cosph0 * sinphi -
  73. P->sinph0 * cosphi * coslam);
  74. else
  75. xy.y = Krho * sinphi;
  76. break;
  77. case S_POLE:
  78. case N_POLE:
  79. lp.phi = fabs(P->p_halfpi - lp.phi);
  80. if (!P->no_cut && (lp.phi - EPS) > HALFPI)
  81. F_ERROR;
  82. if ((lp.phi *= 0.5) > EPS) {
  83. t = tan(lp.phi);
  84. Krho = -2.*(log(cos(lp.phi)) / t + t * P->Cb);
  85. xy.x = Krho * sinlam;
  86. xy.y = Krho * coslam;
  87. if (P->mode == N_POLE)
  88. xy.y = -xy.y;
  89. } else
  90. xy.x = xy.y = 0.;
  91. }
  92. return (xy);
  93. }
  94. FREEUP; if (P) pj_dalloc(P); }
  95. ENTRY0(airy)
  96. double beta;
  97. P->no_cut = pj_param(P->params, "bno_cut").i;
  98. beta = 0.5 * (HALFPI - pj_param(P->params, "rlat_b").f);
  99. if (fabs(beta) < EPS)
  100. P->Cb = -0.5;
  101. else {
  102. P->Cb = 1./tan(beta);
  103. P->Cb *= P->Cb * log(cos(beta));
  104. }
  105. if (fabs(fabs(P->phi0) - HALFPI) < EPS)
  106. if (P->phi0 < 0.) {
  107. P->p_halfpi = -HALFPI;
  108. P->mode = S_POLE;
  109. } else {
  110. P->p_halfpi = HALFPI;
  111. P->mode = N_POLE;
  112. }
  113. else {
  114. if (fabs(P->phi0) < EPS)
  115. P->mode = EQUIT;
  116. else {
  117. P->mode = OBLIQ;
  118. P->sinph0 = sin(P->phi0);
  119. P->cosph0 = cos(P->phi0);
  120. }
  121. }
  122. P->fwd = s_forward;
  123. P->es = 0.;
  124. ENDENTRY(P)