1 #define PJ_LIB__
2 #include <projects.h>
3
4 PROJ_HEAD(stere, "Stereographic") "\n\tAzi, Sph&Ell\n\tlat_ts=";
5 PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Sph&Ell\n\tsouth";
6
7
8 struct pj_opaque {
9 double phits;
10 double sinX1;
11 double cosX1;
12 double akm1;
13 int mode;
14 };
15
16 #define sinph0 P->opaque->sinX1
17 #define cosph0 P->opaque->cosX1
18 #define EPS10 1.e-10
19 #define TOL 1.e-8
20 #define NITER 8
21 #define CONV 1.e-10
22 #define S_POLE 0
23 #define N_POLE 1
24 #define OBLIQ 2
25 #define EQUIT 3
26
ssfn_(double phit,double sinphi,double eccen)27 static double ssfn_ (double phit, double sinphi, double eccen) {
28 sinphi *= eccen;
29 return (tan (.5 * (M_HALFPI + phit)) *
30 pow ((1. - sinphi) / (1. + sinphi), .5 * eccen));
31 }
32
33
e_forward(LP lp,PJ * P)34 static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */
35 XY xy = {0.0,0.0};
36 struct pj_opaque *Q = P->opaque;
37 double coslam, sinlam, sinX = 0.0, cosX = 0.0, X, A, sinphi;
38
39 coslam = cos (lp.lam);
40 sinlam = sin (lp.lam);
41 sinphi = sin (lp.phi);
42 if (Q->mode == OBLIQ || Q->mode == EQUIT) {
43 sinX = sin (X = 2. * atan(ssfn_(lp.phi, sinphi, P->e)) - M_HALFPI);
44 cosX = cos (X);
45 }
46
47 switch (Q->mode) {
48 case OBLIQ:
49 A = Q->akm1 / (Q->cosX1 * (1. + Q->sinX1 * sinX +
50 Q->cosX1 * cosX * coslam));
51 xy.y = A * (Q->cosX1 * sinX - Q->sinX1 * cosX * coslam);
52 goto xmul; /* but why not just xy.x = A * cosX; break; ? */
53
54 case EQUIT:
55 A = 2. * Q->akm1 / (1. + cosX * coslam);
56 xy.y = A * sinX;
57 xmul:
58 xy.x = A * cosX;
59 break;
60
61 case S_POLE:
62 lp.phi = -lp.phi;
63 coslam = - coslam;
64 sinphi = -sinphi;
65 case N_POLE:
66 xy.x = Q->akm1 * pj_tsfn (lp.phi, sinphi, P->e);
67 xy.y = - xy.x * coslam;
68 break;
69 }
70
71 xy.x = xy.x * sinlam;
72 return xy;
73 }
74
75
s_forward(LP lp,PJ * P)76 static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
77 XY xy = {0.0,0.0};
78 struct pj_opaque *Q = P->opaque;
79 double sinphi, cosphi, coslam, sinlam;
80
81 sinphi = sin(lp.phi);
82 cosphi = cos(lp.phi);
83 coslam = cos(lp.lam);
84 sinlam = sin(lp.lam);
85
86 switch (Q->mode) {
87 case EQUIT:
88 xy.y = 1. + cosphi * coslam;
89 goto oblcon;
90 case OBLIQ:
91 xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;
92 oblcon:
93 if (xy.y <= EPS10) F_ERROR;
94 xy.x = (xy.y = Q->akm1 / xy.y) * cosphi * sinlam;
95 xy.y *= (Q->mode == EQUIT) ? sinphi :
96 cosph0 * sinphi - sinph0 * cosphi * coslam;
97 break;
98 case N_POLE:
99 coslam = - coslam;
100 lp.phi = - lp.phi;
101 case S_POLE:
102 if (fabs (lp.phi - M_HALFPI) < TOL) F_ERROR;
103 xy.x = sinlam * ( xy.y = Q->akm1 * tan (M_FORTPI + .5 * lp.phi) );
104 xy.y *= coslam;
105 break;
106 }
107 return xy;
108 }
109
110
e_inverse(XY xy,PJ * P)111 static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */
112 LP lp = {0.0,0.0};
113 struct pj_opaque *Q = P->opaque;
114 double cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0;
115 int i;
116
117 rho = hypot (xy.x, xy.y);
118
119 switch (Q->mode) {
120 case OBLIQ:
121 case EQUIT:
122 cosphi = cos ( tp = 2. * atan2 (rho * Q->cosX1 , Q->akm1) );
123 sinphi = sin (tp);
124 if ( rho == 0.0 )
125 phi_l = asin (cosphi * Q->sinX1);
126 else
127 phi_l = asin (cosphi * Q->sinX1 + (xy.y * sinphi * Q->cosX1 / rho));
128
129 tp = tan (.5 * (M_HALFPI + phi_l));
130 xy.x *= sinphi;
131 xy.y = rho * Q->cosX1 * cosphi - xy.y * Q->sinX1* sinphi;
132 halfpi = M_HALFPI;
133 halfe = .5 * P->e;
134 break;
135 case N_POLE:
136 xy.y = -xy.y;
137 case S_POLE:
138 phi_l = M_HALFPI - 2. * atan (tp = - rho / Q->akm1);
139 halfpi = -M_HALFPI;
140 halfe = -.5 * P->e;
141 break;
142 }
143
144 for (i = NITER; i--; phi_l = lp.phi) {
145 sinphi = P->e * sin(phi_l);
146 lp.phi = 2. * atan (tp * pow ((1.+sinphi)/(1.-sinphi), halfe)) - halfpi;
147 if (fabs (phi_l - lp.phi) < CONV) {
148 if (Q->mode == S_POLE)
149 lp.phi = -lp.phi;
150 lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2 (xy.x, xy.y);
151 return lp;
152 }
153 }
154 I_ERROR;
155 }
156
157
s_inverse(XY xy,PJ * P)158 static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
159 LP lp = {0.0,0.0};
160 struct pj_opaque *Q = P->opaque;
161 double c, rh, sinc, cosc;
162
163 sinc = sin (c = 2. * atan ((rh = hypot (xy.x, xy.y)) / Q->akm1));
164 cosc = cos (c);
165 lp.lam = 0.;
166
167 switch (Q->mode) {
168 case EQUIT:
169 if (fabs (rh) <= EPS10)
170 lp.phi = 0.;
171 else
172 lp.phi = asin (xy.y * sinc / rh);
173 if (cosc != 0. || xy.x != 0.)
174 lp.lam = atan2 (xy.x * sinc, cosc * rh);
175 break;
176 case OBLIQ:
177 if (fabs (rh) <= EPS10)
178 lp.phi = P->phi0;
179 else
180 lp.phi = asin (cosc * sinph0 + xy.y * sinc * cosph0 / rh);
181 if ((c = cosc - sinph0 * sin (lp.phi)) != 0. || xy.x != 0.)
182 lp.lam = atan2 (xy.x * sinc * cosph0, c * rh);
183 break;
184 case N_POLE:
185 xy.y = -xy.y;
186 case S_POLE:
187 if (fabs (rh) <= EPS10)
188 lp.phi = P->phi0;
189 else
190 lp.phi = asin (Q->mode == S_POLE ? - cosc : cosc);
191 lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2 (xy.x, xy.y);
192 break;
193 }
194 return lp;
195 }
196
197
freeup_new(PJ * P)198 static void *freeup_new (PJ *P) { /* Destructor */
199 if (0==P)
200 return 0;
201 pj_dealloc (P->opaque);
202 return pj_dealloc(P);
203 }
204
freeup(PJ * P)205 static void freeup (PJ *P) {
206 freeup_new (P);
207 return;
208 }
209
210
setup(PJ * P)211 static PJ *setup(PJ *P) { /* general initialization */
212 double t;
213 struct pj_opaque *Q = P->opaque;
214
215 if (fabs ((t = fabs (P->phi0)) - M_HALFPI) < EPS10)
216 Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
217 else
218 Q->mode = t > EPS10 ? OBLIQ : EQUIT;
219 Q->phits = fabs (Q->phits);
220
221 if (P->es) {
222 double X;
223
224 switch (Q->mode) {
225 case N_POLE:
226 case S_POLE:
227 if (fabs (Q->phits - M_HALFPI) < EPS10)
228 Q->akm1 = 2. * P->k0 /
229 sqrt (pow (1+P->e,1+P->e) * pow (1-P->e,1-P->e));
230 else {
231 Q->akm1 = cos (Q->phits) /
232 pj_tsfn (Q->phits, t = sin (Q->phits), P->e);
233 t *= P->e;
234 Q->akm1 /= sqrt(1. - t * t);
235 }
236 break;
237 case EQUIT:
238 case OBLIQ:
239 t = sin (P->phi0);
240 X = 2. * atan (ssfn_(P->phi0, t, P->e)) - M_HALFPI;
241 t *= P->e;
242 Q->akm1 = 2. * P->k0 * cos (P->phi0) / sqrt(1. - t * t);
243 Q->sinX1 = sin (X);
244 Q->cosX1 = cos (X);
245 break;
246 }
247 P->inv = e_inverse;
248 P->fwd = e_forward;
249 } else {
250 switch (Q->mode) {
251 case OBLIQ:
252 sinph0 = sin (P->phi0);
253 cosph0 = cos (P->phi0);
254 case EQUIT:
255 Q->akm1 = 2. * P->k0;
256 break;
257 case S_POLE:
258 case N_POLE:
259 Q->akm1 = fabs (Q->phits - M_HALFPI) >= EPS10 ?
260 cos (Q->phits) / tan (M_FORTPI - .5 * Q->phits) :
261 2. * P->k0 ;
262 break;
263 }
264
265 P->inv = s_inverse;
266 P->fwd = s_forward;
267 }
268 return P;
269 }
270
271
PROJECTION(stere)272 PJ *PROJECTION(stere) {
273 struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
274 if (0==Q)
275 return freeup_new (P);
276 P->opaque = Q;
277
278 Q->phits = pj_param (P->ctx, P->params, "tlat_ts").i ?
279 pj_param (P->ctx, P->params, "rlat_ts").f : M_HALFPI;
280
281 return setup(P);
282 }
283
284
PROJECTION(ups)285 PJ *PROJECTION(ups) {
286 struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
287 if (0==Q)
288 return freeup_new (P);
289 P->opaque = Q;
290
291 /* International Ellipsoid */
292 P->phi0 = pj_param(P->ctx, P->params, "bsouth").i ? - M_HALFPI: M_HALFPI;
293 if (!P->es) E_ERROR(-34);
294 P->k0 = .994;
295 P->x0 = 2000000.;
296 P->y0 = 2000000.;
297 Q->phits = M_HALFPI;
298 P->lam0 = 0.;
299
300 return setup(P);
301 }
302
303
304 #ifndef PJ_SELFTEST
pj_stere_selftest(void)305 int pj_stere_selftest (void) {return 0;}
306 #else
307
pj_stere_selftest(void)308 int pj_stere_selftest (void) {
309 double tolerance_lp = 1e-10;
310 double tolerance_xy = 1e-7;
311
312 char e_args[] = {"+proj=stere +ellps=GRS80 +lat_1=0.5 +lat_2=2 +n=0.5"};
313 char s_args[] = {"+proj=stere +a=6400000 +lat_1=0.5 +lat_2=2 +n=0.5"};
314
315 LP fwd_in[] = {
316 { 2, 1},
317 { 2,-1},
318 {-2, 1},
319 {-2,-1}
320 };
321
322 XY e_fwd_expect[] = {
323 { 445289.70910023432, 221221.76694834774},
324 { 445289.70910023432, -221221.76694835056},
325 {-445289.70910023432, 221221.76694834774},
326 {-445289.70910023432, -221221.76694835056},
327 };
328
329 XY s_fwd_expect[] = {
330 { 223407.81025950745, 111737.938996443},
331 { 223407.81025950745, -111737.938996443},
332 {-223407.81025950745, 111737.938996443},
333 {-223407.81025950745, -111737.938996443},
334 };
335
336 XY inv_in[] = {
337 { 200, 100},
338 { 200,-100},
339 {-200, 100},
340 {-200,-100}
341 };
342
343 LP e_inv_expect[] = {
344 { 0.0017966305682022392, 0.00090436947502443507},
345 { 0.0017966305682022392, -0.00090436947502443507},
346 {-0.0017966305682022392, 0.00090436947502443507},
347 {-0.0017966305682022392, -0.00090436947502443507},
348 };
349
350 LP s_inv_expect[] = {
351 { 0.001790493109747395, 0.00089524655465513144},
352 { 0.001790493109747395, -0.00089524655465513144},
353 {-0.001790493109747395, 0.00089524655465513144},
354 {-0.001790493109747395, -0.00089524655465513144},
355 };
356
357 return pj_generic_selftest (e_args, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, s_fwd_expect, inv_in, e_inv_expect, s_inv_expect);
358 }
359
360
361 #endif
362
363
364
365
366
367 #ifndef PJ_SELFTEST
pj_ups_selftest(void)368 int pj_ups_selftest (void) {return 0;}
369 #else
370
pj_ups_selftest(void)371 int pj_ups_selftest (void) {
372 double tolerance_lp = 1e-10;
373 double tolerance_xy = 1e-7;
374
375 char e_args[] = {"+proj=ups +ellps=GRS80 +lat_1=0.5 +lat_2=2 +n=0.5"};
376
377 LP fwd_in[] = {
378 { 2, 1},
379 { 2,-1},
380 {-2, 1},
381 {-2,-1}
382 };
383
384 XY e_fwd_expect[] = {
385 {2433455.5634384668, -10412543.301512826},
386 {2448749.1185681992, -10850493.419804076},
387 {1566544.4365615332, -10412543.301512826},
388 {1551250.8814318008, -10850493.419804076},
389 };
390
391 XY inv_in[] = {
392 { 200, 100},
393 { 200,-100},
394 {-200, 100},
395 {-200,-100}
396 };
397
398 LP e_inv_expect[] = {
399 {-44.998567498074834, 64.9182362867341},
400 {-44.995702709112308, 64.917020250675748},
401 {-45.004297076028529, 64.915804280954518},
402 {-45.001432287066002, 64.914588377560719},
403 };
404
405 return pj_generic_selftest (e_args, 0, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, 0, inv_in, e_inv_expect, 0);
406 }
407
408
409 #endif
410