2
2
** libproj -- library of cartographic projections
4
4
** Copyright (c) 2004 Gerald I. Evenden
5
** Copyright (c) 2012 Martin Raspaud
7
LIBPROJ_ID[] = "$Id: PJ_geos.c 1504 2009-01-06 02:11:57Z warmerdam $";
8
LIBPROJ_ID[] = "$Id: PJ_geos.c 2176 2012-02-27 07:56:32Z warmerdam $";
9
10
** See also (section 4.4.3.2):
10
11
** http://www.eumetsat.int/en/area4/msg/news/us_doc/cgms_03_26.pdf
54
57
if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz) < 0.) F_ERROR;
55
58
/* Calculation based on view angles from satellite.*/
56
59
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));
62
xy.x = P->radius_g_1 * atan(Vy / hypot(Vz, tmp));
63
xy.y = P->radius_g_1 * atan(Vz / tmp);
67
xy.x = P->radius_g_1 * atan(Vy / tmp);
68
xy.y = P->radius_g_1 * atan(Vz / hypot(Vy, tmp));
61
72
FORWARD(e_forward); /* ellipsoid */
75
86
/* Calculation based on view angles from satellite. */
76
87
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));
90
xy.x = P->radius_g_1 * atan (Vy / hypot (Vz, tmp));
91
xy.y = P->radius_g_1 * atan (Vz / tmp);
95
xy.x = P->radius_g_1 * atan (Vy / tmp);
96
xy.y = P->radius_g_1 * atan (Vz / hypot (Vy, tmp));
81
100
INVERSE(s_inverse); /* spheroid */
82
double Vx, Vy, Vz, a, b, c, det, k;
101
double Vx, Vy, Vz, a, b, det, k;
84
103
/* Setting three components of vector from satellite to position.*/
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);
107
Vz = tan (xy.y / (P->radius_g - 1.0));
108
Vy = tan (xy.x / (P->radius_g - 1.0)) * sqrt (1.0 + Vz * Vz);
112
Vy = tan (xy.x / (P->radius_g - 1.0));
113
Vz = tan (xy.y / (P->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy);
88
115
/* Calculation of terms in cubic equation and determinant.*/
89
116
a = Vy * Vy + Vz * Vz + Vx * Vx;
90
117
b = 2 * P->radius_g * Vx;
102
129
INVERSE(e_inverse); /* ellipsoid */
103
double Vx, Vy, Vz, a, b, c, det, k;
130
double Vx, Vy, Vz, a, b, det, k;
105
132
/* Setting three components of vector from satellite to position.*/
107
Vy = tan (xy.x / P->radius_g_1);
108
Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy);
136
Vz = tan (xy.y / P->radius_g_1);
137
Vy = tan (xy.x / P->radius_g_1) * hypot(1.0, Vz);
141
Vy = tan (xy.x / P->radius_g_1);
142
Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy);
109
144
/* Calculation of terms in cubic equation and determinant.*/
110
145
a = Vz / P->radius_p;
111
146
a = Vy * Vy + a * a + Vx * Vx;
125
160
FREEUP; if (P) free(P); }
127
if ((P->h = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30);
162
if ((P->h = pj_param(P->ctx, P->params, "dh").f) <= 0.) E_ERROR(-30);
128
163
if (P->phi0) E_ERROR(-46);
129
P->radius_g = 1. + (P->radius_g_1 = P->h / P->a);
164
P->sweep_axis = pj_param(P->ctx, P->params, "ssweep").s;
165
if (P->sweep_axis == NULL)
169
if (P->sweep_axis[1] != '\0' ||
170
(P->sweep_axis[0] != 'x' &&
171
P->sweep_axis[0] != 'y'))
173
if (P->sweep_axis[0] == 'y')
178
P->radius_g_1 = P->h / P->a;
179
P->radius_g = 1. + P->radius_g_1;
130
180
P->C = P->radius_g * P->radius_g - 1.0;
132
182
P->radius_p = sqrt (P->one_es);