GeographicLib  1.21
TransverseMercatorExact.cpp
Go to the documentation of this file.
00001 /**
00002  * \file TransverseMercatorExact.cpp
00003  * \brief Implementation for GeographicLib::TransverseMercatorExact class
00004  *
00005  * Copyright (c) Charles Karney (2008-2011) <charles@karney.com> and licensed
00006  * under the MIT/X11 License.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  *
00009  * The relevant section of Lee's paper is part V, pp 67&ndash;101,
00010  * <a href="http://dx.doi.org/10.3138/X687-1574-4325-WM62">Conformal
00011  * Projections Based On Jacobian Elliptic Functions</a>.
00012  *
00013  * The method entails using the Thompson Transverse Mercator as an
00014  * intermediate projection.  The projections from the intermediate
00015  * coordinates to [\e phi, \e lam] and [\e x, \e y] are given by elliptic
00016  * functions.  The inverse of these projections are found by Newton's method
00017  * with a suitable starting guess.
00018  *
00019  * This implementation and notation closely follows Lee, with the following
00020  * exceptions:
00021  * <center><table>
00022  * <tr><th>Lee    <th>here    <th>Description
00023  * <tr><td>x/a    <td>xi      <td>Northing (unit Earth)
00024  * <tr><td>y/a    <td>eta     <td>Easting (unit Earth)
00025  * <tr><td>s/a    <td>sigma   <td>xi + i * eta
00026  * <tr><td>y      <td>x       <td>Easting
00027  * <tr><td>x      <td>y       <td>Northing
00028  * <tr><td>k      <td>e       <td>eccentricity
00029  * <tr><td>k^2    <td>mu      <td>elliptic function parameter
00030  * <tr><td>k'^2   <td>mv      <td>elliptic function complementary parameter
00031  * <tr><td>m      <td>k       <td>scale
00032  * <tr><td>zeta   <td>zeta    <td>complex longitude = Mercator = chi in paper
00033  * <tr><td>s      <td>sigma   <td>complex GK = zeta in paper
00034  * </table></center>
00035  *
00036  * Minor alterations have been made in some of Lee's expressions in an
00037  * attempt to control round-off.  For example atanh(sin(phi)) is replaced by
00038  * asinh(tan(phi)) which maintains accuracy near phi = pi/2.  Such changes
00039  * are noted in the code.
00040  **********************************************************************/
00041 
00042 #include <GeographicLib/TransverseMercatorExact.hpp>
00043 
00044 #define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_CPP \
00045   "$Id: 125a2d11919018a84fb0c09ea2e77011a35a4a2d $"
00046 
00047 RCSID_DECL(GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_CPP)
00048 RCSID_DECL(GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP)
00049 
00050 namespace GeographicLib {
00051 
00052   using namespace std;
00053 
00054   const Math::real TransverseMercatorExact::tol_ =
00055     numeric_limits<real>::epsilon();
00056   const Math::real TransverseMercatorExact::tol1_ = real(0.1) * sqrt(tol_);
00057   const Math::real TransverseMercatorExact::tol2_ = real(0.1) * tol_;
00058   const Math::real TransverseMercatorExact::taytol_ = pow(tol_, real(0.6));
00059   // Overflow value s.t. atan(overflow_) = pi/2
00060   const Math::real TransverseMercatorExact::overflow_ = 1 / Math::sq(tol_);
00061 
00062   TransverseMercatorExact::TransverseMercatorExact(real a, real f, real k0,
00063                                                    bool extendp)
00064     : _a(a)
00065     , _f(f <= 1 ? f : 1/f)
00066     , _k0(k0)
00067     , _mu(_f * (2 - _f))        // e^2
00068     , _mv(1 - _mu)              // 1 - e^2
00069     , _e(sqrt(_mu))
00070     , _ep2(_mu / _mv)           // e^2 / (1 - e^2)
00071     , _extendp(extendp)
00072     , _Eu(_mu)
00073     , _Ev(_mv)
00074   {
00075     if (!(Math::isfinite(_a) && _a > 0))
00076       throw GeographicErr("Major radius is not positive");
00077     if (!(_f > 0))
00078       throw GeographicErr("Flattening is not positive");
00079     if (!(_f < 1))
00080       throw GeographicErr("Minor radius is not positive");
00081     if (!(Math::isfinite(_k0) && _k0 > 0))
00082       throw GeographicErr("Scale is not positive");
00083   }
00084 
00085   const TransverseMercatorExact
00086   TransverseMercatorExact::UTM(Constants::WGS84_a<real>(),
00087                                Constants::WGS84_f<real>(),
00088                                Constants::UTM_k0<real>());
00089 
00090   // tau = tan(phi), taup = sinh(psi)
00091   Math::real TransverseMercatorExact::taup(real tau) const throw() {
00092     real
00093       tau1 = Math::hypot(real(1), tau),
00094       sig = sinh( _e * Math::atanh(_e * tau / tau1) );
00095     return Math::hypot(real(1), sig) * tau - sig * tau1;
00096   }
00097 
00098   Math::real TransverseMercatorExact::taupinv(real taup) const throw() {
00099     real
00100       // See comment in TransverseMercator.cpp about the initial guess
00101       tau = taup/_mv,
00102       stol = tol_ * max(real(1), abs(taup));
00103     // min iterations = 1, max iterations = 2; mean = 1.94
00104     for (int i = 0; i < numit_; ++i) {
00105       real
00106         tau1 = Math::hypot(real(1), tau),
00107         sig = sinh( _e * Math::atanh(_e * tau / tau1 ) ),
00108         taupa = Math::hypot(real(1), sig) * tau - sig * tau1,
00109         dtau = (taup - taupa) * (1 + _mv * Math::sq(tau)) /
00110         ( _mv * tau1 * Math::hypot(real(1), taupa) );
00111       tau += dtau;
00112       if (!(abs(dtau) >= stol))
00113         break;
00114     }
00115     return tau;
00116   }
00117 
00118   void TransverseMercatorExact::zeta(real /*u*/, real snu, real cnu, real dnu,
00119                                      real /*v*/, real snv, real cnv, real dnv,
00120                                      real& taup, real& lam) const throw() {
00121     // Lee 54.17 but write
00122     // atanh(snu * dnv) = asinh(snu * dnv / sqrt(cnu^2 + _mv * snu^2 * snv^2))
00123     // atanh(_e * snu / dnv) =
00124     //         asinh(_e * snu / sqrt(_mu * cnu^2 + _mv * cnv^2))
00125     real
00126       d1 = sqrt(Math::sq(cnu) + _mv * Math::sq(snu * snv)),
00127       d2 = sqrt(_mu * Math::sq(cnu) + _mv * Math::sq(cnv)),
00128       t1 = (d1 ? snu * dnv / d1 : (snu < 0 ? -overflow_ : overflow_)),
00129       t2 = (d2 ? sinh( _e * Math::asinh(_e * snu / d2) ) :
00130             (snu < 0 ? -overflow_ : overflow_));
00131     // psi = asinh(t1) - asinh(t2)
00132     // taup = sinh(psi)
00133     taup = t1 * Math::hypot(real(1), t2) - t2 * Math::hypot(real(1), t1);
00134     lam = (d1 != 0 && d2 != 0) ?
00135       atan2(dnu * snv, cnu * cnv) - _e * atan2(_e * cnu * snv, dnu * cnv) :
00136       0;
00137   }
00138 
00139   void TransverseMercatorExact::dwdzeta(real /*u*/,
00140                                         real snu, real cnu, real dnu,
00141                                         real /*v*/,
00142                                         real snv, real cnv, real dnv,
00143                                         real& du, real& dv) const throw() {
00144     // Lee 54.21 but write (1 - dnu^2 * snv^2) = (cnv^2 + _mu * snu^2 * snv^2)
00145     // (see A+S 16.21.4)
00146     real d = _mv * Math::sq(Math::sq(cnv) + _mu * Math::sq(snu * snv));
00147     du =  cnu * dnu * dnv * (Math::sq(cnv) - _mu * Math::sq(snu * snv)) / d;
00148     dv = -snu * snv * cnv * (Math::sq(dnu * dnv) + _mu * Math::sq(cnu)) / d;
00149   }
00150 
00151   // Starting point for zetainv
00152   bool TransverseMercatorExact::zetainv0(real psi, real lam, real& u, real& v)
00153     const throw() {
00154     bool retval = false;
00155     if (psi < -_e * Math::pi<real>()/4 &&
00156         lam > (1 - 2 * _e) * Math::pi<real>()/2 &&
00157         psi < lam - (1 - _e) * Math::pi<real>()/2) {
00158       // N.B. this branch is normally not taken because psi < 0 is converted
00159       // psi > 0 by Forward.
00160       //
00161       // There's a log singularity at w = w0 = Eu.K() + i * Ev.K(),
00162       // corresponding to the south pole, where we have, approximately
00163       //
00164       //   psi = _e + i * pi/2 - _e * atanh(cos(i * (w - w0)/(1 + _mu/2)))
00165       //
00166       // Inverting this gives:
00167       real
00168         psix = 1 - psi / _e,
00169         lamx = (Math::pi<real>()/2 - lam) / _e;
00170       u = Math::asinh(sin(lamx) / Math::hypot(cos(lamx), sinh(psix))) *
00171         (1 + _mu/2);
00172       v = atan2(cos(lamx), sinh(psix)) * (1 + _mu/2);
00173       u = _Eu.K() - u;
00174       v = _Ev.K() - v;
00175     } else if (psi < _e * Math::pi<real>()/2 &&
00176                lam > (1 - 2 * _e) * Math::pi<real>()/2) {
00177       // At w = w0 = i * Ev.K(), we have
00178       //
00179       //     zeta = zeta0 = i * (1 - _e) * pi/2
00180       //     zeta' = zeta'' = 0
00181       //
00182       // including the next term in the Taylor series gives:
00183       //
00184       // zeta = zeta0 - (_mv * _e) / 3 * (w - w0)^3
00185       //
00186       // When inverting this, we map arg(w - w0) = [-90, 0] to
00187       // arg(zeta - zeta0) = [-90, 180]
00188       real
00189         dlam = lam - (1 - _e) * Math::pi<real>()/2,
00190         rad = Math::hypot(psi, dlam),
00191         // atan2(dlam-psi, psi+dlam) + 45d gives arg(zeta - zeta0) in range
00192         // [-135, 225).  Subtracting 180 (since multiplier is negative) makes
00193         // range [-315, 45).  Multiplying by 1/3 (for cube root) gives range
00194         // [-105, 15).  In particular the range [-90, 180] in zeta space maps
00195         // to [-90, 0] in w space as required.
00196         ang = atan2(dlam-psi, psi+dlam) - real(0.75) * Math::pi<real>();
00197       // Error using this guess is about 0.21 * (rad/e)^(5/3)
00198       retval = rad < _e * taytol_;
00199       rad = Math::cbrt(3 / (_mv * _e) * rad);
00200       ang /= 3;
00201       u = rad * cos(ang);
00202       v = rad * sin(ang) + _Ev.K();
00203     } else {
00204       // Use spherical TM, Lee 12.6 -- writing atanh(sin(lam) / cosh(psi)) =
00205       // asinh(sin(lam) / hypot(cos(lam), sinh(psi))).  This takes care of the
00206       // log singularity at zeta = Eu.K() (corresponding to the north pole)
00207       v = Math::asinh(sin(lam) / Math::hypot(cos(lam), sinh(psi)));
00208       u = atan2(sinh(psi), cos(lam));
00209       // But scale to put 90,0 on the right place
00210       u *= _Eu.K() / (Math::pi<real>()/2);
00211       v *= _Eu.K() / (Math::pi<real>()/2);
00212     }
00213     return retval;
00214   }
00215 
00216   // Invert zeta using Newton's method
00217   void TransverseMercatorExact::zetainv(real taup, real lam, real& u, real& v)
00218     const throw()  {
00219     real
00220       psi = Math::asinh(taup),
00221       scal = 1/Math::hypot(real(1), taup);
00222     if (zetainv0(psi, lam, u, v))
00223       return;
00224     real stol2 = tol2_ / Math::sq(max(psi, real(1)));
00225     // min iterations = 2, max iterations = 6; mean = 4.0
00226     for (int i = 0, trip = 0; i < numit_; ++i) {
00227       real snu, cnu, dnu, snv, cnv, dnv;
00228       _Eu.sncndn(u, snu, cnu, dnu);
00229       _Ev.sncndn(v, snv, cnv, dnv);
00230       real tau1, lam1, du1, dv1;
00231       zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau1, lam1);
00232       dwdzeta(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
00233       tau1 -= taup;
00234       lam1 -= lam;
00235       tau1 *= scal;
00236       real
00237         delu = tau1 * du1 - lam1 * dv1,
00238         delv = tau1 * dv1 + lam1 * du1;
00239       u -= delu;
00240       v -= delv;
00241       if (trip)
00242         break;
00243       real delw2 = Math::sq(delu) + Math::sq(delv);
00244       if (!(delw2 >= stol2))
00245         ++trip;
00246     }
00247   }
00248 
00249   void TransverseMercatorExact::sigma(real /*u*/, real snu, real cnu, real dnu,
00250                                       real v, real snv, real cnv, real dnv,
00251                                       real& xi, real& eta) const throw() {
00252     // Lee 55.4 writing
00253     // dnu^2 + dnv^2 - 1 = _mu * cnu^2 + _mv * cnv^2
00254     real d = _mu * Math::sq(cnu) + _mv * Math::sq(cnv);
00255     xi = _Eu.E(snu, cnu, dnu) - _mu * snu * cnu * dnu / d;
00256     eta = v - _Ev.E(snv, cnv, dnv) + _mv * snv * cnv * dnv / d;
00257   }
00258 
00259   void TransverseMercatorExact::dwdsigma(real /*u*/,
00260                                          real snu, real cnu, real dnu,
00261                                          real /*v*/,
00262                                          real snv, real cnv, real dnv,
00263                                          real& du, real& dv) const throw() {
00264     // Reciprocal of 55.9: dw/ds = dn(w)^2/_mv, expanding complex dn(w) using
00265     // A+S 16.21.4
00266     real d = _mv * Math::sq(Math::sq(cnv) + _mu * Math::sq(snu * snv));
00267     real
00268       dnr = dnu * cnv * dnv,
00269       dni = - _mu * snu * cnu * snv;
00270     du = (Math::sq(dnr) - Math::sq(dni)) / d;
00271     dv = 2 * dnr * dni / d;
00272   }
00273 
00274   // Starting point for sigmainv
00275   bool TransverseMercatorExact::sigmainv0(real xi, real eta, real& u, real& v)
00276     const throw() {
00277     bool retval = false;
00278     if (eta > real(1.25) * _Ev.KE() ||
00279         (xi < -real(0.25) * _Eu.E() && xi < eta - _Ev.KE())) {
00280       // sigma as a simple pole at w = w0 = Eu.K() + i * Ev.K() and sigma is
00281       // approximated by
00282       //
00283       // sigma = (Eu.E() + i * Ev.KE()) + 1/(w - w0)
00284       real
00285         x = xi - _Eu.E(),
00286         y = eta - _Ev.KE(),
00287         r2 = Math::sq(x) + Math::sq(y);
00288       u = _Eu.K() + x/r2;
00289       v = _Ev.K() - y/r2;
00290     } else if ((eta > real(0.75) * _Ev.KE() && xi < real(0.25) * _Eu.E())
00291                || eta > _Ev.KE()) {
00292       // At w = w0 = i * Ev.K(), we have
00293       //
00294       //     sigma = sigma0 = i * Ev.KE()
00295       //     sigma' = sigma'' = 0
00296       //
00297       // including the next term in the Taylor series gives:
00298       //
00299       // sigma = sigma0 - _mv / 3 * (w - w0)^3
00300       //
00301       // When inverting this, we map arg(w - w0) = [-pi/2, -pi/6] to
00302       // arg(sigma - sigma0) = [-pi/2, pi/2]
00303       // mapping arg = [-pi/2, -pi/6] to [-pi/2, pi/2]
00304       real
00305         deta = eta - _Ev.KE(),
00306         rad = Math::hypot(xi, deta),
00307         // Map the range [-90, 180] in sigma space to [-90, 0] in w space.  See
00308         // discussion in zetainv0 on the cut for ang.
00309         ang = atan2(deta-xi, xi+deta) - real(0.75) * Math::pi<real>();
00310       // Error using this guess is about 0.068 * rad^(5/3)
00311       retval = rad < 2 * taytol_;
00312       rad = Math::cbrt(3 / _mv * rad);
00313       ang /= 3;
00314       u = rad * cos(ang);
00315       v = rad * sin(ang) + _Ev.K();
00316     } else {
00317       // Else use w = sigma * Eu.K/Eu.E (which is correct in the limit _e -> 0)
00318       u = xi * _Eu.K()/_Eu.E();
00319       v = eta * _Eu.K()/_Eu.E();
00320     }
00321     return retval;
00322   }
00323 
00324   // Invert sigma using Newton's method
00325   void TransverseMercatorExact::sigmainv(real xi, real eta, real& u, real& v)
00326     const throw() {
00327     if (sigmainv0(xi, eta, u, v))
00328       return;
00329     // min iterations = 2, max iterations = 7; mean = 3.9
00330     for (int i = 0, trip = 0; i < numit_; ++i) {
00331       real snu, cnu, dnu, snv, cnv, dnv;
00332       _Eu.sncndn(u, snu, cnu, dnu);
00333       _Ev.sncndn(v, snv, cnv, dnv);
00334       real xi1, eta1, du1, dv1;
00335       sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi1, eta1);
00336       dwdsigma(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
00337       xi1 -= xi;
00338       eta1 -= eta;
00339       real
00340         delu = xi1 * du1 - eta1 * dv1,
00341         delv = xi1 * dv1 + eta1 * du1;
00342       u -= delu;
00343       v -= delv;
00344       if (trip)
00345         break;
00346       real delw2 = Math::sq(delu) + Math::sq(delv);
00347       if (!(delw2 >= tol2_))
00348         ++trip;
00349     }
00350   }
00351 
00352   void TransverseMercatorExact::Scale(real tau, real /*lam*/,
00353                                        real snu, real cnu, real dnu,
00354                                        real snv, real cnv, real dnv,
00355                                        real& gamma, real& k) const throw() {
00356     real sec2 = 1 + Math::sq(tau);    // sec(phi)^2
00357     // Lee 55.12 -- negated for our sign convention.  gamma gives the bearing
00358     // (clockwise from true north) of grid north
00359     gamma = atan2(_mv * snu * snv * cnv, cnu * dnu * dnv);
00360     // Lee 55.13 with nu given by Lee 9.1 -- in sqrt change the numerator
00361     // from
00362     //
00363     //    (1 - snu^2 * dnv^2) to (_mv * snv^2 + cnu^2 * dnv^2)
00364     //
00365     // to maintain accuracy near phi = 90 and change the denomintor from
00366     //
00367     //    (dnu^2 + dnv^2 - 1) to (_mu * cnu^2 + _mv * cnv^2)
00368     //
00369     // to maintain accuracy near phi = 0, lam = 90 * (1 - e).  Similarly
00370     // rewrite sqrt term in 9.1 as
00371     //
00372     //    _mv + _mu * c^2 instead of 1 - _mu * sin(phi)^2
00373     k = sqrt(_mv + _mu / sec2) * sqrt(sec2) *
00374       sqrt( (_mv * Math::sq(snv) + Math::sq(cnu * dnv)) /
00375             (_mu * Math::sq(cnu) + _mv * Math::sq(cnv)) );
00376   }
00377 
00378   void TransverseMercatorExact::Forward(real lon0, real lat, real lon,
00379                                         real& x, real& y, real& gamma, real& k)
00380     const throw() {
00381     // Avoid losing a bit of accuracy in lon (assuming lon0 is an integer)
00382     if (lon - lon0 > 180)
00383       lon -= lon0 + 360;
00384     else if (lon - lon0 <= -180)
00385       lon -= lon0 - 360;
00386     else
00387       lon -= lon0;
00388     // Now lon in (-180, 180]
00389     // Explicitly enforce the parity
00390     int
00391       latsign = !_extendp && lat < 0 ? -1 : 1,
00392       lonsign = !_extendp && lon < 0 ? -1 : 1;
00393     lon *= lonsign;
00394     lat *= latsign;
00395     bool backside = !_extendp && lon > 90;
00396     if (backside) {
00397       if (lat == 0)
00398         latsign = -1;
00399       lon = 180 - lon;
00400     }
00401     real
00402       phi = lat * Math::degree<real>(),
00403       lam = lon * Math::degree<real>(),
00404       tau = tanx(phi);
00405 
00406     // u,v = coordinates for the Thompson TM, Lee 54
00407     real u, v;
00408     if (lat == 90) {
00409       u = _Eu.K();
00410       v = 0;
00411     } else if (lat == 0 && lon == 90 * (1 - _e)) {
00412       u = 0;
00413       v = _Ev.K();
00414     } else
00415       zetainv(taup(tau), lam, u, v);
00416 
00417     real snu, cnu, dnu, snv, cnv, dnv;
00418     _Eu.sncndn(u, snu, cnu, dnu);
00419     _Ev.sncndn(v, snv, cnv, dnv);
00420 
00421     real xi, eta;
00422     sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi, eta);
00423     if (backside)
00424       xi = 2 * _Eu.E() - xi;
00425     y = xi * _a * _k0 * latsign;
00426     x = eta * _a * _k0 * lonsign;
00427 
00428     // Recompute (tau, lam) from (u, v) to improve accuracy of Scale
00429     zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
00430     tau=taupinv(tau);
00431     Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
00432     gamma /= Math::degree<real>();
00433     if (backside)
00434       gamma = 180 - gamma;
00435     gamma *= latsign * lonsign;
00436     k *= _k0;
00437   }
00438 
00439   void TransverseMercatorExact::Reverse(real lon0, real x, real y,
00440                                         real& lat, real& lon,
00441                                         real& gamma, real& k)
00442     const throw() {
00443     // This undoes the steps in Forward.
00444     real
00445       xi = y / (_a * _k0),
00446       eta = x / (_a * _k0);
00447     // Explicitly enforce the parity
00448     int
00449       latsign = !_extendp && y < 0 ? -1 : 1,
00450       lonsign = !_extendp && x < 0 ? -1 : 1;
00451     xi *= latsign;
00452     eta *= lonsign;
00453     bool backside = !_extendp && xi > _Eu.E();
00454     if (backside)
00455       xi = 2 * _Eu.E()- xi;
00456 
00457     // u,v = coordinates for the Thompson TM, Lee 54
00458     real u, v;
00459     if (xi == 0 && eta == _Ev.KE()) {
00460       u = 0;
00461       v = _Ev.K();
00462     } else
00463       sigmainv(xi, eta, u, v);
00464 
00465     real snu, cnu, dnu, snv, cnv, dnv;
00466     _Eu.sncndn(u, snu, cnu, dnu);
00467     _Ev.sncndn(v, snv, cnv, dnv);
00468     real phi, lam, tau;
00469     if (v != 0 || u != _Eu.K()) {
00470       zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
00471       tau = taupinv(tau);
00472       phi = atan(tau);
00473       lat = phi / Math::degree<real>();
00474       lon = lam / Math::degree<real>();
00475     } else {
00476       tau = overflow_;
00477       phi = Math::pi<real>()/2;
00478       lat = 90;
00479       lon = lam = 0;
00480     }
00481     Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
00482     gamma /= Math::degree<real>();
00483     if (backside)
00484       lon = 180 - lon;
00485     lon *= lonsign;
00486     // Avoid losing a bit of accuracy in lon (assuming lon0 is an integer)
00487     if (lon + lon0 >= 180)
00488       lon += lon0 - 360;
00489     else if (lon + lon0 < -180)
00490       lon += lon0 + 360;
00491     else
00492       lon += lon0;
00493     lat *= latsign;
00494     if (backside)
00495       y = 2 * _Eu.E() - y;
00496     y *= _a * _k0 * latsign;
00497     x *= _a * _k0 * lonsign;
00498     if (backside)
00499       gamma = 180 - gamma;
00500     gamma *= latsign * lonsign;
00501     k *= _k0;
00502   }
00503 
00504 } // namespace GeographicLib