1 /*
2 ** libproj -- library of cartographic projections
3 **
4 ** Copyright (c) 2003 Gerald I. Evenden
5 */
6 /*
7 ** Permission is hereby granted, free of charge, to any person obtaining
8 ** a copy of this software and associated documentation files (the
9 ** "Software"), to deal in the Software without restriction, including
10 ** without limitation the rights to use, copy, modify, merge, publish,
11 ** distribute, sublicense, and/or sell copies of the Software, and to
12 ** permit persons to whom the Software is furnished to do so, subject to
13 ** the following conditions:
14 **
15 ** The above copyright notice and this permission notice shall be
16 ** included in all copies or substantial portions of the Software.
17 **
18 ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
19 ** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
20 ** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
21 ** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
22 ** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
23 ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
24 ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
25 */
26 #define PJ_LIB__
27 #include <projects.h>
28
29
30 struct pj_opaque {
31 double phic0;
32 double cosc0, sinc0;
33 double R2;
34 void *en;
35 };
36
37
38 PROJ_HEAD(sterea, "Oblique Stereographic Alternative") "\n\tAzimuthal, Sph&Ell";
39 # define DEL_TOL 1.e-14
40 # define MAX_ITER 10
41
42
43
e_forward(LP lp,PJ * P)44 static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */
45 XY xy = {0.0,0.0};
46 struct pj_opaque *Q = P->opaque;
47 double cosc, sinc, cosl, k;
48
49 lp = pj_gauss(P->ctx, lp, Q->en);
50 sinc = sin(lp.phi);
51 cosc = cos(lp.phi);
52 cosl = cos(lp.lam);
53 k = P->k0 * Q->R2 / (1. + Q->sinc0 * sinc + Q->cosc0 * cosc * cosl);
54 xy.x = k * cosc * sin(lp.lam);
55 xy.y = k * (Q->cosc0 * sinc - Q->sinc0 * cosc * cosl);
56 return xy;
57 }
58
59
e_inverse(XY xy,PJ * P)60 static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */
61 LP lp = {0.0,0.0};
62 struct pj_opaque *Q = P->opaque;
63 double rho, c, sinc, cosc;
64
65 xy.x /= P->k0;
66 xy.y /= P->k0;
67 if ( (rho = hypot (xy.x, xy.y)) ) {
68 c = 2. * atan2 (rho, Q->R2);
69 sinc = sin (c);
70 cosc = cos (c);
71 lp.phi = asin (cosc * Q->sinc0 + xy.y * sinc * Q->cosc0 / rho);
72 lp.lam = atan2 (xy.x * sinc, rho * Q->cosc0 * cosc - xy.y * Q->sinc0 * sinc);
73 } else {
74 lp.phi = Q->phic0;
75 lp.lam = 0.;
76 }
77 return pj_inv_gauss(P->ctx, lp, Q->en);
78 }
79
80
freeup_new(PJ * P)81 static void *freeup_new (PJ *P) { /* Destructor */
82 if (0==P)
83 return 0;
84 if (0==P->opaque)
85 return pj_dealloc (P);
86
87 pj_dealloc (P->opaque->en);
88 pj_dealloc (P->opaque);
89 return pj_dealloc(P);
90 }
91
92
freeup(PJ * P)93 static void freeup (PJ *P) {
94 freeup_new (P);
95 return;
96 }
97
98
PROJECTION(sterea)99 PJ *PROJECTION(sterea) {
100 double R;
101 struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
102
103 if (0==Q)
104 return freeup_new (P);
105 P->opaque = Q;
106
107 Q->en = pj_gauss_ini(P->e, P->phi0, &(Q->phic0), &R);
108 if (0==Q->en)
109 return freeup_new (P);
110
111 Q->sinc0 = sin (Q->phic0);
112 Q->cosc0 = cos (Q->phic0);
113 Q->R2 = 2. * R;
114
115 P->inv = e_inverse;
116 P->fwd = e_forward;
117 return P;
118 }
119
120
121 #ifndef PJ_SELFTEST
pj_sterea_selftest(void)122 int pj_sterea_selftest (void) {return 0;}
123 #else
124
pj_sterea_selftest(void)125 int pj_sterea_selftest (void) {
126 double tolerance_lp = 1e-10;
127 double tolerance_xy = 1e-7;
128
129 char e_args[] = {"+proj=sterea +ellps=GRS80 +lat_1=0.5 +lat_2=2 +n=0.5"};
130 char s_args[] = {"+proj=sterea +a=6400000 +lat_1=0.5 +lat_2=2 +n=0.5"};
131
132 LP fwd_in[] = {
133 { 2, 1},
134 { 2,-1},
135 {-2, 1},
136 {-2,-1}
137 };
138
139 XY e_fwd_expect[] = {
140 { 222644.89410919772, 110611.09187173686},
141 { 222644.89410919772, -110611.09187173827},
142 {-222644.89410919772, 110611.09187173686},
143 {-222644.89410919772, -110611.09187173827},
144 };
145
146 XY s_fwd_expect[] = {
147 { 223407.81025950745, 111737.93899644315},
148 { 223407.81025950745, -111737.93899644315},
149 {-223407.81025950745, 111737.93899644315},
150 {-223407.81025950745, -111737.93899644315},
151 };
152
153 XY inv_in[] = {
154 { 200, 100},
155 { 200,-100},
156 {-200, 100},
157 {-200,-100}
158 };
159
160 LP e_inv_expect[] = {
161 { 0.0017966305682019911, 0.00090436947683099009},
162 { 0.0017966305682019911, -0.00090436947684371233},
163 {-0.0017966305682019911, 0.00090436947683099009},
164 {-0.0017966305682019911, -0.00090436947684371233},
165 };
166
167 LP s_inv_expect[] = {
168 { 0.001790493109747395, 0.00089524655465446378},
169 { 0.001790493109747395, -0.00089524655465446378},
170 {-0.001790493109747395, 0.00089524655465446378},
171 {-0.001790493109747395, -0.00089524655465446378},
172 };
173
174 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);
175 }
176
177
178 #endif
179