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