1# -*- coding: utf-8 -*-

2

3 u'''A pure Python version of I{Karney}'s C++ classes U{Rhumb

4<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Rhumb.html>} and U{RhumbLine

5<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1RhumbLine.html>} from

6 I{GeographicLib version 2.2+} renamed as L{RhumbAux} respectively L{RhumbLineAux}.

7

8 Class L{RhumbLineAux} has been enhanced with methods C{intersection2} and C{nearestOn4} to iteratively

9 find the intersection of two rhumb lines, respectively the nearest point on a rumb line along a

10 geodesic or perpendicular rhumb line from an other point.

11

12 For more details, see the I{2.2} U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/index.html>}

13 documentation, especially the U{Class List<https://GeographicLib.SourceForge.io/C++/doc/annotated.html>},

14 the background information on U{Rhumb lines<https://GeographicLib.SourceForge.io/C++/doc/rhumb.html>},

15 utility U{RhumbSolve<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html>} and U{Online rhumb

16 line calculations<https://GeographicLib.SourceForge.io/cgi-bin/RhumbSolve>}.

17

18 Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2022-2023) and licensed under the MIT/X11

19 License. For more information, see the U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.

20

21@note: Class L{AuxDST} requires package U{numpy<https://PyPI.org/project/numpy>} to be installed,

22 version 1.16 or newer, needed only for I{exact} area calculations C{S12} in L{RhumbAux}

23 and L{RhumbLineAux}.

24'''

25# make sure int/int division yields float quotient

26 from __future__ import division as _; del _ # PYCHOK semicolon

27

28 from pygeodesy.auxilats.auxAngle import AuxMu, AuxPhi, hypot

29 from pygeodesy.auxilats.auxDLat import AuxDLat, _DClenshaw

30# from pygeodesy.auxilats.auxDST import AuxDST # _MODS

31 from pygeodesy.auxilats.auxily import _Dlam, _Dp0Dpsi, _Ufloats

32 from pygeodesy.basics import copysign0, _reverange, _xkwds_get

33 from pygeodesy.constants import EPS_2, MANT_DIG, PI4, isinf, \

34 _0_0, _4_0, _720_0, _log2, _over

35# from pygeodesy.datums import _WGS84 # from .karney

36# from pygeodesy.errors import _xkwds_get # from .basics

37 from pygeodesy.karney import Caps, _polynomial, NN, _WGS84

38# from pygeodesy.fmath import hypot # from .auxilats.auxAngle

39# from pygeodesy.interns import NN # from .karney

40 from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS

41# from pygeodesy.props import Property_RO # from .rhumbBase

42 from pygeodesy.rhumbBase import RhumbBase, RhumbLineBase, Property_RO

43

44 from math import ceil as _ceil, fabs, radians

45

46__all__ = _ALL_LAZY.rhumbaux

47__version__ = '23.10.27'

48

49# DIGITS = (sizeof(real) * 8) bits

50# = (ctypes.sizeof(ctypes.c_double(1.0)) * 8) bits

51# For |n| <= 0.99, actual max for doubles is 2163. This scales

52# as DIGITS and for long doubles (GEOGRAPHICLIB_PRECISION = 3,

53# DIGITS = 64), this becomes 2163 * 64 / 53 = 2612. Round this

54# up to 2^12 = 4096 and scale this by DIGITS//64 if DIGITS > 64.

55#

56# 64 = DIGITS for long double, 6 = 12 - _log2(64)

57_Lbits = 1 << (int(_ceil(_log2(max(MANT_DIG, 64)))) + 6)

58

59

60 class RhumbAux(RhumbBase):

61 '''Class to solve the I{direct} and I{inverse rhumb} problems, based

62 on I{Auxiliary Latitudes} for accuracy near the poles.

63

64 @note: Package U{numpy<https://PyPI.org/project/numpy>} must be

65 installed, version 1.16 or later.

66 '''

67

68 def __init__(self, a_earth=_WGS84, f=None, exact=True, name=NN, **TMorder): # PYCHOK signature

69 '''New C{rhumbaux.RhumbAux}.

70

71 @kwarg a_earth: This rhumb's earth model (L{Datum}, L{Ellipsoid},

72 L{Ellipsoid2}, L{a_f2Tuple}, 2-tuple C{(a, f)}) or

73 the (equatorial) radius (C{meter}, conventionally).

74 @kwarg f: The ellipsoid's flattening (C{scalar}), iff B{C{a_earth}} is

75 C{scalar}, ignored otherwise.

76 @kwarg exact: If C{True}, use the exact expressions for the I{Auxiliary

77 Latitudes}, otherwise use the I{Fourier} series expansion

78 (C{bool}), see also property C{exact}.

79 @kwarg name: Optional name (C{str}).

80 @kwarg TMorder: Optional keyword argument B{C{TMorder}}, see property

81 C{TMorder}.

82

83 @raise ImportError: Package C{numpy} not found or not installed, only

84 required for area C{S12} when C{B{exact} is True}.

85

86 @raise RhumbError: Invalid B{C{a_earth}}, B{C{f}} or B{C{RA_TMorder}}.

87 '''

88 RhumbBase.__init__(self, a_earth, f, exact, name)

89 if TMorder:

90 self.Tmorder = _xkwds_get(TMorder, TMorder=RhumbBase._mTM)

91

92 def areaux(self, **exact):

93 '''Get this ellipsoid's B{C{exact}} surface area (C{meter} I{squared}).

94

95 @kwarg exact: Optional C{exact} (C{bool}), overriding this rhumb's

96 C{exact} setting, if C{True}, use the exact expression

97 for the authalic radius otherwise the I{Taylor} series.

98

99 @return: The (signed?) surface area (C{meter} I{squared}).

100

101 @raise AuxError: If C{B{exact}=False} and C{abs(flattening)} exceeds

102 property C{f_max}.

103

104 @note: The area of a polygon encircling a pole can be found by adding

105 C{areaux / 2} to the sum of C{S12} for each side of the polygon.

106

107 @see: U{The area of rhumb polygons<https://ArXiv.org/pdf/2303.03219.pdf>}

108 and method L{auxilats.AuxLat.AuthalicRadius2}.

109 '''

110 x = _xkwds_get(exact, exact=self.exact)

111 a = (self._c2 * _720_0) if bool(x) is self.exact else \

112 (self._auxD.AuthalicRadius2(exact=x, f_max=self.f_max) * PI4)

113 return a

114

115 @Property_RO

116 def _auxD(self):

117 return AuxDLat(self.ellipsoid)

118

119 @Property_RO

120 def _c2(self): # radians makes _c2 a factor per degree

121 return radians(self._auxD.AuthalicRadius2(exact=self.exact, f_max=self.f_max))

122

123 def _DMu_DPsi(self, Phi1, Phi2, Chi1, Chi2):

124 xD = self._auxD

125 r = xD.DRectifying(Phi1, Phi2) if self.exact else \

126 xD.CRectifying(Chi1, Chi2)

127 if r:

128 r = _over(r, xD.DIsometric(Phi1, Phi2) if self.exact else

129 _Dlam(Chi1.tan, Chi2.tan)) # not Lambertian!

130 return r

131

132 def _Inverse4(self, lon12, r, outmask):

133 '''(INTERNAL) See method C{RhumbBase.Inverse}.

134 '''

135 psi1, Chi1, Phi1 = self._psiChiPhi3(r.lat1)

136 psi2, Chi2, Phi2 = self._psiChiPhi3(r.lat2)

137 psi12 = psi2 - psi1 # radians

138 lam12 = radians(lon12)

139 if (outmask & Caps.DISTANCE):

140 if isinf(psi1) or isinf(psi2): # PYCHOK no cover

141 d = Phi2.toMu(self).toRadians

142 d -= Phi1.toMu(self).toRadians

143 s = fabs(d)

144 else: # dmu/dpsi = dmu/dchi/dpsi/dchi

145 s = hypot(lam12, psi12)

146 if s:

147 s *= self._DMu_DPsi(Phi1, Phi2, Chi1, Chi2)

148 s *= self._rrm

149 a = _over(s, self._mpd)

150 r.set_(a12=copysign0(a, s), s12=s)

151 return lam12, psi12, Chi1, Chi2

152

153 def _latPhi2(self, mu):

154 Mu = AuxMu.fromDegrees(mu)

155 Phi = Mu.toPhi(self)

156 return Phi.toDegrees, Phi

157

158 @Property_RO

159 def _mpd(self): # meter per degree

160 return radians(self._rrm) # == self.ellipsoid._Lpd

161

162 def _psiChiPhi3(self, lat):

163 Phi = AuxPhi.fromDegrees(lat)

164 Chi = Phi.toChi(self)

165 psi = Chi.toLambertianRadians

166 return psi, Chi, Phi

167

168 @Property_RO

169 def _RA(self): # get the coefficients for area calculation

170 return tuple(_RAintegrate(self._auxD) if self.exact else

171 _RAseries(self._auxD))

172

173# _RhumbLine = RhumbLineAux # see further below

174

175 @Property_RO

176 def _rrm(self):

177 return self._auxD.RectifyingRadius(exact=self.exact)

178

179 _mpr = _rrm # meter per radian, see _mpd

180

181 def _S12d(self, Chix, Chiy, lon12): # degrees

182 '''(INTERNAL) Compute the area C{S12} from C{._meanSinXi(Chix, Chiy) * .c2 * lon12}.

183 '''

184 pP, xD = self._RA, self._auxD

185

186 tx, Phix = Chix.tan, Chix.toPhi(self)

187 ty, Phiy = Chiy.tan, Chiy.toPhi(self)

188

189 dD = xD.DParametric(Phix, Phiy) if self.exact else \

190 xD.CParametric(Chix, Chiy)

191 if dD:

192 dD = _over(dD, xD.DIsometric(Phix, Phiy) if self.exact else

193 _Dlam(tx, ty)) # not Lambertian!

194 dD *= _DClenshaw(False, Phix.toBeta(self).normalized,

195 Phiy.toBeta(self).normalized,

196 pP, min(len(pP), xD.ALorder)) # Fsum

197 dD += _Dp0Dpsi(tx, ty)

198 dD *= self._c2 * lon12

199 return float(dD)

200

201

202 class RhumbLineAux(RhumbLineBase):

203 '''Compute one or several points on a single rhumb line.

204

205 Class C{RhumbLineAux} facilitates the determination of points

206 on a single rhumb line. The starting point (C{lat1}, C{lon1})

207 and the azimuth C{azi12} are specified once.

208

209 Method C{RhumbLineAux.Position} returns the location of an

210 other point and optionally the distance C{s12} along and the

211 area C{S12} under the rhumb line.

212

213 Method C{RhumbLineAux.intersection2} finds the intersection

214 between two rhumb lines.

215

216 Method C{RhumbLineAux.nearestOn4} computes the nearest point

217 on and the distance to a rhumb line in different ways.

218 '''

219 _Rhumb = RhumbAux # rhumbaux.RhumbAux

220

221 def __init__(self, rhumb, lat1=0, lon1=0, azi12=None, **caps_name): # PYCHOK signature

222 '''New C{rhumbaux.RhumbLineAux}.

223

224 @arg rhumb: The rhumb reference (C{rhumbaux.RhumbAux}).

225 @kwarg lat1: Latitude of the start point (C{degrees90}).

226 @kwarg lon1: Longitude of the start point (C{degrees180}).

227 @kwarg azi12: Azimuth of this rhumb line (compass C{degrees}).

228 @kwarg caps_name: Optional keyword arguments C{B{name}=NN} and

229 C{B{caps}=0}, a bit-or'ed combination of L{Caps}

230 values specifying the required capabilities. Include

231 C{Caps.LINE_OFF} if updates to the B{C{rhumb}} should

232 I{not} be reflected in this rhumb line.

233 '''

234 RhumbLineBase.__init__(self, rhumb, lat1, lon1, azi12, **caps_name)

235

236 @Property_RO

237 def _Chi1(self):

238 return self._Phi1.toChi(self.rhumb)

239

240 @Property_RO

241 def _mu1(self):

242 '''(INTERNAL) Get the I{rectifying auxiliary} latitude (C{degrees}).

243 '''

244 return self._Phi1.toMu(self.rhumb).toDegrees

245

246 def _mu2lat(self, mu):

247 '''(INTERNAL) Get the inverse I{rectifying auxiliary} latitude (C{degrees}).

248 '''

249 lat, _ = self.rhumb._latPhi2(mu)

250 return lat

251

252 @Property_RO

253 def _Phi1(self):

254 return AuxPhi.fromDegrees(self.lat1)

255

256 def _Position4(self, a12, mu2, *unused): # PYCHOK s12, mu2

257 '''(INTERNAL) See method C{RhumbLineBase._Position}.

258 '''

259 R = self.rhumb

260 lat2, Phi2 = R._latPhi2(mu2)

261 Chi2 = Phi2.toChi(R)

262 Chi1 = self._Chi1

263 lon2 = self._salp * a12

264 if lon2:

265 m = R._DMu_DPsi(self._Phi1, Phi2, Chi1, Chi2)

266 lon2 = _over(lon2, m)

267 return lat2, lon2, Chi1, Chi2

268

269# @Property_RO

270# def _psi1(self):

271# return self._Chi1.toLambertianRadians

272

273 RhumbAux._RhumbLine = RhumbLineAux # PYCHOK see RhumbBase._RhumbLine

274

275

276 def _RAintegrate(auxD):

277 # Compute coefficients by Fourier transform of integrand

278 L = 2

279 fft = _MODS.auxilats.auxDST.AuxDST(L)

280 f = auxD._qIntegrand

281 c_ = fft.transform(f)

282 pP = []

283 _P = pP.append

284 # assert L < _Lbits

285 while L < _Lbits:

286 L = fft.reset(L) * 2

287 c_ = fft.refine(f, c_, _0_0) # sentine[L]

288 # assert len(c_) == L + 1

289 pP[:], k = [], -1

290 for j in range(1, L + 1):

291 # Compute Fourier coefficients of integral

292 p = (c_[j - 1] + c_[j]) / (_4_0 * j)

293 if fabs(p) > EPS_2:

294 k = -1 # run interrupted

295 else:

296 if k < 0:

297 k = 1 # mark as first small value

298 if (j - k) >= ((j + 7) // 8):

299 # run of at least (j - 1) // 8 small values

300 return pP[:j] # break while L loop

301 _P(-p)

302 return pP # no convergence, use pP as-is

303

304

305 def _RAseries(auxD):

306 # Series expansions in n for Fourier coeffients of the integral

307 # @see: U{"Series expansions for computing rhumb areas"

308 # <https:#DOI.org/10.5281/zenodo.7685484>}.

309 d = n = auxD._n

310 i = 0

311 aL = auxD.ALorder

312 Cs = _RACoeffs[aL]

313 # assert len(Cs) == (aL * (aL + 1)) // 2

314 pP = []

315 _P = pP.append

316 _p = _polynomial

317 for m in _reverange(aL): # order

318 j = i + m + 1

319 _P(_p(n, Cs, i, j) * d)

320 d *= n

321 i = j

322 # assert i == len(pP)

323 return pP

324

325

326_f, _u = float, _Ufloats()

327_RACoeffs = { # Rhumb Area Coefficients in matrix Q

328 4: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 4

329 596 / _f(2025), -398 / _f(945), 22 / _f(45), -1 / _f(3),

330 1543 / _f(4725), -118 / _f(315), 1 / _f(5),

331 152 / _f(945), -17 / _f(315),

332 5 / _f(252)),

333 5: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 5

334 -102614 / _f(467775), 596 / _f(2025), -398 / _f(945), 22 / _f(45),

335 -1 / _f(3),

336 -24562 / _f(155925), 1543 / _f(4725), -118 / _f(315), 1 / _f(5),

337 -38068 / _f(155925), 152 / _f(945), -17 / _f(315),

338 -752 / _f(10395), 5 / _f(252),

339 -101 / _f(17325)),

340 6: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 6

341 138734126 / _f(638512875), -102614 / _f(467775), 596 / _f(2025),

342 -398 / _f(945), 22 / _f(45), -1 / _f(3),

343 17749373 / _f(425675250), -24562 / _f(155925), 1543 / _f(4725),

344 -118 / _f(315), 1 / _f(5),

345 1882432 / _f(8513505), -38068 / _f(155925), 152 / _f(945),

346 -17 / _f(315),

347 268864 / _f(2027025), -752 / _f(10395), 5 / _f(252),

348 62464 / _f(2027025), -101 / _f(17325),

349 11537 / _f(4054050)),

350 7: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 7

351 -565017322 / _f(1915538625), 138734126 / _f(638512875),

352 -102614 / _f(467775), 596 / _f(2025), -398 / _f(945), 22 / _f(45),

353 -1 / _f(3),

354 -1969276 / _f(58046625), 17749373 / _f(425675250), -24562 / _f(155925),

355 1543 / _f(4725), -118 / _f(315), 1 / _f(5),

356 -58573784 / _f(638512875), 1882432 / _f(8513505), -38068 / _f(155925),

357 152 / _f(945), -17 / _f(315),

358 -6975184 / _f(42567525), 268864 / _f(2027025), -752 / _f(10395),

359 5 / _f(252),

360 -112832 / _f(1447875), 62464 / _f(2027025), -101 / _f(17325),

361 -4096 / _f(289575), 11537 / _f(4054050),

362 -311 / _f(525525)),

363 8: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 8

364 188270561816 / _f(488462349375), -565017322 / _f(1915538625),

365 138734126 / _f(638512875), -102614 / _f(467775), 596 / _f(2025),

366 -398 / _f(945), 22 / _f(45), -1 / _f(3),

367 2332829602 / _f(23260111875), -1969276 / _f(58046625),

368 17749373 / _f(425675250), -24562 / _f(155925), 1543 / _f(4725),

369 -118 / _f(315), 1 / _f(5),

370 -41570288 / _f(930404475), -58573784 / _f(638512875),

371 1882432 / _f(8513505), -38068 / _f(155925), 152 / _f(945),

372 -17 / _f(315),

373 1538774036 / _f(10854718875), -6975184 / _f(42567525),

374 268864 / _f(2027025), -752 / _f(10395), 5 / _f(252),

375 436821248 / _f(3618239625), -112832 / _f(1447875),

376 62464 / _f(2027025), -101 / _f(17325),

377 3059776 / _f(80405325), -4096 / _f(289575), 11537 / _f(4054050),

378 4193792 / _f(723647925), -311 / _f(525525),

379 1097653 / _f(1929727800))

380}

381 del _f, _u, _Ufloats

382

383__all__ += _ALL_DOCS(Caps, RhumbAux, RhumbLineAux)

384

385 if __name__ == '__main__':

386

387 from pygeodesy import printf, RhumbAux # PYCHOK RhumbAux is __main__.RhumbAux

388 from pygeodesy.basics import _zip

389

390 def _re(fmt, r3, x3):

391 e3 = []

392 for r, x in _zip(r3, x3): # strict=True

393 e = fabs(r - x) / fabs(x)

394 e3.append('%.g' % (e,))

395 printf((fmt % r3) + ' rel errors: ' + ', '.join(e3))

396

397 # <https://GeographicLib.SourceForge.io/cgi-bin/RhumbSolve -p 9> version 2.2

398 rhumb = RhumbAux(exact=True) # WGS84 default

399 printf('# %r\n', rhumb)

400 r = rhumb.Direct8(40.6, -73.8, 51, 5.5e6) # from JFK about NE

401 _re('# JFK NE lat2=%.8f, lon2=%.8f, S12=%.1f', (r.lat2, r.lon2, r.S12), (71.688899882813, 0.2555198244234, 44095641862956.11))

402 r = rhumb.Inverse8(40.6, -73.8, 51.6, -0.5) # JFK to LHR

403 _re('# JFK-LHR azi12=%.8f, s12=%.3f S12=%.1f', (r.azi12, r.s12, r.S12), (77.7683897102557, 5771083.38332803, 37395209100030.39))

404 r = rhumb.Inverse8(40.6, -73.8, 35.8, 140.3) # JFK to Tokyo Narita

405 _re('# JFK-NRT azi12=%.8f, s12=%.3f S12=%.1f', (r.azi12, r.s12, r.S12), (-92.38888798169965, 12782581.067684170, -63760642939072.50))

406

407# % python3 -m pygeodesy.rhumbaux

408

409# RhumbAux(TMorder=6, ellipsoid=Ellipsoid(name='WGS84', a=6378137, b=6356752.31424518, f_=298.25722356, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367449.14582341, L=10001965.72931272, R1=6371008.77141506, R2=6371007.18091847, R3=6371000.79000916, Rbiaxial=6367453.63451633, Rtriaxial=6372797.5559594), exact=True)

410

411# JFK NE lat2=71.68889988, lon2=0.25551982, S12=44095641862956.1 rel errors: 4e-11, 2e-08, 5e-16

412# JFK-LHR azi12=77.76838971, s12=5771083.383 S12=37395209100030.3 rel errors: 3e-12, 5e-15, 6e-16

413# JFK-NRT azi12=-92.38888798, s12=12782581.068 S12=-63760642939072.5 rel errors: 2e-16, 3e-16, 6e-16