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(cass, "Cassini") "\n\tCyl, Sph&Ell";
10
11
12 # define C1 .16666666666666666666
13 # define C2 .00833333333333333333
14 # define C3 .04166666666666666666
15 # define C4 .33333333333333333333
16 # define C5 .06666666666666666666
17
18
19 namespace { // anonymous namespace
20 struct pj_opaque {
21 double *en;
22 double m0;
23 };
24 } // anonymous namespace
25
26
27
cass_e_forward(PJ_LP lp,PJ * P)28 static PJ_XY cass_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
29 double n, t, a1, c, a2, tn;
30 PJ_XY xy = {0.0, 0.0};
31 struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
32
33 n = sin (lp.phi);
34 c = cos (lp.phi);
35 xy.y = pj_mlfn (lp.phi, n, c, Q->en);
36
37 n = 1./sqrt(1. - P->es * n*n);
38 tn = tan(lp.phi);
39 t = tn * tn;
40 a1 = lp.lam * c;
41 c *= P->es * c / (1 - P->es);
42 a2 = a1 * a1;
43
44 xy.x = n * a1 * (1. - a2 * t *
45 (C1 - (8. - t + 8. * c) * a2 * C2));
46 xy.y -= Q->m0 - n * tn * a2 *
47 (.5 + (5. - t + 6. * c) * a2 * C3);
48
49 return xy;
50 }
51
52
cass_s_forward(PJ_LP lp,PJ * P)53 static PJ_XY cass_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */
54 PJ_XY xy = {0.0, 0.0};
55 xy.x = asin (cos (lp.phi) * sin (lp.lam));
56 xy.y = atan2 (tan (lp.phi), cos (lp.lam)) - P->phi0;
57 return xy;
58 }
59
60
cass_e_inverse(PJ_XY xy,PJ * P)61 static PJ_LP cass_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
62 double n, t, r, dd, d2, tn, ph1;
63 PJ_LP lp = {0.0, 0.0};
64 struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
65
66 ph1 = pj_inv_mlfn (P->ctx, Q->m0 + xy.y, P->es, Q->en);
67 tn = tan (ph1);
68 t = tn*tn;
69 n = sin (ph1);
70 r = 1. / (1. - P->es * n * n);
71 n = sqrt (r);
72 r *= (1. - P->es) * n;
73 dd = xy.x / n;
74 d2 = dd * dd;
75 lp.phi = ph1 - (n * tn / r) * d2 *
76 (.5 - (1. + 3. * t) * d2 * C3);
77 lp.lam = dd * (1. + t * d2 *
78 (-C4 + (1. + 3. * t) * d2 * C5)) / cos (ph1);
79 return lp;
80 }
81
82
cass_s_inverse(PJ_XY xy,PJ * P)83 static PJ_LP cass_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
84 PJ_LP lp = {0.0,0.0};
85 double dd;
86 lp.phi = asin(sin(dd = xy.y + P->phi0) * cos(xy.x));
87 lp.lam = atan2(tan(xy.x), cos(dd));
88 return lp;
89 }
90
destructor(PJ * P,int errlev)91 static PJ *destructor (PJ *P, int errlev) { /* Destructor */
92 if (nullptr==P)
93 return nullptr;
94
95 if (nullptr==P->opaque)
96 return pj_default_destructor (P, errlev);
97
98 pj_dealloc (static_cast<struct pj_opaque*>(P->opaque)->en);
99 return pj_default_destructor (P, errlev);
100 }
101
102
103
PROJECTION(cass)104 PJ *PROJECTION(cass) {
105
106 /* Spheroidal? */
107 if (0==P->es) {
108 P->inv = cass_s_inverse;
109 P->fwd = cass_s_forward;
110 return P;
111 }
112
113 /* otherwise it's ellipsoidal */
114 P->opaque = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
115 if (nullptr==P->opaque)
116 return pj_default_destructor (P, ENOMEM);
117 P->destructor = destructor;
118
119 static_cast<struct pj_opaque*>(P->opaque)->en = pj_enfn (P->es);
120 if (nullptr==static_cast<struct pj_opaque*>(P->opaque)->en)
121 return pj_default_destructor (P, ENOMEM);
122
123 static_cast<struct pj_opaque*>(P->opaque)->m0 = pj_mlfn (P->phi0, sin (P->phi0), cos (P->phi0), static_cast<struct pj_opaque*>(P->opaque)->en);
124 P->inv = cass_e_inverse;
125 P->fwd = cass_e_forward;
126
127 return P;
128 }
129