1 /*
2 ** libproj -- library of cartographic projections
3 **
4 ** Copyright (c) 2003, 2006 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 struct pj_opaque {
30 double s0;
31 double A1, A2, A3, A4, A5, A6;
32 double B1, B2, B3, B4, B5, B6, B7, B8;
33 double C1, C2, C3, C4, C5, C6, C7, C8;
34 double D1, D2, D3, D4, D5, D6, D7, D8, D9, D10, D11;
35 void *en;
36 };
37 PROJ_HEAD(rouss, "Roussilhe Stereographic") "\n\tAzi., Ellps.";
38
39
e_forward(LP lp,PJ * P)40 static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */
41 XY xy = {0.0,0.0};
42 struct pj_opaque *Q = P->opaque;
43 double s, al, cp, sp, al2, s2;
44
45 cp = cos(lp.phi);
46 sp = sin(lp.phi);
47 s = proj_mdist(lp.phi, sp, cp, Q->en) - Q->s0;
48 s2 = s * s;
49 al = lp.lam * cp / sqrt(1. - P->es * sp * sp);
50 al2 = al * al;
51 xy.x = P->k0 * al*(1.+s2*(Q->A1+s2*Q->A4)-al2*(Q->A2+s*Q->A3+s2*Q->A5
52 +al2*Q->A6));
53 xy.y = P->k0 * (al2*(Q->B1+al2*Q->B4)+
54 s*(1.+al2*(Q->B3-al2*Q->B6)+s2*(Q->B2+s2*Q->B8)+
55 s*al2*(Q->B5+s*Q->B7)));
56
57 return xy;
58 }
59
60
e_inverse(XY xy,PJ * P)61 static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */
62 LP lp = {0.0,0.0};
63 struct pj_opaque *Q = P->opaque;
64 double s, al, x = xy.x / P->k0, y = xy.y / P->k0, x2, y2;;
65
66 x2 = x * x;
67 y2 = y * y;
68 al = x*(1.-Q->C1*y2+x2*(Q->C2+Q->C3*y-Q->C4*x2+Q->C5*y2-Q->C7*x2*y)
69 +y2*(Q->C6*y2-Q->C8*x2*y));
70 s = Q->s0 + y*(1.+y2*(-Q->D2+Q->D8*y2))+
71 x2*(-Q->D1+y*(-Q->D3+y*(-Q->D5+y*(-Q->D7+y*Q->D11)))+
72 x2*(Q->D4+y*(Q->D6+y*Q->D10)-x2*Q->D9));
73 lp.phi=proj_inv_mdist(P->ctx, s, Q->en);
74 s = sin(lp.phi);
75 lp.lam=al * sqrt(1. - P->es * s * s)/cos(lp.phi);
76
77 return lp;
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 if (P->opaque->en)
88 pj_dealloc (P->opaque->en);
89 pj_dealloc (P->opaque);
90 return pj_dealloc(P);
91 }
92
freeup(PJ * P)93 static void freeup (PJ *P) {
94 freeup_new (P);
95 return;
96 }
97
98
PROJECTION(rouss)99 PJ *PROJECTION(rouss) {
100 double N0, es2, t, t2, R_R0_2, R_R0_4;
101
102 struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
103 if (0==Q)
104 return freeup_new (P);
105 P->opaque = Q;
106
107 if (!((Q->en = proj_mdist_ini(P->es))))
108 E_ERROR_0;
109 es2 = sin(P->phi0);
110 Q->s0 = proj_mdist(P->phi0, es2, cos(P->phi0), Q->en);
111 t = 1. - (es2 = P->es * es2 * es2);
112 N0 = 1./sqrt(t);
113 R_R0_2 = t * t / P->one_es;
114 R_R0_4 = R_R0_2 * R_R0_2;
115 t = tan(P->phi0);
116 t2 = t * t;
117 Q->C1 = Q->A1 = R_R0_2 / 4.;
118 Q->C2 = Q->A2 = R_R0_2 * (2 * t2 - 1. - 2. * es2) / 12.;
119 Q->A3 = R_R0_2 * t * (1. + 4. * t2)/ ( 12. * N0);
120 Q->A4 = R_R0_4 / 24.;
121 Q->A5 = R_R0_4 * ( -1. + t2 * (11. + 12. * t2))/24.;
122 Q->A6 = R_R0_4 * ( -2. + t2 * (11. - 2. * t2))/240.;
123 Q->B1 = t / (2. * N0);
124 Q->B2 = R_R0_2 / 12.;
125 Q->B3 = R_R0_2 * (1. + 2. * t2 - 2. * es2)/4.;
126 Q->B4 = R_R0_2 * t * (2. - t2)/(24. * N0);
127 Q->B5 = R_R0_2 * t * (5. + 4.* t2)/(8. * N0);
128 Q->B6 = R_R0_4 * (-2. + t2 * (-5. + 6. * t2))/48.;
129 Q->B7 = R_R0_4 * (5. + t2 * (19. + 12. * t2))/24.;
130 Q->B8 = R_R0_4 / 120.;
131 Q->C3 = R_R0_2 * t * (1. + t2)/(3. * N0);
132 Q->C4 = R_R0_4 * (-3. + t2 * (34. + 22. * t2))/240.;
133 Q->C5 = R_R0_4 * (4. + t2 * (13. + 12. * t2))/24.;
134 Q->C6 = R_R0_4 / 16.;
135 Q->C7 = R_R0_4 * t * (11. + t2 * (33. + t2 * 16.))/(48. * N0);
136 Q->C8 = R_R0_4 * t * (1. + t2 * 4.)/(36. * N0);
137 Q->D1 = t / (2. * N0);
138 Q->D2 = R_R0_2 / 12.;
139 Q->D3 = R_R0_2 * (2 * t2 + 1. - 2. * es2) / 4.;
140 Q->D4 = R_R0_2 * t * (1. + t2)/(8. * N0);
141 Q->D5 = R_R0_2 * t * (1. + t2 * 2.)/(4. * N0);
142 Q->D6 = R_R0_4 * (1. + t2 * (6. + t2 * 6.))/16.;
143 Q->D7 = R_R0_4 * t2 * (3. + t2 * 4.)/8.;
144 Q->D8 = R_R0_4 / 80.;
145 Q->D9 = R_R0_4 * t * (-21. + t2 * (178. - t2 * 26.))/720.;
146 Q->D10 = R_R0_4 * t * (29. + t2 * (86. + t2 * 48.))/(96. * N0);
147 Q->D11 = R_R0_4 * t * (37. + t2 * 44.)/(96. * N0);
148
149 P->fwd = e_forward;
150 P->inv = e_inverse;
151
152 return P;
153 }
154
155
156 #ifndef PJ_SELFTEST
pj_rouss_selftest(void)157 int pj_rouss_selftest (void) {return 0;}
158 #else
159
pj_rouss_selftest(void)160 int pj_rouss_selftest (void) {
161 double tolerance_lp = 1e-10;
162 double tolerance_xy = 1e-7;
163
164 char e_args[] = {"+proj=rouss +ellps=GRS80 +lat_1=0.5 +lat_2=2"};
165
166 LP fwd_in[] = {
167 { 2, 1},
168 { 2,-1},
169 {-2, 1},
170 {-2,-1}
171 };
172
173 XY e_fwd_expect[] = {
174 { 222644.89413161727, 110611.09186837047},
175 { 222644.89413161727, -110611.09186837047},
176 {-222644.89413161727, 110611.09186837047},
177 {-222644.89413161727, -110611.09186837047},
178 };
179
180 XY inv_in[] = {
181 { 200, 100},
182 { 200,-100},
183 {-200, 100},
184 {-200,-100}
185 };
186
187 LP e_inv_expect[] = {
188 { 0.0017966305682019911, 0.00090436947683699559},
189 { 0.0017966305682019911, -0.00090436947683699559},
190 {-0.0017966305682019911, 0.00090436947683699559},
191 {-0.0017966305682019911, -0.00090436947683699559},
192 };
193
194 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);
195 }
196
197
198 #endif
199