/Proj4/PJ_geos.c
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*/ 6static const char 7LIBPROJ_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 42PROJ_HEAD(geos, "Geostationary Satellite View") "\n\tAzi, Sph&Ell\n\th="; 43 44FORWARD(s_forward); /* spheroid */ 45 double Vx, Vy, Vz, tmp; 46 47/* Calculation of the three components of the vector from satellite to 48** position on earth surface (lon,lat).*/ 49 tmp = cos(lp.phi); 50 Vx = cos (lp.lam) * tmp; 51 Vy = sin (lp.lam) * tmp; 52 Vz = sin (lp.phi); 53/* Check visibility.*/ 54 if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz) < 0.) F_ERROR; 55/* Calculation based on view angles from satellite.*/ 56 tmp = P->radius_g - Vx; 57 xy.x = P->radius_g_1 * atan(Vy / tmp); 58 xy.y = P->radius_g_1 * atan(Vz / hypot(Vy, tmp)); 59 return (xy); 60} 61FORWARD(e_forward); /* ellipsoid */ 62 double r, Vx, Vy, Vz, tmp; 63 64/* Calculation of geocentric latitude. */ 65 lp.phi = atan (P->radius_p2 * tan (lp.phi)); 66/* Calculation of the three components of the vector from satellite to 67** position on earth surface (lon,lat).*/ 68 r = (P->radius_p) / hypot(P->radius_p * cos (lp.phi), sin (lp.phi)); 69 Vx = r * cos (lp.lam) * cos (lp.phi); 70 Vy = r * sin (lp.lam) * cos (lp.phi); 71 Vz = r * sin (lp.phi); 72/* Check visibility. */ 73 if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * P->radius_p_inv2) < 0.) 74 F_ERROR; 75/* Calculation based on view angles from satellite. */ 76 tmp = P->radius_g - Vx; 77 xy.x = P->radius_g_1 * atan (Vy / tmp); 78 xy.y = P->radius_g_1 * atan (Vz / hypot (Vy, tmp)); 79 return (xy); 80} 81INVERSE(s_inverse); /* spheroid */ 82 double Vx, Vy, Vz, a, b, c, det, k; 83 84/* Setting three components of vector from satellite to position.*/ 85 Vx = -1.0; 86 Vy = tan (xy.x / (P->radius_g - 1.0)); 87 Vz = tan (xy.y / (P->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy); 88/* Calculation of terms in cubic equation and determinant.*/ 89 a = Vy * Vy + Vz * Vz + Vx * Vx; 90 b = 2 * P->radius_g * Vx; 91 if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR; 92/* Calculation of three components of vector from satellite to position.*/ 93 k = (-b - sqrt(det)) / (2 * a); 94 Vx = P->radius_g + k * Vx; 95 Vy *= k; 96 Vz *= k; 97/* Calculation of longitude and latitude.*/ 98 lp.lam = atan2 (Vy, Vx); 99 lp.phi = atan (Vz * cos (lp.lam) / Vx); 100 return (lp); 101} 102INVERSE(e_inverse); /* ellipsoid */ 103 double Vx, Vy, Vz, a, b, c, det, k; 104 105/* Setting three components of vector from satellite to position.*/ 106 Vx = -1.0; 107 Vy = tan (xy.x / P->radius_g_1); 108 Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy); 109/* Calculation of terms in cubic equation and determinant.*/ 110 a = Vz / P->radius_p; 111 a = Vy * Vy + a * a + Vx * Vx; 112 b = 2 * P->radius_g * Vx; 113 if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR; 114/* Calculation of three components of vector from satellite to position.*/ 115 k = (-b - sqrt(det)) / (2. * a); 116 Vx = P->radius_g + k * Vx; 117 Vy *= k; 118 Vz *= k; 119/* Calculation of longitude and latitude.*/ 120 lp.lam = atan2 (Vy, Vx); 121 lp.phi = atan (Vz * cos (lp.lam) / Vx); 122 lp.phi = atan (P->radius_p_inv2 * tan (lp.phi)); 123 return (lp); 124} 125FREEUP; if (P) free(P); } 126ENTRY0(geos) 127 if ((P->h = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30); 128 if (P->phi0) E_ERROR(-46); 129 P->radius_g = 1. + (P->radius_g_1 = P->h / P->a); 130 P->C = P->radius_g * P->radius_g - 1.0; 131 if (P->es) { 132 P->radius_p = sqrt (P->one_es); 133 P->radius_p2 = P->one_es; 134 P->radius_p_inv2 = P->rone_es; 135 P->inv = e_inverse; 136 P->fwd = e_forward; 137 } else { 138 P->radius_p = P->radius_p2 = P->radius_p_inv2 = 1.0; 139 P->inv = s_inverse; 140 P->fwd = s_forward; 141 } 142ENDENTRY(P) 143/* 144** $Log: PJ_geos.c,v $ 145** Revision 1.2 2005/02/04 19:27:58 fwarmerdam 146** Added link to reference info. 147** 148** Revision 1.1 2004/10/20 17:04:00 fwarmerdam 149** New 150** 151** Revision 1.2 2004/07/14 18:08:57 gie 152** corrected P->phi_0 to P->phi0 153** 154** Revision 1.1 2004/07/12 17:58:25 gie 155** Initial revision 156** 157*/