/Proj4/PJ_geos.c
C | 157 lines | 91 code | 6 blank | 60 comment | 9 complexity | 8fbd37c899538300c5336885d7f6dd1b MD5 | raw file
- /*
- ** libproj -- library of cartographic projections
- **
- ** Copyright (c) 2004 Gerald I. Evenden
- */
- static const char
- LIBPROJ_ID[] = "$Id: PJ_geos.c,v 1.2 2005/02/04 19:27:58 fwarmerdam Exp $";
- /*
- ** See also (section 4.4.3.2):
- ** http://www.eumetsat.int/en/area4/msg/news/us_doc/cgms_03_26.pdf
- **
- ** Permission is hereby granted, free of charge, to any person obtaining
- ** a copy of this software and associated documentation files (the
- ** "Software"), to deal in the Software without restriction, including
- ** without limitation the rights to use, copy, modify, merge, publish,
- ** distribute, sublicense, and/or sell copies of the Software, and to
- ** permit persons to whom the Software is furnished to do so, subject to
- ** the following conditions:
- **
- ** The above copyright notice and this permission notice shall be
- ** included in all copies or substantial portions of the Software.
- **
- ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
- ** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- ** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- ** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
- ** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
- ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
- ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
- */
- #define PROJ_PARMS__ \
- double h; \
- double radius_p; \
- double radius_p2; \
- double radius_p_inv2; \
- double radius_g; \
- double radius_g_1; \
- double C;
- #define PJ_LIB__
- #include "projects.h"
- PROJ_HEAD(geos, "Geostationary Satellite View") "\n\tAzi, Sph&Ell\n\th=";
- FORWARD(s_forward); /* spheroid */
- double Vx, Vy, Vz, tmp;
- /* Calculation of the three components of the vector from satellite to
- ** position on earth surface (lon,lat).*/
- tmp = cos(lp.phi);
- Vx = cos (lp.lam) * tmp;
- Vy = sin (lp.lam) * tmp;
- Vz = sin (lp.phi);
- /* Check visibility.*/
- if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz) < 0.) F_ERROR;
- /* Calculation based on view angles from satellite.*/
- tmp = P->radius_g - Vx;
- xy.x = P->radius_g_1 * atan(Vy / tmp);
- xy.y = P->radius_g_1 * atan(Vz / hypot(Vy, tmp));
- return (xy);
- }
- FORWARD(e_forward); /* ellipsoid */
- double r, Vx, Vy, Vz, tmp;
- /* Calculation of geocentric latitude. */
- lp.phi = atan (P->radius_p2 * tan (lp.phi));
- /* Calculation of the three components of the vector from satellite to
- ** position on earth surface (lon,lat).*/
- r = (P->radius_p) / hypot(P->radius_p * cos (lp.phi), sin (lp.phi));
- Vx = r * cos (lp.lam) * cos (lp.phi);
- Vy = r * sin (lp.lam) * cos (lp.phi);
- Vz = r * sin (lp.phi);
- /* Check visibility. */
- if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * P->radius_p_inv2) < 0.)
- F_ERROR;
- /* Calculation based on view angles from satellite. */
- tmp = P->radius_g - Vx;
- xy.x = P->radius_g_1 * atan (Vy / tmp);
- xy.y = P->radius_g_1 * atan (Vz / hypot (Vy, tmp));
- return (xy);
- }
- INVERSE(s_inverse); /* spheroid */
- double Vx, Vy, Vz, a, b, c, det, k;
- /* Setting three components of vector from satellite to position.*/
- Vx = -1.0;
- Vy = tan (xy.x / (P->radius_g - 1.0));
- Vz = tan (xy.y / (P->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy);
- /* Calculation of terms in cubic equation and determinant.*/
- a = Vy * Vy + Vz * Vz + Vx * Vx;
- b = 2 * P->radius_g * Vx;
- if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR;
- /* Calculation of three components of vector from satellite to position.*/
- k = (-b - sqrt(det)) / (2 * a);
- Vx = P->radius_g + k * Vx;
- Vy *= k;
- Vz *= k;
- /* Calculation of longitude and latitude.*/
- lp.lam = atan2 (Vy, Vx);
- lp.phi = atan (Vz * cos (lp.lam) / Vx);
- return (lp);
- }
- INVERSE(e_inverse); /* ellipsoid */
- double Vx, Vy, Vz, a, b, c, det, k;
- /* Setting three components of vector from satellite to position.*/
- Vx = -1.0;
- Vy = tan (xy.x / P->radius_g_1);
- Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy);
- /* Calculation of terms in cubic equation and determinant.*/
- a = Vz / P->radius_p;
- a = Vy * Vy + a * a + Vx * Vx;
- b = 2 * P->radius_g * Vx;
- if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR;
- /* Calculation of three components of vector from satellite to position.*/
- k = (-b - sqrt(det)) / (2. * a);
- Vx = P->radius_g + k * Vx;
- Vy *= k;
- Vz *= k;
- /* Calculation of longitude and latitude.*/
- lp.lam = atan2 (Vy, Vx);
- lp.phi = atan (Vz * cos (lp.lam) / Vx);
- lp.phi = atan (P->radius_p_inv2 * tan (lp.phi));
- return (lp);
- }
- FREEUP; if (P) free(P); }
- ENTRY0(geos)
- if ((P->h = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30);
- if (P->phi0) E_ERROR(-46);
- P->radius_g = 1. + (P->radius_g_1 = P->h / P->a);
- P->C = P->radius_g * P->radius_g - 1.0;
- if (P->es) {
- P->radius_p = sqrt (P->one_es);
- P->radius_p2 = P->one_es;
- P->radius_p_inv2 = P->rone_es;
- P->inv = e_inverse;
- P->fwd = e_forward;
- } else {
- P->radius_p = P->radius_p2 = P->radius_p_inv2 = 1.0;
- P->inv = s_inverse;
- P->fwd = s_forward;
- }
- ENDENTRY(P)
- /*
- ** $Log: PJ_geos.c,v $
- ** Revision 1.2 2005/02/04 19:27:58 fwarmerdam
- ** Added link to reference info.
- **
- ** Revision 1.1 2004/10/20 17:04:00 fwarmerdam
- ** New
- **
- ** Revision 1.2 2004/07/14 18:08:57 gie
- ** corrected P->phi_0 to P->phi0
- **
- ** Revision 1.1 2004/07/12 17:58:25 gie
- ** Initial revision
- **
- */