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
128
129
130
131
132
133
134
135
136
137
|
/* Determine latitude angle phi-2. */
#include <math.h>
#include <limits>
#include <algorithm>
#include "proj.h"
#include "proj_internal.h"
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()),
tol = rooteps / 10, // the convergence criterion for Newton's method
tmax = 2 / rooteps; // the large arg limit is exact for tau > tmax
const double
e2m = 1 - e * e,
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),
sig = sinh( e * atanh(e * tau / tau1) ),
taupa = sqrt(1 + sig * sig) * tau - sig * tau1,
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)
* 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));
}
|
