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