1 #define PJ_LIB__
2
3 #include <errno.h>
4 #include <math.h>
5
6 #include "proj.h"
7 #include "proj_internal.h"
8
9 PROJ_HEAD(poly, "Polyconic (American)")
10 "\n\tConic, Sph&Ell";
11
12 namespace { // anonymous namespace
13 struct pj_opaque {
14 double ml0; \
15 double *en;
16 };
17 } // anonymous namespace
18
19 #define TOL 1e-10
20 #define CONV 1e-10
21 #define N_ITER 10
22 #define I_ITER 20
23 #define ITOL 1.e-12
24
25
poly_e_forward(PJ_LP lp,PJ * P)26 static PJ_XY poly_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
27 PJ_XY xy = {0.0,0.0};
28 struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
29 double ms, sp, cp;
30
31 if (fabs(lp.phi) <= TOL) {
32 xy.x = lp.lam;
33 xy.y = -Q->ml0;
34 } else {
35 sp = sin(lp.phi);
36 cp = cos(lp.phi);
37 ms = fabs(cp) > TOL ? pj_msfn(sp, cp, P->es) / sp : 0.;
38 lp.lam *= sp;
39 xy.x = ms * sin(lp.lam);
40 xy.y = (pj_mlfn(lp.phi, sp, cp, Q->en) - Q->ml0) + ms * (1. - cos(lp.lam));
41 }
42
43 return xy;
44 }
45
46
poly_s_forward(PJ_LP lp,PJ * P)47 static PJ_XY poly_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */
48 PJ_XY xy = {0.0,0.0};
49 struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
50
51 if (fabs(lp.phi) <= TOL) {
52 xy.x = lp.lam;
53 xy.y = Q->ml0;
54 } else {
55 const double cot = 1. / tan(lp.phi);
56 const double E = lp.lam * sin(lp.phi);
57 xy.x = sin(E) * cot;
58 xy.y = lp.phi - P->phi0 + cot * (1. - cos(E));
59 }
60
61 return xy;
62 }
63
64
poly_e_inverse(PJ_XY xy,PJ * P)65 static PJ_LP poly_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
66 PJ_LP lp = {0.0,0.0};
67 struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
68
69 xy.y += Q->ml0;
70 if (fabs(xy.y) <= TOL) {
71 lp.lam = xy.x;
72 lp.phi = 0.;
73 } else {
74 int i;
75
76 const double r = xy.y * xy.y + xy.x * xy.x;
77 lp.phi = xy.y;
78 for (i = I_ITER; i ; --i) {
79 const double sp = sin(lp.phi);
80 const double cp = cos(lp.phi);
81 const double s2ph = sp * cp;
82 if (fabs(cp) < ITOL) {
83 proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
84 return lp;
85 }
86 double mlp = sqrt(1. - P->es * sp * sp);
87 const double c = sp * mlp / cp;
88 const double ml = pj_mlfn(lp.phi, sp, cp, Q->en);
89 const double mlb = ml * ml + r;
90 mlp = P->one_es / (mlp * mlp * mlp);
91 const double dPhi =
92 ( ml + ml + c * mlb - 2. * xy.y * (c * ml + 1.) ) / (
93 P->es * s2ph * (mlb - 2. * xy.y * ml) / c +
94 2.* (xy.y - ml) * (c * mlp - 1. / s2ph) - mlp - mlp );
95 lp.phi += dPhi;
96 if (fabs(dPhi) <= ITOL)
97 break;
98 }
99 if (!i) {
100 proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
101 return lp;
102 }
103 const double c = sin(lp.phi);
104 lp.lam = asin(xy.x * tan(lp.phi) * sqrt(1. - P->es * c * c)) / sin(lp.phi);
105 }
106
107 return lp;
108 }
109
110
poly_s_inverse(PJ_XY xy,PJ * P)111 static PJ_LP poly_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
112 PJ_LP lp = {0.0,0.0};
113
114 if (fabs(xy.y = P->phi0 + xy.y) <= TOL) {
115 lp.lam = xy.x;
116 lp.phi = 0.;
117 } else {
118 lp.phi = xy.y;
119 const double B = xy.x * xy.x + xy.y * xy.y;
120 int i = N_ITER;
121 while(true) {
122 const double tp = tan(lp.phi);
123 const double dphi = (xy.y * (lp.phi * tp + 1.) - lp.phi -
124 .5 * ( lp.phi * lp.phi + B) * tp) /
125 ((lp.phi - xy.y) / tp - 1.);
126 lp.phi -= dphi;
127 if( !(fabs(dphi) > CONV) )
128 break;
129 --i;
130 if( i == 0 ) {
131 proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
132 return lp;
133 }
134 }
135 lp.lam = asin(xy.x * tan(lp.phi)) / sin(lp.phi);
136 }
137
138 return lp;
139 }
140
141
destructor(PJ * P,int errlev)142 static PJ *destructor(PJ *P, int errlev) {
143 if (nullptr==P)
144 return nullptr;
145
146 if (nullptr==P->opaque)
147 return pj_default_destructor (P, errlev);
148
149 if (static_cast<struct pj_opaque*>(P->opaque)->en)
150 pj_dealloc (static_cast<struct pj_opaque*>(P->opaque)->en);
151
152 return pj_default_destructor(P, errlev);
153 }
154
155
PROJECTION(poly)156 PJ *PROJECTION(poly) {
157 struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
158 if (nullptr==Q)
159 return pj_default_destructor (P, ENOMEM);
160
161 P->opaque = Q;
162 P->destructor = destructor;
163
164 if (P->es != 0.0) {
165 if (!(Q->en = pj_enfn(P->es)))
166 return pj_default_destructor (P, ENOMEM);
167 Q->ml0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), Q->en);
168 P->inv = poly_e_inverse;
169 P->fwd = poly_e_forward;
170 } else {
171 Q->ml0 = -P->phi0;
172 P->inv = poly_s_inverse;
173 P->fwd = poly_s_forward;
174 }
175
176 return P;
177 }
178