00001 /** 00002 * \file Ellipsoid.hpp 00003 * \brief Header for GeographicLib::Ellipsoid class 00004 * 00005 * Copyright (c) Charles Karney (2012-2014) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_ELLIPSOID_HPP) 00011 #define GEOGRAPHICLIB_ELLIPSOID_HPP 1 00012 00013 #include <GeographicLib/Constants.hpp> 00014 #include <GeographicLib/TransverseMercator.hpp> 00015 #include <GeographicLib/EllipticFunction.hpp> 00016 #include <GeographicLib/AlbersEqualArea.hpp> 00017 00018 namespace GeographicLib { 00019 00020 /** 00021 * \brief Properties of an ellipsoid 00022 * 00023 * This class returns various properties of the ellipsoid and converts 00024 * between various types of latitudes. The latitude conversions are also 00025 * possible using the various projections supported by %GeographicLib; but 00026 * Ellipsoid provides more direct access (sometimes using private functions 00027 * of the projection classes). Ellipsoid::RectifyingLatitude, 00028 * Ellipsoid::InverseRectifyingLatitude, and Ellipsoid::MeridianDistance 00029 * provide functionality which can be provided by the Geodesic class. 00030 * However Geodesic uses a series approximation (valid for abs \e f < 1/150), 00031 * whereas Ellipsoid computes these quantities using EllipticFunction which 00032 * provides accurate results even when \e f is large. Use of this class 00033 * should be limited to −3 < \e f < 3/4 (i.e., 1/4 < b/a < 4). 00034 * 00035 * Example of use: 00036 * \include example-Ellipsoid.cpp 00037 **********************************************************************/ 00038 00039 class GEOGRAPHICLIB_EXPORT Ellipsoid { 00040 private: 00041 typedef Math::real real; 00042 static const int numit_ = 10; 00043 real stol_; 00044 real _a, _f, _f1, _f12, _e2, _e, _e12, _n, _b; 00045 TransverseMercator _tm; 00046 EllipticFunction _ell; 00047 AlbersEqualArea _au; 00048 static inline real tand(real x) { 00049 using std::abs; using std::tan; 00050 return 00051 abs(x) == real(90) ? (x < 0 ? 00052 - TransverseMercator::overflow() 00053 : TransverseMercator::overflow()) : 00054 tan(x * Math::degree()); 00055 } 00056 static inline real atand(real x) 00057 { using std::atan; return atan(x) / Math::degree(); } 00058 00059 // These are the alpha and beta coefficients in the Krueger series from 00060 // TransverseMercator. Thy are used by RhumbSolve to compute 00061 // (psi2-psi1)/(mu2-mu1). 00062 const Math::real* ConformalToRectifyingCoeffs() const { return _tm._alp; } 00063 const Math::real* RectifyingToConformalCoeffs() const { return _tm._bet; } 00064 friend class Rhumb; friend class RhumbLine; 00065 public: 00066 /** \name Constructor 00067 **********************************************************************/ 00068 ///@{ 00069 00070 /** 00071 * Constructor for a ellipsoid with 00072 * 00073 * @param[in] a equatorial radius (meters). 00074 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00075 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set 00076 * flattening to 1/\e f. 00077 * @exception GeographicErr if \e a or (1 − \e f) \e a is not 00078 * positive. 00079 **********************************************************************/ 00080 Ellipsoid(real a, real f); 00081 ///@} 00082 00083 /** \name %Ellipsoid dimensions. 00084 **********************************************************************/ 00085 ///@{ 00086 00087 /** 00088 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00089 * the value used in the constructor. 00090 **********************************************************************/ 00091 Math::real MajorRadius() const { return _a; } 00092 00093 /** 00094 * @return \e b the polar semi-axis (meters). 00095 **********************************************************************/ 00096 Math::real MinorRadius() const { return _b; } 00097 00098 /** 00099 * @return \e L the distance between the equator and a pole along a 00100 * meridian (meters). For a sphere \e L = (π/2) \e a. The radius 00101 * of a sphere with the same meridian length is \e L / (π/2). 00102 **********************************************************************/ 00103 Math::real QuarterMeridian() const; 00104 00105 /** 00106 * @return \e A the total area of the ellipsoid (meters<sup>2</sup>). For 00107 * a sphere \e A = 4π <i>a</i><sup>2</sup>. The radius of a sphere 00108 * with the same area is sqrt(\e A / (4π)). 00109 **********************************************************************/ 00110 Math::real Area() const; 00111 00112 /** 00113 * @return \e V the total volume of the ellipsoid (meters<sup>3</sup>). 00114 * For a sphere \e V = (4π / 3) <i>a</i><sup>3</sup>. The radius of 00115 * a sphere with the same volume is cbrt(\e V / (4π/3)). 00116 **********************************************************************/ 00117 Math::real Volume() const 00118 { return (4 * Math::pi()) * Math::sq(_a) * _b / 3; } 00119 ///@} 00120 00121 /** \name %Ellipsoid shape 00122 **********************************************************************/ 00123 ///@{ 00124 00125 /** 00126 * @return \e f = (\e a − \e b) / \e a, the flattening of the 00127 * ellipsoid. This is the value used in the constructor. This is zero, 00128 * positive, or negative for a sphere, oblate ellipsoid, or prolate 00129 * ellipsoid. 00130 **********************************************************************/ 00131 Math::real Flattening() const { return _f; } 00132 00133 /** 00134 * @return \e f ' = (\e a − \e b) / \e b, the second flattening of 00135 * the ellipsoid. This is zero, positive, or negative for a sphere, 00136 * oblate ellipsoid, or prolate ellipsoid. 00137 **********************************************************************/ 00138 Math::real SecondFlattening() const { return _f / (1 - _f); } 00139 00140 /** 00141 * @return \e n = (\e a − \e b) / (\e a + \e b), the third flattening 00142 * of the ellipsoid. This is zero, positive, or negative for a sphere, 00143 * oblate ellipsoid, or prolate ellipsoid. 00144 **********************************************************************/ 00145 Math::real ThirdFlattening() const { return _n; } 00146 00147 /** 00148 * @return <i>e</i><sup>2</sup> = (<i>a</i><sup>2</sup> − 00149 * <i>b</i><sup>2</sup>) / <i>a</i><sup>2</sup>, the eccentricity squared 00150 * of the ellipsoid. This is zero, positive, or negative for a sphere, 00151 * oblate ellipsoid, or prolate ellipsoid. 00152 **********************************************************************/ 00153 Math::real EccentricitySq() const { return _e2; } 00154 00155 /** 00156 * @return <i>e'</i> <sup>2</sup> = (<i>a</i><sup>2</sup> − 00157 * <i>b</i><sup>2</sup>) / <i>b</i><sup>2</sup>, the second eccentricity 00158 * squared of the ellipsoid. This is zero, positive, or negative for a 00159 * sphere, oblate ellipsoid, or prolate ellipsoid. 00160 **********************************************************************/ 00161 Math::real SecondEccentricitySq() const { return _e12; } 00162 00163 /** 00164 * @return <i>e''</i> <sup>2</sup> = (<i>a</i><sup>2</sup> − 00165 * <i>b</i><sup>2</sup>) / (<i>a</i><sup>2</sup> + <i>b</i><sup>2</sup>), 00166 * the third eccentricity squared of the ellipsoid. This is zero, 00167 * positive, or negative for a sphere, oblate ellipsoid, or prolate 00168 * ellipsoid. 00169 **********************************************************************/ 00170 Math::real ThirdEccentricitySq() const { return _e2 / (2 - _e2); } 00171 ///@} 00172 00173 /** \name Latitude conversion. 00174 **********************************************************************/ 00175 ///@{ 00176 00177 /** 00178 * @param[in] phi the geographic latitude (degrees). 00179 * @return β the parametric latitude (degrees). 00180 * 00181 * The geographic latitude, φ, is the angle beween the equatorial 00182 * plane and a vector normal to the surface of the ellipsoid. 00183 * 00184 * The parametric latitude (also called the reduced latitude), β, 00185 * allows the cartesian coordinated of a meridian to be expressed 00186 * conveniently in parametric form as 00187 * - \e R = \e a cos β 00188 * - \e Z = \e b sin β 00189 * . 00190 * where \e a and \e b are the equatorial radius and the polar semi-axis. 00191 * For a sphere β = φ. 00192 * 00193 * φ must lie in the range [−90°, 90°]; the 00194 * result is undefined if this condition does not hold. The returned value 00195 * β lies in [−90°, 90°]. 00196 **********************************************************************/ 00197 Math::real ParametricLatitude(real phi) const; 00198 00199 /** 00200 * @param[in] beta the parametric latitude (degrees). 00201 * @return φ the geographic latitude (degrees). 00202 * 00203 * β must lie in the range [−90°, 90°]; the 00204 * result is undefined if this condition does not hold. The returned value 00205 * φ lies in [−90°, 90°]. 00206 **********************************************************************/ 00207 Math::real InverseParametricLatitude(real beta) const; 00208 00209 /** 00210 * @param[in] phi the geographic latitude (degrees). 00211 * @return θ the geocentric latitude (degrees). 00212 * 00213 * The geocentric latitude, θ, is the angle beween the equatorial 00214 * plane and a line between the center of the ellipsoid and a point on the 00215 * ellipsoid. For a sphere θ = φ. 00216 * 00217 * φ must lie in the range [−90°, 90°]; the 00218 * result is undefined if this condition does not hold. The returned value 00219 * θ lies in [−90°, 90°]. 00220 **********************************************************************/ 00221 Math::real GeocentricLatitude(real phi) const; 00222 00223 /** 00224 * @param[in] theta the geocentric latitude (degrees). 00225 * @return φ the geographic latitude (degrees). 00226 * 00227 * θ must lie in the range [−90°, 90°]; the 00228 * result is undefined if this condition does not hold. The returned value 00229 * φ lies in [−90°, 90°]. 00230 **********************************************************************/ 00231 Math::real InverseGeocentricLatitude(real theta) const; 00232 00233 /** 00234 * @param[in] phi the geographic latitude (degrees). 00235 * @return μ the rectifying latitude (degrees). 00236 * 00237 * The rectifying latitude, μ, has the property that the distance along 00238 * a meridian of the ellipsoid between two points with rectifying latitudes 00239 * μ<sub>1</sub> and μ<sub>2</sub> is equal to 00240 * (μ<sub>2</sub> - μ<sub>1</sub>) \e L / 90°, 00241 * where \e L = QuarterMeridian(). For a sphere μ = φ. 00242 * 00243 * φ must lie in the range [−90°, 90°]; the 00244 * result is undefined if this condition does not hold. The returned value 00245 * μ lies in [−90°, 90°]. 00246 **********************************************************************/ 00247 Math::real RectifyingLatitude(real phi) const; 00248 00249 /** 00250 * @param[in] mu the rectifying latitude (degrees). 00251 * @return φ the geographic latitude (degrees). 00252 * 00253 * μ must lie in the range [−90°, 90°]; the 00254 * result is undefined if this condition does not hold. The returned value 00255 * φ lies in [−90°, 90°]. 00256 **********************************************************************/ 00257 Math::real InverseRectifyingLatitude(real mu) const; 00258 00259 /** 00260 * @param[in] phi the geographic latitude (degrees). 00261 * @return ξ the authalic latitude (degrees). 00262 * 00263 * The authalic latitude, ξ, has the property that the area of the 00264 * ellipsoid between two circles with authalic latitudes 00265 * ξ<sub>1</sub> and ξ<sub>2</sub> is equal to (sin 00266 * ξ<sub>2</sub> - sin ξ<sub>1</sub>) \e A / 2, where \e A 00267 * = Area(). For a sphere ξ = φ. 00268 * 00269 * φ must lie in the range [−90°, 90°]; the 00270 * result is undefined if this condition does not hold. The returned value 00271 * ξ lies in [−90°, 90°]. 00272 **********************************************************************/ 00273 Math::real AuthalicLatitude(real phi) const; 00274 00275 /** 00276 * @param[in] xi the authalic latitude (degrees). 00277 * @return φ the geographic latitude (degrees). 00278 * 00279 * ξ must lie in the range [−90°, 90°]; the 00280 * result is undefined if this condition does not hold. The returned value 00281 * φ lies in [−90°, 90°]. 00282 **********************************************************************/ 00283 Math::real InverseAuthalicLatitude(real xi) const; 00284 00285 /** 00286 * @param[in] phi the geographic latitude (degrees). 00287 * @return χ the conformal latitude (degrees). 00288 * 00289 * The conformal latitude, χ, gives the mapping of the ellipsoid to a 00290 * sphere which which is conformal (angles are preserved) and in which the 00291 * equator of the ellipsoid maps to the equator of the sphere. For a 00292 * sphere χ = φ. 00293 * 00294 * φ must lie in the range [−90°, 90°]; the 00295 * result is undefined if this condition does not hold. The returned value 00296 * χ lies in [−90°, 90°]. 00297 **********************************************************************/ 00298 Math::real ConformalLatitude(real phi) const; 00299 00300 /** 00301 * @param[in] chi the conformal latitude (degrees). 00302 * @return φ the geographic latitude (degrees). 00303 * 00304 * χ must lie in the range [−90°, 90°]; the 00305 * result is undefined if this condition does not hold. The returned value 00306 * φ lies in [−90°, 90°]. 00307 **********************************************************************/ 00308 Math::real InverseConformalLatitude(real chi) const; 00309 00310 /** 00311 * @param[in] phi the geographic latitude (degrees). 00312 * @return ψ the isometric latitude (degrees). 00313 * 00314 * The isometric latitude gives the mapping of the ellipsoid to a plane 00315 * which which is conformal (angles are preserved) and in which the equator 00316 * of the ellipsoid maps to a straight line of constant scale; this mapping 00317 * defines the Mercator projection. For a sphere ψ = 00318 * sinh<sup>−1</sup> tan φ. 00319 * 00320 * φ must lie in the range [−90°, 90°]; the result is 00321 * undefined if this condition does not hold. The value returned for φ 00322 * = ±90° is some (positive or negative) large but finite value, 00323 * such that InverseIsometricLatitude returns the original value of φ. 00324 **********************************************************************/ 00325 Math::real IsometricLatitude(real phi) const; 00326 00327 /** 00328 * @param[in] psi the isometric latitude (degrees). 00329 * @return φ the geographic latitude (degrees). 00330 * 00331 * The returned value φ lies in [−90°, 90°]. For a 00332 * sphere φ = tan<sup>−1</sup> sinh ψ. 00333 **********************************************************************/ 00334 Math::real InverseIsometricLatitude(real psi) const; 00335 ///@} 00336 00337 /** \name Other quantities. 00338 **********************************************************************/ 00339 ///@{ 00340 00341 /** 00342 * @param[in] phi the geographic latitude (degrees). 00343 * @return \e R = \e a cos β the radius of a circle of latitude 00344 * φ (meters). \e R (π/180°) gives meters per degree 00345 * longitude measured along a circle of latitude. 00346 * 00347 * φ must lie in the range [−90°, 90°]; the 00348 * result is undefined if this condition does not hold. 00349 **********************************************************************/ 00350 Math::real CircleRadius(real phi) const; 00351 00352 /** 00353 * @param[in] phi the geographic latitude (degrees). 00354 * @return \e Z = \e b sin β the distance of a circle of latitude 00355 * φ from the equator measured parallel to the ellipsoid axis 00356 * (meters). 00357 * 00358 * φ must lie in the range [−90°, 90°]; the 00359 * result is undefined if this condition does not hold. 00360 **********************************************************************/ 00361 Math::real CircleHeight(real phi) const; 00362 00363 /** 00364 * @param[in] phi the geographic latitude (degrees). 00365 * @return \e s the distance along a meridian 00366 * between the equator and a point of latitude φ (meters). \e s is 00367 * given by \e s = μ \e L / 90°, where \e L = 00368 * QuarterMeridian()). 00369 * 00370 * φ must lie in the range [−90°, 90°]; the 00371 * result is undefined if this condition does not hold. 00372 **********************************************************************/ 00373 Math::real MeridianDistance(real phi) const; 00374 00375 /** 00376 * @param[in] phi the geographic latitude (degrees). 00377 * @return ρ the meridional radius of curvature of the ellipsoid at 00378 * latitude φ (meters); this is the curvature of the meridian. \e 00379 * rho is given by ρ = (180°/π) d\e s / dφ, 00380 * where \e s = MeridianDistance(); thus ρ (π/180°) 00381 * gives meters per degree latitude measured along a meridian. 00382 * 00383 * φ must lie in the range [−90°, 90°]; the 00384 * result is undefined if this condition does not hold. 00385 **********************************************************************/ 00386 Math::real MeridionalCurvatureRadius(real phi) const; 00387 00388 /** 00389 * @param[in] phi the geographic latitude (degrees). 00390 * @return ν the transverse radius of curvature of the ellipsoid at 00391 * latitude φ (meters); this is the curvature of a curve on the 00392 * ellipsoid which also lies in a plane perpendicular to the ellipsoid 00393 * and to the meridian. ν is related to \e R = CircleRadius() by \e 00394 * R = ν cos φ. 00395 * 00396 * φ must lie in the range [−90°, 90°]; the 00397 * result is undefined if this condition does not hold. 00398 **********************************************************************/ 00399 Math::real TransverseCurvatureRadius(real phi) const; 00400 00401 /** 00402 * @param[in] phi the geographic latitude (degrees). 00403 * @param[in] azi the angle between the meridian and the normal section 00404 * (degrees). 00405 * @return the radius of curvature of the ellipsoid in the normal 00406 * section at latitude φ inclined at an angle \e azi to the 00407 * meridian (meters). 00408 * 00409 * φ must lie in the range [−90°, 90°] and \e 00410 * azi must lie in the range [−540°, 540°); the 00411 * result is undefined if either of conditions does not hold. 00412 **********************************************************************/ 00413 Math::real NormalCurvatureRadius(real phi, real azi) const; 00414 ///@} 00415 00416 /** \name Eccentricity conversions. 00417 **********************************************************************/ 00418 ///@{ 00419 00420 /** 00421 * @param[in] fp = \e f ' = (\e a − \e b) / \e b, the second 00422 * flattening. 00423 * @return \e f = (\e a − \e b) / \e a, the flattening. 00424 * 00425 * \e f ' should lie in (−1, ∞). 00426 * The returned value \e f lies in (−∞, 1). 00427 **********************************************************************/ 00428 static Math::real SecondFlatteningToFlattening(real fp) 00429 { return fp / (1 + fp); } 00430 00431 /** 00432 * @param[in] f = (\e a − \e b) / \e a, the flattening. 00433 * @return \e f ' = (\e a − \e b) / \e b, the second flattening. 00434 * 00435 * \e f should lie in (−∞, 1). 00436 * The returned value \e f ' lies in (−1, ∞). 00437 **********************************************************************/ 00438 static Math::real FlatteningToSecondFlattening(real f) 00439 { return f / (1 - f); } 00440 00441 /** 00442 * @param[in] n = (\e a − \e b) / (\e a + \e b), the third 00443 * flattening. 00444 * @return \e f = (\e a − \e b) / \e a, the flattening. 00445 * 00446 * \e n should lie in (−1, 1). 00447 * The returned value \e f lies in (−∞, 1). 00448 **********************************************************************/ 00449 static Math::real ThirdFlatteningToFlattening(real n) 00450 { return 2 * n / (1 + n); } 00451 00452 /** 00453 * @param[in] f = (\e a − \e b) / \e a, the flattening. 00454 * @return \e n = (\e a − \e b) / (\e a + \e b), the third 00455 * flattening. 00456 * 00457 * \e f should lie in (−∞, 1). 00458 * The returned value \e n lies in (−1, 1). 00459 **********************************************************************/ 00460 static Math::real FlatteningToThirdFlattening(real f) 00461 { return f / (2 - f); } 00462 00463 /** 00464 * @param[in] e2 = <i>e</i><sup>2</sup> = (<i>a</i><sup>2</sup> − 00465 * <i>b</i><sup>2</sup>) / <i>a</i><sup>2</sup>, the eccentricity 00466 * squared. 00467 * @return \e f = (\e a − \e b) / \e a, the flattening. 00468 * 00469 * <i>e</i><sup>2</sup> should lie in (−∞, 1). 00470 * The returned value \e f lies in (−∞, 1). 00471 **********************************************************************/ 00472 static Math::real EccentricitySqToFlattening(real e2) 00473 { using std::sqrt; return e2 / (sqrt(1 - e2) + 1); } 00474 00475 /** 00476 * @param[in] f = (\e a − \e b) / \e a, the flattening. 00477 * @return <i>e</i><sup>2</sup> = (<i>a</i><sup>2</sup> − 00478 * <i>b</i><sup>2</sup>) / <i>a</i><sup>2</sup>, the eccentricity 00479 * squared. 00480 * 00481 * \e f should lie in (−∞, 1). 00482 * The returned value <i>e</i><sup>2</sup> lies in (−∞, 1). 00483 **********************************************************************/ 00484 static Math::real FlatteningToEccentricitySq(real f) 00485 { return f * (2 - f); } 00486 00487 /** 00488 * @param[in] ep2 = <i>e'</i> <sup>2</sup> = (<i>a</i><sup>2</sup> − 00489 * <i>b</i><sup>2</sup>) / <i>b</i><sup>2</sup>, the second eccentricity 00490 * squared. 00491 * @return \e f = (\e a − \e b) / \e a, the flattening. 00492 * 00493 * <i>e'</i> <sup>2</sup> should lie in (−1, ∞). 00494 * The returned value \e f lies in (−∞, 1). 00495 **********************************************************************/ 00496 static Math::real SecondEccentricitySqToFlattening(real ep2) 00497 { using std::sqrt; return ep2 / (sqrt(1 + ep2) + 1 + ep2); } 00498 00499 /** 00500 * @param[in] f = (\e a − \e b) / \e a, the flattening. 00501 * @return <i>e'</i> <sup>2</sup> = (<i>a</i><sup>2</sup> − 00502 * <i>b</i><sup>2</sup>) / <i>b</i><sup>2</sup>, the second eccentricity 00503 * squared. 00504 * 00505 * \e f should lie in (−∞, 1). 00506 * The returned value <i>e'</i> <sup>2</sup> lies in (−1, ∞). 00507 **********************************************************************/ 00508 static Math::real FlatteningToSecondEccentricitySq(real f) 00509 { return f * (2 - f) / Math::sq(1 - f); } 00510 00511 /** 00512 * @param[in] epp2 = <i>e''</i> <sup>2</sup> = (<i>a</i><sup>2</sup> 00513 * − <i>b</i><sup>2</sup>) / (<i>a</i><sup>2</sup> + 00514 * <i>b</i><sup>2</sup>), the third eccentricity squared. 00515 * @return \e f = (\e a − \e b) / \e a, the flattening. 00516 * 00517 * <i>e''</i> <sup>2</sup> should lie in (−1, 1). 00518 * The returned value \e f lies in (−∞, 1). 00519 **********************************************************************/ 00520 static Math::real ThirdEccentricitySqToFlattening(real epp2) 00521 { return 2 * epp2 / (sqrt((1 - epp2) * (1 + epp2)) + 1 + epp2); } 00522 00523 /** 00524 * @param[in] f = (\e a − \e b) / \e a, the flattening. 00525 * @return <i>e''</i> <sup>2</sup> = (<i>a</i><sup>2</sup> − 00526 * <i>b</i><sup>2</sup>) / (<i>a</i><sup>2</sup> + <i>b</i><sup>2</sup>), 00527 * the third eccentricity squared. 00528 * 00529 * \e f should lie in (−∞, 1). 00530 * The returned value <i>e''</i> <sup>2</sup> lies in (−1, 1). 00531 **********************************************************************/ 00532 static Math::real FlatteningToThirdEccentricitySq(real f) 00533 { return f * (2 - f) / (1 + Math::sq(1 - f)); } 00534 00535 ///@} 00536 00537 /** 00538 * A global instantiation of Ellipsoid with the parameters for the WGS84 00539 * ellipsoid. 00540 **********************************************************************/ 00541 static const Ellipsoid& WGS84(); 00542 }; 00543 00544 } // namespace GeographicLib 00545 00546 #endif // GEOGRAPHICLIB_ELLIPSOID_HPP