GeographicLib  1.38
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Rhumb.cpp
Go to the documentation of this file.
1 /**
2  * \file Rhumb.cpp
3  * \brief Implementation for GeographicLib::Rhumb and GeographicLib::RhumbLine
4  * classes
5  *
6  * Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
7  * the MIT/X11 License. For more information, see
8  * http://geographiclib.sourceforge.net/
9  **********************************************************************/
10 
11 #include <algorithm>
12 #include <GeographicLib/Rhumb.hpp>
13 
14 namespace GeographicLib {
15 
16  using namespace std;
17 
18  const Rhumb& Rhumb::WGS84() {
19  static const Rhumb wgs84(Constants::WGS84_a(), Constants::WGS84_f(), false);
20  return wgs84;
21  }
22 
23  void Rhumb::Inverse(real lat1, real lon1, real lat2, real lon2,
24  real& s12, real& azi12) const {
25  real
27  psi1 = _ell.IsometricLatitude(lat1),
28  psi2 = _ell.IsometricLatitude(lat2),
29  psi12 = psi2 - psi1,
30  h = Math::hypot(lon12, psi12);
31  azi12 = 0 - atan2(-lon12, psi12) / Math::degree();
32  real dmudpsi = DIsometricToRectifying(psi2 * Math::degree(),
33  psi1 * Math::degree());
34  s12 = h * dmudpsi * _ell.QuarterMeridian() / 90;
35  }
36 
37  RhumbLine Rhumb::Line(real lat1, real lon1, real azi12) const
38  { return RhumbLine(*this, lat1, lon1, azi12, _exact); }
39 
40  void Rhumb::Direct(real lat1, real lon1, real azi12, real s12,
41  real& lat2, real& lon2) const
42  { Line(lat1, lon1, azi12).Position(s12, lat2, lon2); }
43 
44  Math::real Rhumb::DE(real x, real y) const {
45  const EllipticFunction& ei = _ell._ell;
46  real d = x - y;
47  if (x * y <= 0)
48  return d ? (ei.E(x) - ei.E(y)) / d : 1;
49  // See DLMF: Eqs (19.11.2) and (19.11.4) letting
50  // theta -> x, phi -> -y, psi -> z
51  //
52  // (E(x) - E(y)) / d = E(z)/d - k2 * sin(x) * sin(y) * sin(z)/d
53  //
54  // tan(z/2) = (sin(x)*Delta(y) - sin(y)*Delta(x)) / (cos(x) + cos(y))
55  // = d * Dsin(x,y) * (sin(x) + sin(y))/(cos(x) + cos(y)) /
56  // (sin(x)*Delta(y) + sin(y)*Delta(x))
57  // = t = d * Dt
58  // sin(z) = 2*t/(1+t^2); cos(z) = (1-t^2)/(1+t^2)
59  // Alt (this only works for |z| <= pi/2 -- however, this conditions holds
60  // if x*y > 0):
61  // sin(z) = d * Dsin(x,y) * (sin(x) + sin(y))/
62  // (sin(x)*cos(y)*Delta(y) + sin(y)*cos(x)*Delta(x))
63  // cos(z) = sqrt((1-sin(z))*(1+sin(z)))
64  real sx = sin(x), sy = sin(y), cx = cos(x), cy = cos(y);
65  real Dt = Dsin(x, y) * (sx + sy) /
66  ((cx + cy) * (sx * ei.Delta(sy, cy) + sy * ei.Delta(sx, cx))),
67  t = d * Dt, Dsz = 2 * Dt / (1 + t*t),
68  sz = d * Dsz, cz = (1 - t) * (1 + t) / (1 + t*t);
69  return ((sz ? ei.E(sz, cz, ei.Delta(sz, cz)) / sz : 1)
70  - ei.k2() * sx * sy) * Dsz;
71  }
72 
73  Math::real Rhumb::DRectifying(real latx, real laty) const {
74  real
75  phix = latx * Math::degree(), tbetx = _ell._f1 * tano(phix),
76  phiy = laty * Math::degree(), tbety = _ell._f1 * tano(phiy);
77  return (Math::pi()/2) * _ell._b * _ell._f1 * DE(atan(tbetx), atan(tbety))
78  * Dtan(phix, phiy) * Datan(tbetx, tbety) / _ell.QuarterMeridian();
79  }
80 
81  Math::real Rhumb::DIsometric(real latx, real laty) const {
82  real
83  phix = latx * Math::degree(), tx = tano(phix),
84  phiy = laty * Math::degree(), ty = tano(phiy);
85  return Dasinh(tx, ty) * Dtan(phix, phiy)
86  - Deatanhe(sin(phix), sin(phiy)) * Dsin(phix, phiy);
87  }
88 
89  Math::real Rhumb::SinSeries(real x, real y, const real c[], int n) {
90  // N.B. n >= 0 and c[] has n+1 elements 0..n, of which c[0] is ignored.
91  //
92  // Use Clenshaw summation to evaluate
93  // m = (g(x) + g(y)) / 2 -- mean value
94  // s = (g(x) - g(y)) / (x - y) -- average slope
95  // where
96  // g(x) = sum(c[j]*sin(2*j*x), j = 1..n)
97  //
98  // This function returns only s and m is discarded.
99  //
100  // Write
101  // t = [m; s]
102  // t = sum(c[j] * f[j](x,y), j = 1..n)
103  // where
104  // f[j](x,y) = [ (sin(2*j*x)+sin(2*j*y))/2 ]
105  // [ (sin(2*j*x)-sin(2*j*y))/d ]
106  //
107  // = [ sin(j*p)*cos(j*d) ]
108  // [ (2/d)*sin(j*d)*cos(j*p) ]
109  // and
110  // p = x+y, d = x-y
111  //
112  // f[j+1](x,y) = A * f[j](x,y) - f[j-1](x,y)
113  //
114  // A = [ 2*cos(p)*cos(d) -sin(p)*sin(d)*d]
115  // [ -4*sin(p)*sin(d)/d 2*cos(p)*cos(d) ]
116  //
117  // Let b[n+1] = b[n+2] = [0 0; 0 0]
118  // b[j] = A * b[j+1] - b[j+2] + c[j] * I for j = n..1
119  // t = b[1] * f[1](x,y)
120  real p = x + y, d = x - y,
121  cp = cos(p), cd = cos(d),
122  sp = sin(p), sd = d ? sin(d)/d : 1,
123  m = 2 * cp * cd, s = sp * sd;
124  // 2x2 matrices stored in row-major order
125  const real a[4] = {m, -s * d * d, -4 * s, m};
126  real ba[4] = {0, 0, 0, 0};
127  real bb[4] = {0, 0, 0, 0};
128  real* b0 = ba;
129  real* b1 = bb;
130  if (n > 0) b0[0] = b0[3] = c[n];
131  for (int j = n - 1; j > 0; --j) { // j = n-1 .. 1
132  std::swap(b0, b1);
133  // b0 = A * b1 - b0 + c[j] * I
134  b0[0] = a[0] * b1[0] + a[1] * b1[2] - b0[0] + c[j];
135  b0[1] = a[0] * b1[1] + a[1] * b1[3] - b0[1];
136  b0[2] = a[2] * b1[0] + a[3] * b1[2] - b0[2];
137  b0[3] = a[2] * b1[1] + a[3] * b1[3] - b0[3] + c[j];
138  }
139  b1[0] = sp * cd; b1[2] = 2 * sd * cp;
140  // Here is the (unused) expression for m
141  // m = b0[0] * b1[0] + b0[1] * b1[2];
142  s = b0[2] * b1[0] + b0[3] * b1[2];
143  return s;
144  }
145 
146  Math::real Rhumb::DConformalToRectifying(real chix, real chiy) const {
147  return 1 + SinSeries(chix, chiy,
148  _ell.ConformalToRectifyingCoeffs(), tm_maxord);
149  }
150 
151  Math::real Rhumb::DRectifyingToConformal(real mux, real muy) const {
152  return 1 - SinSeries(mux, muy,
153  _ell.RectifyingToConformalCoeffs(), tm_maxord);
154  }
155 
156  Math::real Rhumb::DIsometricToRectifying(real psix, real psiy) const {
157  if (_exact) {
158  real
159  latx = _ell.InverseIsometricLatitude(psix/Math::degree()),
160  laty = _ell.InverseIsometricLatitude(psiy/Math::degree());
161  return DRectifying(latx, laty) / DIsometric(latx, laty);
162  } else
163  return DConformalToRectifying(gd(psix), gd(psiy)) * Dgd(psix, psiy);
164  }
165 
166  Math::real Rhumb::DRectifyingToIsometric(real mux, real muy) const {
167  real
168  latx = _ell.InverseRectifyingLatitude(mux/Math::degree()),
169  laty = _ell.InverseRectifyingLatitude(muy/Math::degree());
170  return _exact ?
171  DIsometric(latx, laty) / DRectifying(latx, laty) :
172  Dgdinv(_ell.ConformalLatitude(latx) * Math::degree(),
173  _ell.ConformalLatitude(laty) * Math::degree()) *
174  DRectifyingToConformal(mux, muy);
175  }
176 
177  RhumbLine::RhumbLine(const Rhumb& rh, real lat1, real lon1, real azi12,
178  bool exact)
179  : _rh(rh)
180  , _exact(exact)
181  , _lat1(lat1)
182  , _lon1(Math::AngNormalize(lon1))
183  , _azi12(Math::AngNormalize(azi12))
184  {
185  real alp12 = azi12 * Math::degree();
186  _salp = azi12 == -180 ? 0 : sin(alp12);
187  _calp = abs(azi12) == 90 ? 0 : cos(alp12);
188  _mu1 = _rh._ell.RectifyingLatitude(lat1);
189  _psi1 = _rh._ell.IsometricLatitude(lat1);
190  _r1 = _rh._ell.CircleRadius(lat1);
191  }
192 
193  void RhumbLine::Position(real s12, real& lat2, real& lon2) const {
194  real
195  mu12 = s12 * _calp * 90 / _rh._ell.QuarterMeridian(),
196  mu2 = _mu1 + mu12;
197  if (abs(mu2) <= 90) {
198  if (_calp) {
199  lat2 = _rh._ell.InverseRectifyingLatitude(mu2);
200  real psi12 = _rh.DRectifyingToIsometric( mu2 * Math::degree(),
201  _mu1 * Math::degree()) * mu12;
202  lon2 = _salp * psi12 / _calp;
203  } else {
204  lat2 = _lat1;
205  lon2 = _salp * s12 / (_r1 * Math::degree());
206  }
207  lon2 = Math::AngNormalize2(_lon1 + lon2);
208  } else {
209  // Reduce to the interval [-180, 180)
210  mu2 = Math::AngNormalize2(mu2);
211  // Deal with points on the anti-meridian
212  if (abs(mu2) > 90) mu2 = Math::AngNormalize(180 - mu2);
213  lat2 = _rh._ell.InverseRectifyingLatitude(mu2);
214  lon2 = Math::NaN();
215  }
216  }
217 
218 } // namespace GeographicLib
static T AngNormalize(T x)
Definition: Math.hpp:399
static T NaN()
Definition: Math.hpp:460
Math::real InverseRectifyingLatitude(real mu) const
Definition: Ellipsoid.cpp:65
static T pi()
Definition: Math.hpp:213
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
void Position(real s12, real &lat2, real &lon2) const
Definition: Rhumb.cpp:193
Header for GeographicLib::Rhumb and GeographicLib::RhumbLine classes.
Elliptic integrals and functions.
Math::real Delta(real sn, real cn) const
static T hypot(T x, T y)
Definition: Math.hpp:254
static const Rhumb & WGS84()
Definition: Rhumb.cpp:18
RhumbLine Line(real lat1, real lon1, real azi12) const
Definition: Rhumb.cpp:37
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
Math::real QuarterMeridian() const
Definition: Ellipsoid.cpp:37
static T degree()
Definition: Math.hpp:227
static T AngDiff(T x, T y)
Definition: Math.hpp:429
Solve of the direct and inverse rhumb problems.
Definition: Rhumb.hpp:48
void Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi12) const
Definition: Rhumb.cpp:23
Find a sequence of points on a single rhumb line.
Definition: Rhumb.hpp:275
void Direct(real lat1, real lon1, real azi12, real s12, real &lat2, real &lon2) const
Definition: Rhumb.cpp:40
static T AngNormalize2(T x)
Definition: Math.hpp:411