svn-gvsig-desktop / trunk / libraries / libjni-proj4 / src / PJ_mod_ster.c @ 7098
History | View | Annotate | Download (4.8 KB)
1 |
#ifndef lint
|
---|---|
2 |
static const char SCCSID[]="@(#)PJ_mod_ster.c 4.1 94/02/15 GIE REL"; |
3 |
#endif
|
4 |
/* based upon Snyder and Linck, USGS-NMD */
|
5 |
#define PROJ_PARMS__ \
|
6 |
COMPLEX *zcoeff; \ |
7 |
double cchio, schio; \
|
8 |
int n;
|
9 |
#define PJ_LIB__
|
10 |
#include <projects.h> |
11 |
PROJ_HEAD(mil_os, "Miller Oblated Stereographic") "\n\tAzi(mod)"; |
12 |
PROJ_HEAD(lee_os, "Lee Oblated Stereographic") "\n\tAzi(mod)"; |
13 |
PROJ_HEAD(gs48, "Mod. Stererographics of 48 U.S.") "\n\tAzi(mod)"; |
14 |
PROJ_HEAD(alsk, "Mod. Stererographics of Alaska") "\n\tAzi(mod)"; |
15 |
PROJ_HEAD(gs50, "Mod. Stererographics of 50 U.S.") "\n\tAzi(mod)"; |
16 |
#define EPSLN 1e-10 |
17 |
|
18 |
FORWARD(e_forward); /* ellipsoid */
|
19 |
double sinlon, coslon, esphi, chi, schi, cchi, s;
|
20 |
COMPLEX p; |
21 |
|
22 |
sinlon = sin(lp.lam); |
23 |
coslon = cos(lp.lam); |
24 |
esphi = P->e * sin(lp.phi); |
25 |
chi = 2. * atan(tan((HALFPI + lp.phi) * .5) * |
26 |
pow((1. - esphi) / (1. + esphi), P->e * .5)) - HALFPI; |
27 |
schi = sin(chi); |
28 |
cchi = cos(chi); |
29 |
s = 2. / (1. + P->schio * schi + P->cchio * cchi * coslon); |
30 |
p.r = s * cchi * sinlon; |
31 |
p.i = s * (P->cchio * schi - P->schio * cchi * coslon); |
32 |
p = pj_zpoly1(p, P->zcoeff, P->n); |
33 |
xy.x = p.r; |
34 |
xy.y = p.i; |
35 |
return xy;
|
36 |
} |
37 |
INVERSE(e_inverse); /* ellipsoid */
|
38 |
int nn;
|
39 |
COMPLEX p, fxy, fpxy, dp; |
40 |
double den, rh, z, sinz, cosz, chi, phi, dphi, esphi;
|
41 |
|
42 |
p.r = xy.x; |
43 |
p.i = xy.y; |
44 |
for (nn = 20; nn ;--nn) { |
45 |
fxy = pj_zpolyd1(p, P->zcoeff, P->n, &fpxy); |
46 |
fxy.r -= xy.x; |
47 |
fxy.i -= xy.y; |
48 |
den = fpxy.r * fpxy.r + fpxy.i * fpxy.i; |
49 |
dp.r = -(fxy.r * fpxy.r + fxy.i * fpxy.i) / den; |
50 |
dp.i = -(fxy.i * fpxy.r - fxy.r * fpxy.i) / den; |
51 |
p.r += dp.r; |
52 |
p.i += dp.i; |
53 |
if ((fabs(dp.r) + fabs(dp.i)) <= EPSLN)
|
54 |
break;
|
55 |
} |
56 |
if (nn) {
|
57 |
rh = hypot(p.r, p.i); |
58 |
z = 2. * atan(.5 * rh); |
59 |
sinz = sin(z); |
60 |
cosz = cos(z); |
61 |
lp.lam = P->lam0; |
62 |
if (fabs(rh) <= EPSLN) {
|
63 |
lp.phi = P->phi0; |
64 |
return lp;
|
65 |
} |
66 |
chi = aasin(cosz * P->schio + p.i * sinz * P->cchio / rh); |
67 |
phi = chi; |
68 |
for (nn = 20; nn ;--nn) { |
69 |
esphi = P->e * sin(phi); |
70 |
dphi = 2. * atan(tan((HALFPI + chi) * .5) * |
71 |
pow((1. + esphi) / (1. - esphi), P->e * .5)) - HALFPI - phi; |
72 |
phi += dphi; |
73 |
if (fabs(dphi) <= EPSLN)
|
74 |
break;
|
75 |
} |
76 |
} |
77 |
if (nn) {
|
78 |
lp.phi = phi; |
79 |
lp.lam = atan2(p.r * sinz, rh * P->cchio * cosz - p.i * |
80 |
P->schio * sinz); |
81 |
} else
|
82 |
lp.lam = lp.phi = HUGE_VAL; |
83 |
return lp;
|
84 |
} |
85 |
FREEUP; if (P) pj_dalloc(P); }
|
86 |
static PJ *
|
87 |
setup(PJ *P) { /* general initialization */
|
88 |
double esphi, chio;
|
89 |
|
90 |
if (P->es) {
|
91 |
esphi = P->e * sin(P->phi0); |
92 |
chio = 2. * atan(tan((HALFPI + P->phi0) * .5) * |
93 |
pow((1. - esphi) / (1. + esphi), P->e * .5)) - HALFPI; |
94 |
} else
|
95 |
chio = P->phi0; |
96 |
P->schio = sin(chio); |
97 |
P->cchio = cos(chio); |
98 |
P->inv = e_inverse; P->fwd = e_forward; |
99 |
return P;
|
100 |
} |
101 |
ENTRY0(mil_os) |
102 |
static COMPLEX /* Miller Oblated Stereographic */ |
103 |
AB[] = { |
104 |
{0.924500, 0.}, |
105 |
{0., 0.}, |
106 |
{0.019430, 0.} |
107 |
}; |
108 |
|
109 |
P->n = 2;
|
110 |
P->lam0 = DEG_TO_RAD * 20.; |
111 |
P->phi0 = DEG_TO_RAD * 18.; |
112 |
P->zcoeff = AB; |
113 |
P->es = 0.;
|
114 |
ENDENTRY(setup(P)) |
115 |
ENTRY0(lee_os) |
116 |
static COMPLEX /* Lee Oblated Stereographic */ |
117 |
AB[] = { |
118 |
{0.721316, 0.}, |
119 |
{0., 0.}, |
120 |
{-0.0088162, -0.00617325} |
121 |
}; |
122 |
|
123 |
P->n = 2;
|
124 |
P->lam0 = DEG_TO_RAD * -165.; |
125 |
P->phi0 = DEG_TO_RAD * -10.; |
126 |
P->zcoeff = AB; |
127 |
P->es = 0.;
|
128 |
ENDENTRY(setup(P)) |
129 |
ENTRY0(gs48) |
130 |
static COMPLEX /* 48 United States */ |
131 |
AB[] = { |
132 |
{0.98879, 0.}, |
133 |
{0., 0.}, |
134 |
{-0.050909, 0.}, |
135 |
{0., 0.}, |
136 |
{0.075528, 0.} |
137 |
}; |
138 |
|
139 |
P->n = 4;
|
140 |
P->lam0 = DEG_TO_RAD * -96.; |
141 |
P->phi0 = DEG_TO_RAD * -39.; |
142 |
P->zcoeff = AB; |
143 |
P->es = 0.;
|
144 |
P->a = 6370997.; |
145 |
ENDENTRY(setup(P)) |
146 |
ENTRY0(alsk) |
147 |
static COMPLEX
|
148 |
ABe[] = { /* Alaska ellipsoid */
|
149 |
{.9945303, 0.}, |
150 |
{.0052083, -.0027404}, |
151 |
{.0072721, .0048181}, |
152 |
{-.0151089, -.1932526}, |
153 |
{.0642675, -.1381226}, |
154 |
{.3582802, -.2884586}}, |
155 |
ABs[] = { /* Alaska sphere */
|
156 |
{.9972523, 0.}, |
157 |
{.0052513, -.0041175}, |
158 |
{.0074606, .0048125}, |
159 |
{-.0153783, -.1968253}, |
160 |
{.0636871, -.1408027}, |
161 |
{.3660976, -.2937382} |
162 |
}; |
163 |
|
164 |
P->n = 5;
|
165 |
P->lam0 = DEG_TO_RAD * -152.; |
166 |
P->phi0 = DEG_TO_RAD * 64.; |
167 |
if (P->es) { /* fixed ellipsoid/sphere */ |
168 |
P->zcoeff = ABe; |
169 |
P->a = 6378206.4; |
170 |
P->e = sqrt(P->es = 0.00676866); |
171 |
} else {
|
172 |
P->zcoeff = ABs; |
173 |
P->a = 6370997.; |
174 |
} |
175 |
ENDENTRY(setup(P)) |
176 |
ENTRY0(gs50) |
177 |
static COMPLEX
|
178 |
ABe[] = { /* GS50 ellipsoid */
|
179 |
{.9827497, 0.}, |
180 |
{.0210669, .0053804}, |
181 |
{-.1031415, -.0571664}, |
182 |
{-.0323337, -.0322847}, |
183 |
{.0502303, .1211983}, |
184 |
{.0251805, .0895678}, |
185 |
{-.0012315, -.1416121}, |
186 |
{.0072202, -.1317091}, |
187 |
{-.0194029, .0759677}, |
188 |
{-.0210072, .0834037} |
189 |
}, |
190 |
ABs[] = { /* GS50 sphere */
|
191 |
{.9842990, 0.}, |
192 |
{.0211642, .0037608}, |
193 |
{-.1036018, -.0575102}, |
194 |
{-.0329095, -.0320119}, |
195 |
{.0499471, .1223335}, |
196 |
{.0260460, .0899805}, |
197 |
{.0007388, -.1435792}, |
198 |
{.0075848, -.1334108}, |
199 |
{-.0216473, .0776645}, |
200 |
{-.0225161, .0853673} |
201 |
}; |
202 |
|
203 |
P->n = 9;
|
204 |
P->lam0 = DEG_TO_RAD * -120.; |
205 |
P->phi0 = DEG_TO_RAD * 45.; |
206 |
if (P->es) { /* fixed ellipsoid/sphere */ |
207 |
P->zcoeff = ABe; |
208 |
P->a = 6378206.4; |
209 |
P->e = sqrt(P->es = 0.00676866); |
210 |
} else {
|
211 |
P->zcoeff = ABs; |
212 |
P->a = 6370997.; |
213 |
} |
214 |
ENDENTRY(setup(P)) |