aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/phi2.cpp175
-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.cpp26
-rw-r--r--src/projections/tobmerc.cpp19
-rw-r--r--src/strerrno.cpp2
-rw-r--r--src/tsfn.cpp29
8 files changed, 168 insertions, 114 deletions
diff --git a/src/phi2.cpp b/src/phi2.cpp
index b81456b0..eb6d5c82 100644
--- a/src/phi2.cpp
+++ b/src/phi2.cpp
@@ -1,6 +1,7 @@
/* Determine latitude angle phi-2. */
#include <math.h>
+#include <limits>
#include "proj.h"
#include "proj_internal.h"
@@ -8,61 +9,125 @@
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
+ constexpr double tol = sqrt(std::numeric_limits<double>::epsilon()) / 10;
+ constexpr double tmax = 2 / sqrt(std::numeric_limits<double>::epsilon());
+ double
+ e2m = 1 - e * e,
+ tau = fabs(taup) > 70 ? taup * exp(e * atanh(e)) : taup / e2m,
+ stol = tol * std::max(1.0, fabs(taup));
+ if (!(fabs(tau) < tmax)) return tau; // handles +/-inf and nan and e = 1
+ // if (e2m < 0) return std::numeric_limits<double>::quiet_NaN();
+ int i = numit;
+ for (; i; --i) {
+ double tau1 = sqrt(1 + tau * tau),
+ sig = sinh( e * atanh(e * tau / tau1) ),
+ taupa = sqrt(1 + sig * sig) * tau - sig * tau1,
+ dtau = (taup - taupa) * (1 + e2m * tau * tau) /
+ ( e2m * sqrt(1 + tau * tau) * sqrt(1 + taupa * taupa) );
+ tau += dtau;
+ if (!(fabs(dtau) >= stol))
+ 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)
+ * 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))
+ *
+ * OLD: This routine inverts this relation using the iterative scheme given
+ * by Snyder (1987), Eqs. (7-9) - (7-11).
+ *
+ * NEW: This routine writes converts t = exp(-psi) to
+ *
+ * tau' = 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..b83b8072 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( abs(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( abs(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..472cb43f 100644
--- a/src/projections/merc.cpp
+++ b/src/projections/merc.cpp
@@ -10,45 +10,25 @@
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));
+ xy.y = P->k0 * (asinh(tan(lp.phi)) - P->e * atanh(P->e * sin(lp.phi)));
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..fe8f29ed 100644
--- a/src/tsfn.cpp
+++ b/src/tsfn.cpp
@@ -4,14 +4,25 @@
#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)
+ * 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))
+ ***************************************************************************/
- /* 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 );
}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 22:05:54 +0000