1

2# -*- coding: utf-8 -*-

3

4 u'''A pure Python version of I{Karney}'s I{Auxiliary Latitudes}, C++ classes U{Rhumb

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

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

7 I{GeographicLib version 2.2+} renamed to L{RhumbAux} respectively L{RhumbLineAux}.

8

9 Class L{RhumbLineAux} has been enhanced with methods C{Intersecant2}, C{Intersection} and C{PlumbTo}

10 to iteratively find the intersection of a rhumb line and a circle or an other rhumb line, respectively

11 a perpendicular geodesic or other rhumb line.

12

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

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

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

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

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

18

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

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

21

22@note: C{S12} area calculations in classes L{RhumbAux} and L{RhumbLineAux} depend on class L{AuxDST} which

23 requires U{numpy<https://PyPI.org/project/numpy>} to be installed, version 1.16 or newer.

24

25@note: Windows reserves file names U{AUX, COM[#], CON, LPT[#], NUL, PRN<https://learn.Microsoft.com/en-us/

26 windows/win32/fileio/naming-a-file#naming-conventions>} with and without extension.

27'''

28# make sure int/int division yields float quotient

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

30

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

32 from pygeodesy.auxilats.auxDLat import AuxDLat, _DClenshaw

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

34 from pygeodesy.auxilats.auxily import _Dlam, _Dp0Dpsi

35 from pygeodesy.auxilats._CX_Rs import _Rdict, _Rtuple

36 from pygeodesy.basics import copysign0, _reverange, _xkwds_get1

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

38 _0_0, _4_0, _720_0, _log2, _over

39# from pygeodesy.datums import _WGS84 # from .rhumb.bases

40# from pygeodesy.errors import _xkwds_get1 # from .basics

41 from pygeodesy.karney import Caps, _polynomial

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

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

44# from pygeodesy.props import Property_RO # from .rhumb.bases

45 from pygeodesy.rhumb.bases import RhumbBase, RhumbLineBase, \

46 Property_RO, _WGS84

47

48 from math import ceil as _ceil, fabs, radians

49

50__all__ = _ALL_LAZY.rhumb_aux_

51__version__ = '24.09.02'

52

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

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

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

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

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

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

59#

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

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

62

63

64 class RhumbAux(RhumbBase):

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

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

67

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

69 installed, version 1.16 or later.

70 '''

71

72 def __init__(self, a_earth=_WGS84, f=None, exact=True, **TMorder_name): # PYCHOK signature

73 '''New C{RhumbAux}.

74

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

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

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

78 @kwarg f: The ellipsoid's flattening (C{scalar}), required if B{C{a_earth}}

79 is C{scalar}, ignored otherwise.

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

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

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

83 @kwarg TMorder_name: Optional C{B{name}=NN} (C{str}) and optional

84 keyword argument C{B{TMorder}=6}, the order of the

85 L{KTransverseMercator}, see property C{TMorder}.

86

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

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

89

90 @raise RhumbError: Invalid B{C{a_earth}}, B{C{f}} or B{C{TMorder}}.

91 '''

92 RhumbBase.__init__(self, a_earth, f, exact, TMorder_name)

93

94 def areaux(self, **exact):

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

96

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

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

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

100

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

102

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

104 property C{f_max}.

105

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

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

108

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

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

111 '''

112 x = _xkwds_get1(exact, exact=self.exact)

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

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

115 return a

116

117 @Property_RO

118 def _auxD(self):

119 return AuxDLat(self.ellipsoid)

120

121 @Property_RO

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

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

124

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

126 xD = self._auxD

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

128 xD.CRectifying(Chi1, Chi2)

129 if r:

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

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

132 return r

133

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

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

136 '''

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

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

139 psi12 = psi2 - psi1 # radians

140 lam12 = radians(lon12)

141 if (outmask & Caps.DISTANCE):

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

143 s = fabs(Phi2.toMu(self).toRadians -

144 Phi1.toMu(self).toRadians)

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

146 s = hypot(lam12, psi12)

147 if s:

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

149 s *= self._rrm

150 a = _over(s, self._mpd)

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

152 return lam12, psi12, Chi1, Chi2

153

154 def _latPhi2(self, mu):

155 Mu = AuxMu.fromDegrees(mu)

156 Phi = Mu.toPhi(self)

157 return Phi.toDegrees, Phi

158

159 @Property_RO

160 def _mpd(self): # meter per degree

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

162

163 def _psiChiPhi3(self, lat):

164 Phi = AuxPhi.fromDegrees(lat)

165 Chi = Phi.toChi(self)

166 psi = Chi.toLambertianRadians

167 return psi, Chi, Phi

168

169 @Property_RO

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

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

172 _RAseries(self._auxD))

173

174# _RhumbLine = RhumbLineAux # see further below

175

176 @Property_RO

177 def _rrm(self):

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

179

180 _mpr = _rrm # meter per radian, see _mpd

181

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

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

184 '''

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

186

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

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

189

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

191 xD.CParametric(Chix, Chiy)

192 if dD:

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

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

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

196 Phiy.toBeta(self).normalized,

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

198 dD += _Dp0Dpsi(tx, ty)

199 dD *= self._c2 * lon12

200 return float(dD)

201

202

203 class RhumbLineAux(RhumbLineBase):

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

205

206 Class C{RhumbLineAux} facilitates the determination of points

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

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

209 '''

210 _Rhumb = RhumbAux # rhumb.aux_.RhumbAux

211

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

213 '''New C{RhumbLineAux}.

214

215 @arg rhumb: The rhumb reference (L{RhumbAux}).

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

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

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

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

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

221 values specifying the required capabilities. Include

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

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

224 '''

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

226

227 @Property_RO

228 def _Chi1(self):

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

230

231 @Property_RO

232 def _mu1(self):

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

234 '''

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

236

237 def _mu2lat(self, mu):

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

239 '''

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

241 return lat

242

243 @Property_RO

244 def _Phi1(self):

245 return AuxPhi.fromDegrees(self.lat1)

246

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

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

249 '''

250 R = self.rhumb

251 lat2, Phi2 = R._latPhi2(mu2)

252 Chi2 = Phi2.toChi(R)

253 Chi1 = self._Chi1

254 lon2 = self._salp * a12

255 if lon2:

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

257 lon2 = _over(lon2, m)

258 return lat2, lon2, Chi1, Chi2

259

260# @Property_RO

261# def _psi1(self):

262# return self._Chi1.toLambertianRadians

263

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

265

266

267 def _RAintegrate(auxD):

268 # Compute coefficients by Fourier transform of integrand

269 L = 2

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

271 f = auxD._qIntegrand

272 c_ = fft.transform(f)

273 pP = []

274 _P = pP.append

275 # assert L < _Lbits

276 while L < _Lbits:

277 L = fft.reset(L) * 2

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

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

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

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

282 # Compute Fourier coefficients of integral

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

284 if fabs(p) > EPS_2:

285 k = -1 # run interrupted

286 else:

287 if k < 0:

288 k = 1 # mark as first small value

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

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

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

292 _P(-p)

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

294

295

296 def _RAseries(auxD):

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

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

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

300 d = n = auxD._n

301 i = 0

302 aL = auxD.ALorder

303 Cs = _RACoeffs[aL]

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

305 pP = []

306 _p = _polynomial

307 for m in _reverange(aL): # order

308 j = i + m + 1

309 pP.append(_p(n, Cs, i, j) * d)

310 d *= n

311 i = j

312 # assert i == len(pP)

313 return pP

314

315

316_RACoeffs = _Rdict(110, # Rhumb Area Coefficients in matrix Q

317 _Rtuple(4, 10, # GEOGRAPHICLIB_RHUMBAREA_ORDER == 4

318 '596/2025, -398/945, 22/45, -1/3',

319 '1543/4725, -118/315, 1/5',

320 '152/945, -17/315',

321 '5/252'),

322 _Rtuple(5, 15, # GEOGRAPHICLIB_RHUMBAREA_ORDER == 5

323 '-102614/467775, 596/2025, -398/945, 22/45, -1/3',

324 '-24562/155925, 1543/4725, -118/315, 1/5',

325 '-38068/155925, 152/945, -17/315',

326 '-752/10395, 5/252',

327 '-101/17325'),

328 _Rtuple(6, 21, # GEOGRAPHICLIB_RHUMBAREA_ORDER == 6

329 '138734126/638512875, -102614/467775, 596/2025, -398/945, 22/45, -1/3',

330 '17749373/425675250, -24562/155925, 1543/4725, -118/315, 1/5',

331 '1882432/8513505, -38068/155925, 152/945, -17/315',

332 '268864/2027025, -752/10395, 5/252',

333 '62464/2027025, -101/17325',

334 '11537/4054050'),

335 _Rtuple(7, 28, # GEOGRAPHICLIB_RHUMBAREA_ORDER == 7

336 '-565017322/1915538625, 138734126/638512875, -102614/467775, 596/2025, -398/945, 22/45, -1/3',

337 '-1969276/58046625, 17749373/425675250, -24562/155925, 1543/4725, -118/315, 1/5',

338 '-58573784/638512875, 1882432/8513505, -38068/155925, 152/945, -17/315',

339 '-6975184/42567525, 268864/2027025, -752/10395, 5/252',

340 '-112832/1447875, 62464/2027025, -101/17325',

341 '-4096/289575, 11537/4054050',

342 '-311/525525'),

343 _Rtuple(8, 36, # GEOGRAPHICLIB_RHUMBAREA_ORDER == 8

344 '188270561816/488462349375, -565017322/1915538625, 138734126/638512875, -102614/467775, 596/2025, -398/945, 22/45, -1/3',

345 '2332829602/23260111875, -1969276/58046625, 17749373/425675250, -24562/155925, 1543/4725, -118/315, 1/5',

346 '-41570288/930404475, -58573784/638512875, 1882432/8513505, -38068/155925, 152/945, -17/315',

347 '1538774036/10854718875, -6975184/42567525, 268864/2027025, -752/10395, 5/252',

348 '436821248/3618239625, -112832/1447875, 62464/2027025, -101/17325',

349 '3059776/80405325, -4096/289575, 11537/4054050',

350 '4193792/723647925, -311/525525',

351 '1097653/1929727800')

352)

353 del _Rdict, _Rtuple

354

355__all__ += _ALL_DOCS(Caps)

356

357# **) MIT License

358#

359# Copyright (C) 2023-2024 -- mrJean1 at Gmail -- All Rights Reserved.

360#

361# Permission is hereby granted, free of charge, to any person obtaining a

362# copy of this software and associated documentation files (the "Software"),

363# to deal in the Software without restriction, including without limitation

364# the rights to use, copy, modify, merge, publish, distribute, sublicense,

365# and/or sell copies of the Software, and to permit persons to whom the

366# Software is furnished to do so, subject to the following conditions:

367#

368# The above copyright notice and this permission notice shall be included

369# in all copies or substantial portions of the Software.

370#

371# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

372# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

373# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

374# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

375# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

376# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

377# OTHER DEALINGS IN THE SOFTWARE.