00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <algorithm>
00012 #include <GeographicLib/Rhumb.hpp>
00013
00014 namespace GeographicLib {
00015
00016 using namespace std;
00017
00018 const Rhumb& Rhumb::WGS84() {
00019 static const Rhumb wgs84(Constants::WGS84_a(), Constants::WGS84_f(), false);
00020 return wgs84;
00021 }
00022
00023 void Rhumb::Inverse(real lat1, real lon1, real lat2, real lon2,
00024 real& s12, real& azi12) const {
00025 real
00026 lon12 = Math::AngDiff(Math::AngNormalize(lon1), Math::AngNormalize(lon2)),
00027 psi1 = _ell.IsometricLatitude(lat1),
00028 psi2 = _ell.IsometricLatitude(lat2),
00029 psi12 = psi2 - psi1,
00030 h = Math::hypot(lon12, psi12);
00031 azi12 = 0 - atan2(-lon12, psi12) / Math::degree();
00032 real dmudpsi = DIsometricToRectifying(psi2 * Math::degree(),
00033 psi1 * Math::degree());
00034 s12 = h * dmudpsi * _ell.QuarterMeridian() / 90;
00035 }
00036
00037 RhumbLine Rhumb::Line(real lat1, real lon1, real azi12) const
00038 { return RhumbLine(*this, lat1, lon1, azi12, _exact); }
00039
00040 void Rhumb::Direct(real lat1, real lon1, real azi12, real s12,
00041 real& lat2, real& lon2) const
00042 { Line(lat1, lon1, azi12).Position(s12, lat2, lon2); }
00043
00044 Math::real Rhumb::DE(real x, real y) const {
00045 const EllipticFunction& ei = _ell._ell;
00046 real d = x - y;
00047 if (x * y <= 0)
00048 return d ? (ei.E(x) - ei.E(y)) / d : 1;
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064 real sx = sin(x), sy = sin(y), cx = cos(x), cy = cos(y);
00065 real Dt = Dsin(x, y) * (sx + sy) /
00066 ((cx + cy) * (sx * ei.Delta(sy, cy) + sy * ei.Delta(sx, cx))),
00067 t = d * Dt, Dsz = 2 * Dt / (1 + t*t),
00068 sz = d * Dsz, cz = (1 - t) * (1 + t) / (1 + t*t);
00069 return ((sz ? ei.E(sz, cz, ei.Delta(sz, cz)) / sz : 1)
00070 - ei.k2() * sx * sy) * Dsz;
00071 }
00072
00073 Math::real Rhumb::DRectifying(real latx, real laty) const {
00074 real
00075 phix = latx * Math::degree(), tbetx = _ell._f1 * tano(phix),
00076 phiy = laty * Math::degree(), tbety = _ell._f1 * tano(phiy);
00077 return (Math::pi()/2) * _ell._b * _ell._f1 * DE(atan(tbetx), atan(tbety))
00078 * Dtan(phix, phiy) * Datan(tbetx, tbety) / _ell.QuarterMeridian();
00079 }
00080
00081 Math::real Rhumb::DIsometric(real latx, real laty) const {
00082 real
00083 phix = latx * Math::degree(), tx = tano(phix),
00084 phiy = laty * Math::degree(), ty = tano(phiy);
00085 return Dasinh(tx, ty) * Dtan(phix, phiy)
00086 - Deatanhe(sin(phix), sin(phiy)) * Dsin(phix, phiy);
00087 }
00088
00089 Math::real Rhumb::SinSeries(real x, real y, const real c[], int n) {
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120 real p = x + y, d = x - y,
00121 cp = cos(p), cd = cos(d),
00122 sp = sin(p), sd = d ? sin(d)/d : 1,
00123 m = 2 * cp * cd, s = sp * sd;
00124
00125 const real a[4] = {m, -s * d * d, -4 * s, m};
00126 real ba[4] = {0, 0, 0, 0};
00127 real bb[4] = {0, 0, 0, 0};
00128 real* b0 = ba;
00129 real* b1 = bb;
00130 if (n > 0) b0[0] = b0[3] = c[n];
00131 for (int j = n - 1; j > 0; --j) {
00132 std::swap(b0, b1);
00133
00134 b0[0] = a[0] * b1[0] + a[1] * b1[2] - b0[0] + c[j];
00135 b0[1] = a[0] * b1[1] + a[1] * b1[3] - b0[1];
00136 b0[2] = a[2] * b1[0] + a[3] * b1[2] - b0[2];
00137 b0[3] = a[2] * b1[1] + a[3] * b1[3] - b0[3] + c[j];
00138 }
00139 b1[0] = sp * cd; b1[2] = 2 * sd * cp;
00140
00141
00142 s = b0[2] * b1[0] + b0[3] * b1[2];
00143 return s;
00144 }
00145
00146 Math::real Rhumb::DConformalToRectifying(real chix, real chiy) const {
00147 return 1 + SinSeries(chix, chiy,
00148 _ell.ConformalToRectifyingCoeffs(), tm_maxord);
00149 }
00150
00151 Math::real Rhumb::DRectifyingToConformal(real mux, real muy) const {
00152 return 1 - SinSeries(mux, muy,
00153 _ell.RectifyingToConformalCoeffs(), tm_maxord);
00154 }
00155
00156 Math::real Rhumb::DIsometricToRectifying(real psix, real psiy) const {
00157 if (_exact) {
00158 real
00159 latx = _ell.InverseIsometricLatitude(psix/Math::degree()),
00160 laty = _ell.InverseIsometricLatitude(psiy/Math::degree());
00161 return DRectifying(latx, laty) / DIsometric(latx, laty);
00162 } else
00163 return DConformalToRectifying(gd(psix), gd(psiy)) * Dgd(psix, psiy);
00164 }
00165
00166 Math::real Rhumb::DRectifyingToIsometric(real mux, real muy) const {
00167 real
00168 latx = _ell.InverseRectifyingLatitude(mux/Math::degree()),
00169 laty = _ell.InverseRectifyingLatitude(muy/Math::degree());
00170 return _exact ?
00171 DIsometric(latx, laty) / DRectifying(latx, laty) :
00172 Dgdinv(_ell.ConformalLatitude(latx) * Math::degree(),
00173 _ell.ConformalLatitude(laty) * Math::degree()) *
00174 DRectifyingToConformal(mux, muy);
00175 }
00176
00177 RhumbLine::RhumbLine(const Rhumb& rh, real lat1, real lon1, real azi12,
00178 bool exact)
00179 : _rh(rh)
00180 , _exact(exact)
00181 , _lat1(lat1)
00182 , _lon1(Math::AngNormalize(lon1))
00183 , _azi12(Math::AngNormalize(azi12))
00184 {
00185 real alp12 = azi12 * Math::degree();
00186 _salp = azi12 == -180 ? 0 : sin(alp12);
00187 _calp = abs(azi12) == 90 ? 0 : cos(alp12);
00188 _mu1 = _rh._ell.RectifyingLatitude(lat1);
00189 _psi1 = _rh._ell.IsometricLatitude(lat1);
00190 _r1 = _rh._ell.CircleRadius(lat1);
00191 }
00192
00193 void RhumbLine::Position(real s12, real& lat2, real& lon2) const {
00194 real
00195 mu12 = s12 * _calp * 90 / _rh._ell.QuarterMeridian(),
00196 mu2 = _mu1 + mu12;
00197 if (abs(mu2) <= 90) {
00198 if (_calp) {
00199 lat2 = _rh._ell.InverseRectifyingLatitude(mu2);
00200 real psi12 = _rh.DRectifyingToIsometric( mu2 * Math::degree(),
00201 _mu1 * Math::degree()) * mu12;
00202 lon2 = _salp * psi12 / _calp;
00203 } else {
00204 lat2 = _lat1;
00205 lon2 = _salp * s12 / (_r1 * Math::degree());
00206 }
00207 lon2 = Math::AngNormalize2(_lon1 + lon2);
00208 } else {
00209
00210 mu2 = Math::AngNormalize2(mu2);
00211
00212 if (abs(mu2) > 90) mu2 = Math::AngNormalize(180 - mu2);
00213 lat2 = _rh._ell.InverseRectifyingLatitude(mu2);
00214 lon2 = Math::NaN();
00215 }
00216 }
00217
00218 }