GeographicLib  1.42
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 (2014-2015) <charles@karney.com> and licensed
7  * under 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  Rhumb::Rhumb(real a, real f, bool exact)
19  : _ell(a, f)
20  , _exact(exact)
21  , _c2(_ell.Area() / 720)
22  {
23  real n = _ell._n, nx = n;
24  switch (maxpow_) {
25  case 4:
26  _R[1] = nx*(n*(n*(1772*n-5340)+6930)-4725)/14175;
27  nx *= n;
28  _R[2] = nx*((1590-1747*n)*n-630)/4725;
29  nx *= n;
30  _R[3] = nx*(104*n-31)/315;
31  nx *= n;
32  _R[4] = -41*nx/420;
33  break;
34  case 5:
35  _R[1] = nx*(n*(n*(n*(41662*n+58476)-176220)+228690)-155925)/467775;
36  nx *= n;
37  _R[2] = nx*(n*(n*(18118*n-57651)+52470)-20790)/155925;
38  nx *= n;
39  _R[3] = nx*((17160-23011*n)*n-5115)/51975;
40  nx *= n;
41  _R[4] = nx*(5480*n-1353)/13860;
42  nx *= n;
43  _R[5] = -668*nx/5775;
44  break;
45  case 6:
46  _R[1] = nx*(n*(n*(n*((56868630-114456994*n)*n+79819740)-240540300)+
47  312161850)-212837625)/638512875;
48  nx *= n;
49  _R[2] = nx*(n*(n*(n*(51304574*n+24731070)-78693615)+71621550)-28378350)/
50  212837625;
51  nx *= n;
52  _R[3] = nx*(n*(n*(1554472*n-6282003)+4684680)-1396395)/14189175;
53  nx *= n;
54  _R[4] = nx*((3205800-4913956*n)*n-791505)/8108100;
55  nx *= n;
56  _R[5] = nx*(1092376*n-234468)/2027025;
57  nx *= n;
58  _R[6] = -313076*nx/2027025;
59  break;
60  case 7:
61  _R[1] = nx*(n*(n*(n*(n*(n*(258618446*n-343370982)+170605890)+239459220)-
62  721620900)+936485550)-638512875)/1915538625;
63  nx *= n;
64  _R[2] = nx*(n*(n*(n*((153913722-248174686*n)*n+74193210)-236080845)+
65  214864650)-85135050)/638512875;
66  nx *= n;
67  _R[3] = nx*(n*(n*(n*(114450437*n+23317080)-94230045)+70270200)-20945925)/
68  212837625;
69  nx *= n;
70  _R[4] = nx*(n*(n*(15445736*n-103193076)+67321800)-16621605)/170270100;
71  nx *= n;
72  _R[5] = nx*((16385640-27766753*n)*n-3517020)/30405375;
73  nx *= n;
74  _R[6] = nx*(4892722*n-939228)/6081075;
75  nx *= n;
76  _R[7] = -3189007*nx/14189175;
77  break;
78  case 8:
79  _R[1] = nx*(n*(n*(n*(n*(n*((65947703730LL-13691187484LL*n)*n-
80  87559600410LL)+43504501950LL)+61062101100LL)-
81  184013329500LL)+238803815250LL)-162820783125LL)/
82  488462349375LL;
83  nx *= n;
84  _R[2] = nx*(n*(n*(n*(n*(n*(30802104839LL*n-63284544930LL)+39247999110LL)+
85  18919268550LL)-60200615475LL)+54790485750LL)-
86  21709437750LL)/162820783125LL;
87  nx *= n;
88  _R[3] = nx*(n*(n*(n*((5836972287LL-8934064508LL*n)*n+1189171080)-
89  4805732295LL)+3583780200LL)-1068242175)/10854718875LL;
90  nx *= n;
91  _R[4] = nx*(n*(n*(n*(50072287748LL*n+3938662680LL)-26314234380LL)+
92  17167059000LL)-4238509275LL)/43418875500LL;
93  nx *= n;
94  _R[5] = nx*(n*(n*(359094172*n-9912730821LL)+5849673480LL)-1255576140)/
95  10854718875LL;
96  nx *= n;
97  _R[6] = nx*((8733508770LL-16053944387LL*n)*n-1676521980)/10854718875LL;
98  nx *= n;
99  _R[7] = nx*(930092876*n-162639357)/723647925;
100  nx *= n;
101  _R[8] = -673429061*nx/1929727800;
102  break;
103  default:
104  GEOGRAPHICLIB_STATIC_ASSERT(maxpow_ >= 4 && maxpow_ <= 8,
105  "Bad value of maxpow_");
106  }
107  }
108 
109  const Rhumb& Rhumb::WGS84() {
110  static const Rhumb wgs84(Constants::WGS84_a(), Constants::WGS84_f(), false);
111  return wgs84;
112  }
113 
114  void Rhumb::GenInverse(real lat1, real lon1, real lat2, real lon2,
115  unsigned outmask,
116  real& s12, real& azi12, real& S12) const {
117  real
119  psi1 = _ell.IsometricLatitude(lat1),
120  psi2 = _ell.IsometricLatitude(lat2),
121  psi12 = psi2 - psi1,
122  h = Math::hypot(lon12, psi12);
123  if (outmask & AZIMUTH)
124  azi12 = Math::atan2d(lon12, psi12);
125  if (outmask & DISTANCE) {
126  real dmudpsi = DIsometricToRectifying(psi2, psi1);
127  s12 = h * dmudpsi * _ell.QuarterMeridian() / 90;
128  }
129  if (outmask & AREA)
130  S12 = _c2 * lon12 *
131  MeanSinXi(psi2 * Math::degree(), psi1 * Math::degree());
132  }
133 
134  RhumbLine Rhumb::Line(real lat1, real lon1, real azi12) const
135  { return RhumbLine(*this, lat1, lon1, azi12, _exact); }
136 
137  void Rhumb::GenDirect(real lat1, real lon1, real azi12, real s12,
138  unsigned outmask,
139  real& lat2, real& lon2, real& S12) const
140  { Line(lat1, lon1, azi12).GenPosition(s12, outmask, lat2, lon2, S12); }
141 
142  Math::real Rhumb::DE(real x, real y) const {
143  const EllipticFunction& ei = _ell._ell;
144  real d = x - y;
145  if (x * y <= 0)
146  return d ? (ei.E(x) - ei.E(y)) / d : 1;
147  // See DLMF: Eqs (19.11.2) and (19.11.4) letting
148  // theta -> x, phi -> -y, psi -> z
149  //
150  // (E(x) - E(y)) / d = E(z)/d - k2 * sin(x) * sin(y) * sin(z)/d
151  //
152  // tan(z/2) = (sin(x)*Delta(y) - sin(y)*Delta(x)) / (cos(x) + cos(y))
153  // = d * Dsin(x,y) * (sin(x) + sin(y))/(cos(x) + cos(y)) /
154  // (sin(x)*Delta(y) + sin(y)*Delta(x))
155  // = t = d * Dt
156  // sin(z) = 2*t/(1+t^2); cos(z) = (1-t^2)/(1+t^2)
157  // Alt (this only works for |z| <= pi/2 -- however, this conditions holds
158  // if x*y > 0):
159  // sin(z) = d * Dsin(x,y) * (sin(x) + sin(y))/
160  // (sin(x)*cos(y)*Delta(y) + sin(y)*cos(x)*Delta(x))
161  // cos(z) = sqrt((1-sin(z))*(1+sin(z)))
162  real sx = sin(x), sy = sin(y), cx = cos(x), cy = cos(y);
163  real Dt = Dsin(x, y) * (sx + sy) /
164  ((cx + cy) * (sx * ei.Delta(sy, cy) + sy * ei.Delta(sx, cx))),
165  t = d * Dt, Dsz = 2 * Dt / (1 + t*t),
166  sz = d * Dsz, cz = (1 - t) * (1 + t) / (1 + t*t);
167  return ((sz ? ei.E(sz, cz, ei.Delta(sz, cz)) / sz : 1)
168  - ei.k2() * sx * sy) * Dsz;
169  }
170 
171  Math::real Rhumb::DRectifying(real latx, real laty) const {
172  real
173  tbetx = _ell._f1 * Math::tand(latx),
174  tbety = _ell._f1 * Math::tand(laty);
175  return (Math::pi()/2) * _ell._b * _ell._f1 * DE(atan(tbetx), atan(tbety))
176  * Dtan(latx, laty) * Datan(tbetx, tbety) / _ell.QuarterMeridian();
177  }
178 
179  Math::real Rhumb::DIsometric(real latx, real laty) const {
180  real
181  phix = latx * Math::degree(), tx = Math::tand(latx),
182  phiy = laty * Math::degree(), ty = Math::tand(laty);
183  return Dasinh(tx, ty) * Dtan(latx, laty)
184  - Deatanhe(sin(phix), sin(phiy)) * Dsin(phix, phiy);
185  }
186 
187  Math::real Rhumb::SinCosSeries(bool sinp,
188  real x, real y, const real c[], int n) {
189  // N.B. n >= 0 and c[] has n+1 elements 0..n, of which c[0] is ignored.
190  //
191  // Use Clenshaw summation to evaluate
192  // m = (g(x) + g(y)) / 2 -- mean value
193  // s = (g(x) - g(y)) / (x - y) -- average slope
194  // where
195  // g(x) = sum(c[j]*SC(2*j*x), j = 1..n)
196  // SC = sinp ? sin : cos
197  // CS = sinp ? cos : sin
198  //
199  // This function returns only s; m is discarded.
200  //
201  // Write
202  // t = [m; s]
203  // t = sum(c[j] * f[j](x,y), j = 1..n)
204  // where
205  // f[j](x,y) = [ (SC(2*j*x)+SC(2*j*y))/2 ]
206  // [ (SC(2*j*x)-SC(2*j*y))/d ]
207  //
208  // = [ cos(j*d)*SC(j*p) ]
209  // [ +/-(2/d)*sin(j*d)*CS(j*p) ]
210  // (+/- = sinp ? + : -) and
211  // p = x+y, d = x-y
212  //
213  // f[j+1](x,y) = A * f[j](x,y) - f[j-1](x,y)
214  //
215  // A = [ 2*cos(p)*cos(d) -sin(p)*sin(d)*d]
216  // [ -4*sin(p)*sin(d)/d 2*cos(p)*cos(d) ]
217  //
218  // Let b[n+1] = b[n+2] = [0 0; 0 0]
219  // b[j] = A * b[j+1] - b[j+2] + c[j] * I for j = n..1
220  // t = (c[0] * I - b[2]) * f[0](x,y) + b[1] * f[1](x,y)
221  // c[0] is not accessed for s = t[2]
222  real p = x + y, d = x - y,
223  cp = cos(p), cd = cos(d),
224  sp = sin(p), sd = d ? sin(d)/d : 1,
225  m = 2 * cp * cd, s = sp * sd;
226  // 2x2 matrices stored in row-major order
227  const real a[4] = {m, -s * d * d, -4 * s, m};
228  real ba[4] = {0, 0, 0, 0};
229  real bb[4] = {0, 0, 0, 0};
230  real* b1 = ba;
231  real* b2 = bb;
232  if (n > 0) b1[0] = b1[3] = c[n];
233  for (int j = n - 1; j > 0; --j) { // j = n-1 .. 1
234  std::swap(b1, b2);
235  // b1 = A * b2 - b1 + c[j] * I
236  b1[0] = a[0] * b2[0] + a[1] * b2[2] - b1[0] + c[j];
237  b1[1] = a[0] * b2[1] + a[1] * b2[3] - b1[1];
238  b1[2] = a[2] * b2[0] + a[3] * b2[2] - b1[2];
239  b1[3] = a[2] * b2[1] + a[3] * b2[3] - b1[3] + c[j];
240  }
241  // Here are the full expressions for m and s
242  // m = (c[0] - b2[0]) * f01 - b2[1] * f02 + b1[0] * f11 + b1[1] * f12;
243  // s = - b2[2] * f01 + (c[0] - b2[3]) * f02 + b1[2] * f11 + b1[3] * f12;
244  if (sinp) {
245  // real f01 = 0, f02 = 0;
246  real f11 = cd * sp, f12 = 2 * sd * cp;
247  // m = b1[0] * f11 + b1[1] * f12;
248  s = b1[2] * f11 + b1[3] * f12;
249  } else {
250  // real f01 = 1, f02 = 0;
251  real f11 = cd * cp, f12 = - 2 * sd * sp;
252  // m = c[0] - b2[0] + b1[0] * f11 + b1[1] * f12;
253  s = - b2[2] + b1[2] * f11 + b1[3] * f12;
254  }
255  return s;
256  }
257 
258  Math::real Rhumb::DConformalToRectifying(real chix, real chiy) const {
259  return 1 + SinCosSeries(true, chix, chiy,
260  _ell.ConformalToRectifyingCoeffs(), tm_maxord);
261  }
262 
263  Math::real Rhumb::DRectifyingToConformal(real mux, real muy) const {
264  return 1 - SinCosSeries(true, mux, muy,
265  _ell.RectifyingToConformalCoeffs(), tm_maxord);
266  }
267 
268  Math::real Rhumb::DIsometricToRectifying(real psix, real psiy) const {
269  if (_exact) {
270  real
271  latx = _ell.InverseIsometricLatitude(psix),
272  laty = _ell.InverseIsometricLatitude(psiy);
273  return DRectifying(latx, laty) / DIsometric(latx, laty);
274  } else {
275  psix *= Math::degree();
276  psiy *= Math::degree();
277  return DConformalToRectifying(gd(psix), gd(psiy)) * Dgd(psix, psiy);
278  }
279  }
280 
281  Math::real Rhumb::DRectifyingToIsometric(real mux, real muy) const {
282  real
283  latx = _ell.InverseRectifyingLatitude(mux/Math::degree()),
284  laty = _ell.InverseRectifyingLatitude(muy/Math::degree());
285  return _exact ?
286  DIsometric(latx, laty) / DRectifying(latx, laty) :
287  Dgdinv(Math::taupf(Math::tand(latx), _ell._es),
288  Math::taupf(Math::tand(laty), _ell._es)) *
289  DRectifyingToConformal(mux, muy);
290  }
291 
292  Math::real Rhumb::MeanSinXi(real psix, real psiy) const {
293  return Dlog(cosh(psix), cosh(psiy)) * Dcosh(psix, psiy)
294  + SinCosSeries(false, gd(psix), gd(psiy), _R, maxpow_) * Dgd(psix, psiy);
295  }
296 
297  RhumbLine::RhumbLine(const Rhumb& rh, real lat1, real lon1, real azi12,
298  bool exact)
299  : _rh(rh)
300  , _exact(exact)
301  , _lat1(lat1)
302  , _lon1(lon1)
303  , _azi12(Math::AngNormalize(azi12))
304  {
305  real alp12 = _azi12 * Math::degree();
306  _salp = _azi12 == -180 ? 0 : sin(alp12);
307  _calp = abs(_azi12) == 90 ? 0 : cos(alp12);
308  _mu1 = _rh._ell.RectifyingLatitude(lat1);
309  _psi1 = _rh._ell.IsometricLatitude(lat1);
310  _r1 = _rh._ell.CircleRadius(lat1);
311  }
312 
313  void RhumbLine::GenPosition(real s12, unsigned outmask,
314  real& lat2, real& lon2, real& S12) const {
315  real
316  mu12 = s12 * _calp * 90 / _rh._ell.QuarterMeridian(),
317  mu2 = _mu1 + mu12;
318  real psi2, lat2x, lon2x;
319  if (abs(mu2) <= 90) {
320  if (_calp) {
321  lat2x = _rh._ell.InverseRectifyingLatitude(mu2);
322  real psi12 = _rh.DRectifyingToIsometric( mu2 * Math::degree(),
323  _mu1 * Math::degree()) * mu12;
324  lon2x = _salp * psi12 / _calp;
325  psi2 = _psi1 + psi12;
326  } else {
327  lat2x = _lat1;
328  lon2x = _salp * s12 / (_r1 * Math::degree());
329  psi2 = _psi1;
330  }
331  if (outmask & AREA)
332  S12 = _rh._c2 * lon2x *
333  _rh.MeanSinXi(_psi1 * Math::degree(), psi2 * Math::degree());
334  lon2x = outmask & LONG_NOWRAP ? _lon1 + lon2x :
335  Math::AngNormalize2(Math::AngNormalize(_lon1) + lon2x);
336  } else {
337  // Reduce to the interval [-180, 180)
338  mu2 = Math::AngNormalize2(mu2);
339  // Deal with points on the anti-meridian
340  if (abs(mu2) > 90) mu2 = Math::AngNormalize(180 - mu2);
341  lat2x = _rh._ell.InverseRectifyingLatitude(mu2);
342  lon2x = Math::NaN();
343  if (outmask & AREA)
344  S12 = Math::NaN();
345  }
346  if (outmask & LATITUDE) lat2 = lat2x;
347  if (outmask & LONGITUDE) lon2 = lon2x;
348  }
349 
350 } // namespace GeographicLib
static T AngNormalize(T x)
Definition: Math.hpp:428
static T NaN()
Definition: Math.hpp:611
Math::real InverseRectifyingLatitude(real mu) const
Definition: Ellipsoid.cpp:65
Math::real IsometricLatitude(real phi) const
Definition: Ellipsoid.cpp:84
static T pi()
Definition: Math.hpp:214
GeographicLib::Math::real real
Definition: GeodSolve.cpp:32
Rhumb(real a, real f, bool exact=true)
Definition: Rhumb.cpp:18
void GenPosition(real s12, unsigned outmask, real &lat2, real &lon2, real &S12) const
Definition: Rhumb.cpp:313
Header for GeographicLib::Rhumb and GeographicLib::RhumbLine classes.
Elliptic integrals and functions.
Math::real Delta(real sn, real cn) const
friend class RhumbLine
Definition: Rhumb.hpp:69
static T hypot(T x, T y)
Definition: Math.hpp:255
static const Rhumb & WGS84()
Definition: Rhumb.cpp:109
RhumbLine Line(real lat1, real lon1, real azi12) const
Definition: Rhumb.cpp:134
static T atan2d(T y, T x)
Definition: Math.hpp:534
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
Math::real QuarterMeridian() const
Definition: Ellipsoid.cpp:37
static T degree()
Definition: Math.hpp:228
static T AngDiff(T x, T y)
Definition: Math.hpp:458
Math::real InverseIsometricLatitude(real psi) const
Definition: Ellipsoid.cpp:87
static T tand(T x)
Definition: Math.hpp:500
Solve of the direct and inverse rhumb problems.
Definition: Rhumb.hpp:66
Find a sequence of points on a single rhumb line.
Definition: Rhumb.hpp:442
static T AngNormalize2(T x)
Definition: Math.hpp:440