Make tmerc an alias for etmerc. (#1234)

* Make tmerc an alias for etmerc

This switches the algorithm used in tmerc to the Poder/Engsager
tmerc algorithm. The original tmerc algorithm of Evenden/Snyder
origin can still be accessed by adding the +approx flag when
initializing a tmerc projection. The +approx flag can also
be used when initializing UTM projections, in which case the
Evenden/Snyder algorithm is used as well.

If a tmerc projection is instantiated on a spherical earth
the Evenden/Snyder algorithm is used as well since the
Poder/Engsager algorithm is only defined on the ellipsoid.

+proj=etmerc can still be instantiated for backwards compatibility
reasons.

Co-authored-by: Kristian Evers <kristianevers@gmail.com>
Co-authored-by: Even Rouault <even.rouault@spatialys.com>

19 files changed, 497 insertions, 485 deletions

diff --git a/docs/source/operations/projections/etmerc.rst b/docs/source/operations/projections/etmerc.rst
deleted file mode 100644
index 2a007130..00000000
--- a/docs/source/operations/projections/etmerc.rst
+++ /dev/null
@@ -1,32 +0,0 @@
-.. _etmerc:
-
-********************************************************************************
-Extended Transverse Mercator
-********************************************************************************
-
-.. figure:: ./images/etmerc.png
-   :width: 500 px
-   :align: center
-   :alt:   Extended Transverse Mercator
-
-   proj-string: ``+proj=etmerc +lon_0=-40``
-
-Parameters
-################################################################################
-
-.. note:: All parameters for the projection are optional.
-
-.. include:: ../options/lon_0.rst
-
-.. include:: ../options/lat_0.rst
-
-.. include:: ../options/ellps.rst
-
-.. include:: ../options/R.rst
-
-.. include:: ../options/k_0.rst
-
-.. include:: ../options/x_0.rst
-
-.. include:: ../options/y_0.rst
-
diff --git a/docs/source/operations/projections/images/etmerc.png b/docs/source/operations/projections/images/etmerc.png
deleted file mode 100644
index 35bce083..00000000
--- a/docs/source/operations/projections/images/etmerc.png
+++ /dev/null
diff --git a/docs/source/operations/projections/index.rst b/docs/source/operations/projections/index.rst
index 572deb1a..e6ec9076 100644
--- a/docs/source/operations/projections/index.rst
+++ b/docs/source/operations/projections/index.rst
@@ -43,7 +43,6 @@ Projections map the spherical 3D space to a flat 2D space.
    eqdc
    eqearth
    euler
-   etmerc
    fahey
    fouc
    fouc_s
diff --git a/docs/source/operations/projections/tmerc.rst b/docs/source/operations/projections/tmerc.rst
index 6d1c145d..fd600001 100644
--- a/docs/source/operations/projections/tmerc.rst
+++ b/docs/source/operations/projections/tmerc.rst
@@ -68,13 +68,19 @@ Example using Gauss-Kruger on Germany area (aka EPSG:31467) ::
 Example using Gauss Boaga on Italy area (EPSG:3004) ::
 
     $ echo 15 42 | proj +proj=tmerc +lat_0=0 +lon_0=15 +k_0=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m +no_defs
-    2520000.00	4649858.60 
+    2520000.00	4649858.60
 
 Parameters
 ################################################################################
 
 .. note:: All parameters for the projection are optional.
 
+.. option:: +approx
+
+    .. versionadded:: 6.0.0
+
+    Use faster, less accurate algorithm for the Transverse Mercator.
+
 .. include:: ../options/lon_0.rst
 
 .. include:: ../options/lat_0.rst
diff --git a/docs/source/operations/projections/utm.rst b/docs/source/operations/projections/utm.rst
index ac65aa1a..8759f4eb 100644
--- a/docs/source/operations/projections/utm.rst
+++ b/docs/source/operations/projections/utm.rst
@@ -66,6 +66,12 @@ Required
 
     Add this flag when using the UTM on the southern hemisphere.
 
+.. option:: +approx
+
+    .. versionadded:: 6.0.0
+
+    Use faster, less accurate algorithm for the Transverse Mercator.
+
 Optional
 -------------------------------------------------------------------------------
 
diff --git a/include/proj/io.hpp b/include/proj/io.hpp
index aedf0457..091efcb7 100644
--- a/include/proj/io.hpp
+++ b/include/proj/io.hpp
@@ -365,7 +365,7 @@ class PROJ_GCC_DLL PROJStringFormatter {
     PROJ_DLL ~PROJStringFormatter();
     //! @endcond
 
-    PROJ_DLL void setUseETMercForTMerc(bool flag);
+    PROJ_DLL void setUseApproxTMerc(bool flag);
 
     PROJ_DLL const std::string &toString() const;
 
@@ -378,7 +378,7 @@ class PROJ_GCC_DLL PROJStringFormatter {
     PROJ_DLL void startInversion();
     PROJ_DLL void stopInversion();
     PROJ_INTERNAL bool isInverted() const;
-    PROJ_INTERNAL bool getUseETMercForTMerc(bool &settingSetOut) const;
+    PROJ_INTERNAL bool getUseApproxTMerc() const;
     PROJ_INTERNAL void setCoordinateOperationOptimizations(bool enable);
 
     PROJ_DLL void
diff --git a/src/Makefile.am b/src/Makefile.am
index 2af53e30..698baf97 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -122,7 +122,6 @@ libproj_la_SOURCES = \
 	projections/wag7.cpp \
 	projections/lcca.cpp \
 	projections/geos.cpp \
-	projections/etmerc.cpp \
 	projections/boggs.cpp \
 	projections/collg.cpp \
 	projections/comill.cpp \
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 6fafa2c8..240d11b3 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -1207,9 +1207,8 @@ const char *proj_as_wkt(PJ_CONTEXT *ctx, const PJ *obj, PJ_WKT_TYPE type,
  * @param type PROJ String version.
  * @param options NULL-terminated list of strings with "KEY=VALUE" format. or
  * NULL.
- * The currently recognized option is USE_ETMERC=YES to use
- * +proj=etmerc instead of +proj=tmerc (or USE_ETMERC=NO to disable implicit
- * use of etmerc by utm conversions)
+ * The currently recognized option is USE_APPROX_TMERC=YES to add the +approx
+ * flag to +proj=tmerc or +proj=utm
  * @return a string, or NULL in case of error.
  */
 const char *proj_as_proj_string(PJ_CONTEXT *ctx, const PJ *obj,
@@ -1243,10 +1242,8 @@ const char *proj_as_proj_string(PJ_CONTEXT *ctx, const PJ *obj,
     try {
         auto formatter = PROJStringFormatter::create(convention, dbContext);
         if (options != nullptr && options[0] != nullptr) {
-            if (ci_equal(options[0], "USE_ETMERC=YES")) {
-                formatter->setUseETMercForTMerc(true);
-            } else if (ci_equal(options[0], "USE_ETMERC=NO")) {
-                formatter->setUseETMercForTMerc(false);
+            if (ci_equal(options[0], "USE_APPROX_TMERC=YES")) {
+                formatter->setUseApproxTMerc(true);
             }
         }
         obj->lastPROJString = exportable->exportToPROJString(formatter.get());
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 8a10bc5a..d36d901f 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -5381,13 +5381,11 @@ bool Conversion::addWKTExtensionNode(io::WKTFormatter *formatter) const {
         const auto &l_method = method();
         const auto &methodName = l_method->nameStr();
         const int methodEPSGCode = l_method->getEPSGCode();
-        int zone = 0;
-        bool north = true;
-        if (l_method->getPrivate()->projMethodOverride_ == "etmerc" &&
-            !isUTM(zone, north)) {
+        if (l_method->getPrivate()->projMethodOverride_ == "tmerc approx" ||
+            l_method->getPrivate()->projMethodOverride_ == "utm approx") {
             auto projFormatter = io::PROJStringFormatter::create();
             projFormatter->setCRSExport(true);
-            projFormatter->setUseETMercForTMerc(true);
+            projFormatter->setUseApproxTMerc(true);
             formatter->startNode(io::WKTConstants::EXTENSION, false);
             formatter->addQuotedString("PROJ4");
             _exportToPROJString(projFormatter.get());
@@ -5479,17 +5477,21 @@ void Conversion::_exportToPROJString(
     const auto &convName = nameStr();
     bool bConversionDone = false;
     bool bEllipsoidParametersDone = false;
-    bool useETMerc = false;
+    bool useApprox = false;
     if (methodEPSGCode == EPSG_CODE_METHOD_TRANSVERSE_MERCATOR) {
         // Check for UTM
         int zone = 0;
         bool north = true;
-        bool etMercSettingSet = false;
-        useETMerc = formatter->getUseETMercForTMerc(etMercSettingSet) ||
-                    l_method->getPrivate()->projMethodOverride_ == "etmerc";
-        if (isUTM(zone, north) && !(etMercSettingSet && !useETMerc)) {
+        useApprox =
+            formatter->getUseApproxTMerc() ||
+            l_method->getPrivate()->projMethodOverride_ == "tmerc approx" ||
+            l_method->getPrivate()->projMethodOverride_ == "utm approx";
+        if (isUTM(zone, north)) {
             bConversionDone = true;
             formatter->addStep("utm");
+            if( useApprox ) {
+                formatter->addParam("approx");
+            }
             formatter->addParam("zone", zone);
             if (!north) {
                 formatter->addParam("south");
@@ -5662,7 +5664,10 @@ void Conversion::_exportToPROJString(
     if (!bConversionDone) {
         const MethodMapping *mapping = getMapping(l_method.get());
         if (mapping && mapping->proj_name_main) {
-            formatter->addStep(useETMerc ? "etmerc" : mapping->proj_name_main);
+            formatter->addStep(mapping->proj_name_main);
+            if (useApprox) {
+                formatter->addParam("approx");
+            }
             if (mapping->proj_name_aux) {
                 if (internal::starts_with(mapping->proj_name_aux, "axis=")) {
                     bAxisSpecFound = true;
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 431c75af..60c28201 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -4746,8 +4746,7 @@ struct PROJStringFormatter::Private {
     bool omitProjLongLatIfPossible_ = false;
     bool omitZUnitConversion_ = false;
     DatabaseContextPtr dbContext_{};
-    bool useETMercForTMerc_ = false;
-    bool useETMercForTMercSet_ = false;
+    bool useApproxTMerc_ = false;
     bool addNoDefs_ = true;
     bool coordOperationOptimizations_ = false;
     bool crsExport_ = false;
@@ -4803,11 +4802,9 @@ PROJStringFormatter::create(Convention conventionIn,
 
 // ---------------------------------------------------------------------------
 
-/** \brief Set whether Extended Transverse Mercator (etmerc) should be used
- * instead of tmerc */
-void PROJStringFormatter::setUseETMercForTMerc(bool flag) {
-    d->useETMercForTMerc_ = flag;
-    d->useETMercForTMercSet_ = true;
+/** \brief Set whether approximate Transverse Mercator or UTM should be used */
+void PROJStringFormatter::setUseApproxTMerc(bool flag) {
+    d->useApproxTMerc_ = flag;
 }
 
 // ---------------------------------------------------------------------------
@@ -5256,9 +5253,8 @@ PROJStringFormatter::Convention PROJStringFormatter::convention() const {
 
 // ---------------------------------------------------------------------------
 
-bool PROJStringFormatter::getUseETMercForTMerc(bool &settingSetOut) const {
-    settingSetOut = d->useETMercForTMercSet_;
-    return d->useETMercForTMerc_;
+bool PROJStringFormatter::getUseApproxTMerc() const {
+    return d->useApproxTMerc_;
 }
 
 // ---------------------------------------------------------------------------
@@ -7081,8 +7077,10 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                                 : UnitOfMeasure::NONE)));
         }
 
-        if (step.name == "etmerc") {
-            methodMap.set("proj_method", "etmerc");
+        if (step.name == "tmerc" && hasParamValue(step, "approx")) {
+            methodMap.set("proj_method", "tmerc approx");
+        } else if (step.name == "utm" && hasParamValue(step, "approx")) {
+            methodMap.set("proj_method", "utm approx");
         }
 
         conv = Conversion::create(mapWithUnknownName, methodMap, parameters,
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index cd380b24..f28a1d68 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -117,7 +117,6 @@ SET(SRC_LIBPROJ_PROJECTIONS
         projections/wag7.cpp
         projections/lcca.cpp
         projections/geos.cpp
-        projections/etmerc.cpp
         projections/boggs.cpp
         projections/collg.cpp
         projections/comill.cpp
diff --git a/src/projections/etmerc.cpp b/src/projections/etmerc.cpp
deleted file mode 100644
index e75bc168..00000000
--- a/src/projections/etmerc.cpp
+++ /dev/null
@@ -1,362 +0,0 @@
-/*
-** libproj -- library of cartographic projections
-**
-** Copyright (c) 2008   Gerald I. Evenden
-*/
-
-/*
-** Permission is hereby granted, free of charge, to any person obtaining
-** a copy of this software and associated documentation files (the
-** "Software"), to deal in the Software without restriction, including
-** without limitation the rights to use, copy, modify, merge, publish,
-** distribute, sublicense, and/or sell copies of the Software, and to
-** permit persons to whom the Software is furnished to do so, subject to
-** the following conditions:
-**
-** The above copyright notice and this permission notice shall be
-** included in all copies or substantial portions of the Software.
-**
-** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
-** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
-** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
-** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-
-/* The code in this file is largly based upon procedures:
- *
- * Written by: Knud Poder and Karsten Engsager
- *
- * Based on math from: R.Koenig and K.H. Weise, "Mathematische
- * Grundlagen der hoeheren Geodaesie und Kartographie,
- * Springer-Verlag, Berlin/Goettingen" Heidelberg, 1951.
- *
- * Modified and used here by permission of Reference Networks
- * Division, Kort og Matrikelstyrelsen (KMS), Copenhagen, Denmark
- *
-*/
-
-#define PJ_LIB__
-
-#include <errno.h>
-
-#include "proj.h"
-#include "proj_internal.h"
-#include "proj_math.h"
-
-
-namespace { // anonymous namespace
-struct pj_opaque {
-    double    Qn;    /* Merid. quad., scaled to the projection */ \
-    double    Zb;    /* Radius vector in polar coord. systems  */ \
-    double    cgb[6]; /* Constants for Gauss -> Geo lat */ \
-    double    cbg[6]; /* Constants for Geo lat -> Gauss */ \
-    double    utg[6]; /* Constants for transv. merc. -> geo */ \
-    double    gtu[6]; /* Constants for geo -> transv. merc. */
-};
-} // anonymous namespace
-
-PROJ_HEAD(etmerc, "Extended Transverse Mercator")
-    "\n\tCyl, Sph\n\tlat_ts=(0)\nlat_0=(0)";
-PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)")
-    "\n\tCyl, Sph\n\tzone= south";
-
-#define PROJ_ETMERC_ORDER 6
-
-#ifdef _GNU_SOURCE
-    inline
-#endif
-static double gatg(double *p1, int len_p1, double B) {
-    double *p;
-    double h = 0, h1, h2 = 0, cos_2B;
-
-    cos_2B = 2*cos(2*B);
-    p = p1 + len_p1;
-    h1 = *--p;
-    while (p - p1) {
-        h = -h2 + cos_2B*h1 + *--p;
-        h2 = h1;
-        h1 = h;
-    }
-    return (B + h*sin(2*B));
-}
-
-/* Complex Clenshaw summation */
-#ifdef _GNU_SOURCE
-    inline
-#endif
-static double clenS(double *a, int size, double arg_r, double arg_i, double *R, double *I) {
-    double      *p, r, i, hr, hr1, hr2, hi, hi1, hi2;
-    double      sin_arg_r, cos_arg_r, sinh_arg_i, cosh_arg_i;
-
-    /* arguments */
-    p = a + size;
-#ifdef _GNU_SOURCE
-    sincos(arg_r, &sin_arg_r, &cos_arg_r);
-#else
-    sin_arg_r  = sin(arg_r);
-    cos_arg_r  = cos(arg_r);
-#endif
-    sinh_arg_i = sinh(arg_i);
-    cosh_arg_i = cosh(arg_i);
-    r          =  2*cos_arg_r*cosh_arg_i;
-    i          = -2*sin_arg_r*sinh_arg_i;
-
-    /* summation loop */
-    hi1 = hr1 = hi = 0;
-    hr = *--p;
-    for (; a - p;) {
-        hr2 = hr1;
-        hi2 = hi1;
-        hr1 = hr;
-        hi1 = hi;
-        hr  = -hr2 + r*hr1 - i*hi1 + *--p;
-        hi  = -hi2 + i*hr1 + r*hi1;
-    }
-
-    r   = sin_arg_r*cosh_arg_i;
-    i   = cos_arg_r*sinh_arg_i;
-    *R  = r*hr - i*hi;
-    *I  = r*hi + i*hr;
-    return *R;
-}
-
-
-/* Real Clenshaw summation */
-static double clens(double *a, int size, double arg_r) {
-    double      *p, r, hr, hr1, hr2, cos_arg_r;
-
-    p = a + size;
-    cos_arg_r  = cos(arg_r);
-    r          =  2*cos_arg_r;
-
-    /* summation loop */
-    hr1 = 0;
-    hr = *--p;
-    for (; a - p;) {
-        hr2 = hr1;
-        hr1 = hr;
-        hr  = -hr2 + r*hr1 + *--p;
-    }
-    return sin (arg_r)*hr;
-}
-
-
-
-static PJ_XY e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
-    PJ_XY xy = {0.0,0.0};
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-    double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
-    double Cn = lp.phi, Ce = lp.lam;
-
-    /* ell. LAT, LNG -> Gaussian LAT, LNG */
-    Cn  = gatg (Q->cbg, PROJ_ETMERC_ORDER, Cn);
-    /* Gaussian LAT, LNG -> compl. sph. LAT */
-#ifdef _GNU_SOURCE
-    sincos (Cn, &sin_Cn, &cos_Cn);
-    sincos (Ce, &sin_Ce, &cos_Ce);
-#else
-    sin_Cn = sin (Cn);
-    cos_Cn = cos (Cn);
-    sin_Ce = sin (Ce);
-    cos_Ce = cos (Ce);
-#endif
-
-    Cn     = atan2 (sin_Cn, cos_Ce*cos_Cn);
-    Ce     = atan2 (sin_Ce*cos_Cn,  hypot (sin_Cn, cos_Cn*cos_Ce));
-
-    /* compl. sph. N, E -> ell. norm. N, E */
-    Ce  = asinh ( tan (Ce) );     /* Replaces: Ce  = log(tan(FORTPI + Ce*0.5)); */
-    Cn += clenS (Q->gtu, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
-    Ce += dCe;
-    if (fabs (Ce) <= 2.623395162778) {
-        xy.y  = Q->Qn * Cn + Q->Zb;  /* Northing */
-        xy.x  = Q->Qn * Ce;          /* Easting  */
-    } else
-        xy.x = xy.y = HUGE_VAL;
-    return xy;
-}
-
-
-
-static PJ_LP e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
-    PJ_LP lp = {0.0,0.0};
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-    double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
-    double Cn = xy.y, Ce = xy.x;
-
-    /* normalize N, E */
-    Cn = (Cn - Q->Zb)/Q->Qn;
-    Ce = Ce/Q->Qn;
-
-    if (fabs(Ce) <= 2.623395162778) { /* 150 degrees */
-        /* norm. N, E -> compl. sph. LAT, LNG */
-        Cn += clenS(Q->utg, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
-        Ce += dCe;
-        Ce = atan (sinh (Ce)); /* Replaces: Ce = 2*(atan(exp(Ce)) - FORTPI); */
-        /* compl. sph. LAT -> Gaussian LAT, LNG */
-#ifdef _GNU_SOURCE
-        sincos (Cn, &sin_Cn, &cos_Cn);
-        sincos (Ce, &sin_Ce, &cos_Ce);
-#else
-        sin_Cn = sin (Cn);
-        cos_Cn = cos (Cn);
-        sin_Ce = sin (Ce);
-        cos_Ce = cos (Ce);
-#endif
-        Ce     = atan2 (sin_Ce, cos_Ce*cos_Cn);
-        Cn     = atan2 (sin_Cn*cos_Ce,  hypot (sin_Ce, cos_Ce*cos_Cn));
-        /* Gaussian LAT, LNG -> ell. LAT, LNG */
-        lp.phi = gatg (Q->cgb,  PROJ_ETMERC_ORDER, Cn);
-        lp.lam = Ce;
-    }
-    else
-        lp.phi = lp.lam = HUGE_VAL;
-    return lp;
-}
-
-
-static PJ *setup(PJ *P) { /* general initialization */
-    double f, n, np, Z;
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-
-    if (P->es <= 0) {
-        return pj_default_destructor(P, PJD_ERR_ELLIPSOID_USE_REQUIRED);
-    }
-
-    /* flattening */
-    f = P->es / (1 + sqrt (1 -  P->es)); /* Replaces: f = 1 - sqrt(1-P->es); */
-
-    /* third flattening */
-    np = n = f/(2 - f);
-
-    /* COEF. OF TRIG SERIES GEO <-> GAUSS */
-    /* cgb := Gaussian -> Geodetic, KW p190 - 191 (61) - (62) */
-    /* cbg := Geodetic -> Gaussian, KW p186 - 187 (51) - (52) */
-    /* PROJ_ETMERC_ORDER = 6th degree : Engsager and Poder: ICC2007 */
-
-    Q->cgb[0] = n*( 2 + n*(-2/3.0  + n*(-2      + n*(116/45.0 + n*(26/45.0 +
-                n*(-2854/675.0 ))))));
-    Q->cbg[0] = n*(-2 + n*( 2/3.0  + n*( 4/3.0  + n*(-82/45.0 + n*(32/45.0 +
-                n*( 4642/4725.0))))));
-    np     *= n;
-    Q->cgb[1] = np*(7/3.0 + n*( -8/5.0  + n*(-227/45.0 + n*(2704/315.0 +
-                n*( 2323/945.0)))));
-    Q->cbg[1] = np*(5/3.0 + n*(-16/15.0 + n*( -13/9.0  + n*( 904/315.0 +
-                n*(-1522/945.0)))));
-    np     *= n;
-    /* n^5 coeff corrected from 1262/105 -> -1262/105 */
-    Q->cgb[2] = np*( 56/15.0  + n*(-136/35.0 + n*(-1262/105.0 +
-                n*( 73814/2835.0))));
-    Q->cbg[2] = np*(-26/15.0  + n*(  34/21.0 + n*(    8/5.0   +
-                n*(-12686/2835.0))));
-    np     *= n;
-    /* n^5 coeff corrected from 322/35 -> 332/35 */
-    Q->cgb[3] = np*(4279/630.0 + n*(-332/35.0 + n*(-399572/14175.0)));
-    Q->cbg[3] = np*(1237/630.0 + n*( -12/5.0  + n*( -24832/14175.0)));
-    np     *= n;
-    Q->cgb[4] = np*(4174/315.0 + n*(-144838/6237.0 ));
-    Q->cbg[4] = np*(-734/315.0 + n*( 109598/31185.0));
-    np     *= n;
-    Q->cgb[5] = np*(601676/22275.0 );
-    Q->cbg[5] = np*(444337/155925.0);
-
-    /* Constants of the projections */
-    /* Transverse Mercator (UTM, ITM, etc) */
-    np = n*n;
-    /* Norm. mer. quad, K&W p.50 (96), p.19 (38b), p.5 (2) */
-    Q->Qn = P->k0/(1 + n) * (1 + np*(1/4.0 + np*(1/64.0 + np/256.0)));
-    /* coef of trig series */
-    /* utg := ell. N, E -> sph. N, E,  KW p194 (65) */
-    /* gtu := sph. N, E -> ell. N, E,  KW p196 (69) */
-    Q->utg[0] = n*(-0.5  + n*( 2/3.0 + n*(-37/96.0 + n*( 1/360.0 +
-                n*(  81/512.0 + n*(-96199/604800.0))))));
-    Q->gtu[0] = n*( 0.5  + n*(-2/3.0 + n*(  5/16.0 + n*(41/180.0 +
-                n*(-127/288.0 + n*(  7891/37800.0 ))))));
-    Q->utg[1] = np*(-1/48.0 + n*(-1/15.0 + n*(437/1440.0 + n*(-46/105.0 +
-                n*( 1118711/3870720.0)))));
-    Q->gtu[1] = np*(13/48.0 + n*(-3/5.0  + n*(557/1440.0 + n*(281/630.0 +
-                n*(-1983433/1935360.0)))));
-    np      *= n;
-    Q->utg[2] = np*(-17/480.0 + n*(  37/840.0 + n*(  209/4480.0  +
-                n*( -5569/90720.0 ))));
-    Q->gtu[2] = np*( 61/240.0 + n*(-103/140.0 + n*(15061/26880.0 +
-                n*(167603/181440.0))));
-    np      *= n;
-    Q->utg[3] = np*(-4397/161280.0 + n*(  11/504.0 + n*( 830251/7257600.0)));
-    Q->gtu[3] = np*(49561/161280.0 + n*(-179/168.0 + n*(6601661/7257600.0)));
-    np     *= n;
-    Q->utg[4] = np*(-4583/161280.0 + n*(  108847/3991680.0));
-    Q->gtu[4] = np*(34729/80640.0  + n*(-3418889/1995840.0));
-    np     *= n;
-    Q->utg[5] = np*(-20648693/638668800.0);
-    Q->gtu[5] = np*(212378941/319334400.0);
-
-    /* Gaussian latitude value of the origin latitude */
-    Z = gatg (Q->cbg, PROJ_ETMERC_ORDER, P->phi0);
-
-    /* Origin northing minus true northing at the origin latitude */
-    /* i.e. true northing = N - P->Zb                         */
-    Q->Zb  = - Q->Qn*(Z + clens(Q->gtu, PROJ_ETMERC_ORDER, 2*Z));
-    P->inv = e_inverse;
-    P->fwd = e_forward;
-    return P;
-}
-
-
-
-PJ *PROJECTION(etmerc) {
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
-    if (nullptr==Q)
-        return pj_default_destructor (P, ENOMEM);
-    P->opaque = Q;
-   return setup (P);
-}
-
-
-
-/* utm uses etmerc for the underlying projection */
-
-
-PJ *PROJECTION(utm) {
-    long zone;
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
-    if (nullptr==Q)
-        return pj_default_destructor (P, ENOMEM);
-    P->opaque = Q;
-
-    if (P->es == 0.0) {
-        proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED);
-        return pj_default_destructor(P, ENOMEM);
-    }
-    if (P->lam0 < -1000.0 || P->lam0 > 1000.0) {
-        return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE);
-    }
-
-    P->y0 = pj_param (P->ctx, P->params, "bsouth").i ? 10000000. : 0.;
-    P->x0 = 500000.;
-    if (pj_param (P->ctx, P->params, "tzone").i) /* zone input ? */
-    {
-        zone = pj_param(P->ctx, P->params, "izone").i;
-        if (zone > 0 && zone <= 60)
-            --zone;
-        else {
-            return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE);
-        }
-    }
-    else /* nearest central meridian input */
-    {
-        zone = lround((floor ((adjlon (P->lam0) + M_PI) * 30. / M_PI)));
-        if (zone < 0)
-            zone = 0;
-        else if (zone >= 60)
-            zone = 59;
-    }
-    P->lam0 = (zone + .5) * M_PI / 30. - M_PI;
-    P->k0 = 0.9996;
-    P->phi0 = 0.;
-
-    return setup (P);
-}
diff --git a/src/projections/tmerc.cpp b/src/projections/tmerc.cpp
index d1938116..c91c5174 100644
--- a/src/projections/tmerc.cpp
+++ b/src/projections/tmerc.cpp
@@ -1,3 +1,15 @@
+/*
+*                   Transverse Mercator implementations
+*
+* In this file two transverse mercator implementations are found. One of Gerald
+* Evenden/John Snyder origin and one of Knud Poder/Karsten Engsager origin. The
+* former is regarded as "approximate" in the following and the latter is "exact".
+* This word choice has been made to distinguish between the two algorithms, where
+* the Evenden/Snyder implementation is the faster, less accurate implementation
+* and the Poder/Engsager algorithm is a slightly slower, but more accurate
+* implementation.
+*/
+
 #define PJ_LIB__
 
 #include <errno.h>
@@ -5,18 +17,32 @@
 
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 
-PROJ_HEAD(tmerc, "Transverse Mercator") "\n\tCyl, Sph&Ell";
 
+PROJ_HEAD(tmerc, "Transverse Mercator") "\n\tCyl, Sph&Ell\n\tapprox";
+PROJ_HEAD(etmerc, "Extended Transverse Mercator") "\n\tCyl, Sph";
+PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)") "\n\tCyl, Sph\n\tzone= south approx";
 
 namespace { // anonymous namespace
-struct pj_opaque {
+struct pj_opaque_approx {
     double  esp;
     double  ml0;
     double  *en;
 };
+
+struct pj_opaque_exact {
+    double    Qn;     /* Merid. quad., scaled to the projection */
+    double    Zb;     /* Radius vector in polar coord. systems  */
+    double    cgb[6]; /* Constants for Gauss -> Geo lat */
+    double    cbg[6]; /* Constants for Geo lat -> Gauss */
+    double    utg[6]; /* Constants for transv. merc. -> geo */
+    double    gtu[6]; /* Constants for geo -> transv. merc. */
+};
+
 } // anonymous namespace
 
+/* Constants for "approximate" transverse mercator */
 #define EPS10   1.e-10
 #define FC1 1.
 #define FC2 .5
@@ -27,10 +53,18 @@ struct pj_opaque {
 #define FC7 .02380952380952380952
 #define FC8 .01785714285714285714
 
+/* Constant for "exact" transverse mercator */
+#define PROJ_ETMERC_ORDER 6
 
-static PJ_XY e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
+/*****************************************************************************/
+//
+//                  Approximate Transverse Mercator functions
+//
+/*****************************************************************************/
+
+static PJ_XY approx_e_fwd (PJ_LP lp, PJ *P) {
     PJ_XY xy = {0.0, 0.0};
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+    struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(P->opaque);
     double al, als, n, cosphi, sinphi, t;
 
     /*
@@ -70,7 +104,7 @@ static PJ_XY e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
 }
 
 
-static PJ_XY s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
+static PJ_XY approx_s_fwd (PJ_LP lp, PJ *P) {
     PJ_XY xy = {0.0,0.0};
     double b, cosphi;
 
@@ -95,7 +129,7 @@ static PJ_XY s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
         return xy;
     }
 
-    xy.x = static_cast<struct pj_opaque*>(P->opaque)->ml0 * log ((1. + b) / (1. - b));
+    xy.x = static_cast<struct pj_opaque_approx*>(P->opaque)->ml0 * log ((1. + b) / (1. - b));
     xy.y = cosphi * cos (lp.lam) / sqrt (1. - b * b);
 
     b = fabs ( xy.y );
@@ -110,14 +144,14 @@ static PJ_XY s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
 
     if (lp.phi < 0.)
         xy.y = -xy.y;
-    xy.y = static_cast<struct pj_opaque*>(P->opaque)->esp * (xy.y - P->phi0);
+    xy.y = static_cast<struct pj_opaque_approx*>(P->opaque)->esp * (xy.y - P->phi0);
     return xy;
 }
 
 
-static PJ_LP e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+static PJ_LP approx_e_inv (PJ_XY xy, PJ *P) {
     PJ_LP lp = {0.0,0.0};
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+    struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(P->opaque);
     double n, con, cosphi, d, ds, sinphi, t;
 
     lp.phi = pj_inv_mlfn(P->ctx, Q->ml0 + xy.y / P->k0, P->es, Q->en);
@@ -149,13 +183,13 @@ static PJ_LP e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
 }
 
 
-static PJ_LP s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
+static PJ_LP approx_s_inv (PJ_XY xy, PJ *P) {
     PJ_LP lp = {0.0, 0.0};
     double h, g;
 
-    h = exp(xy.x / static_cast<struct pj_opaque*>(P->opaque)->esp);
+    h = exp(xy.x / static_cast<struct pj_opaque_approx*>(P->opaque)->esp);
     g = .5 * (h - 1. / h);
-    h = cos (P->phi0 + xy.y / static_cast<struct pj_opaque*>(P->opaque)->esp);
+    h = cos (P->phi0 + xy.y / static_cast<struct pj_opaque_approx*>(P->opaque)->esp);
     lp.phi = asin(sqrt((1. - h * h) / (1. + g * g)));
 
     /* Make sure that phi is on the correct hemisphere when false northing is used */
@@ -166,45 +200,386 @@ static PJ_LP s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
 }
 
 
-static PJ *destructor(PJ *P, int errlev) {                       /* Destructor */
+static PJ *destructor_approx(PJ *P, int errlev) {
     if (nullptr==P)
         return nullptr;
 
     if (nullptr==P->opaque)
         return pj_default_destructor(P, errlev);
 
-    pj_dealloc (static_cast<struct pj_opaque*>(P->opaque)->en);
+    pj_dealloc (static_cast<struct pj_opaque_approx*>(P->opaque)->en);
     return pj_default_destructor(P, errlev);
 }
 
 
-static PJ *setup(PJ *P) {                   /* general initialization */
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+static PJ *setup_approx(PJ *P) {
+    struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(P->opaque);
+
+    P->destructor = destructor_approx;
+
     if (P->es != 0.0) {
         if (!(Q->en = pj_enfn(P->es)))
             return pj_default_destructor(P, ENOMEM);
 
         Q->ml0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), Q->en);
         Q->esp = P->es / (1. - P->es);
-        P->inv = e_inverse;
-        P->fwd = e_forward;
+        P->inv = approx_e_inv;
+        P->fwd = approx_e_fwd;
     } else {
         Q->esp = P->k0;
         Q->ml0 = .5 * Q->esp;
-        P->inv = s_inverse;
-        P->fwd = s_forward;
+        P->inv = approx_s_inv;
+        P->fwd = approx_s_fwd;
     }
     return P;
 }
 
 
+
+/*****************************************************************************/
+//
+//                  Exact Transverse Mercator functions
+//
+//
+// The code in this file is largly based upon procedures:
+//
+// Written by: Knud Poder and Karsten Engsager
+//
+// Based on math from: R.Koenig and K.H. Weise, "Mathematische
+// Grundlagen der hoeheren Geodaesie und Kartographie,
+// Springer-Verlag, Berlin/Goettingen" Heidelberg, 1951.
+//
+// Modified and used here by permission of Reference Networks
+// Division, Kort og Matrikelstyrelsen (KMS), Copenhagen, Denmark
+//
+/*****************************************************************************/
+
+/* Helper functios for "exact" transverse mercator */
+#ifdef _GNU_SOURCE
+    inline
+#endif
+static double gatg(double *p1, int len_p1, double B) {
+    double *p;
+    double h = 0, h1, h2 = 0, cos_2B;
+
+    cos_2B = 2*cos(2*B);
+    p = p1 + len_p1;
+    h1 = *--p;
+    while (p - p1) {
+        h = -h2 + cos_2B*h1 + *--p;
+        h2 = h1;
+        h1 = h;
+    }
+    return (B + h*sin(2*B));
+}
+
+/* Complex Clenshaw summation */
+#ifdef _GNU_SOURCE
+    inline
+#endif
+static double clenS(double *a, int size, double arg_r, double arg_i, double *R, double *I) {
+    double      *p, r, i, hr, hr1, hr2, hi, hi1, hi2;
+    double      sin_arg_r, cos_arg_r, sinh_arg_i, cosh_arg_i;
+
+    /* arguments */
+    p = a + size;
+#ifdef _GNU_SOURCE
+    sincos(arg_r, &sin_arg_r, &cos_arg_r);
+#else
+    sin_arg_r  = sin(arg_r);
+    cos_arg_r  = cos(arg_r);
+#endif
+    sinh_arg_i = sinh(arg_i);
+    cosh_arg_i = cosh(arg_i);
+    r          =  2*cos_arg_r*cosh_arg_i;
+    i          = -2*sin_arg_r*sinh_arg_i;
+
+    /* summation loop */
+    hi1 = hr1 = hi = 0;
+    hr = *--p;
+    for (; a - p;) {
+        hr2 = hr1;
+        hi2 = hi1;
+        hr1 = hr;
+        hi1 = hi;
+        hr  = -hr2 + r*hr1 - i*hi1 + *--p;
+        hi  = -hi2 + i*hr1 + r*hi1;
+    }
+
+    r   = sin_arg_r*cosh_arg_i;
+    i   = cos_arg_r*sinh_arg_i;
+    *R  = r*hr - i*hi;
+    *I  = r*hi + i*hr;
+    return *R;
+}
+
+
+/* Real Clenshaw summation */
+static double clens(double *a, int size, double arg_r) {
+    double      *p, r, hr, hr1, hr2, cos_arg_r;
+
+    p = a + size;
+    cos_arg_r  = cos(arg_r);
+    r          =  2*cos_arg_r;
+
+    /* summation loop */
+    hr1 = 0;
+    hr = *--p;
+    for (; a - p;) {
+        hr2 = hr1;
+        hr1 = hr;
+        hr  = -hr2 + r*hr1 + *--p;
+    }
+    return sin (arg_r)*hr;
+}
+
+/* Ellipsoidal, forward */
+static PJ_XY exact_e_fwd (PJ_LP lp, PJ *P) {
+    PJ_XY xy = {0.0,0.0};
+    struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(P->opaque);
+    double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
+    double Cn = lp.phi, Ce = lp.lam;
+
+    /* ell. LAT, LNG -> Gaussian LAT, LNG */
+    Cn  = gatg (Q->cbg, PROJ_ETMERC_ORDER, Cn);
+    /* Gaussian LAT, LNG -> compl. sph. LAT */
+#ifdef _GNU_SOURCE
+    sincos (Cn, &sin_Cn, &cos_Cn);
+    sincos (Ce, &sin_Ce, &cos_Ce);
+#else
+    sin_Cn = sin (Cn);
+    cos_Cn = cos (Cn);
+    sin_Ce = sin (Ce);
+    cos_Ce = cos (Ce);
+#endif
+
+    Cn     = atan2 (sin_Cn, cos_Ce*cos_Cn);
+    Ce     = atan2 (sin_Ce*cos_Cn,  hypot (sin_Cn, cos_Cn*cos_Ce));
+
+    /* compl. sph. N, E -> ell. norm. N, E */
+    Ce  = asinh ( tan (Ce) );     /* Replaces: Ce  = log(tan(FORTPI + Ce*0.5)); */
+    Cn += clenS (Q->gtu, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
+    Ce += dCe;
+    if (fabs (Ce) <= 2.623395162778) {
+        xy.y  = Q->Qn * Cn + Q->Zb;  /* Northing */
+        xy.x  = Q->Qn * Ce;          /* Easting  */
+    } else
+        xy.x = xy.y = HUGE_VAL;
+    return xy;
+}
+
+
+/* Ellipsoidal, inverse */
+static PJ_LP exact_e_inv (PJ_XY xy, PJ *P) {
+    PJ_LP lp = {0.0,0.0};
+    struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(P->opaque);
+    double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
+    double Cn = xy.y, Ce = xy.x;
+
+    /* normalize N, E */
+    Cn = (Cn - Q->Zb)/Q->Qn;
+    Ce = Ce/Q->Qn;
+
+    if (fabs(Ce) <= 2.623395162778) { /* 150 degrees */
+        /* norm. N, E -> compl. sph. LAT, LNG */
+        Cn += clenS(Q->utg, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
+        Ce += dCe;
+        Ce = atan (sinh (Ce)); /* Replaces: Ce = 2*(atan(exp(Ce)) - FORTPI); */
+        /* compl. sph. LAT -> Gaussian LAT, LNG */
+#ifdef _GNU_SOURCE
+        sincos (Cn, &sin_Cn, &cos_Cn);
+        sincos (Ce, &sin_Ce, &cos_Ce);
+#else
+        sin_Cn = sin (Cn);
+        cos_Cn = cos (Cn);
+        sin_Ce = sin (Ce);
+        cos_Ce = cos (Ce);
+#endif
+        Ce     = atan2 (sin_Ce, cos_Ce*cos_Cn);
+        Cn     = atan2 (sin_Cn*cos_Ce,  hypot (sin_Ce, cos_Ce*cos_Cn));
+        /* Gaussian LAT, LNG -> ell. LAT, LNG */
+        lp.phi = gatg (Q->cgb,  PROJ_ETMERC_ORDER, Cn);
+        lp.lam = Ce;
+    }
+    else
+        lp.phi = lp.lam = HUGE_VAL;
+    return lp;
+}
+
+static PJ *setup_exact(PJ *P) {
+    double f, n, np, Z;
+    struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(P->opaque);
+
+    if (P->es <= 0) {
+        return pj_default_destructor(P, PJD_ERR_ELLIPSOID_USE_REQUIRED);
+    }
+
+    /* flattening */
+    f = P->es / (1 + sqrt (1 -  P->es)); /* Replaces: f = 1 - sqrt(1-P->es); */
+
+    /* third flattening */
+    np = n = f/(2 - f);
+
+    /* COEF. OF TRIG SERIES GEO <-> GAUSS */
+    /* cgb := Gaussian -> Geodetic, KW p190 - 191 (61) - (62) */
+    /* cbg := Geodetic -> Gaussian, KW p186 - 187 (51) - (52) */
+    /* PROJ_ETMERC_ORDER = 6th degree : Engsager and Poder: ICC2007 */
+
+    Q->cgb[0] = n*( 2 + n*(-2/3.0  + n*(-2      + n*(116/45.0 + n*(26/45.0 +
+                n*(-2854/675.0 ))))));
+    Q->cbg[0] = n*(-2 + n*( 2/3.0  + n*( 4/3.0  + n*(-82/45.0 + n*(32/45.0 +
+                n*( 4642/4725.0))))));
+    np     *= n;
+    Q->cgb[1] = np*(7/3.0 + n*( -8/5.0  + n*(-227/45.0 + n*(2704/315.0 +
+                n*( 2323/945.0)))));
+    Q->cbg[1] = np*(5/3.0 + n*(-16/15.0 + n*( -13/9.0  + n*( 904/315.0 +
+                n*(-1522/945.0)))));
+    np     *= n;
+    /* n^5 coeff corrected from 1262/105 -> -1262/105 */
+    Q->cgb[2] = np*( 56/15.0  + n*(-136/35.0 + n*(-1262/105.0 +
+                n*( 73814/2835.0))));
+    Q->cbg[2] = np*(-26/15.0  + n*(  34/21.0 + n*(    8/5.0   +
+                n*(-12686/2835.0))));
+    np     *= n;
+    /* n^5 coeff corrected from 322/35 -> 332/35 */
+    Q->cgb[3] = np*(4279/630.0 + n*(-332/35.0 + n*(-399572/14175.0)));
+    Q->cbg[3] = np*(1237/630.0 + n*( -12/5.0  + n*( -24832/14175.0)));
+    np     *= n;
+    Q->cgb[4] = np*(4174/315.0 + n*(-144838/6237.0 ));
+    Q->cbg[4] = np*(-734/315.0 + n*( 109598/31185.0));
+    np     *= n;
+    Q->cgb[5] = np*(601676/22275.0 );
+    Q->cbg[5] = np*(444337/155925.0);
+
+    /* Constants of the projections */
+    /* Transverse Mercator (UTM, ITM, etc) */
+    np = n*n;
+    /* Norm. mer. quad, K&W p.50 (96), p.19 (38b), p.5 (2) */
+    Q->Qn = P->k0/(1 + n) * (1 + np*(1/4.0 + np*(1/64.0 + np/256.0)));
+    /* coef of trig series */
+    /* utg := ell. N, E -> sph. N, E,  KW p194 (65) */
+    /* gtu := sph. N, E -> ell. N, E,  KW p196 (69) */
+    Q->utg[0] = n*(-0.5  + n*( 2/3.0 + n*(-37/96.0 + n*( 1/360.0 +
+                n*(  81/512.0 + n*(-96199/604800.0))))));
+    Q->gtu[0] = n*( 0.5  + n*(-2/3.0 + n*(  5/16.0 + n*(41/180.0 +
+                n*(-127/288.0 + n*(  7891/37800.0 ))))));
+    Q->utg[1] = np*(-1/48.0 + n*(-1/15.0 + n*(437/1440.0 + n*(-46/105.0 +
+                n*( 1118711/3870720.0)))));
+    Q->gtu[1] = np*(13/48.0 + n*(-3/5.0  + n*(557/1440.0 + n*(281/630.0 +
+                n*(-1983433/1935360.0)))));
+    np      *= n;
+    Q->utg[2] = np*(-17/480.0 + n*(  37/840.0 + n*(  209/4480.0  +
+                n*( -5569/90720.0 ))));
+    Q->gtu[2] = np*( 61/240.0 + n*(-103/140.0 + n*(15061/26880.0 +
+                n*(167603/181440.0))));
+    np      *= n;
+    Q->utg[3] = np*(-4397/161280.0 + n*(  11/504.0 + n*( 830251/7257600.0)));
+    Q->gtu[3] = np*(49561/161280.0 + n*(-179/168.0 + n*(6601661/7257600.0)));
+    np     *= n;
+    Q->utg[4] = np*(-4583/161280.0 + n*(  108847/3991680.0));
+    Q->gtu[4] = np*(34729/80640.0  + n*(-3418889/1995840.0));
+    np     *= n;
+    Q->utg[5] = np*(-20648693/638668800.0);
+    Q->gtu[5] = np*(212378941/319334400.0);
+
+    /* Gaussian latitude value of the origin latitude */
+    Z = gatg (Q->cbg, PROJ_ETMERC_ORDER, P->phi0);
+
+    /* Origin northing minus true northing at the origin latitude */
+    /* i.e. true northing = N - P->Zb                         */
+    Q->Zb  = - Q->Qn*(Z + clens(Q->gtu, PROJ_ETMERC_ORDER, 2*Z));
+    P->inv = exact_e_inv;
+    P->fwd = exact_e_fwd;
+    return P;
+}
+
+
+
+
+/*****************************************************************************/
+//
+//                                Operation Setups
+//
+/*****************************************************************************/
+
 PJ *PROJECTION(tmerc) {
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
+    /* exact transverse mercator only exists in ellipsoidal form, */
+    /* use approximate version if +a sphere is requested          */
+    if (pj_param (P->ctx, P->params, "bapprox").i || P->es <= 0) {
+        struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(pj_calloc (1, sizeof (struct pj_opaque_approx)));
+        if (nullptr==Q)
+            return pj_default_destructor (P, ENOMEM);
+
+        P->opaque = Q;
+
+        return setup_approx(P);
+    } else {
+        struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(pj_calloc (1, sizeof (struct pj_opaque_exact)));
+        if (nullptr==Q)
+            return pj_default_destructor (P, ENOMEM);
+        P->opaque = Q;
+       return setup_exact (P);
+    }
+}
+
+
+PJ *PROJECTION(etmerc) {
+    struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(pj_calloc (1, sizeof (struct pj_opaque_exact)));
     if (nullptr==Q)
         return pj_default_destructor (P, ENOMEM);
-
     P->opaque = Q;
-    P->destructor = destructor;
+   return setup_exact (P);
+}
+
 
-    return setup(P);
+/* UTM uses the Poder/Engsager implementation for the underlying projection      */
+/* UNLESS +approx is set in which case the Evenden/Snyder implemenation is used. */
+PJ *PROJECTION(utm) {
+    long zone;
+    if (P->es == 0.0) {
+        proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED);
+        return pj_default_destructor(P, ENOMEM);
+    }
+    if (P->lam0 < -1000.0 || P->lam0 > 1000.0) {
+        return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE);
+    }
+
+    P->y0 = pj_param (P->ctx, P->params, "bsouth").i ? 10000000. : 0.;
+    P->x0 = 500000.;
+    if (pj_param (P->ctx, P->params, "tzone").i) /* zone input ? */
+    {
+        zone = pj_param(P->ctx, P->params, "izone").i;
+        if (zone > 0 && zone <= 60)
+            --zone;
+        else {
+            return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE);
+        }
+    }
+    else /* nearest central meridian input */
+    {
+        zone = lround((floor ((adjlon (P->lam0) + M_PI) * 30. / M_PI)));
+        if (zone < 0)
+            zone = 0;
+        else if (zone >= 60)
+            zone = 59;
+    }
+    P->lam0 = (zone + .5) * M_PI / 30. - M_PI;
+    P->k0 = 0.9996;
+    P->phi0 = 0.;
+
+    if (pj_param(P->ctx, P->params, "bapprox").i) {
+        struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(pj_calloc (1, sizeof (struct pj_opaque_approx)));
+        if (nullptr==Q)
+            return pj_default_destructor (P, ENOMEM);
+        P->opaque = Q;
+
+        return setup_approx(P);
+    } else {
+        struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(pj_calloc (1, sizeof (struct pj_opaque_exact)));
+        if (nullptr==Q)
+            return pj_default_destructor (P, ENOMEM);
+        P->opaque = Q;
+
+        return setup_exact(P);
+    }
 }
diff --git a/test/cli/testvarious b/test/cli/testvarious
index 7ec50bb3..c1fa61df 100755
--- a/test/cli/testvarious
+++ b/test/cli/testvarious
@@ -239,7 +239,7 @@ EOF
 echo "##############################################################" >> ${OUT}
 echo "Test transverse mercator (#97)" >> ${OUT}
 #
-$EXE +proj=tmerc +k=0.998 +lon_0=-20 +datum=WGS84 +x_0=10000 +y_0=20000 \
+$EXE +proj=tmerc +approx +k=0.998 +lon_0=-20 +datum=WGS84 +x_0=10000 +y_0=20000 \
  +to +proj=latlong +datum=WGS84 \
  -E >>${OUT} <<EOF
 10000 20000
@@ -253,7 +253,7 @@ echo "##############################################################" >> ${OUT}
 echo "Test transverse mercator inverse (#97)" >> ${OUT}
 #
 $EXE +proj=latlong +datum=WGS84 \
- +to +proj=tmerc +k=0.998 +lon_0=-20 +datum=WGS84 +x_0=10000 +y_0=20000 \
+ +to +proj=tmerc +approx +k=0.998 +lon_0=-20 +datum=WGS84 +x_0=10000 +y_0=20000 \
  -E >>${OUT} <<EOF
 0dN 0.000
 15d22'16.108"W	17d52'53.478"N 0.000
diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist
index 148d413d..72e95634 100644
--- a/test/cli/tv_out.dist
+++ b/test/cli/tv_out.dist
@@ -352,7 +352,7 @@ Test inverse handling
 10 20	-1384841.19	7581707.88 0.00
 ##############################################################
 Test MGI datum gives expected results (#207)
-16.33 48.20	595710.3732102	5357598.4645755 -44.4951085
+16.33 48.20	595710.3731286	5357598.4645652 -44.4951085
 ##############################################################
 Test omerc sensitivity with locations 90d from origin(#114)
 56.958381652832 72.8798	-9985.16336453	-227.67701050 0.00000000
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index 82222210..930ca71f 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -5218,6 +5218,16 @@ accept  -200 -100
 expect  -7.490535682 -0.000901935
 
 
+operation +proj=utm +zone=32
+tolerance 0.001 mm
+accept	12 56 0 2000
+expect  687071.43910944  6210141.32674801    0.00000000     2000.0000
+
+operation +proj=utm +zone=32 +approx
+tolerance 0.001 mm
+accept	12 56 0 2000
+expect  687071.43911000  6210141.32675053    0.00000000     2000.0000
+
 ===============================================================================
 van der Grinten (I)
 	Misc Sph
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 8ef426ac..d535e412 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -507,36 +507,21 @@ TEST_F(CApi, proj_as_proj_string_incompatible_WKT1) {
 
 // ---------------------------------------------------------------------------
 
-TEST_F(CApi, proj_as_proj_string_etmerc_option_yes) {
+TEST_F(CApi, proj_as_proj_string_approx_tmerc_option_yes) {
     auto obj = proj_create(m_ctxt, "+proj=tmerc +type=crs");
     ObjectKeeper keeper(obj);
     ASSERT_NE(obj, nullptr);
 
-    const char *options[] = {"USE_ETMERC=YES", nullptr};
+    const char *options[] = {"USE_APPROX_TMERC=YES", nullptr};
     auto str = proj_as_proj_string(m_ctxt, obj, PJ_PROJ_4, options);
     ASSERT_NE(str, nullptr);
     EXPECT_EQ(str,
-              std::string("+proj=etmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 "
+              std::string("+proj=tmerc +approx +lat_0=0 +lon_0=0 +k=1 +x_0=0 "
                           "+y_0=0 +datum=WGS84 +units=m +no_defs +type=crs"));
 }
 
 // ---------------------------------------------------------------------------
 
-TEST_F(CApi, proj_as_proj_string_etmerc_option_no) {
-    auto obj = proj_create(m_ctxt, "+proj=utm +zone=31 +type=crs");
-    ObjectKeeper keeper(obj);
-    ASSERT_NE(obj, nullptr);
-
-    const char *options[] = {"USE_ETMERC=NO", nullptr};
-    auto str = proj_as_proj_string(m_ctxt, obj, PJ_PROJ_4, options);
-    ASSERT_NE(str, nullptr);
-    EXPECT_EQ(str, std::string("+proj=tmerc +lat_0=0 +lon_0=3 +k=0.9996 "
-                               "+x_0=500000 +y_0=0 +datum=WGS84 +units=m "
-                               "+no_defs +type=crs"));
-}
-
-// ---------------------------------------------------------------------------
-
 TEST_F(CApi, proj_crs_create_bound_crs_to_WGS84) {
     auto crs = proj_create_from_database(m_ctxt, "EPSG", "3844",
                                          PJ_CATEGORY_CRS, false, nullptr);
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 69ef9073..3cb50352 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -7654,12 +7654,39 @@ TEST(io, projparse_etmerc) {
         crs->exportToPROJString(
             PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
                 .get()),
-        "+proj=etmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 "
+        "+proj=tmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 "
         "+datum=WGS84 +units=m +no_defs +type=crs");
 
     auto wkt1 = crs->exportToWKT(
         WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
-    EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\",\"+proj=etmerc +lat_0=0 "
+    EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\"") == std::string::npos)
+        << wkt1;
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_tmerc_approx) {
+    auto obj =
+        PROJStringParser().createFromPROJString("+proj=tmerc +approx +type=crs");
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    auto wkt2 = crs->exportToWKT(
+        &WKTFormatter::create()->simulCurNodeHasId().setMultiLine(false));
+    EXPECT_TRUE(
+        wkt2.find("METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]]") !=
+        std::string::npos)
+        << wkt2;
+
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        "+proj=tmerc +approx +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 "
+        "+datum=WGS84 +units=m +no_defs +type=crs");
+
+    auto wkt1 = crs->exportToWKT(
+        WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
+    EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\",\"+proj=tmerc +approx +lat_0=0 "
                           "+lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m "
                           "+no_defs\"]") != std::string::npos)
         << wkt1;
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 6c1ecf1c..586226e2 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -1237,9 +1237,9 @@ TEST(operation, tmerc_export) {
 
     {
         auto formatter = PROJStringFormatter::create();
-        formatter->setUseETMercForTMerc(true);
+        formatter->setUseApproxTMerc(true);
         EXPECT_EQ(conv->exportToPROJString(formatter.get()),
-                  "+proj=etmerc +lat_0=1 +lon_0=2 +k=3 +x_0=4 +y_0=5");
+                  "+proj=tmerc +approx +lat_0=1 +lon_0=2 +k=3 +x_0=4 +y_0=5");
     }
 
     EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()),
@@ -6061,21 +6061,21 @@ TEST(operation, compoundCRS_to_compoundCRS_with_vertical_transform) {
               "+ellps=WGS84");
     {
         auto formatter = PROJStringFormatter::create();
-        formatter->setUseETMercForTMerc(true);
+        formatter->setUseApproxTMerc(true);
         EXPECT_EQ(op->exportToPROJString(formatter.get()),
-                  "+proj=pipeline +step +inv +proj=etmerc +lat_0=1 +lon_0=2 "
+                  "+proj=pipeline +step +inv +proj=tmerc +approx +lat_0=1 +lon_0=2 "
                   "+k=3 +x_0=4 +y_0=5 +ellps=WGS84 +step "
                   "+proj=vgridshift +grids=bla.gtx +multiplier=0.001 +step "
-                  "+proj=utm +zone=32 "
+                  "+proj=utm +approx +zone=32 "
                   "+ellps=WGS84");
     }
     {
         auto formatter = PROJStringFormatter::create();
-        formatter->setUseETMercForTMerc(true);
+        formatter->setUseApproxTMerc(true);
         EXPECT_EQ(op->inverse()->exportToPROJString(formatter.get()),
-                  "+proj=pipeline +step +inv +proj=utm +zone=32 +ellps=WGS84 "
+                  "+proj=pipeline +step +inv +proj=utm +approx +zone=32 +ellps=WGS84 "
                   "+step +inv +proj=vgridshift +grids=bla.gtx "
-                  "+multiplier=0.001 +step +proj=etmerc +lat_0=1 +lon_0=2 "
+                  "+multiplier=0.001 +step +proj=tmerc +approx +lat_0=1 +lon_0=2 "
                   "+k=3 +x_0=4 +y_0=5 +ellps=WGS84");
     }