/** * \file Ellipsoid.hpp * \brief Header for GeographicLib::Ellipsoid class * * Copyright (c) Charles Karney (2012-2020) and licensed * under the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ **********************************************************************/ #if !defined(GEOGRAPHICLIB_ELLIPSOID_HPP) #define GEOGRAPHICLIB_ELLIPSOID_HPP 1 #include #include #include #include namespace GeographicLib { /** * \brief Properties of an ellipsoid * * This class returns various properties of the ellipsoid and converts * between various types of latitudes. The latitude conversions are also * possible using the various projections supported by %GeographicLib; but * Ellipsoid provides more direct access (sometimes using private functions * of the projection classes). Ellipsoid::RectifyingLatitude, * Ellipsoid::InverseRectifyingLatitude, and Ellipsoid::MeridianDistance * provide functionality which can be provided by the Geodesic class. * However Geodesic uses a series approximation (valid for abs \e f < 1/150), * whereas Ellipsoid computes these quantities using EllipticFunction which * provides accurate results even when \e f is large. Use of this class * should be limited to −3 < \e f < 3/4 (i.e., 1/4 < b/a < 4). * * Example of use: * \include example-Ellipsoid.cpp **********************************************************************/ class GEOGRAPHICLIB_EXPORT Ellipsoid { private: typedef Math::real real; static const int numit_ = 10; real stol_; real _a, _f, _f1, _f12, _e2, _es, _e12, _n, _b; TransverseMercator _tm; EllipticFunction _ell; AlbersEqualArea _au; // These are the alpha and beta coefficients in the Krueger series from // TransverseMercator. Thy are used by RhumbSolve to compute // (psi2-psi1)/(mu2-mu1). const Math::real* ConformalToRectifyingCoeffs() const { return _tm._alp; } const Math::real* RectifyingToConformalCoeffs() const { return _tm._bet; } friend class Rhumb; friend class RhumbLine; public: /** \name Constructor **********************************************************************/ ///@{ /** * Constructor for a ellipsoid with * * @param[in] a equatorial radius (meters). * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. * Negative \e f gives a prolate ellipsoid. * @exception GeographicErr if \e a or (1 − \e f) \e a is not * positive. **********************************************************************/ Ellipsoid(real a, real f); ///@} /** \name %Ellipsoid dimensions. **********************************************************************/ ///@{ /** * @return \e a the equatorial radius of the ellipsoid (meters). This is * the value used in the constructor. **********************************************************************/ Math::real EquatorialRadius() const { return _a; } /** * @return \e b the polar semi-axis (meters). **********************************************************************/ Math::real MinorRadius() const { return _b; } /** * @return \e L the distance between the equator and a pole along a * meridian (meters). For a sphere \e L = (π/2) \e a. The radius * of a sphere with the same meridian length is \e L / (π/2). **********************************************************************/ Math::real QuarterMeridian() const; /** * @return \e A the total area of the ellipsoid (meters²). For * a sphere \e A = 4π a². The radius of a sphere * with the same area is sqrt(\e A / (4π)). **********************************************************************/ Math::real Area() const; /** * @return \e V the total volume of the ellipsoid (meters³). * For a sphere \e V = (4π / 3) a³. The radius of * a sphere with the same volume is cbrt(\e V / (4π/3)). **********************************************************************/ Math::real Volume() const { return (4 * Math::pi()) * Math::sq(_a) * _b / 3; } /** * \deprecated An old name for EquatorialRadius(). **********************************************************************/ GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()") Math::real MajorRadius() const { return EquatorialRadius(); } ///@} /** \name %Ellipsoid shape **********************************************************************/ ///@{ /** * @return \e f = (\e a − \e b) / \e a, the flattening of the * ellipsoid. This is the value used in the constructor. This is zero, * positive, or negative for a sphere, oblate ellipsoid, or prolate * ellipsoid. **********************************************************************/ Math::real Flattening() const { return _f; } /** * @return \e f ' = (\e a − \e b) / \e b, the second flattening of * the ellipsoid. This is zero, positive, or negative for a sphere, * oblate ellipsoid, or prolate ellipsoid. **********************************************************************/ Math::real SecondFlattening() const { return _f / (1 - _f); } /** * @return \e n = (\e a − \e b) / (\e a + \e b), the third flattening * of the ellipsoid. This is zero, positive, or negative for a sphere, * oblate ellipsoid, or prolate ellipsoid. **********************************************************************/ Math::real ThirdFlattening() const { return _n; } /** * @return e² = (a² − * b²) / a², the eccentricity squared * of the ellipsoid. This is zero, positive, or negative for a sphere, * oblate ellipsoid, or prolate ellipsoid. **********************************************************************/ Math::real EccentricitySq() const { return _e2; } /** * @return e' ² = (a² − * b²) / b², the second eccentricity * squared of the ellipsoid. This is zero, positive, or negative for a * sphere, oblate ellipsoid, or prolate ellipsoid. **********************************************************************/ Math::real SecondEccentricitySq() const { return _e12; } /** * @return e'' ² = (a² − * b²) / (a² + b²), * the third eccentricity squared of the ellipsoid. This is zero, * positive, or negative for a sphere, oblate ellipsoid, or prolate * ellipsoid. **********************************************************************/ Math::real ThirdEccentricitySq() const { return _e2 / (2 - _e2); } ///@} /** \name Latitude conversion. **********************************************************************/ ///@{ /** * @param[in] phi the geographic latitude (degrees). * @return β the parametric latitude (degrees). * * The geographic latitude, φ, is the angle between the equatorial * plane and a vector normal to the surface of the ellipsoid. * * The parametric latitude (also called the reduced latitude), β, * allows the cartesian coordinated of a meridian to be expressed * conveniently in parametric form as * - \e R = \e a cos β * - \e Z = \e b sin β * . * where \e a and \e b are the equatorial radius and the polar semi-axis. * For a sphere β = φ. * * φ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. The returned value * β lies in [−90°, 90°]. **********************************************************************/ Math::real ParametricLatitude(real phi) const; /** * @param[in] beta the parametric latitude (degrees). * @return φ the geographic latitude (degrees). * * β must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. The returned value * φ lies in [−90°, 90°]. **********************************************************************/ Math::real InverseParametricLatitude(real beta) const; /** * @param[in] phi the geographic latitude (degrees). * @return θ the geocentric latitude (degrees). * * The geocentric latitude, θ, is the angle between the equatorial * plane and a line between the center of the ellipsoid and a point on the * ellipsoid. For a sphere θ = φ. * * φ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. The returned value * θ lies in [−90°, 90°]. **********************************************************************/ Math::real GeocentricLatitude(real phi) const; /** * @param[in] theta the geocentric latitude (degrees). * @return φ the geographic latitude (degrees). * * θ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. The returned value * φ lies in [−90°, 90°]. **********************************************************************/ Math::real InverseGeocentricLatitude(real theta) const; /** * @param[in] phi the geographic latitude (degrees). * @return μ the rectifying latitude (degrees). * * The rectifying latitude, μ, has the property that the distance along * a meridian of the ellipsoid between two points with rectifying latitudes * μ₁ and μ₂ is equal to * (μ₂ - μ₁) \e L / 90°, * where \e L = QuarterMeridian(). For a sphere μ = φ. * * φ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. The returned value * μ lies in [−90°, 90°]. **********************************************************************/ Math::real RectifyingLatitude(real phi) const; /** * @param[in] mu the rectifying latitude (degrees). * @return φ the geographic latitude (degrees). * * μ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. The returned value * φ lies in [−90°, 90°]. **********************************************************************/ Math::real InverseRectifyingLatitude(real mu) const; /** * @param[in] phi the geographic latitude (degrees). * @return ξ the authalic latitude (degrees). * * The authalic latitude, ξ, has the property that the area of the * ellipsoid between two circles with authalic latitudes * ξ₁ and ξ₂ is equal to (sin * ξ₂ - sin ξ₁) \e A / 2, where \e A * = Area(). For a sphere ξ = φ. * * φ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. The returned value * ξ lies in [−90°, 90°]. **********************************************************************/ Math::real AuthalicLatitude(real phi) const; /** * @param[in] xi the authalic latitude (degrees). * @return φ the geographic latitude (degrees). * * ξ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. The returned value * φ lies in [−90°, 90°]. **********************************************************************/ Math::real InverseAuthalicLatitude(real xi) const; /** * @param[in] phi the geographic latitude (degrees). * @return χ the conformal latitude (degrees). * * The conformal latitude, χ, gives the mapping of the ellipsoid to a * sphere which which is conformal (angles are preserved) and in which the * equator of the ellipsoid maps to the equator of the sphere. For a * sphere χ = φ. * * φ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. The returned value * χ lies in [−90°, 90°]. **********************************************************************/ Math::real ConformalLatitude(real phi) const; /** * @param[in] chi the conformal latitude (degrees). * @return φ the geographic latitude (degrees). * * χ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. The returned value * φ lies in [−90°, 90°]. **********************************************************************/ Math::real InverseConformalLatitude(real chi) const; /** * @param[in] phi the geographic latitude (degrees). * @return ψ the isometric latitude (degrees). * * The isometric latitude gives the mapping of the ellipsoid to a plane * which which is conformal (angles are preserved) and in which the equator * of the ellipsoid maps to a straight line of constant scale; this mapping * defines the Mercator projection. For a sphere ψ = * sinh⁻¹ tan φ. * * φ must lie in the range [−90°, 90°]; the result is * undefined if this condition does not hold. The value returned for φ * = ±90° is some (positive or negative) large but finite value, * such that InverseIsometricLatitude returns the original value of φ. **********************************************************************/ Math::real IsometricLatitude(real phi) const; /** * @param[in] psi the isometric latitude (degrees). * @return φ the geographic latitude (degrees). * * The returned value φ lies in [−90°, 90°]. For a * sphere φ = tan⁻¹ sinh ψ. **********************************************************************/ Math::real InverseIsometricLatitude(real psi) const; ///@} /** \name Other quantities. **********************************************************************/ ///@{ /** * @param[in] phi the geographic latitude (degrees). * @return \e R = \e a cos β the radius of a circle of latitude * φ (meters). \e R (π/180°) gives meters per degree * longitude measured along a circle of latitude. * * φ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. **********************************************************************/ Math::real CircleRadius(real phi) const; /** * @param[in] phi the geographic latitude (degrees). * @return \e Z = \e b sin β the distance of a circle of latitude * φ from the equator measured parallel to the ellipsoid axis * (meters). * * φ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. **********************************************************************/ Math::real CircleHeight(real phi) const; /** * @param[in] phi the geographic latitude (degrees). * @return \e s the distance along a meridian * between the equator and a point of latitude φ (meters). \e s is * given by \e s = μ \e L / 90°, where \e L = * QuarterMeridian()). * * φ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. **********************************************************************/ Math::real MeridianDistance(real phi) const; /** * @param[in] phi the geographic latitude (degrees). * @return ρ the meridional radius of curvature of the ellipsoid at * latitude φ (meters); this is the curvature of the meridian. \e * rho is given by ρ = (180°/π) d\e s / dφ, * where \e s = MeridianDistance(); thus ρ (π/180°) * gives meters per degree latitude measured along a meridian. * * φ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. **********************************************************************/ Math::real MeridionalCurvatureRadius(real phi) const; /** * @param[in] phi the geographic latitude (degrees). * @return ν the transverse radius of curvature of the ellipsoid at * latitude φ (meters); this is the curvature of a curve on the * ellipsoid which also lies in a plane perpendicular to the ellipsoid * and to the meridian. ν is related to \e R = CircleRadius() by \e * R = ν cos φ. * * φ must lie in the range [−90°, 90°]; the * result is undefined if this condition does not hold. **********************************************************************/ Math::real TransverseCurvatureRadius(real phi) const; /** * @param[in] phi the geographic latitude (degrees). * @param[in] azi the angle between the meridian and the normal section * (degrees). * @return the radius of curvature of the ellipsoid in the normal * section at latitude φ inclined at an angle \e azi to the * meridian (meters). * * φ must lie in the range [−90°, 90°]; the result is * undefined this condition does not hold. **********************************************************************/ Math::real NormalCurvatureRadius(real phi, real azi) const; ///@} /** \name Eccentricity conversions. **********************************************************************/ ///@{ /** * @param[in] fp = \e f ' = (\e a − \e b) / \e b, the second * flattening. * @return \e f = (\e a − \e b) / \e a, the flattening. * * \e f ' should lie in (−1, ∞). * The returned value \e f lies in (−∞, 1). **********************************************************************/ static Math::real SecondFlatteningToFlattening(real fp) { return fp / (1 + fp); } /** * @param[in] f = (\e a − \e b) / \e a, the flattening. * @return \e f ' = (\e a − \e b) / \e b, the second flattening. * * \e f should lie in (−∞, 1). * The returned value \e f ' lies in (−1, ∞). **********************************************************************/ static Math::real FlatteningToSecondFlattening(real f) { return f / (1 - f); } /** * @param[in] n = (\e a − \e b) / (\e a + \e b), the third * flattening. * @return \e f = (\e a − \e b) / \e a, the flattening. * * \e n should lie in (−1, 1). * The returned value \e f lies in (−∞, 1). **********************************************************************/ static Math::real ThirdFlatteningToFlattening(real n) { return 2 * n / (1 + n); } /** * @param[in] f = (\e a − \e b) / \e a, the flattening. * @return \e n = (\e a − \e b) / (\e a + \e b), the third * flattening. * * \e f should lie in (−∞, 1). * The returned value \e n lies in (−1, 1). **********************************************************************/ static Math::real FlatteningToThirdFlattening(real f) { return f / (2 - f); } /** * @param[in] e2 = e² = (a² − * b²) / a², the eccentricity * squared. * @return \e f = (\e a − \e b) / \e a, the flattening. * * e² should lie in (−∞, 1). * The returned value \e f lies in (−∞, 1). **********************************************************************/ static Math::real EccentricitySqToFlattening(real e2) { using std::sqrt; return e2 / (sqrt(1 - e2) + 1); } /** * @param[in] f = (\e a − \e b) / \e a, the flattening. * @return e² = (a² − * b²) / a², the eccentricity * squared. * * \e f should lie in (−∞, 1). * The returned value e² lies in (−∞, 1). **********************************************************************/ static Math::real FlatteningToEccentricitySq(real f) { return f * (2 - f); } /** * @param[in] ep2 = e' ² = (a² − * b²) / b², the second eccentricity * squared. * @return \e f = (\e a − \e b) / \e a, the flattening. * * e' ² should lie in (−1, ∞). * The returned value \e f lies in (−∞, 1). **********************************************************************/ static Math::real SecondEccentricitySqToFlattening(real ep2) { using std::sqrt; return ep2 / (sqrt(1 + ep2) + 1 + ep2); } /** * @param[in] f = (\e a − \e b) / \e a, the flattening. * @return e' ² = (a² − * b²) / b², the second eccentricity * squared. * * \e f should lie in (−∞, 1). * The returned value e' ² lies in (−1, ∞). **********************************************************************/ static Math::real FlatteningToSecondEccentricitySq(real f) { return f * (2 - f) / Math::sq(1 - f); } /** * @param[in] epp2 = e'' ² = (a² * − b²) / (a² + * b²), the third eccentricity squared. * @return \e f = (\e a − \e b) / \e a, the flattening. * * e'' ² should lie in (−1, 1). * The returned value \e f lies in (−∞, 1). **********************************************************************/ static Math::real ThirdEccentricitySqToFlattening(real epp2) { using std::sqrt; return 2 * epp2 / (sqrt((1 - epp2) * (1 + epp2)) + 1 + epp2); } /** * @param[in] f = (\e a − \e b) / \e a, the flattening. * @return e'' ² = (a² − * b²) / (a² + b²), * the third eccentricity squared. * * \e f should lie in (−∞, 1). * The returned value e'' ² lies in (−1, 1). **********************************************************************/ static Math::real FlatteningToThirdEccentricitySq(real f) { return f * (2 - f) / (1 + Math::sq(1 - f)); } ///@} /** * A global instantiation of Ellipsoid with the parameters for the WGS84 * ellipsoid. **********************************************************************/ static const Ellipsoid& WGS84(); }; } // namespace GeographicLib #endif // GEOGRAPHICLIB_ELLIPSOID_HPP