Document the default poder_engsager algorithm for tmerc. (#2379)

* Document the default poder_engsager algorithm for tmerc.

* Give exact expression for phi' in terms of phi

* Aadd another datapoint on range of applicability + explanation for
complex numbers.

* Update tmerc figure with one reflecting poder/engsager algo.

Courtesy of @hobu.

Co-authored-by: Charles Karney <charles.karney@sri.com>

3 files changed, 238 insertions, 42 deletions

diff --git a/docs/source/operations/projections/images/tmerc.png b/docs/source/operations/projections/images/tmerc.png
index 29d4776e..8ef73986 100644
--- a/docs/source/operations/projections/images/tmerc.png
+++ b/docs/source/operations/projections/images/tmerc.png
diff --git a/docs/source/operations/projections/tmerc.rst b/docs/source/operations/projections/tmerc.rst
index cf428cff..3163262f 100644
--- a/docs/source/operations/projections/tmerc.rst
+++ b/docs/source/operations/projections/tmerc.rst
@@ -11,7 +11,7 @@ The transverse Mercator projection in its various forms is the most widely used
 +---------------------+----------------------------------------------------------+
 | **Available forms** | Forward and inverse, spherical and ellipsoidal           |
 +---------------------+----------------------------------------------------------+
-| **Defined area**    | Global, but reasonably accurate only within 15 degrees   |
+| **Defined area**    | Global, with full accuracy within 3900 km                |
 |                     | of the central meridian                                  |
 +---------------------+----------------------------------------------------------+
 | **Alias**           | tmerc                                                    |
@@ -42,21 +42,21 @@ In the post-war years, these concepts were extended into the Universal Transvers
 
 The following table gives special cases of the Transverse Mercator projection.
 
-+-------------------------------------+-----------------------------------------------------+--------------------------------+------------------------------------------+--------------+
-| Projection Name                     | Areas                                               | Central meridian               | Zone width                               | Scale Factor |
-+=====================================+=====================================================+================================+==========================================+==============+
-| Transverse Mercator                 | World wide                                          | Various                        | less than 6°                             | Various      |
-+-------------------------------------+-----------------------------------------------------+--------------------------------+------------------------------------------+--------------+
-| Transverse Mercator south oriented  | Southern Africa                                     | 2° intervals E of 11°E         | 2°                                       | 1.000        |
-+-------------------------------------+-----------------------------------------------------+--------------------------------+------------------------------------------+--------------+
-| UTM North hemisphere                | World wide equator to 84°N                          | 6° intervals E & W of 3° E & W | Always 6°                                | 0.9996       |
-+-------------------------------------+-----------------------------------------------------+--------------------------------+------------------------------------------+--------------+
-| UTM South hemisphere                | World wide north of 80°S to equator                 | 6° intervals E & W of 3° E & W | Always 6°                                | 0.9996       |
-+-------------------------------------+-----------------------------------------------------+--------------------------------+------------------------------------------+--------------+
-| Gauss-Kruger                        | Former USSR, Yugoslavia, Germany, S. America, China | Various, according to area     | Usually less than 6°, often less than 4° | 1.0000       |
-+-------------------------------------+-----------------------------------------------------+--------------------------------+------------------------------------------+--------------+
-| Gauss Boaga                         | Italy                                               | Various, according to area     | 6°                                       | 0.9996       |
-+-------------------------------------+-----------------------------------------------------+--------------------------------+------------------------------------------+--------------+
++-------------------------------------+-----------------------------------------------------+--------------------------------+-------------------------------------------+--------------+
+| Projection Name                     | Areas                                               | Central meridian               | Zone width                                | Scale Factor |
++=====================================+=====================================================+================================+===========================================+==============+
+| Transverse Mercator                 | World wide                                          | Various                        | less than 1000 km                         | Various      |
++-------------------------------------+-----------------------------------------------------+--------------------------------+-------------------------------------------+--------------+
+| Transverse Mercator south oriented  | Southern Africa                                     | 2° intervals E of 11°E         | 2°                                        | 1.000        |
++-------------------------------------+-----------------------------------------------------+--------------------------------+-------------------------------------------+--------------+
+| UTM North hemisphere                | World wide equator to 84°N                          | 6° intervals E & W of 3° E & W | Usually 6°, wider for Norway and Svalbard | 0.9996       |
++-------------------------------------+-----------------------------------------------------+--------------------------------+-------------------------------------------+--------------+
+| UTM South hemisphere                | World wide north of 80°S to equator                 | 6° intervals E & W of 3° E & W | Always 6°                                 | 0.9996       |
++-------------------------------------+-----------------------------------------------------+--------------------------------+-------------------------------------------+--------------+
+| Gauss-Kruger                        | Former USSR, Yugoslavia, Germany, S. America, China | Various, according to area     | Usually less than 6°, often less than 4°  | 1.0000       |
++-------------------------------------+-----------------------------------------------------+--------------------------------+-------------------------------------------+--------------+
+| Gauss Boaga                         | Italy                                               | Various, according to area     | 6°                                        | 0.9996       |
++-------------------------------------+-----------------------------------------------------+--------------------------------+-------------------------------------------+--------------+
 
 
 
@@ -79,16 +79,16 @@ Parameters
 
     .. versionadded:: 6.0.0
 
-    Use the algorithm described in section "Ellipsoidal Form" below.
-    It is faster than the default algorithm, but also diverges faster
-    as the distance from the central meridian increases.
+    Use the Evenden-Snyder algorithm described below under "Legacy
+    ellipsoidal form".  It is faster than the default algorithm, but is
+    less accurate and diverges beyond 3° from the central meridian.
 
 .. option:: +algo=auto/evenden_snyder/poder_engsager
 
     .. versionadded:: 7.1
 
     Selects the algorithm to use. The hardcoded value and the one defined in
-    :ref:`proj-ini` default to ``poder_engsager``, that is the most precise
+    :ref:`proj-ini` default to ``poder_engsager``; that is the most precise
     one.
 
     When using auto, a heuristics based on the input coordinate to transform
@@ -115,60 +115,222 @@ Parameters
 Mathematical definition
 #######################
 
-The formulas describing the Transverse Mercator below are quoted from Evenden's [Evenden2005]_.
+The formulation given here for the Transverse Mercator projection is due
+to Krüger :cite:`Krueger1912` who gave the series expansions accurate to
+:math:`n^4`, where :math:`n = (a-b)/(a+b)` is the third flattening.
+These series were extended to sixth order by Engsager and Poder in
+:cite:`Poder1998` and :cite:`Engsager2007`.  This gives full
+double-precision accuracy withing 3900 km of the central meridian (about
+57% of the surface of the earth) :cite:`Karney2011tm`.  The error is
+less than 0.1 mm within 7000 km of the central meridian (about 89% of
+the surface of the earth).
+
+This formulation consists of three steps: a conformal projection from
+the ellipsoid to a sphere, the spherical transverse Mercator
+projection, rectifying this projection to give constant scale on the
+central meridian.
+
+The scale on the central meridian is :math:`k_0` and is set by ``+k_0``.
+
+Option :option:`+lon_0` sets the central meridian; in the formulation
+below :math:`\lambda` is the longitude relative to the central meridian.
+
+Options :option:`+lat_0`, :option:`+x_0`, and :option:`+y_0` serve to
+translate the projected coordinates so that at :math:`(\phi, \lambda) =
+(\phi_0, \lambda_0)`, the projected coordinates are :math:`(x,y) =
+(x_0,y_0)`.  To simplify the formulas below, these options are set to
+zero (their default values).
+
+Because the projection is conformal, the formulation is most
+conveniently given in terms of complex numbers.  In particular, the
+unscaled projected coordinates :math:`\eta` (proportional to the
+easting, :math:`x`) and :math:`\xi` (proportional to the northing,
+:math:`y`) are combined into the single complex quantity :math:`\zeta =
+\xi + i\eta`, where :math:`i=\sqrt{-1}`.  Then any analytic function
+:math:`f(\zeta)` defines a conformal mapping (this follows from the
+Cauchy-Riemann conditions).
 
-:math:`\phi_0` is the latitude of origin that match the center of the map. It can be set with ``+lat_0``.
+Spherical form
+**************
 
-:math:`k_0` is the scale factor at the natural origin (on the central meridian). It can be set with ``+k_0``.
+Because the full (ellipsoidal) projection includes the spherical
+projection as one of the components, we present the spherial form first
+with the coordinates tagged with primes, :math:`\phi'`,
+:math:`\lambda'`, :math:`\zeta' = \xi' + i\eta'`, :math:`x'`,
+:math:`y'`, so that they can be distinguished from the corresponding
+ellipsoidal coordinates (without the primes).  The projected coordinates
+for the sphere are given by
 
-:math:`M(\phi)` is the meridional distance.
+.. math::
 
-Spherical form
-**************
+   x' = k_0 R \eta';\qquad y' = k_0 R \xi'
 
 Forward projection
 ==================
 
 .. math::
 
-   B = \cos \phi \sin \lambda
+   \xi' = \tan^{-1}\biggl(\frac{\tan\phi'}{\cos\lambda'}\biggr)
 
 .. math::
 
-   x = \frac{k_0}{2} \ln(\frac{1+B}{1-B})
+   \eta' = \sinh^{-1}\biggl(\frac{\sin\lambda'}
+   {\sqrt{\tan^2\phi' + \cos^2\lambda'}}\biggr)
+
+
+Inverse projection
+==================
 
 .. math::
 
-   y = k_0 ( \arctan(\frac{\tan(\phi)}{\cos \lambda}) - \phi_0)
+  \phi' = \tan^{-1}\biggl(\frac{\sin\xi'}
+  {\sqrt{\sinh^2\eta' + \cos^2\xi'}}\biggr)
 
+.. math::
+
+  \lambda' = \tan^{-1}\biggl(\frac{\sinh\eta'}{\cos\xi'}\biggr)
 
-Inverse projection
-==================
+
+Ellipsoidal form
+****************
+
+The projected coordinates are given by
 
 .. math::
 
-  D = \frac{y}{k_0} + \phi_0
+   \zeta = \xi + i\eta;\qquad x = k_0 A \eta;\qquad y = k_0 A \xi
 
 .. math::
 
-  x' = \frac{x}{k_0}
+   A = \frac a{1+n}\biggl(1 + \frac14 n^2 + \frac1{64} n^4 +
+   \frac1{256}n^6\biggr)
+
+The series for conversion between ellipsoidal and spherical geographic
+coordinates and ellipsoidal and spherical projected coordinates are
+given in matrix notation where :math:`\mathbf S(\theta)` and
+:math:`\mathbf N` are the row and column vectors of length 6
 
 .. math::
 
-  \phi = \arcsin(\frac{\sin D}{\cosh x'})
+   \mathbf S(\theta) = \begin{pmatrix}
+   \sin  2\theta &
+   \sin  4\theta &
+   \sin  6\theta &
+   \sin  8\theta &
+   \sin 10\theta &
+   \sin 12\theta
+   \end{pmatrix}
 
 .. math::
 
-  \lambda = \arctan(\frac{\sinh x'}{\cos D})
+   \mathbf N = \begin{pmatrix}
+   n \\ n^2 \\ n^3\\ n^4 \\ n^5 \\ n^6
+   \end{pmatrix}
 
+and :math:`\mathsf C_{\alpha,\beta}` are upper triangular
+:math:`6\times6` matrices.
 
-Ellipsoidal form
-****************
+Relation between geographic coordinates
+=======================================
+
+.. math::
+
+  \lambda' = \lambda
+
+.. math::
+
+   \phi' = \tan^{-1}\sinh\bigl(\sinh^{-1}\tan\phi
+   - e \tanh^{-1}(e\sin\phi)\bigr)
+
+Instead of using this analytical formula for :math:`\phi'`, the
+conversions between :math:`\phi` and :math:`\phi'` use the series
+approximations:
+
+.. math::
+
+  \phi' = \phi + \mathbf S(\phi) \cdot \mathsf C_{\chi,\phi} \cdot \mathbf N
+
+.. math::
+
+  \phi = \phi' + \mathbf S(\phi') \cdot \mathsf C_{\phi,\chi} \cdot \mathbf N
+
+.. math::
+
+  \mathsf C_{\chi,\phi} = \begin{pmatrix}
+  -2& \frac{2}{3}& \frac{4}{3}& -\frac{82}{45}& \frac{32}{45}& \frac{4642}{4725} \\
+  & \frac{5}{3}& -\frac{16}{15}& -\frac{13}{9}& \frac{904}{315}& -\frac{1522}{945} \\
+  & & -\frac{26}{15}& \frac{34}{21}& \frac{8}{5}& -\frac{12686}{2835} \\
+  & & & \frac{1237}{630}& -\frac{12}{5}& -\frac{24832}{14175} \\
+  & & & & -\frac{734}{315}& \frac{109598}{31185} \\
+  & & & & & \frac{444337}{155925} \\
+  \end{pmatrix}
+
+.. math::
+
+  \mathsf C_{\phi,\chi} = \begin{pmatrix}
+  2& -\frac{2}{3}& -2& \frac{116}{45}& \frac{26}{45}& -\frac{2854}{675} \\
+  & \frac{7}{3}& -\frac{8}{5}& -\frac{227}{45}& \frac{2704}{315}& \frac{2323}{945} \\
+  & & \frac{56}{15}& -\frac{136}{35}& -\frac{1262}{105}& \frac{73814}{2835} \\
+  & & & \frac{4279}{630}& -\frac{332}{35}& -\frac{399572}{14175} \\
+  & & & & \frac{4174}{315}& -\frac{144838}{6237} \\
+  & & & & & \frac{601676}{22275} \\
+  \end{pmatrix}
+
+Here :math:`\phi'` is the conformal latitude (sometimes denoted by
+:math:`\chi`) and :math:`\mathsf C_{\chi,\phi}` and :math:`\mathsf
+C_{\phi,\chi}` are the coefficients in the trigonometric series for
+converting between :math:`\phi` and :math:`\chi`.
+
+Relation between projected coordinates
+======================================
+
+.. math::
+
+  \zeta = \zeta' + \mathbf S(\zeta') \cdot \mathsf C_{\mu,\chi} \cdot \mathbf N
+
+.. math::
+
+  \zeta' = \zeta + \mathbf S(\zeta) \cdot \mathsf C_{\chi,\mu} \cdot \mathbf N
+
+.. math::
+
+  \mathsf C_{\mu,\chi} = \begin{pmatrix}
+  \frac{1}{2}& -\frac{2}{3}& \frac{5}{16}& \frac{41}{180}& -\frac{127}{288}& \frac{7891}{37800} \\
+  & \frac{13}{48}& -\frac{3}{5}& \frac{557}{1440}& \frac{281}{630}& -\frac{1983433}{1935360} \\
+  & & \frac{61}{240}& -\frac{103}{140}& \frac{15061}{26880}& \frac{167603}{181440} \\
+  & & & \frac{49561}{161280}& -\frac{179}{168}& \frac{6601661}{7257600} \\
+  & & & & \frac{34729}{80640}& -\frac{3418889}{1995840} \\
+  & & & & & \frac{212378941}{319334400} \\
+  \end{pmatrix}
+
+.. math::
+
+  \mathsf C_{\chi,\mu} = \begin{pmatrix}
+  -\frac{1}{2}& \frac{2}{3}& -\frac{37}{96}& \frac{1}{360}& \frac{81}{512}& -\frac{96199}{604800} \\
+  & -\frac{1}{48}& -\frac{1}{15}& \frac{437}{1440}& -\frac{46}{105}& \frac{1118711}{3870720} \\
+  & & -\frac{17}{480}& \frac{37}{840}& \frac{209}{4480}& -\frac{5569}{90720} \\
+  & & & -\frac{4397}{161280}& \frac{11}{504}& \frac{830251}{7257600} \\
+  & & & & -\frac{4583}{161280}& \frac{108847}{3991680} \\
+  & & & & & -\frac{20648693}{638668800} \\
+  \end{pmatrix}
+
+On the central meridian (:math:`\lambda = \lambda' = 0`), :math:`\zeta'
+= \phi'` is the conformal latitude :math:`\chi` and :math:`\zeta` plays
+the role of the rectifying latitude (sometimes denoted by :math:`\mu`).
+:math:`\mathsf C_{\mu,\chi}` and :math:`\mathsf C_{\chi,\mu}` are the
+coefficients in the trigonometric series for converting between
+:math:`\chi` and :math:`\mu`.
+
+Legacy ellipsoidal form
+***********************
 
 The formulas below describe the algorithm used when giving the
 :option:`+approx` option. They are originally from :cite:`Snyder1987`,
-but here quoted from :cite:`Evenden1995`.
-The default algorithm is given by Poder and Engsager in :cite:`Poder1998`
+but here quoted from :cite:`Evenden1995` and :cite:`Evenden2005`.  These
+are less accurate that the formulation above and are only valid within
+about 5 degrees of the central meridian.  Here :math:`M(\phi)` is the
+meridional distance.
+
 
 Forward projection
 ==================
@@ -209,7 +371,7 @@ Inverse projection
 
 .. math::
 
-  \phi_1 = M^-1(y)
+  \phi_1 = M^{-1}(y)
 
 .. math::
 
@@ -245,5 +407,4 @@ Inverse projection
 Further reading
 ###############
 
-#. `Wikipedia <https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system>`_
-#. `EPSG, POSC literature pertaining to Coordinate Conversions and Transformations including Formulas  <http://www.ihsenergy.com/epsg/guid7.pdf>`_
+#. `Wikipedia <https://en.wikipedia.org/wiki/Transverse_Mercator_projection>`_
diff --git a/docs/source/references.bib b/docs/source/references.bib
index 13853764..cda58402 100644
--- a/docs/source/references.bib
+++ b/docs/source/references.bib
@@ -84,6 +84,16 @@
   Url                      = {http://www.calcofi.org/publications/calcofireports/v20/Vol_20_Eber___Hewitt.pdf}
 }
 
+@InProceedings{Engsager2007,
+  author =	 {Karsten E. Engsager and Knud Poder},
+  title =	 {A highly accurate world wide algorithm for the
+		  transverse {M}ercator mapping (almost)},
+  booktitle =	 {Proc. XXIII Intl. Cartographic Conf. (ICC2007), Moscow},
+  pages =	 {2.1.2},
+  year =	 2007,
+  month =	 aug
+}
+
 @Manual{Evenden2005,
   Title                    = {libproj4: A Comprehensive Library of Cartographic Projection Functions (Preliminary Draft)},
   Author                   = {Gerald I. Evenden},
@@ -197,6 +207,20 @@
   Doi                      = {10.1007/s00190-012-0578-z}
 }
 
+@Article{Karney2011tm,
+  author =	 {Charles F. F. Karney},
+  title =	 {Transverse {M}ercator with an accuracy of a few
+		  nanometers},
+  journal =	 {J. Geod.},
+  year =	 2011,
+  volume =	 85,
+  number =	 8,
+  pages =	 {475-485},
+  month =	 aug,
+  eprint =	 {1002.1417},
+  doi =		 {10.1007/s00190-011-0445-3}
+}
+
 @Article{Karney2011,
   Title                    = {Geodesics on an ellipsoid of revolution},
   Author                   = {Charles F. F. Karney},
@@ -218,6 +242,17 @@
   Volume                   = {71}
 }
 
+@TechReport{Krueger1912,
+  author =	 {Johann Heinrich Louis Kr\"uger},
+  title =	 {Konforme {A}bbildung des {E}rdellipsoids in der {E}bene},
+  institution =	 {Royal Prussian Geodetic Institute},
+  address =	 {Potsdam},
+  year =	 1912,
+  number =	 52,
+  type =	 {New Series},
+  doi =		 {10.2312/GFZ.b103-krueger28}
+}
+
 @InProceedings{LambersKolb2012,
   Title                    = {Ellipsoidal Cube Maps for Accurate Rendering of Planetary-Scale Terrain Data},
   Author                   = {Lambers, Martin and Kolb, Andreas},