root / org.gvsig.jcrs / libjni-proj4 / src / PJ_krovak.c @ 36
History | View | Annotate | Download (7.81 KB)
1 |
/******************************************************************************
|
---|---|
2 |
* $Id: PJ_krovak.c,v 1.4 2002/12/15 22:31:04 warmerda Exp $
|
3 |
*
|
4 |
* Project: PROJ.4
|
5 |
* Purpose: Implementation of the krovak (Krovak) projection.
|
6 |
* Definition: http://www.ihsenergy.com/epsg/guid7.html#1.4.3
|
7 |
* Author: Thomas Flemming, tf@ttqv.com
|
8 |
*
|
9 |
******************************************************************************
|
10 |
* Copyright (c) 2001, Thomas Flemming, tf@ttqv.com
|
11 |
*
|
12 |
* Permission is hereby granted, free of charge, to any person obtaining a
|
13 |
* copy of this software and associated documentation files (the "Software"),
|
14 |
* to deal in the Software without restriction, including without limitation
|
15 |
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
|
16 |
* and/or sell copies of the Software, and to permit persons to whom the
|
17 |
* Software is furnished to do so, subject to the following conditions:
|
18 |
*
|
19 |
* The above copyright notice and this permission notice shall be included
|
20 |
* in all copies or substantial portions of the Software.
|
21 |
*
|
22 |
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
23 |
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
|
24 |
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
25 |
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
|
26 |
* BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
|
27 |
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
|
28 |
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
29 |
* SOFTWARE.
|
30 |
******************************************************************************
|
31 |
*
|
32 |
* $Log: PJ_krovak.c,v $
|
33 |
* Revision 1.4 2002/12/15 22:31:04 warmerda
|
34 |
* handle lon_0, k, and prime meridian properly
|
35 |
*
|
36 |
* Revision 1.3 2002/12/15 00:13:30 warmerda
|
37 |
* lat_0 may now be set by user, but still defaults to 49d30N
|
38 |
*
|
39 |
* Revision 1.2 2002/12/14 19:35:21 warmerda
|
40 |
* updated headers
|
41 |
*
|
42 |
*/
|
43 |
|
44 |
#define PROJ_PARMS__ \
|
45 |
double C_x;
|
46 |
#define PJ_LIB__
|
47 |
|
48 |
#include <projects.h> |
49 |
#include <string.h> |
50 |
#include <stdio.h> |
51 |
|
52 |
PJ_CVSID("$Id: PJ_krovak.c,v 1.4 2002/12/15 22:31:04 warmerda Exp $");
|
53 |
|
54 |
PROJ_HEAD(krovak, "Krovak") "\n\tPCyl., Sph."; |
55 |
|
56 |
/**
|
57 |
NOTES: According to EPSG the full Krovak projection method should have
|
58 |
the following parameters. Within PROJ.4 the azimuth, and pseudo
|
59 |
standard parallel are hardcoded in the algorithm and can't be
|
60 |
altered from outside. The others all have defaults to match the
|
61 |
common usage with Krovak projection.
|
62 |
|
63 |
lat_0 = latitude of centre of the projection
|
64 |
|
65 |
lon_0 = longitude of centre of the projection
|
66 |
|
67 |
** = azimuth (true) of the centre line passing through the centre of the projection
|
68 |
|
69 |
** = latitude of pseudo standard parallel
|
70 |
|
71 |
k = scale factor on the pseudo standard parallel
|
72 |
|
73 |
x_0 = False Easting of the centre of the projection at the apex of the cone
|
74 |
|
75 |
y_0 = False Northing of the centre of the projection at the apex of the cone
|
76 |
|
77 |
**/
|
78 |
|
79 |
|
80 |
|
81 |
FORWARD(s_forward); /* spheroid */
|
82 |
/* calculate xy from lat/lon */
|
83 |
|
84 |
char errmess[255]; |
85 |
char tmp[16]; |
86 |
|
87 |
/* Constants, identical to inverse transform function */
|
88 |
double s45, s90, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
|
89 |
double gfi, u, fi0, lon17, lamdd, deltav, s, d, eps, ro;
|
90 |
|
91 |
|
92 |
s45 = 0.785398163397448; /* 45? */ |
93 |
s90 = 2 * s45;
|
94 |
fi0 = P->phi0; /* Latitude of projection centre 49? 30' */
|
95 |
|
96 |
/* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
|
97 |
be set to 1 here.
|
98 |
Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
|
99 |
e2=0.006674372230614;
|
100 |
*/
|
101 |
a = 1; /* 6377397.155; */ |
102 |
/* e2 = P->es;*/ /* 0.006674372230614; */ |
103 |
e2 = 0.006674372230614; |
104 |
e = sqrt(e2); |
105 |
|
106 |
alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2)); |
107 |
|
108 |
uq = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */ |
109 |
u0 = asin(sin(fi0) / alfa); |
110 |
g = pow( (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2. ); |
111 |
|
112 |
k = tan( u0 / 2. + s45) / pow (tan(fi0 / 2. + s45) , alfa) * g; |
113 |
|
114 |
k1 = P->k0; |
115 |
n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2)); |
116 |
s0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78? 30'00" N */ |
117 |
n = sin(s0); |
118 |
ro0 = k1 * n0 / tan(s0); |
119 |
ad = s90 - uq; |
120 |
|
121 |
/* Transformation */
|
122 |
|
123 |
gfi =pow ( ((1. + e * sin(lp.phi)) /
|
124 |
(1. - e * sin(lp.phi))) , (alfa * e / 2.)); |
125 |
|
126 |
u= 2. * (atan(k * pow( tan(lp.phi / 2. + s45), alfa) / gfi)-s45); |
127 |
|
128 |
deltav = - lp.lam * alfa; |
129 |
|
130 |
s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav)); |
131 |
d = asin(cos(u) * sin(deltav) / cos(s)); |
132 |
eps = n * d; |
133 |
ro = ro0 * pow(tan(s0 / 2. + s45) , n) / pow(tan(s / 2. + s45) , n) ; |
134 |
|
135 |
/* x and y are reverted! */
|
136 |
xy.y = ro * cos(eps) / a; |
137 |
xy.x = ro * sin(eps) / a; |
138 |
|
139 |
#ifdef DEBUG
|
140 |
strcpy(errmess,"a: ");
|
141 |
strcpy(tmp," ");
|
142 |
ltoa((long)(a*1000000000),tmp,10); |
143 |
strcat(errmess,tmp); |
144 |
strcat(errmess,"e2: ");
|
145 |
strcpy(tmp," ");
|
146 |
ltoa((long)(e2*1000000000),tmp,10); |
147 |
strcat(errmess,tmp); |
148 |
|
149 |
MessageBox(NULL, errmess, NULL, 0); |
150 |
#endif
|
151 |
|
152 |
|
153 |
return (xy);
|
154 |
} |
155 |
|
156 |
|
157 |
|
158 |
INVERSE(s_inverse); /* spheroid */
|
159 |
/* calculate lat/lon from xy */
|
160 |
|
161 |
/* Constants, identisch wie in der Umkehrfunktion */
|
162 |
double s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
|
163 |
double u, l24, lamdd, deltav, s, d, eps, ro, fi1, xy0, lon17;
|
164 |
int ok;
|
165 |
|
166 |
s45 = 0.785398163397448; /* 45? */ |
167 |
s90 = 2 * s45;
|
168 |
fi0 = P->phi0; /* Latitude of projection centre 49? 30' */
|
169 |
|
170 |
|
171 |
/* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
|
172 |
be set to 1 here.
|
173 |
Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
|
174 |
e2=0.006674372230614;
|
175 |
*/
|
176 |
a = 1; /* 6377397.155; */ |
177 |
/* e2 = P->es; */ /* 0.006674372230614; */ |
178 |
e2 = 0.006674372230614; |
179 |
e = sqrt(e2); |
180 |
|
181 |
alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2)); |
182 |
uq = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */ |
183 |
u0 = asin(sin(fi0) / alfa); |
184 |
g = pow( (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2. ); |
185 |
|
186 |
k = tan( u0 / 2. + s45) / pow (tan(fi0 / 2. + s45) , alfa) * g; |
187 |
|
188 |
k1 = P->k0; |
189 |
n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2)); |
190 |
s0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78? 30'00" N */ |
191 |
n = sin(s0); |
192 |
ro0 = k1 * n0 / tan(s0); |
193 |
ad = s90 - uq; |
194 |
|
195 |
|
196 |
/* Transformation */
|
197 |
/* revert y, x*/
|
198 |
xy0=xy.x; |
199 |
xy.x=xy.y; |
200 |
xy.y=xy0; |
201 |
|
202 |
ro = sqrt(xy.x * xy.x + xy.y * xy.y); |
203 |
eps = atan2(xy.y, xy.x); |
204 |
d = eps / sin(s0); |
205 |
s = 2. * (atan( pow(ro0 / ro, 1. / n) * tan(s0 / 2. + s45)) - s45); |
206 |
|
207 |
u = asin(cos(ad) * sin(s) - sin(ad) * cos(s) * cos(d)); |
208 |
deltav = asin(cos(s) * sin(d) / cos(u)); |
209 |
|
210 |
lp.lam = P->lam0 - deltav / alfa; |
211 |
|
212 |
/* ITERATION FOR lp.phi */
|
213 |
fi1 = u; |
214 |
|
215 |
ok = 0;
|
216 |
do
|
217 |
{ |
218 |
lp.phi = 2. * ( atan( pow( k, -1. / alfa) * |
219 |
pow( tan(u / 2. + s45) , 1. / alfa) * |
220 |
pow( (1. + e * sin(fi1)) / (1. - e * sin(fi1)) , e / 2.) |
221 |
) - s45); |
222 |
|
223 |
if (fabs(fi1 - lp.phi) < 0.000000000000001) ok=1; |
224 |
fi1 = lp.phi; |
225 |
|
226 |
} |
227 |
while (ok==0); |
228 |
|
229 |
lp.lam -= P->lam0; |
230 |
|
231 |
return (lp);
|
232 |
} |
233 |
|
234 |
FREEUP; if (P) pj_dalloc(P); }
|
235 |
|
236 |
ENTRY0(krovak) |
237 |
double ts;
|
238 |
/* read some Parameters,
|
239 |
* here Latitude Truescale */
|
240 |
|
241 |
ts = pj_param(P->params, "rlat_ts").f;
|
242 |
P->C_x = ts; |
243 |
|
244 |
/* we want Bessel as fixed ellipsoid */
|
245 |
P->a = 6377397.155; |
246 |
P->e = sqrt(P->es = 0.006674372230614); |
247 |
|
248 |
/* if latitude of projection center is not set, use 49d30'N */
|
249 |
if (!pj_param(P->params, "tlat_0").i) |
250 |
P->phi0 = 0.863937979737193; |
251 |
|
252 |
/* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */
|
253 |
/* that will correspond to using longitudes relative to greenwich */
|
254 |
/* as input and output, instead of lat/long relative to Ferro */
|
255 |
if (!pj_param(P->params, "tlon_0").i) |
256 |
P->lam0 = 0.7417649320975901 - 0.308341501185665; |
257 |
; |
258 |
|
259 |
/* if scale not set default to 0.9999 */
|
260 |
if (!pj_param(P->params, "tk").i) |
261 |
P->k0 = 0.9999; |
262 |
|
263 |
/* always the same */
|
264 |
P->inv = s_inverse; |
265 |
P->fwd = s_forward; |
266 |
|
267 |
ENDENTRY(P) |
268 |
|