root / org.gvsig.jcrs / libjni-proj4 / src / PJ_geos.c @ 36
History | View | Annotate | Download (4.94 KB)
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.1 2004/10/20 17:04:00 fwarmerdam Exp $";
|
8 |
/*
|
9 |
** Permission is hereby granted, free of charge, to any person obtaining
|
10 |
** a copy of this software and associated documentation files (the
|
11 |
** "Software"), to deal in the Software without restriction, including
|
12 |
** without limitation the rights to use, copy, modify, merge, publish,
|
13 |
** distribute, sublicense, and/or sell copies of the Software, and to
|
14 |
** permit persons to whom the Software is furnished to do so, subject to
|
15 |
** the following conditions:
|
16 |
**
|
17 |
** The above copyright notice and this permission notice shall be
|
18 |
** included in all copies or substantial portions of the Software.
|
19 |
**
|
20 |
** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
21 |
** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
|
22 |
** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
|
23 |
** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
|
24 |
** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
|
25 |
** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
|
26 |
** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
|
27 |
*/
|
28 |
#define PROJ_PARMS__ \
|
29 |
double h; \
|
30 |
double radius_p; \
|
31 |
double radius_p2; \
|
32 |
double radius_p_inv2; \
|
33 |
double radius_g; \
|
34 |
double radius_g_1; \
|
35 |
double C;
|
36 |
#define PJ_LIB__
|
37 |
#include <projects.h> |
38 |
|
39 |
PROJ_HEAD(geos, "Geostationary Satellite View") "\n\tAzi, Sph&Ell\n\th="; |
40 |
|
41 |
FORWARD(s_forward); /* spheroid */
|
42 |
double Vx, Vy, Vz, tmp;
|
43 |
|
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 |
|
61 |
/* Calculation of geocentric latitude. */
|
62 |
lp.phi = atan (P->radius_p2 * tan (lp.phi)); |
63 |
/* Calculation of the three components of the vector from satellite to
|
64 |
** position on earth surface (lon,lat).*/
|
65 |
r = (P->radius_p) / hypot(P->radius_p * cos (lp.phi), sin (lp.phi)); |
66 |
Vx = r * cos (lp.lam) * cos (lp.phi); |
67 |
Vy = r * sin (lp.lam) * cos (lp.phi); |
68 |
Vz = r * sin (lp.phi); |
69 |
/* Check visibility. */
|
70 |
if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * P->radius_p_inv2) < 0.) |
71 |
F_ERROR; |
72 |
/* Calculation based on view angles from satellite. */
|
73 |
tmp = P->radius_g - Vx; |
74 |
xy.x = P->radius_g_1 * atan (Vy / tmp); |
75 |
xy.y = P->radius_g_1 * atan (Vz / hypot (Vy, tmp)); |
76 |
return (xy);
|
77 |
} |
78 |
INVERSE(s_inverse); /* spheroid */
|
79 |
double Vx, Vy, Vz, a, b, c, det, k;
|
80 |
|
81 |
/* Setting three components of vector from satellite to position.*/
|
82 |
Vx = -1.0; |
83 |
Vy = tan (xy.x / (P->radius_g - 1.0)); |
84 |
Vz = tan (xy.y / (P->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy); |
85 |
/* Calculation of terms in cubic equation and determinant.*/
|
86 |
a = Vy * Vy + Vz * Vz + Vx * Vx; |
87 |
b = 2 * P->radius_g * Vx;
|
88 |
if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR; |
89 |
/* Calculation of three components of vector from satellite to position.*/
|
90 |
k = (-b - sqrt(det)) / (2 * a);
|
91 |
Vx = P->radius_g + k * Vx; |
92 |
Vy *= k; |
93 |
Vz *= k; |
94 |
/* Calculation of longitude and latitude.*/
|
95 |
lp.lam = atan2 (Vy, Vx); |
96 |
lp.phi = atan (Vz * cos (lp.lam) / Vx); |
97 |
return (lp);
|
98 |
} |
99 |
INVERSE(e_inverse); /* ellipsoid */
|
100 |
double Vx, Vy, Vz, a, b, c, det, k;
|
101 |
|
102 |
/* Setting three components of vector from satellite to position.*/
|
103 |
Vx = -1.0; |
104 |
Vy = tan (xy.x / P->radius_g_1); |
105 |
Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy); |
106 |
/* Calculation of terms in cubic equation and determinant.*/
|
107 |
a = Vz / P->radius_p; |
108 |
a = Vy * Vy + a * a + Vx * Vx; |
109 |
b = 2 * P->radius_g * Vx;
|
110 |
if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR; |
111 |
/* Calculation of three components of vector from satellite to position.*/
|
112 |
k = (-b - sqrt(det)) / (2. * a);
|
113 |
Vx = P->radius_g + k * Vx; |
114 |
Vy *= k; |
115 |
Vz *= k; |
116 |
/* Calculation of longitude and latitude.*/
|
117 |
lp.lam = atan2 (Vy, Vx); |
118 |
lp.phi = atan (Vz * cos (lp.lam) / Vx); |
119 |
lp.phi = atan (P->radius_p_inv2 * tan (lp.phi)); |
120 |
return (lp);
|
121 |
} |
122 |
FREEUP; if (P) free(P); }
|
123 |
ENTRY0(geos) |
124 |
if ((P->h = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30); |
125 |
if (P->phi0) E_ERROR(-46); |
126 |
P->radius_g = 1. + (P->radius_g_1 = P->h / P->a);
|
127 |
P->C = P->radius_g * P->radius_g - 1.0; |
128 |
if (P->es) {
|
129 |
P->radius_p = sqrt (P->one_es); |
130 |
P->radius_p2 = P->one_es; |
131 |
P->radius_p_inv2 = P->rone_es; |
132 |
P->inv = e_inverse; |
133 |
P->fwd = e_forward; |
134 |
} else {
|
135 |
P->radius_p = P->radius_p2 = P->radius_p_inv2 = 1.0; |
136 |
P->inv = s_inverse; |
137 |
P->fwd = s_forward; |
138 |
} |
139 |
ENDENTRY(P) |
140 |
/*
|
141 |
** $Log: PJ_geos.c,v $
|
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 |
*/
|