aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--docs/source/operations/projections/merc.rst30
-rw-r--r--src/apps/gie.cpp2
-rw-r--r--src/phi2.cpp180
-rw-r--r--src/proj_internal.h3
-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
-rw-r--r--src/strerrno.cpp2
-rw-r--r--src/tsfn.cpp32
-rw-r--r--test/gie/builtins.gie116
-rw-r--r--test/unit/pj_phi2_test.cpp87
12 files changed, 344 insertions, 185 deletions
diff --git a/docs/source/operations/projections/merc.rst b/docs/source/operations/projections/merc.rst
index 7b6e13da..2107b7b2 100644
--- a/docs/source/operations/projections/merc.rst
+++ b/docs/source/operations/projections/merc.rst
@@ -158,7 +158,35 @@ Inverse projection
\lambda = \frac{x}{k_0 a}; \quad \psi = \frac{y}{k_0 a}
The latitude :math:`\phi` is found by inverting the equation for
-:math:`\psi` iteratively.
+:math:`\psi`. This follows the method given by :cite:`Karney2011tm`.
+Start by introducing the conformal latitude
+
+.. math::
+ \chi = \tan^{-1}\sinh\psi
+
+The tangents of the latitudes :math:`\tau = \tan\phi` and :math:`\tau' =
+\tan\chi = \sinh\psi` are related by
+
+.. math::
+ \tau' = \tau \sqrt{1 + \sigma^2} - \sigma \sqrt{1 + \tau^2}
+
+where
+
+.. math::
+ \sigma = \sinh\bigl(e \tanh^{-1}(e \tau/\sqrt{1 + \tau^2}) \bigr)
+
+This is obtained by taking the :math:`\sinh` of the equation for
+:math:`\psi` and using the multiple argument formula. The equation for
+:math:`\tau'` can be solved to give :math:`\tau` using Newton's method
+using :math:`\tau = \tau'/(1 - e^2)` as an initial guess and with the
+needed derivative given by
+
+..math::
+ \frac{d\tau'}{d\tau} = \frac{1 - e^2}{1 + (1 - e^2)\tau^2}
+ \sqrt{1 + \tau'^2} \sqrt{1 + \tau^2}
+
+This converges after no more than 2 iterations. Finally set
+:math:`\phi=\tan^{-1}\tau`.
Further reading
###############
diff --git a/src/apps/gie.cpp b/src/apps/gie.cpp
index 8940afde..2fe854fa 100644
--- a/src/apps/gie.cpp
+++ b/src/apps/gie.cpp
@@ -1124,7 +1124,7 @@ static const struct errno_vs_err_const lookup[] = {
{"pjd_err_invalid_x_or_y" , -15},
{"pjd_err_wrong_format_dms_value" , -16},
{"pjd_err_non_conv_inv_meri_dist" , -17},
- {"pjd_err_non_con_inv_phi2" , -18},
+ {"pjd_err_non_conv_sinhpsi2tanphi" , -18},
{"pjd_err_acos_asin_arg_too_large" , -19},
{"pjd_err_tolerance_condition" , -20},
{"pjd_err_conic_lat_equal" , -21},
diff --git a/src/phi2.cpp b/src/phi2.cpp
index b81456b0..1c48d67f 100644
--- a/src/phi2.cpp
+++ b/src/phi2.cpp
@@ -1,68 +1,134 @@
/* Determine latitude angle phi-2. */
#include <math.h>
+#include <limits>
+#include <algorithm>
#include "proj.h"
#include "proj_internal.h"
-static const double TOL = 1.0e-10;
-static const int N_ITER = 15;
+double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) {
+ /****************************************************************************
+ * Convert tau' = sinh(psi) = tan(chi) to tau = tan(phi). The code is taken
+ * from GeographicLib::Math::tauf(taup, e).
+ *
+ * Here
+ * phi = geographic latitude (radians)
+ * psi is the isometric latitude
+ * psi = asinh(tan(phi)) - e * atanh(e * sin(phi))
+ * = asinh(tan(chi))
+ * chi is the conformal latitude
+ *
+ * The representation of latitudes via their tangents, tan(phi) and tan(chi),
+ * maintains full *relative* accuracy close to latitude = 0 and +/- pi/2.
+ * This is sometimes important, e.g., to compute the scale of the transverse
+ * Mercator projection which involves cos(phi)/cos(chi) tan(phi)
+ *
+ * From Karney (2011), Eq. 7,
+ *
+ * tau' = sinh(psi) = sinh(asinh(tan(phi)) - e * atanh(e * sin(phi)))
+ * = tan(phi) * cosh(e * atanh(e * sin(phi))) -
+ * sec(phi) * sinh(e * atanh(e * sin(phi)))
+ * = tau * sqrt(1 + sigma^2) - sqrt(1 + tau^2) * sigma
+ * where
+ * sigma = sinh(e * atanh( e * tau / sqrt(1 + tau^2) ))
+ *
+ * For e small,
+ *
+ * tau' = (1 - e^2) * tau
+ *
+ * The relation tau'(tau) can therefore by reliably inverted by Newton's
+ * method with
+ *
+ * tau = tau' / (1 - e^2)
+ *
+ * as an initial guess. Newton's method requires dtau'/dtau. Noting that
+ *
+ * dsigma/dtau = e^2 * sqrt(1 + sigma^2) /
+ * (sqrt(1 + tau^2) * (1 + (1 - e^2) * tau^2))
+ * d(sqrt(1 + tau^2))/dtau = tau / sqrt(1 + tau^2)
+ *
+ * we have
+ *
+ * dtau'/dtau = (1 - e^2) * sqrt(1 + tau'^2) * sqrt(1 + tau^2) /
+ * (1 + (1 - e^2) * tau^2)
+ *
+ * This works fine unless tau^2 and tau'^2 overflows. This may be partially
+ * cured by writing, e.g., sqrt(1 + tau^2) as hypot(1, tau). However, nan
+ * will still be generated with tau' = inf, since (inf - inf) will appear in
+ * the Newton iteration.
+ *
+ * If we note that for sufficiently large |tau|, i.e., |tau| >= 2/sqrt(eps),
+ * sqrt(1 + tau^2) = |tau| and
+ *
+ * tau' = exp(- e * atanh(e)) * tau
+ *
+ * So
+ *
+ * tau = exp(e * atanh(e)) * tau'
+ *
+ * can be returned unless |tau| >= 2/sqrt(eps); this then avoids overflow
+ * problems for large tau' and returns the correct result for tau' = +/-inf
+ * and nan.
+ *
+ * Newton's method usually take 2 iterations to converge to double precision
+ * accuracy (for WGS84 flattening). However only 1 iteration is needed for
+ * |chi| < 3.35 deg. In addition, only 1 iteration is needed for |chi| >
+ * 89.18 deg (tau' > 70), if tau = exp(e * atanh(e)) * tau' is used as the
+ * starting guess.
+ ****************************************************************************/
+
+ constexpr int numit = 5;
+ // min iterations = 1, max iterations = 2; mean = 1.954
+ static const double rooteps = sqrt(std::numeric_limits<double>::epsilon());
+ static const double tol = rooteps / 10; // the criterion for Newton's method
+ static const double tmax = 2 / rooteps; // threshold for large arg limit exact
+ const double e2m = 1 - e * e;
+ const double stol = tol * std::max(1.0, fabs(taup));
+ // The initial guess. 70 corresponds to chi = 89.18 deg (see above)
+ double tau = fabs(taup) > 70 ? taup * exp(e * atanh(e)) : taup / e2m;
+ if (!(fabs(tau) < tmax)) // handles +/-inf and nan and e = 1
+ return tau;
+ // If we need to deal with e > 1, then we could include:
+ // if (e2m < 0) return std::numeric_limits<double>::quiet_NaN();
+ int i = numit;
+ for (; i; --i) {
+ double tau1 = sqrt(1 + tau * tau);
+ double sig = sinh( e * atanh(e * tau / tau1) );
+ double taupa = sqrt(1 + sig * sig) * tau - sig * tau1;
+ double dtau = ( (taup - taupa) * (1 + e2m * (tau * tau)) /
+ (e2m * tau1 * sqrt(1 + taupa * taupa)) );
+ tau += dtau;
+ if (!(fabs(dtau) >= stol)) // backwards test to allow nans to succeed.
+ break;
+ }
+ if (i == 0)
+ pj_ctx_set_errno(ctx, PJD_ERR_NON_CONV_SINHPSI2TANPHI);
+ return tau;
+}
/*****************************************************************************/
double pj_phi2(projCtx ctx, const double ts0, const double e) {
-/******************************************************************************
-Determine latitude angle phi-2.
-Inputs:
- ts = exp(-psi) where psi is the isometric latitude (dimensionless)
- e = eccentricity of the ellipsoid (dimensionless)
-Output:
- phi = geographic latitude (radians)
-Here isometric latitude is defined by
- psi = log( tan(pi/4 + phi/2) *
- ( (1 - e*sin(phi)) / (1 + e*sin(phi)) )^(e/2) )
- = asinh(tan(phi)) - e * atanh(e * sin(phi))
-This routine inverts this relation using the iterative scheme given
-by Snyder (1987), Eqs. (7-9) - (7-11)
-*******************************************************************************/
- const double eccnth = .5 * e;
- double ts = ts0;
-#ifdef no_longer_used_original_convergence_on_exact_dphi
- double Phi = M_HALFPI - 2. * atan(ts);
-#endif
- int i = N_ITER;
-
- for(;;) {
- /*
- * sin(Phi) = sin(PI/2 - 2* atan(ts))
- * = cos(2*atan(ts))
- * = 2*cos^2(atan(ts)) - 1
- * = 2 / (1 + ts^2) - 1
- * = (1 - ts^2) / (1 + ts^2)
- */
- const double sinPhi = (1 - ts * ts) / (1 + ts * ts);
- const double con = e * sinPhi;
- double old_ts = ts;
- ts = ts0 * pow((1. - con) / (1. + con), eccnth);
-#ifdef no_longer_used_original_convergence_on_exact_dphi
- /* The convergence criterion is nominally on exact dphi */
- const double newPhi = M_HALFPI - 2. * atan(ts);
- const double dphi = newPhi - Phi;
- Phi = newPhi;
-#else
- /* If we don't immediately update Phi from this, we can
- * change the conversion criterion to save us computing atan() at each step.
- * Particularly we can observe that:
- * |atan(ts) - atan(old_ts)| <= |ts - old_ts|
- * So if |ts - old_ts| matches our convergence criterion, we're good.
- */
- const double dphi = 2 * (ts - old_ts);
-#endif
- if (fabs(dphi) > TOL && --i) {
- continue;
- }
- break;
- }
- if (i <= 0)
- pj_ctx_set_errno(ctx, PJD_ERR_NON_CON_INV_PHI2);
- return M_HALFPI - 2. * atan(ts);
+ /****************************************************************************
+ * Determine latitude angle phi-2.
+ * Inputs:
+ * ts = exp(-psi) where psi is the isometric latitude (dimensionless)
+ * this variable is defined in Snyder (1987), Eq. (7-10)
+ * e = eccentricity of the ellipsoid (dimensionless)
+ * Output:
+ * phi = geographic latitude (radians)
+ * Here isometric latitude is defined by
+ * psi = log( tan(pi/4 + phi/2) *
+ * ( (1 - e*sin(phi)) / (1 + e*sin(phi)) )^(e/2) )
+ * = asinh(tan(phi)) - e * atanh(e * sin(phi))
+ * = asinh(tan(chi))
+ * chi = conformal latitude
+ *
+ * This routine converts t = exp(-psi) to
+ *
+ * tau' = tan(chi) = sinh(psi) = (1/t - t)/2
+ *
+ * returns atan(sinpsi2tanphi(tau'))
+ ***************************************************************************/
+ return atan(pj_sinhpsi2tanphi(ctx, (1/ts0 - ts0) / 2, e));
}
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 79b1da10..203765a3 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -626,7 +626,7 @@ struct FACTORS {
#define PJD_ERR_INVALID_X_OR_Y -15
#define PJD_ERR_WRONG_FORMAT_DMS_VALUE -16
#define PJD_ERR_NON_CONV_INV_MERI_DIST -17
-#define PJD_ERR_NON_CON_INV_PHI2 -18
+#define PJD_ERR_NON_CONV_SINHPSI2TANPHI -18
#define PJD_ERR_ACOS_ASIN_ARG_TOO_LARGE -19
#define PJD_ERR_TOLERANCE_CONDITION -20
#define PJD_ERR_CONIC_LAT_EQUAL -21
@@ -846,6 +846,7 @@ double pj_qsfn(double, double, double);
double pj_tsfn(double, double, double);
double pj_msfn(double, double, double);
double PROJ_DLL pj_phi2(projCtx_t *, const double, const double);
+double pj_sinhpsi2tanphi(projCtx_t *, const double, const double);
double pj_qsfn_(double, PJ *);
double *pj_authset(double);
double pj_authlat(double, double *);
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;
}
diff --git a/src/strerrno.cpp b/src/strerrno.cpp
index 5ae0d7e1..1c8673d0 100644
--- a/src/strerrno.cpp
+++ b/src/strerrno.cpp
@@ -27,7 +27,7 @@ pj_err_list[] = {
"invalid x or y", /* -15 */
"improperly formed DMS value", /* -16 */
"non-convergent inverse meridional dist", /* -17 */
- "non-convergent inverse phi2", /* -18 */
+ "non-convergent sinh(psi) to tan(phi)", /* -18 */
"acos/asin: |arg| >1.+1e-14", /* -19 */
"tolerance condition error", /* -20 */
"conic lat_1 = -lat_2", /* -21 */
diff --git a/src/tsfn.cpp b/src/tsfn.cpp
index 32da09f2..8ed258d6 100644
--- a/src/tsfn.cpp
+++ b/src/tsfn.cpp
@@ -4,14 +4,28 @@
#include "proj_internal.h"
double pj_tsfn(double phi, double sinphi, double e) {
- double denominator;
- sinphi *= e;
+ /****************************************************************************
+ * Determine function ts(phi) defined in Snyder (1987), Eq. (7-10)
+ * Inputs:
+ * phi = geographic latitude (radians)
+ * e = eccentricity of the ellipsoid (dimensionless)
+ * Output:
+ * ts = exp(-psi) where psi is the isometric latitude (dimensionless)
+ * = 1 / (tan(chi) + sec(chi))
+ * Here isometric latitude is defined by
+ * psi = log( tan(pi/4 + phi/2) *
+ * ( (1 - e*sin(phi)) / (1 + e*sin(phi)) )^(e/2) )
+ * = asinh(tan(phi)) - e * atanh(e * sin(phi))
+ * = asinh(tan(chi))
+ * chi = conformal latitude
+ ***************************************************************************/
- /* avoid zero division, fail gracefully */
- denominator = 1.0 + sinphi;
- if (denominator == 0.0)
- return HUGE_VAL;
-
- return (tan (.5 * (M_HALFPI - phi)) /
- pow((1. - sinphi) / (denominator), .5 * e));
+ double cosphi = cos(phi);
+ // exp(-asinh(tan(phi))) = 1 / (tan(phi) + sec(phi))
+ // = cos(phi) / (1 + sin(phi)) good for phi > 0
+ // = (1 - sin(phi)) / cos(phi) good for phi < 0
+ return exp(e * atanh(e * sinphi)) *
+ ( sinphi > 0 ?
+ cosphi / (1 + sinphi) :
+ (1 - sinphi) / cosphi );
}
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index cfce5041..def30206 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -1423,9 +1423,9 @@ expect -0.001790143 -0.000895247
===============================================================================
-------------------------------------------------------------------------------
-operation +proj=etmerc +ellps=GRS80 +zone=30
+operation +proj=etmerc +ellps=GRS80
-------------------------------------------------------------------------------
-tolerance 0.1 mm
+tolerance 50 nm
accept 2 1
expect 222650.796797586 110642.229411933
accept 2 -1
@@ -1434,17 +1434,28 @@ accept -2 1
expect -222650.796797586 110642.229411933
accept -2 -1
expect -222650.796797586 -110642.229411933
+# near pole
+accept 30 89.9999
+expect 5.584698978 10001956.056248082
+# 3900 km from central meridian
+accept 44.69 35.37
+expect 4168136.489446198 4985511.302287407
direction inverse
accept 200 100
-expect 0.001796631 0.000904369
+expect 0.00179663056816 0.00090436947663
accept 200 -100
-expect 0.001796631 -0.000904369
+expect 0.00179663056816 -0.00090436947663
accept -200 100
-expect -0.001796631 0.000904369
+expect -0.00179663056816 0.00090436947663
accept -200 -100
-expect -0.001796631 -0.000904369
-
+expect -0.00179663056816 -0.00090436947663
+# near pole
+accept 6 1.0001e7
+expect 0.35596960759234 89.99135362646302
+# 3900 km from central meridian
+accept 4168136.489446198 4985511.302287407
+expect 44.69 35.37
===============================================================================
# Fahey
@@ -3355,30 +3366,54 @@ expect -0.001953415 -0.000820580
-------------------------------------------------------------------------------
operation +proj=merc +ellps=GRS80
-------------------------------------------------------------------------------
-tolerance 0.1 mm
+tolerance 0 m
+accept 0 0
+expect 0 0
+tolerance 50 nm
accept 2 1
-expect 222638.981586547 110579.965218250
+expect 222638.981586547 110579.965218249
accept 2 -1
expect 222638.981586547 -110579.965218249
accept -2 1
-expect -222638.981586547 110579.965218250
+expect -222638.981586547 110579.965218249
accept -2 -1
expect -222638.981586547 -110579.965218249
+# inflate tolerance by scale (k = 5.7e15)
+tolerance 3e8
+accept 0 89.99999999999999
+expect 0 235805185.015130176
+accept 0 -89.99999999999999
+expect 0 -235805185.015130176
direction inverse
+tolerance 0 m
+accept 0 0
+expect 0 0
+tolerance 50 nm
accept 200 100
-expect 0.001796631 0.000904369
+expect 0.00179663056824 0.00090436947704
accept 200 -100
-expect 0.001796631 -0.000904369
+expect 0.00179663056824 -0.00090436947704
accept -200 100
-expect -0.001796631 0.000904369
+expect -0.00179663056824 0.00090436947704
accept -200 -100
-expect -0.001796631 -0.000904369
+expect -0.00179663056824 -0.00090436947704
+accept 0 235805185.015130176
+expect 0 89.99999999999999
+accept 0 -235805185.015130176
+expect 0 -89.99999999999999
+accept 0 1e10
+expect 0 90
+accept 0 -1e10
+expect 0 -90
-------------------------------------------------------------------------------
operation +proj=merc +R=6400000
-------------------------------------------------------------------------------
-tolerance 0.1 mm
+tolerance 0 m
+accept 0 0
+expect 0 0
+tolerance 50 nm
accept 2 1
expect 223402.144255274 111706.743574944
accept 2 -1
@@ -3389,25 +3424,32 @@ accept -2 -1
expect -223402.144255274 -111706.743574944
direction inverse
+tolerance 0 m
+accept 0 0
+expect 0 0
+tolerance 50 nm
accept 200 100
-expect 0.001790493 0.000895247
+expect 0.00179049310978 0.00089524655486
accept 200 -100
-expect 0.001790493 -0.000895247
+expect 0.00179049310978 -0.00089524655486
accept -200 100
-expect -0.001790493 0.000895247
+expect -0.00179049310978 0.00089524655486
accept -200 -100
-expect -0.001790493 -0.000895247
-
+expect -0.00179049310978 -0.00089524655486
-------------------------------------------------------------------------------
operation +proj=merc +R=1
-------------------------------------------------------------------------------
# Test the numerical stability of the inverse spherical Mercator
-------------------------------------------------------------------------------
-tolerance 1e-15 m
-accept 0 1e-15
+tolerance 1e-17 m
+accept 0 57.295779513e-15
expect 0 1e-15
+direction inverse
+accept 0 1e-15
+expect 0 57.295779513e-15
+
===============================================================================
# Miller Oblated Stereographic
@@ -5658,25 +5700,37 @@ expect -0.001790143 0.511651393
-------------------------------------------------------------------------------
operation +proj=tmerc +ellps=GRS80
-------------------------------------------------------------------------------
-tolerance 0.1 mm
+tolerance 50 nm
accept 2 1
-expect 222650.796795778 110642.229411927
+expect 222650.796797586 110642.229411933
accept 2 -1
-expect 222650.796795778 -110642.229411927
+expect 222650.796797586 -110642.229411933
accept -2 1
-expect -222650.796795778 110642.229411927
+expect -222650.796797586 110642.229411933
accept -2 -1
-expect -222650.796795778 -110642.229411927
+expect -222650.796797586 -110642.229411933
+# near pole
+accept 30 89.9999
+expect 5.584698978 10001956.056248082
+# 3900 km from central meridian
+accept 44.69 35.37
+expect 4168136.489446198 4985511.302287407
direction inverse
accept 200 100
-expect 0.001796631 0.000904369
+expect 0.00179663056816 0.00090436947663
accept 200 -100
-expect 0.001796631 -0.000904369
+expect 0.00179663056816 -0.00090436947663
accept -200 100
-expect -0.001796631 0.000904369
+expect -0.00179663056816 0.00090436947663
accept -200 -100
-expect -0.001796631 -0.000904369
+expect -0.00179663056816 -0.00090436947663
+# near pole
+accept 6 1.0001e7
+expect 0.35596960759234 89.99135362646302
+# 3900 km from central meridian
+accept 4168136.489446198 4985511.302287407
+expect 44.69 35.37
-------------------------------------------------------------------------------
operation +proj=tmerc +R=6400000
diff --git a/test/unit/pj_phi2_test.cpp b/test/unit/pj_phi2_test.cpp
index c4db6e52..fdedae98 100644
--- a/test/unit/pj_phi2_test.cpp
+++ b/test/unit/pj_phi2_test.cpp
@@ -39,47 +39,62 @@ namespace {
TEST(PjPhi2Test, Basic) {
projCtx ctx = pj_get_default_ctx();
- EXPECT_DOUBLE_EQ(M_PI_2, pj_phi2(ctx, 0.0, 0.0));
-
- EXPECT_NEAR(0.0, pj_phi2(ctx, 1.0, 0.0), 1e-16);
- EXPECT_DOUBLE_EQ(M_PI_2, pj_phi2(ctx, 0.0, 1.0));
- EXPECT_DOUBLE_EQ(M_PI, pj_phi2(ctx, -1.0, 0.0));
- EXPECT_DOUBLE_EQ(M_PI_2, pj_phi2(ctx, 0.0, -1.0));
-
- EXPECT_NEAR(0.0, pj_phi2(ctx, 1.0, 1.0), 1e-16);
- EXPECT_DOUBLE_EQ(M_PI, pj_phi2(ctx, -1.0, -1.0));
-
- // TODO(schwehr): M_PI_4, M_PI_2, M_PI, M_E
- // https://www.gnu.org/software/libc/manual/html_node/Mathematical-Constants.html
-
- EXPECT_DOUBLE_EQ(-0.95445818456292697, pj_phi2(ctx, M_PI, 0.0));
- EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, M_PI)));
- EXPECT_DOUBLE_EQ(4.0960508381527205, pj_phi2(ctx, -M_PI, 0.0));
- EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, -M_PI)));
-
- EXPECT_TRUE(std::isnan(pj_phi2(ctx, M_PI, M_PI)));
- EXPECT_TRUE(std::isnan(pj_phi2(ctx, -M_PI, -M_PI)));
+ // Expectation is that only sane values of e (and nan is here reckoned to
+ // be sane) are passed to pj_phi2. Thus the return value with other values
+ // of e is "implementation dependent".
+
+ constexpr auto inf = std::numeric_limits<double>::infinity();
+ constexpr auto nan = std::numeric_limits<double>::quiet_NaN();
+
+ // Strict equality is demanded here.
+ EXPECT_EQ( M_PI_2, pj_phi2(ctx, +0.0, 0.0));
+ EXPECT_EQ( 0.0 , pj_phi2(ctx, 1.0, 0.0));
+ EXPECT_EQ(-M_PI_2, pj_phi2(ctx, inf, 0.0));
+ // We don't expect pj_phi2 to be called with negative ts (since ts =
+ // exp(-psi)). However, in the current implementation it is odd in ts.
+ // N.B. ts = +0.0 and ts = -0.0 return different results.
+ EXPECT_EQ(-M_PI_2, pj_phi2(ctx, -0.0, 0.0));
+ EXPECT_EQ( 0.0 , pj_phi2(ctx, -1.0, 0.0));
+ EXPECT_EQ(+M_PI_2, pj_phi2(ctx, -inf, 0.0));
+
+ constexpr double e = 0.2;
+ EXPECT_EQ( M_PI_2, pj_phi2(ctx, +0.0, e));
+ EXPECT_EQ( 0.0 , pj_phi2(ctx, 1.0, e));
+ EXPECT_EQ(-M_PI_2, pj_phi2(ctx, inf, e));
+ EXPECT_EQ(-M_PI_2, pj_phi2(ctx, -0.0, e));
+ EXPECT_EQ( 0.0 , pj_phi2(ctx, -1.0, e));
+ EXPECT_EQ(+M_PI_2, pj_phi2(ctx, -inf, e));
+
+ EXPECT_TRUE(std::isnan(pj_phi2(ctx, nan, 0.0)));
+ EXPECT_TRUE(std::isnan(pj_phi2(ctx, nan, e )));
+ EXPECT_TRUE(std::isnan(pj_phi2(ctx, +0.0, nan)));
+ EXPECT_TRUE(std::isnan(pj_phi2(ctx, 1.0, nan)));
+ EXPECT_TRUE(std::isnan(pj_phi2(ctx, inf, nan)));
+ EXPECT_TRUE(std::isnan(pj_phi2(ctx, -0.0, nan)));
+ EXPECT_TRUE(std::isnan(pj_phi2(ctx, -1.0, nan)));
+ EXPECT_TRUE(std::isnan(pj_phi2(ctx, -inf, nan)));
+ EXPECT_TRUE(std::isnan(pj_phi2(ctx, nan, nan)));
+
+ EXPECT_DOUBLE_EQ( M_PI/3, pj_phi2(ctx, 1/(sqrt(3.0)+2), 0.0));
+ EXPECT_DOUBLE_EQ( M_PI/4, pj_phi2(ctx, 1/(sqrt(2.0)+1), 0.0));
+ EXPECT_DOUBLE_EQ( M_PI/6, pj_phi2(ctx, 1/ sqrt(3.0) , 0.0));
+ EXPECT_DOUBLE_EQ(-M_PI/3, pj_phi2(ctx, sqrt(3.0)+2 , 0.0));
+ EXPECT_DOUBLE_EQ(-M_PI/4, pj_phi2(ctx, sqrt(2.0)+1 , 0.0));
+ EXPECT_DOUBLE_EQ(-M_PI/6, pj_phi2(ctx, sqrt(3.0) , 0.0));
+
+ // Generated with exp(e * atanh(e * sin(phi))) / (tan(phi) + sec(phi))
+ EXPECT_DOUBLE_EQ( M_PI/3, pj_phi2(ctx, 0.27749174377027023413, e));
+ EXPECT_DOUBLE_EQ( M_PI/4, pj_phi2(ctx, 0.42617788119104192995, e));
+ EXPECT_DOUBLE_EQ( M_PI/6, pj_phi2(ctx, 0.58905302448626726064, e));
+ EXPECT_DOUBLE_EQ(-M_PI/3, pj_phi2(ctx, 3.6037108218537833089, e));
+ EXPECT_DOUBLE_EQ(-M_PI/4, pj_phi2(ctx, 2.3464380582241712935, e));
+ EXPECT_DOUBLE_EQ(-M_PI/6, pj_phi2(ctx, 1.6976400399134411849, e));
}
-TEST(PjPhi2Test, AvoidUndefinedBehavior) {
- auto ctx = pj_get_default_ctx();
+} // namespace
- const auto nan = std::numeric_limits<double>::quiet_NaN();
- EXPECT_TRUE(std::isnan(pj_phi2(ctx, nan, 0.0)));
- EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, nan)));
- EXPECT_TRUE(std::isnan(pj_phi2(ctx, nan, nan)));
- // We do not really care about the values that follow.
- const auto inf = std::numeric_limits<double>::infinity();
- EXPECT_DOUBLE_EQ(-M_PI_2, pj_phi2(ctx, inf, 0.0));
- EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, inf)));
- EXPECT_DOUBLE_EQ(4.7123889803846897, pj_phi2(ctx, -inf, 0.0));
- EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, -inf)));
- EXPECT_TRUE(std::isnan(pj_phi2(ctx, inf, inf)));
- EXPECT_TRUE(std::isnan(pj_phi2(ctx, -inf, -inf)));
-}
-} // namespace
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 18:14:48 +0000