/** * \file AuxLatitude.hpp * \brief Header for the GeographicLib::AuxLatitude and GeographicLib::AuxAngle * classes. * * \note This is just sample code. It is not part of GeographicLib itself. * * This file is an implementation of the methods described in * - C. F. F. Karney, * On auxiliary latitudes, * Technical Report, SRI International, December 2022. * https://arxiv.org/abs/2212.05818 * . * Copyright (c) Charles Karney (2022) and licensed under * the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ **********************************************************************/ #if !defined(AUXLATITUDE_HPP) #define AUXLATITUDE_HPP 1 #include #if !defined(GEOGRAPHICLIB_AUXLATITUDE_ORDER) /** * The order of the series approximation used in AuxLatitude. * GEOGRAPHICLIB_AUXLATITUDE_ORDER can be set to one of [4, 6, 8]. **********************************************************************/ # define GEOGRAPHICLIB_AUXLATITUDE_ORDER 6 #endif namespace GeographicLib { /** * \brief An accurate representation of angles. * * \note This is just sample code. It is not part of GeographicLib itself. * * This class is an implementation of the methods described in * - C. F. F. Karney, * On auxiliary latitudes, * Technical Report, SRI International, December 2022. * https://arxiv.org/abs/2212.05818 * * An angle is represented be the \e y and \e x coordinates of a point in the * 2d plane. The two coordinates are proportional to the sine and cosine of * the angle. This allows angles close to the cardinal points to be * represented accurately. Only angles in [−180°, 180°] can be * represented. (A possible extension would be to keep count of the number * of turns.) * * @tparam T the floating-point type to use for real numbers. **********************************************************************/ template class AuxAngle { public: /** * The floating-point type for real numbers. This just connects to the * template parameters for the class. **********************************************************************/ typedef T real; /** * The constructor. * * @param[in] y the \e y coordinate. * @param[in] x the \e x coordinate. * * \note the \e y coordinate is specified \e first. * \warning either \e x or \e y can be infinite, but not both. * * The defaults (\e x = 1 and \e y = 0) are such that * + no arguments gives an angle of 0; * + 1 argument specifies the tangent of the angle. **********************************************************************/ AuxAngle(real y = 0, real x = 1) : _y(y), _x(x) {} /** * @return the \e y component. This is the sine of the angle if the * AuxAngle has been normalized. **********************************************************************/ real y() const { return _y; } /** * @return the \e x component. This is the cosine of the angle if the * AuxAngle has been normalized. **********************************************************************/ real x() const { return _x; } /** * @return a reference to the \e y component. This allows this component * to be altered. **********************************************************************/ real& y() { return _y; } /** * @return a reference to the \e x component. This allows this component * to be altered. **********************************************************************/ real& x() { return _x; } /** * @return the AuxAngle converted to the conventional angle measured in * degrees. **********************************************************************/ real degrees() const; /** * @return the AuxAngle converted to the conventional angle measured in * radians. **********************************************************************/ real radians() const; /** * @return the tangent of the angle. **********************************************************************/ real tan() const { return _y / _x; } /** * @return a new normalized AuxAngle with the point lying on the unit * circle and the \e y and \e x components are equal to the sine and * cosine of the angle. **********************************************************************/ AuxAngle normalized() const; /** * Normalize the AuxAngle in place so that the \e y and \e x components are * equal to the sine and cosine of the angle. **********************************************************************/ void normalize() { *this = normalized(); } /** * Set the quadrant for the AuxAngle. * * @param[in] p the AuxAngle from which the quadrant information is taken. * @return the new AuxAngle in the same quadrant as \e p. **********************************************************************/ AuxAngle copyquadrant(const AuxAngle& p) const; /** * Add an AuxAngle. * * @param[in] p the AuxAngle to be added. * @return a reference to the new AuxAngle. * * The addition is done in place, altering the current AuxAngle. * * \warning Neither *this nor \e p should have an infinite component. If * necessary, invoke AuxAngle::normalize on these angles first. **********************************************************************/ AuxAngle& operator+=(const AuxAngle& p); /** * Convert degrees to an AuxAngle. * * @param[in] d the angle measured in degrees. * @return the corresponding AuxAngle. * * This allows a new AuxAngle to be initialized as an angle in degrees with * @code * AuxAngle phi = AuxAngle::degrees(d); * @endcode * This is the so-called "named constructor" idiom. **********************************************************************/ static AuxAngle degrees(real d); /** * Convert radians to an AuxAngle. * * @param[in] r the angle measured in radians. * @return the corresponding AuxAngle. * * This allows a new AuxAngle to be initialized as an angle in radians with * @code * AuxAngle phi = AuxAngle::radians(r); * @endcode * This is the so-called "named constructor" idiom. **********************************************************************/ static AuxAngle radians(real r); /** * @return a "NaN" AuxAngle. **********************************************************************/ static AuxAngle NaN(); /** * Compute the absolute error in another angle. * * @tparam T1 the floating-point type of the other angle. * @param[in] p the other angle * @return the absolute error between p and *this considered as angles in * radians. **********************************************************************/ template real AbsError(const AuxAngle& p) const; /** * Compute the relative error in another angle. * * @tparam T1 the floating-point type of the other angle. * @param[in] p the other angle * @return the relative error between p.tan() and this->tan(). **********************************************************************/ template real RelError(const AuxAngle& p) const; private: real _y, _x; }; /// \cond SKIP template inline AuxAngle AuxAngle::degrees(real d) { real y, x; Math::sincosd(d, y, x); return AuxAngle(y, x); } template inline AuxAngle AuxAngle::radians(real r) { using std::sin; using std::cos; return AuxAngle(sin(r), cos(r)); } template inline T AuxAngle::degrees() const { return Math::atan2d(_y, _x); } template inline T AuxAngle::radians() const { using std::atan2; return atan2(_y, _x); } template template inline T AuxAngle::AbsError(const AuxAngle& p) const { using std::fabs; return fabs((AuxAngle(-T(p.y()), T(p.x())) += *this).radians()); } template template inline T AuxAngle::RelError(const AuxAngle& p) const { using std::fabs; return fabs((T(p.y()) / T(p.x()) - tan()) / tan()); } /// \endcond /** * \brief Conversions between auxiliary latitudes. * * \note This is just sample code. It is not part of GeographicLib itself. * * This class is an implementation of the methods described in * - C. F. F. Karney, * On auxiliary latitudes, * Technical Report, SRI International, December 2022. * https://arxiv.org/abs/2212.05818 * * The provides accurate conversions between geographic (\e phi, φ), * parametric (\e beta, β), geocentric (\e theta, θ), rectifying * (\e mu, μ), conformal (\e chi, χ), and authalic (\e xi, ξ) * latitudes for an ellipsoid of revolution. A latitude is represented by an * AuxAngle in order to maintain precision close to the poles. * * The class implements two methods for the conversion: * - Direct evaluation of the defining equations, the \e exact method. These * equations are formulated so as to preserve relative accuracy of the * tangent of the latitude, ensuring high accuracy near the equator and the * poles. Newton's method is used for those conversions that can't be * expressed in closed form. * - Expansions in powers of &e n, the third flattening, the \e series * method. This delivers full accuracy for abs(\e f) ≤ 1/150. Here, \e * f is the flattening of the ellipsoid. * * The series method is the preferred method of conversion for any conversion * involving μ, χ, or ξ, with abs(\e f) ≤ 1/150. The equations * for the conversions between φ, β, and θ are sufficiently * simple that the exact method should be used for such conversions and also * for conversions with with abs(\e f) > 1/150. * * @tparam T the floating-point type to use. * * Example of use: * \include example-AuxLatitude.cpp * * For more information on this projection, see \ref auxlat. **********************************************************************/ template class AuxLatitude { public: /** * The floating-point type for real numbers. This just connects to the * template parameters for the class. **********************************************************************/ typedef T real; /** * The type used to represent angles. **********************************************************************/ typedef AuxAngle angle; /** * The different auxiliary latitudes. **********************************************************************/ enum aux { /** * Geographic latitude, \e phi, φ * @hideinitializer **********************************************************************/ GEOGRAPHIC = 0, /** * Parametric latitude, \e beta, β * @hideinitializer **********************************************************************/ PARAMETRIC = 1, /** * %Geocentric latitude, \e theta, θ * @hideinitializer **********************************************************************/ GEOCENTRIC = 2, /** * Rectifying latitude, \e mu, μ * @hideinitializer **********************************************************************/ RECTIFYING = 3, /** * Conformal latitude, \e chi, χ * @hideinitializer **********************************************************************/ CONFORMAL = 4, /** * Authalic latitude, \e xi, ξ * @hideinitializer **********************************************************************/ AUTHALIC = 5, /** * The total number of auxiliary latitudes * @hideinitializer **********************************************************************/ AUXNUMBER = 6, /** * An alias for GEOGRAPHIC * @hideinitializer **********************************************************************/ COMMON = GEOGRAPHIC, /** * An alias for GEOGRAPHIC * @hideinitializer **********************************************************************/ GEODETIC = GEOGRAPHIC, /** * An alias for PARAMETRIC * @hideinitializer **********************************************************************/ REDUCED = PARAMETRIC, }; /** * Constructor * * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. * Negative \e f gives a prolate ellipsoid. * * \note the constructor does not precompute the coefficients for the * Fourier series for the series conversions. These are computed and saved * when first needed. **********************************************************************/ AuxLatitude(real f); /** * Constructor * * @param[in] a equatorial radius. * @param[in] b polar semi-axis. * * \note the constructor does not precompute the coefficients for the * Fourier series for the series conversions. These are computed and saved * when first needed. **********************************************************************/ AuxLatitude(real a, real b); /** * Convert between any two auxiliary latitudes. * * @param[in] auxin an AuxLatitude::aux indicating the type of * auxiliary latitude \e zeta. * @param[in] auxout an AuxLatitude::aux indicating the type of * auxiliary latitude \e eta. * @param[in] zeta the input auxiliary latitude. * @param[in] series if true use the Taylor series instead of the exact * equations [default false]. * @return the output auxiliary latitude \e eta. * * With \e series = true, the Fourier coefficients for a specific \e auxin * and \e auxout are computed and saved on the first call; the saved * coefficients are used on subsequent calls. The series method is * accurate for abs(\e f) ≤ 1/150; for other \e f, the exact method * should be used **********************************************************************/ angle Convert(int auxin, int auxout, const angle& zeta, bool series = false) const; /** * Convert geographic latitude to an auxiliary latitude \e eta. * * @param[in] auxout an AuxLatitude::aux indicating the auxiliary * latitude returned. * @param[in] phi the geographic latitude. * @param[out] diff optional pointer to the derivative d tan(\e eta) / d * tan(\e phi). * @return the auxiliary latitude \e eta. * * This uses the exact equations. **********************************************************************/ angle ToAuxiliary(int auxout, const angle& phi, real* diff = nullptr) const; /** * Convert an auxiliary latitude \e zeta to geographic latitude. * * @param[in] auxin an AuxLatitude::aux indicating the type of * auxiliary latitude \e zeta. * @param[in] zeta the input auxiliary latitude. * @param[out] niter optional pointer to the number of iterations. * @return the geographic latitude \e phi. * * This uses the exact equations. **********************************************************************/ angle FromAuxiliary(int auxin, const angle& zeta, int* niter = nullptr) const; /** * @return \e f, the flattening of the ellipsoid. **********************************************************************/ real Flattening() const { return _f; } /** * The order of the series expansions. This is set at compile time to * either 4, 6, or 8, by the preprocessor macro * GEOGRAPHICLIB_AUXLATITUDE_ORDER. * @hideinitializer **********************************************************************/ static const int Lmax = GEOGRAPHICLIB_AUXLATITUDE_ORDER; private: /** * Convert geographic latitude to parametric latitude * * @param[in] phi geographic latitude. * @param[out] diff optional pointer to the derivative d tan(\e beta) / d * tan(\e phi). * @return \e beta, the parametric latitude **********************************************************************/ angle Parametric(const angle& phi, real* diff = nullptr) const; /** * Convert geographic latitude to geocentric latitude * * @param[in] phi geographic latitude. * @param[out] diff optional pointer to the derivative d tan(\e theta) / d * tan(\e phi). * @return \e theta, the geocentric latitude. **********************************************************************/ angle Geocentric(const angle& phi, real* diff = nullptr) const; /** * Convert geographic latitude to rectifying latitude * * @param[in] phi geographic latitude. * @param[out] diff optional pointer to the derivative d tan(\e mu) / d * tan(\e phi). * @return \e mu, the rectifying latitude. **********************************************************************/ angle Rectifying(const angle& phi, real* diff = nullptr) const; /** * Convert geographic latitude to conformal latitude * * @param[in] phi geographic latitude. * @param[out] diff optional pointer to the derivative d tan(\e chi) / d * tan(\e phi). * @return \e chi, the conformal latitude. **********************************************************************/ angle Conformal(const angle& phi, real* diff = nullptr) const; /** * Convert geographic latitude to authalic latitude * * @param[in] phi geographic latitude. * @param[out] diff optional pointer to the derivative d tan(\e xi) / d * tan(\e phi). * @return \e xi, the authalic latitude. **********************************************************************/ angle Authalic(const angle& phi, real* diff = nullptr) const; // Maximum number of iterations for Newton's method static const int numit_ = 1000; real tol_, bmin_, bmax_; // Static consts for Newton's method // Ellipsoid parameters real _f, _fm1, _e2, _e2m1, _e12, _e12p1, _n, _e, _e1, _n2, _q; // To hold computed Fourier coefficients mutable real _c[Lmax * AUXNUMBER * AUXNUMBER]; // 1d index into AUXNUMBER x AUXNUMBER data static int ind(int auxout, int auxin) { return (auxout >= 0 && auxout < AUXNUMBER && auxin >= 0 && auxin < AUXNUMBER) ? AUXNUMBER * auxout + auxin : -1; } // the function sqrt(1 + tphi^2), convert tan to sec static real sc(real tphi) { using std::hypot; return hypot(real(1), tphi); } // the function tphi / sqrt(1 + tphi^2), convert tan to sin static real sn(real tphi) { using std::isfinite; using std::isnan; using std::copysign; return isfinite(tphi) || isnan(tphi) ? tphi / sc(tphi) : copysign(real(1), tphi); } // The symmetric elliptic integral RD static real RD(real x, real y, real z); // The symmetric elliptic integral RF static real RF(real x, real y, real z); // the function atanh(e * sphi)/e; works for e^2 = 0 and e^2 < 0 real atanhee(real tphi) const; // the function atanh(e * sphi)/e + sphi / (1 - (e * sphi)^2); real q(real tphi) const; // The divided difference of (q(1) - q(sphi)) / (1 - sphi) real Dq(real tphi) const; // Populate [_c[Lmax * k], _c[Lmax * (k + 1)]) void fillcoeff(int auxin, int auxout, int k) const; // Clenshaw applied to sum(c[k] * sin( (2*k+2) * zeta), i, 0, K-1) // if alt, use the Reinsch optimizations static real Clenshaw(real szeta, real czeta, const real c[], int K, bool alt = true); }; } // namespace GeographicLib #endif // AUXLATITUDE_HPP