PageRenderTime 34ms CodeModel.GetById 23ms RepoModel.GetById 0ms app.codeStats 0ms

/Proj4/PJ_geos.c

http://github.com/route-me/route-me
C | 157 lines | 91 code | 6 blank | 60 comment | 9 complexity | 8fbd37c899538300c5336885d7f6dd1b MD5 | raw file
  1. /*
  2. ** libproj -- library of cartographic projections
  3. **
  4. ** Copyright (c) 2004 Gerald I. Evenden
  5. */
  6. static const char
  7. LIBPROJ_ID[] = "$Id: PJ_geos.c,v 1.2 2005/02/04 19:27:58 fwarmerdam Exp $";
  8. /*
  9. ** See also (section 4.4.3.2):
  10. ** http://www.eumetsat.int/en/area4/msg/news/us_doc/cgms_03_26.pdf
  11. **
  12. ** Permission is hereby granted, free of charge, to any person obtaining
  13. ** a copy of this software and associated documentation files (the
  14. ** "Software"), to deal in the Software without restriction, including
  15. ** without limitation the rights to use, copy, modify, merge, publish,
  16. ** distribute, sublicense, and/or sell copies of the Software, and to
  17. ** permit persons to whom the Software is furnished to do so, subject to
  18. ** the following conditions:
  19. **
  20. ** The above copyright notice and this permission notice shall be
  21. ** included in all copies or substantial portions of the Software.
  22. **
  23. ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  24. ** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  25. ** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
  26. ** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
  27. ** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
  28. ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
  29. ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  30. */
  31. #define PROJ_PARMS__ \
  32. double h; \
  33. double radius_p; \
  34. double radius_p2; \
  35. double radius_p_inv2; \
  36. double radius_g; \
  37. double radius_g_1; \
  38. double C;
  39. #define PJ_LIB__
  40. #include "projects.h"
  41. PROJ_HEAD(geos, "Geostationary Satellite View") "\n\tAzi, Sph&Ell\n\th=";
  42. FORWARD(s_forward); /* spheroid */
  43. double Vx, Vy, Vz, tmp;
  44. /* Calculation of the three components of the vector from satellite to
  45. ** position on earth surface (lon,lat).*/
  46. tmp = cos(lp.phi);
  47. Vx = cos (lp.lam) * tmp;
  48. Vy = sin (lp.lam) * tmp;
  49. Vz = sin (lp.phi);
  50. /* Check visibility.*/
  51. if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz) < 0.) F_ERROR;
  52. /* Calculation based on view angles from satellite.*/
  53. tmp = P->radius_g - Vx;
  54. xy.x = P->radius_g_1 * atan(Vy / tmp);
  55. xy.y = P->radius_g_1 * atan(Vz / hypot(Vy, tmp));
  56. return (xy);
  57. }
  58. FORWARD(e_forward); /* ellipsoid */
  59. double r, Vx, Vy, Vz, tmp;
  60. /* Calculation of geocentric latitude. */
  61. lp.phi = atan (P->radius_p2 * tan (lp.phi));
  62. /* Calculation of the three components of the vector from satellite to
  63. ** position on earth surface (lon,lat).*/
  64. r = (P->radius_p) / hypot(P->radius_p * cos (lp.phi), sin (lp.phi));
  65. Vx = r * cos (lp.lam) * cos (lp.phi);
  66. Vy = r * sin (lp.lam) * cos (lp.phi);
  67. Vz = r * sin (lp.phi);
  68. /* Check visibility. */
  69. if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * P->radius_p_inv2) < 0.)
  70. F_ERROR;
  71. /* Calculation based on view angles from satellite. */
  72. tmp = P->radius_g - Vx;
  73. xy.x = P->radius_g_1 * atan (Vy / tmp);
  74. xy.y = P->radius_g_1 * atan (Vz / hypot (Vy, tmp));
  75. return (xy);
  76. }
  77. INVERSE(s_inverse); /* spheroid */
  78. double Vx, Vy, Vz, a, b, c, det, k;
  79. /* Setting three components of vector from satellite to position.*/
  80. Vx = -1.0;
  81. Vy = tan (xy.x / (P->radius_g - 1.0));
  82. Vz = tan (xy.y / (P->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy);
  83. /* Calculation of terms in cubic equation and determinant.*/
  84. a = Vy * Vy + Vz * Vz + Vx * Vx;
  85. b = 2 * P->radius_g * Vx;
  86. if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR;
  87. /* Calculation of three components of vector from satellite to position.*/
  88. k = (-b - sqrt(det)) / (2 * a);
  89. Vx = P->radius_g + k * Vx;
  90. Vy *= k;
  91. Vz *= k;
  92. /* Calculation of longitude and latitude.*/
  93. lp.lam = atan2 (Vy, Vx);
  94. lp.phi = atan (Vz * cos (lp.lam) / Vx);
  95. return (lp);
  96. }
  97. INVERSE(e_inverse); /* ellipsoid */
  98. double Vx, Vy, Vz, a, b, c, det, k;
  99. /* Setting three components of vector from satellite to position.*/
  100. Vx = -1.0;
  101. Vy = tan (xy.x / P->radius_g_1);
  102. Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy);
  103. /* Calculation of terms in cubic equation and determinant.*/
  104. a = Vz / P->radius_p;
  105. a = Vy * Vy + a * a + Vx * Vx;
  106. b = 2 * P->radius_g * Vx;
  107. if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR;
  108. /* Calculation of three components of vector from satellite to position.*/
  109. k = (-b - sqrt(det)) / (2. * a);
  110. Vx = P->radius_g + k * Vx;
  111. Vy *= k;
  112. Vz *= k;
  113. /* Calculation of longitude and latitude.*/
  114. lp.lam = atan2 (Vy, Vx);
  115. lp.phi = atan (Vz * cos (lp.lam) / Vx);
  116. lp.phi = atan (P->radius_p_inv2 * tan (lp.phi));
  117. return (lp);
  118. }
  119. FREEUP; if (P) free(P); }
  120. ENTRY0(geos)
  121. if ((P->h = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30);
  122. if (P->phi0) E_ERROR(-46);
  123. P->radius_g = 1. + (P->radius_g_1 = P->h / P->a);
  124. P->C = P->radius_g * P->radius_g - 1.0;
  125. if (P->es) {
  126. P->radius_p = sqrt (P->one_es);
  127. P->radius_p2 = P->one_es;
  128. P->radius_p_inv2 = P->rone_es;
  129. P->inv = e_inverse;
  130. P->fwd = e_forward;
  131. } else {
  132. P->radius_p = P->radius_p2 = P->radius_p_inv2 = 1.0;
  133. P->inv = s_inverse;
  134. P->fwd = s_forward;
  135. }
  136. ENDENTRY(P)
  137. /*
  138. ** $Log: PJ_geos.c,v $
  139. ** Revision 1.2 2005/02/04 19:27:58 fwarmerdam
  140. ** Added link to reference info.
  141. **
  142. ** Revision 1.1 2004/10/20 17:04:00 fwarmerdam
  143. ** New
  144. **
  145. ** Revision 1.2 2004/07/14 18:08:57 gie
  146. ** corrected P->phi_0 to P->phi0
  147. **
  148. ** Revision 1.1 2004/07/12 17:58:25 gie
  149. ** Initial revision
  150. **
  151. */