1 #define PJ_LIB__
2
3 #include <float.h>
4 #include <math.h>
5
6 #include "proj.h"
7 #include "proj_internal.h"
8 #include <math.h>
9
10 PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts=";
11 PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Ell\n\t";
12
13 #define EPS10 1.e-10
logtanpfpim1(double x)14 static double logtanpfpim1(double x) { /* log(tan(x/2 + M_FORTPI)) */
15 if (fabs(x) <= DBL_EPSILON) {
16 /* tan(M_FORTPI + .5 * x) can be approximated by 1.0 + x */
17 return log1p(x);
18 }
19 return log(tan(M_FORTPI + .5 * x));
20 }
21
merc_e_forward(PJ_LP lp,PJ * P)22 static PJ_XY merc_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
23 PJ_XY xy = {0.0,0.0};
24 if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
25 proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
26 return xy;
27 }
28 xy.x = P->k0 * lp.lam;
29 xy.y = - P->k0 * log(pj_tsfn(lp.phi, sin(lp.phi), P->e));
30 return xy;
31 }
32
33
merc_s_forward(PJ_LP lp,PJ * P)34 static PJ_XY merc_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */
35 PJ_XY xy = {0.0,0.0};
36 if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
37 proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
38 return xy;
39 }
40 xy.x = P->k0 * lp.lam;
41 xy.y = P->k0 * logtanpfpim1(lp.phi);
42 return xy;
43 }
44
45
merc_e_inverse(PJ_XY xy,PJ * P)46 static PJ_LP merc_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
47 PJ_LP lp = {0.0,0.0};
48 if ((lp.phi = pj_phi2(P->ctx, exp(- xy.y / P->k0), P->e)) == HUGE_VAL) {
49 proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
50 return lp;
51 }
52 lp.lam = xy.x / P->k0;
53 return lp;
54 }
55
56
merc_s_inverse(PJ_XY xy,PJ * P)57 static PJ_LP merc_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
58 PJ_LP lp = {0.0,0.0};
59 lp.phi = atan(sinh(xy.y / P->k0));
60 lp.lam = xy.x / P->k0;
61 return lp;
62 }
63
64
PROJECTION(merc)65 PJ *PROJECTION(merc) {
66 double phits=0.0;
67 int is_phits;
68
69 if( (is_phits = pj_param(P->ctx, P->params, "tlat_ts").i) ) {
70 phits = fabs(pj_param(P->ctx, P->params, "rlat_ts").f);
71 if (phits >= M_HALFPI)
72 return pj_default_destructor(P, PJD_ERR_LAT_TS_LARGER_THAN_90);
73 }
74
75 if (P->es != 0.0) { /* ellipsoid */
76 if (is_phits)
77 P->k0 = pj_msfn(sin(phits), cos(phits), P->es);
78 P->inv = merc_e_inverse;
79 P->fwd = merc_e_forward;
80 }
81
82 else { /* sphere */
83 if (is_phits)
84 P->k0 = cos(phits);
85 P->inv = merc_s_inverse;
86 P->fwd = merc_s_forward;
87 }
88
89 return P;
90 }
91
PROJECTION(webmerc)92 PJ *PROJECTION(webmerc) {
93
94 /* Overriding k_0 with fixed parameter */
95 P->k0 = 1.0;
96
97 P->inv = merc_s_inverse;
98 P->fwd = merc_s_forward;
99 return P;
100 }
101