/Proj4/PJ_geos.c

  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*/