aboutsummaryrefslogtreecommitdiff
path: root/src/projections
diff options
context:
space:
mode:
Diffstat (limited to 'src/projections')
-rw-r--r--src/projections/gstmerc.cpp18
-rw-r--r--src/projections/lcc.cpp10
-rw-r--r--src/projections/merc.cpp30
-rw-r--r--src/projections/tobmerc.cpp19
4 files changed, 29 insertions, 48 deletions
diff --git a/src/projections/gstmerc.cpp b/src/projections/gstmerc.cpp
index 808d9ef7..50814bb5 100644
--- a/src/projections/gstmerc.cpp
+++ b/src/projections/gstmerc.cpp
@@ -28,9 +28,9 @@ static PJ_XY gstmerc_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forw
double L, Ls, sinLs1, Ls1;
L = Q->n1*lp.lam;
- Ls = Q->c + Q->n1 * log(pj_tsfn(-1.0 * lp.phi, -1.0 * sin(lp.phi), P->e));
+ Ls = Q->c + Q->n1 * log(pj_tsfn(-lp.phi, -sin(lp.phi), P->e));
sinLs1 = sin(L) / cosh(Ls);
- Ls1 = log(pj_tsfn(-1.0 * asin(sinLs1), 0.0, 0.0));
+ Ls1 = log(pj_tsfn(-asin(sinLs1), -sinLs1, 0.0));
xy.x = (Q->XS + Q->n2*Ls1) * P->ra;
xy.y = (Q->YS + Q->n2*atan(sinh(Ls) / cos(L))) * P->ra;
@@ -45,9 +45,9 @@ static PJ_LP gstmerc_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inve
L = atan(sinh((xy.x * P->a - Q->XS) / Q->n2) / cos((xy.y * P->a - Q->YS) / Q->n2));
sinC = sin((xy.y * P->a - Q->YS) / Q->n2) / cosh((xy.x * P->a - Q->XS) / Q->n2);
- LC = log(pj_tsfn(-1.0 * asin(sinC), 0.0, 0.0));
+ LC = log(pj_tsfn(-asin(sinC), -sinC, 0.0));
lp.lam = L / Q->n1;
- lp.phi = -1.0 * pj_phi2(P->ctx, exp((LC - Q->c) / Q->n1), P->e);
+ lp.phi = -pj_phi2(P->ctx, exp((LC - Q->c) / Q->n1), P->e);
return lp;
}
@@ -60,13 +60,13 @@ PJ *PROJECTION(gstmerc) {
P->opaque = Q;
Q->lamc = P->lam0;
- Q->n1 = sqrt(1.0 + P->es * pow(cos(P->phi0), 4.0) / (1.0 - P->es));
+ Q->n1 = sqrt(1 + P->es * pow(cos(P->phi0), 4.0) / (1 - P->es));
Q->phic = asin(sin(P->phi0) / Q->n1);
- Q->c = log(pj_tsfn(-1.0 * Q->phic, 0.0, 0.0))
- - Q->n1 * log(pj_tsfn(-1.0 * P->phi0, -1.0 * sin(P->phi0), P->e));
- Q->n2 = P->k0 * P->a * sqrt(1.0 - P->es) / (1.0 - P->es * sin(P->phi0) * sin(P->phi0));
+ Q->c = log(pj_tsfn(-Q->phic, -sin(P->phi0) / Q->n1, 0.0))
+ - Q->n1 * log(pj_tsfn(-P->phi0, -sin(P->phi0), P->e));
+ Q->n2 = P->k0 * P->a * sqrt(1 - P->es) / (1 - P->es * sin(P->phi0) * sin(P->phi0));
Q->XS = 0;
- Q->YS = -1.0 * Q->n2 * Q->phic;
+ Q->YS = -Q->n2 * Q->phic;
P->inv = gstmerc_s_inverse;
P->fwd = gstmerc_s_forward;
diff --git a/src/projections/lcc.cpp b/src/projections/lcc.cpp
index 91ffc511..46378ce4 100644
--- a/src/projections/lcc.cpp
+++ b/src/projections/lcc.cpp
@@ -106,10 +106,10 @@ PJ *PROJECTION(lcc) {
double ml1, m1;
m1 = pj_msfn(sinphi, cosphi, P->es);
- ml1 = pj_tsfn(Q->phi1, sinphi, P->e);
- if( ml1 == 0 ) {
+ if( fabs(Q->phi1) == M_HALFPI ) {
return pj_default_destructor(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90);
}
+ ml1 = pj_tsfn(Q->phi1, sinphi, P->e);
if (secant) { /* secant cone */
sinphi = sin(Q->phi2);
Q->n = log(m1 / pj_msfn(sinphi, cos(Q->phi2), P->es));
@@ -117,10 +117,10 @@ PJ *PROJECTION(lcc) {
// Not quite, but es is very close to 1...
return pj_default_destructor(P, PJD_ERR_INVALID_ECCENTRICITY);
}
- const double ml2 = pj_tsfn(Q->phi2, sinphi, P->e);
- if( ml2 == 0 ) {
- return pj_default_destructor(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90);
+ if( fabs(Q->phi2) == M_HALFPI ) {
+ return pj_default_destructor(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90);
}
+ const double ml2 = pj_tsfn(Q->phi2, sinphi, P->e);
const double denom = log(ml1 / ml2);
if( denom == 0 ) {
// Not quite, but es is very close to 1...
diff --git a/src/projections/merc.cpp b/src/projections/merc.cpp
index a77d7517..3a0ed7b4 100644
--- a/src/projections/merc.cpp
+++ b/src/projections/merc.cpp
@@ -10,45 +10,29 @@
PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts=";
PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Ell\n\t";
-#define EPS10 1.e-10
-static double logtanpfpim1(double x) { /* log(tan(x/2 + M_FORTPI)) */
- if (fabs(x) <= DBL_EPSILON) {
- /* tan(M_FORTPI + .5 * x) can be approximated by 1.0 + x */
- return log1p(x);
- }
- return log(tan(M_FORTPI + .5 * x));
-}
-
static PJ_XY merc_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
PJ_XY xy = {0.0,0.0};
- if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
- proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
- return xy;
- }
xy.x = P->k0 * lp.lam;
- xy.y = - P->k0 * log(pj_tsfn(lp.phi, sin(lp.phi), P->e));
+ // Instead of calling tan and sin, call sin and cos which the compiler
+ // optimizes to a single call to sincos.
+ double sphi = sin(lp.phi);
+ double cphi = cos(lp.phi);
+ xy.y = P->k0 * (asinh(sphi/cphi) - P->e * atanh(P->e * sphi));
return xy;
}
static PJ_XY merc_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */
PJ_XY xy = {0.0,0.0};
- if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
- proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
- return xy;
-}
xy.x = P->k0 * lp.lam;
- xy.y = P->k0 * logtanpfpim1(lp.phi);
+ xy.y = P->k0 * asinh(tan(lp.phi));
return xy;
}
static PJ_LP merc_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
PJ_LP lp = {0.0,0.0};
- if ((lp.phi = pj_phi2(P->ctx, exp(- xy.y / P->k0), P->e)) == HUGE_VAL) {
- proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
- return lp;
-}
+ lp.phi = atan(pj_sinhpsi2tanphi(P->ctx, sinh(xy.y / P->k0), P->e));
lp.lam = xy.x / P->k0;
return lp;
}
diff --git a/src/projections/tobmerc.cpp b/src/projections/tobmerc.cpp
index a1616036..f05a9b6b 100644
--- a/src/projections/tobmerc.cpp
+++ b/src/projections/tobmerc.cpp
@@ -9,27 +9,24 @@
PROJ_HEAD(tobmerc, "Tobler-Mercator") "\n\tCyl, Sph";
-#define EPS10 1.e-10
-static double logtanpfpim1(double x) { /* log(tan(x/2 + M_FORTPI)) */
- if (fabs(x) <= DBL_EPSILON) {
- /* tan(M_FORTPI + .5 * x) can be approximated by 1.0 + x */
- return log1p(x);
- }
- return log(tan(M_FORTPI + .5 * x));
-}
-
static PJ_XY tobmerc_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */
PJ_XY xy = {0.0, 0.0};
double cosphi;
- if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
+ if (fabs(lp.phi) >= M_HALFPI) {
+ // builtins.gie tests "Test expected failure at the poles:". However
+ // given that M_HALFPI is strictly less than pi/2 in double precision,
+ // it's not clear why shouldn't just return a large result for xy.y (and
+ // it's not even that large, merely 38.025...). Even if the logic was
+ // such that phi was strictly equal to pi/2, allowing xy.y = inf would be
+ // a reasonable result.
proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
return xy;
}
cosphi = cos(lp.phi);
xy.x = P->k0 * lp.lam * cosphi * cosphi;
- xy.y = P->k0 * logtanpfpim1(lp.phi);
+ xy.y = P->k0 * asinh(tan(lp.phi));
return xy;
}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 19:20:22 +0000