00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <GeographicLib/Ellipsoid.hpp>
00011
00012 namespace GeographicLib {
00013
00014 using namespace std;
00015
00016 Ellipsoid::Ellipsoid(real a, real f)
00017 : stol_(real(0.01) * sqrt(numeric_limits<real>::epsilon()))
00018 , _a(a)
00019 , _f(f <= 1 ? f : 1/f)
00020 , _f1(1 - _f)
00021 , _f12(Math::sq(_f1))
00022 , _e2(_f * (2 - _f))
00023 , _e(sqrt(abs(_e2)))
00024 , _e12(_e2 / (1 - _e2))
00025 , _n(_f / (2 - _f))
00026 , _b(_a * _f1)
00027 , _tm(_a, _f, real(1))
00028 , _ell(-_e12)
00029 , _au(_a, _f, real(0), real(1), real(0), real(1), real(1))
00030 {}
00031
00032 const Ellipsoid& Ellipsoid::WGS84() {
00033 static const Ellipsoid wgs84(Constants::WGS84_a(), Constants::WGS84_f());
00034 return wgs84;
00035 }
00036
00037 Math::real Ellipsoid::QuarterMeridian() const
00038 { return _b * _ell.E(); }
00039
00040 Math::real Ellipsoid::Area() const {
00041 return 4 * Math::pi() *
00042 ((Math::sq(_a) + Math::sq(_b) *
00043 (_e2 == 0 ? 1 :
00044 (_e2 > 0 ? Math::atanh(sqrt(_e2)) : atan(sqrt(-_e2))) /
00045 sqrt(abs(_e2))))/2);
00046 }
00047
00048 Math::real Ellipsoid::ParametricLatitude(real phi) const
00049 { return atand(_f1 * tand(phi)); }
00050
00051 Math::real Ellipsoid::InverseParametricLatitude(real beta) const
00052 { return atand(tand(beta) / _f1); }
00053
00054 Math::real Ellipsoid::GeocentricLatitude(real phi) const
00055 { return atand(_f12 * tand(phi)); }
00056
00057 Math::real Ellipsoid::InverseGeocentricLatitude(real theta) const
00058 { return atand(tand(theta) / _f12); }
00059
00060 Math::real Ellipsoid::RectifyingLatitude(real phi) const {
00061 return abs(phi) == 90 ? phi:
00062 90 * MeridianDistance(phi) / QuarterMeridian();
00063 }
00064
00065 Math::real Ellipsoid::InverseRectifyingLatitude(real mu) const {
00066 if (abs(mu) == 90)
00067 return mu;
00068 return InverseParametricLatitude(_ell.Einv(mu * _ell.E() / 90) /
00069 Math::degree());
00070 }
00071
00072 Math::real Ellipsoid::AuthalicLatitude(real phi) const
00073 { return atand(_au.txif(tand(phi))); }
00074
00075 Math::real Ellipsoid::InverseAuthalicLatitude(real xi) const
00076 { return atand(_au.tphif(tand(xi))); }
00077
00078 Math::real Ellipsoid::ConformalLatitude(real phi) const
00079 { return atand(_tm.taupf(tand(phi))); }
00080
00081 Math::real Ellipsoid::InverseConformalLatitude(real chi) const
00082 { return atand(_tm.tauf(tand(chi))); }
00083
00084 Math::real Ellipsoid::IsometricLatitude(real phi) const
00085 { return Math::asinh(_tm.taupf(tand(phi))) / Math::degree(); }
00086
00087 Math::real Ellipsoid::InverseIsometricLatitude(real psi) const
00088 { return atand(_tm.tauf(sinh(psi * Math::degree()))); }
00089
00090 Math::real Ellipsoid::CircleRadius(real phi) const {
00091 return abs(phi) == 90 ? 0 :
00092
00093 _a / Math::hypot(real(1), _f1 * tand(phi));
00094 }
00095
00096 Math::real Ellipsoid::CircleHeight(real phi) const {
00097 real tbeta = _f1 * tand(phi);
00098
00099 return _b * tbeta / Math::hypot(real(1), _f1 * tand(phi));
00100 }
00101
00102 Math::real Ellipsoid::MeridianDistance(real phi) const
00103 { return _b * _ell.Ed( ParametricLatitude(phi) ); }
00104
00105 Math::real Ellipsoid::MeridionalCurvatureRadius(real phi) const {
00106 real v = 1 - _e2 * Math::sq(sin(phi * Math::degree()));
00107 return _a * (1 - _e2) / (v * sqrt(v));
00108 }
00109
00110 Math::real Ellipsoid::TransverseCurvatureRadius(real phi) const {
00111 real v = 1 - _e2 * Math::sq(sin(phi * Math::degree()));
00112 return _a / sqrt(v);
00113 }
00114
00115 Math::real Ellipsoid::NormalCurvatureRadius(real phi, real azi)
00116 const {
00117 real
00118 alpha = azi * Math::degree(),
00119 v = 1 - _e2 * Math::sq(sin(phi * Math::degree()));
00120 return _a / (sqrt(v) *
00121 (Math::sq(cos(alpha)) * v / (1 - _e2) + Math::sq(sin(alpha))));
00122 }
00123
00124 }