Coverage for pygeodesy/rhumb/ekx.py: 98%

231 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-06-10 14:08 -0400

1 

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

3 

4u'''A pure Python version of I{Karney}'s I{elliptic functions}, I{Krüger series expansion}, C++ 

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

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

7from I{GeographicLib version 2.0}, kept for backward compatibility. 

8 

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

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

11a perpendicular geodesic or other rhumb line. 

12 

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

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

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

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

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

18 

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

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

21''' 

22# make sure int/int division yields float quotient 

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

24 

25from pygeodesy.basics import copysign0, neg 

26from pygeodesy.constants import PI_2, _0_0s, _0_0, _0_5, _1_0, \ 

27 _2_0, _4_0, _720_0, _over, _1_over 

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

29# from pygeodesy.deprecated import RhumbOrder2Tuple # _MODS 

30from pygeodesy.errors import RhumbError, _xkwds_pop2, _Xorder 

31from pygeodesy.fmath import hypot, hypot1 

32# from pygeodesy.fsums import fsum1f_ # _MODS 

33# from pygeodesy.karney import Caps # from .rhumb.bases 

34from pygeodesy.ktm import KTransverseMercator, _Xs, \ 

35 _AlpCoeffs, _BetCoeffs # PYCHOK used! 

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

37from pygeodesy.props import deprecated_method, Property, Property_RO, \ 

38 property_RO 

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

40 Caps, _update_all_rls, _WGS84 

41from pygeodesy.utily import atan1, sincos2_ 

42 

43from math import asinh, atan, cos, cosh, radians, sin, sinh, sqrt, tan 

44 

45__all__ = _ALL_LAZY.rhumb_ekx 

46__version__ = '24.05.30' 

47 

48 

49class Rhumb(RhumbBase): 

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

51 I{elliptic functions} or I{Krüger series expansion} 

52 

53 @see: The U{Detailed Description<https://GeographicLib.SourceForge.io/C++/doc/ 

54 classGeographicLib_1_1Rhumb.html>} of I{Karney}'s C++ C{Rhumb Class}. 

55 ''' 

56 _mRA = 6 # see .RAorder 

57 

58 def __init__(self, a_earth=_WGS84, f=None, exact=True, **RA_TMorder_name): 

59 '''New C{Rhumb}. 

60 

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

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

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

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

65 C{scalar}, ignored otherwise. 

66 @kwarg exact: If C{True}, use an addition theorem for elliptic integrals 

67 to compute I{Divided differences}, otherwise use the I{Krüger} 

68 series expansion (C{bool} or C{None}), see also properties 

69 C{exact} and C{TMorder}. 

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

71 @kwarg RA_TMorder_name: Optional C{B{name}=NN} (C{str}) and optional keyword 

72 arguments B{C{RAorder}=6} and B{C{TMorder}=6} to set the respective 

73 C{order}, see properties C{RAorder} and C{TMorder}. 

74 

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

76 ''' 

77 if RA_TMorder_name: 

78 M = self._mRA 

79 m, kwds = _xkwds_pop2(RA_TMorder_name, RAorder=M) 

80 if m != M: 

81 self.RAorder = m 

82 else: 

83 kwds = {} 

84 RhumbBase.__init__(self, a_earth, f, exact, kwds) 

85 

86 @Property_RO 

87 def _A2(self): # Conformal2RectifyingCoeffs 

88 m = self.TMorder 

89 return _Xs(_AlpCoeffs, m, self.ellipsoid), m 

90 

91 @Property_RO 

92 def _B2(self): # Rectifying2ConformalCoeffs 

93 m = self.TMorder 

94 return _Xs(_BetCoeffs, m, self.ellipsoid), m 

95 

96 def _DConformal2Rectifying(self, x, y): # radians 

97 return _1_0 + (_sincosSeries(True, x, y, *self._A2) if self.f else _0_0) 

98 

99 @deprecated_method 

100 def Direct7(self, lat1, lon1, azi12, s12, outmask=Caps.LATITUDE_LONGITUDE_AREA): # PYCHOK no cover 

101 '''DEPRECATED, use method L{Rhumb.Direct8}. 

102 

103 @return: A I{DEPRECATED} L{Rhumb7Tuple}. 

104 ''' 

105 return self.Direct8(lat1, lon1, azi12, s12, outmask=outmask)._to7Tuple() 

106 

107 def _DIsometrict(self, phix, phiy, tphix, tphiy, _Dtan_phix_phiy): 

108 E = self.ellipsoid 

109 return _Dtan_phix_phiy * _Dasinh(tphix, tphiy) - \ 

110 _Dsin(phix, phiy) * _DeatanhE(sin(phix), sin(phiy), E) 

111 

112 def _DIsometric2Rectifyingd(self, psix, psiy): # degrees 

113 if self.exact: 

114 E = self.ellipsoid 

115 phix, phiy, tphix, tphiy = _Eaux4(E.auxIsometric, psix, psiy) 

116 t = _Dtant(phix - phiy, tphix, tphiy) 

117 r = _over(self._DRectifyingt( tphix, tphiy, t), 

118 self._DIsometrict(phix, phiy, tphix, tphiy, t)) 

119 else: 

120 x, y = radians(psix), radians(psiy) 

121 r = self._DConformal2Rectifying(_gd(x), _gd(y)) * _Dgd(x, y) 

122 return r 

123 

124 def _DRectifyingt(self, tphix, tphiy, _Dtan_phix_phiy): 

125 E = self.ellipsoid 

126 tbetx = E.f1 * tphix 

127 tbety = E.f1 * tphiy 

128 return (E.f1 * _Dtan_phix_phiy * E.b * PI_2 

129 * _DfEt( tbetx, tbety, self._eF) 

130 * _Datan(tbetx, tbety)) / E.L 

131 

132 def _DRectifying2Conformal(self, x, y): # radians 

133 return _1_0 - (_sincosSeries(True, x, y, *self._B2) if self.f else _0_0) 

134 

135 def _DRectifying2Isometricd(self, mux, muy): # degrees 

136 E = self.ellipsoid 

137 phix, phiy, tphix, tphiy = _Eaux4(E.auxRectifying, mux, muy) 

138 if self.exact: 

139 t = _Dtant(phix - phiy, tphix, tphiy) 

140 r = _over(self._DIsometrict(phix, phiy, tphix, tphiy, t), 

141 self._DRectifyingt( tphix, tphiy, t)) 

142 else: 

143 r = self._DRectifying2Conformal(radians(mux), radians(muy)) * \ 

144 _Dgdinv(E.es_taupf(tphix), E.es_taupf(tphiy)) 

145 return r 

146 

147 @Property_RO 

148 def _eF(self): 

149 '''(INTERNAL) Get the ellipsoid's elliptic function. 

150 ''' 

151 # .k2 = 0.006739496742276434 

152 return self.ellipsoid._elliptic_e12 # _MODS.elliptic.Elliptic(-self.ellipsoid._e12) 

153 

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

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

156 ''' 

157 E, Cs = self.ellipsoid, Caps 

158 psi1 = E.auxIsometric(r.lat1) 

159 psi2 = E.auxIsometric(r.lat2) 

160 psi12 = psi2 - psi1 # degrees 

161 if (outmask & Cs.DISTANCE): 

162 a = s = hypot(lon12, psi12) 

163 if a: 

164 a *= self._DIsometric2Rectifyingd(psi2, psi1) 

165 s = self._mpd * a # == E._Lpd 

166 a = copysign0(a, s) 

167 r.set_(a12=a, s12=s) 

168 

169 if ((outmask | self._debug) & Cs._DEBUG_INVERSE): # PYCHOK no cover 

170 r.set_(a=E.a, f=E.f, f1=E.f1, L=E.L, 

171 b=E.b, e=E.e, e2=E.e2, k2=self._eF.k2, 

172 lon12=lon12, psi1=psi1, exact=self.exact, 

173 psi12=psi12, psi2=psi2) 

174 return lon12, psi12, psi1, psi2 

175 

176 @deprecated_method 

177 def Inverse7(self, lat1, lon1, azi12, s12, outmask=Caps.AZIMUTH_DISTANCE_AREA): # PYCHOK no cover 

178 '''DEPRECATED, use method L{Rhumb.Inverse8}. 

179 

180 @return: A I{DEPRECATED} L{Rhumb7Tuple}. 

181 ''' 

182 return self.Inverse8(lat1, lon1, azi12, s12, outmask=outmask)._to7Tuple() 

183 

184 @Property_RO 

185 def _mpd(self): # meter per degree 

186 return self.ellipsoid._Lpd 

187 

188 @Property_RO 

189 def _mpr(self): # meter per radian 

190 return self.ellipsoid._Lpr # degrees(._Lpd) 

191 

192 @deprecated_method 

193 def orders(self, RAorder=6, TMorder=6): # PYCHOK no cover 

194 '''DEPRECATED, use properties C{RAorder} and/or C{TMorder}. 

195 

196 Get and set the I{RAorder} and/or I{TMorder}. 

197 

198 @kwarg RAorder: I{Rhumb Area} order (C{int}, 4, 5, 6, 7 

199 or 8). 

200 @kwarg TMorder: I{Transverse Mercator} order (C{int}, 4, 

201 5, 6, 7 or 8). 

202 

203 @return: DEPRECATED L{RhumbOrder2Tuple}C{(RAorder, TMorder)} 

204 with the previous C{RAorder} and C{TMorder} setting. 

205 ''' 

206 t = _MODS.deprecated.classes.RhumbOrder2Tuple(self.RAorder, self.TMorder) 

207 if RAorder != t.RAorder: # PYCHOK attr 

208 self.RAorder = RAorder 

209 if TMorder != t.TMorder: # PYCHOK attr 

210 self.TMorder = TMorder 

211 return t 

212 

213 @Property_RO 

214 def _RA2(self): 

215 # for WGS84: (0, -0.0005583633519275459, -3.743803759172812e-07, -4.633682270824446e-10, 

216 # RAorder 6: -7.709197397676237e-13, -1.5323287106694307e-15, -3.462875359099873e-18) 

217 m = self.RAorder 

218 return _Xs(_RACoeffs, m, self.ellipsoid, RA=True), m 

219 

220 @Property 

221 def RAorder(self): 

222 '''Get the I{Rhumb Area} order (C{int}, 4, 5, 6, 7 or 8). 

223 ''' 

224 return self._mRA 

225 

226 @RAorder.setter # PYCHOK setter! 

227 def RAorder(self, order): 

228 '''Set the I{Rhumb Area} order (C{int}, 4, 5, 6, 7 or 8). 

229 ''' 

230 m = _Xorder(_RACoeffs, RhumbError, RAorder=order) 

231 if self._mRA != m: 

232 _update_all_rls(self) 

233 self._mRA = m 

234 

235# _RhumbLine = RhumbLine # see further below 

236 

237 def _S12d(self, psi1, psi2, lon12): # degrees 

238 '''(INTERNAL) Compute the area C{S12}. 

239 ''' 

240 S = (self.ellipsoid.areax if self.exact else 

241 self.ellipsoid.area) * lon12 / _720_0 

242 if S: 

243 x, y = radians(psi1), radians(psi2) # _meanSinXi(x, y) 

244 s = _Dlog(cosh(x), cosh(y)) * _Dcosh(x, y) 

245 if self.f: 

246 s += _sincosSeries(False, _gd(x), _gd(y), *self._RA2) * _Dgd(x, y) 

247 S *= s 

248 return S 

249 

250 

251class RhumbLine(RhumbLineBase): 

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

253 

254 Class C{RhumbLine} facilitates the determination of points on 

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

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

257 ''' 

258 _Rhumb = Rhumb # rhumb.ekx.Rhumb 

259 

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

261 '''New C{RhumbLine}. 

262 

263 @arg rhumb: The rhumb reference (L{Rhumb}). 

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

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

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

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

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

269 values specifying the required capabilities. Include 

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

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

272 ''' 

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

274 

275 @Property_RO 

276 def _dpm12(self): # PYCHOK no cover 

277 '''(INTERNAL) Get this rhumb line's parallel I{circle radius} (C{degree per meter}). 

278 ''' 

279 r = self._salp 

280 if r: 

281 r = _over(r, self.ellipsoid.circle4(self.lat1).radius) 

282 return r 

283 

284 @Property_RO 

285 def _mu1(self): 

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

287 ''' 

288 return self.ellipsoid.auxRectifying(self.lat1) 

289 

290 def _mu2lat(self, mu): 

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

292 ''' 

293 return self.ellipsoid.auxRectifying(mu, inverse=True) 

294 

295 def _Position4(self, unused, mu2, s12, mu12): 

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

297 ''' 

298 psi1 = psi2 = self._psi1 

299 if mu12: # self._calp != 0 

300 lat2 = self._mu2lat(mu2) 

301 psi12 = self.rhumb._DRectifying2Isometricd(mu2, self._mu1) * mu12 

302 lon2 = self._talp * psi12 

303 psi2 += psi12 

304 else: # meridional 

305 lat2 = self.lat1 

306 lon2 = self._dpm12 * s12 

307 return lat2, lon2, psi1, psi2 

308 

309 @Property_RO 

310 def _psi1(self): 

311 '''(INTERNAL) Get the I{isometric auxiliary} latitude C{psi} (C{degrees}). 

312 ''' 

313 return self.ellipsoid.auxIsometric(self.lat1) 

314 

315 @property_RO 

316 def RAorder(self): 

317 '''Get this rhumb line's I{Rhumb Area} order (C{int}, 4, 5, 6, 7 or 8). 

318 ''' 

319 return self.rhumb.RAorder 

320 

321Rhumb._RhumbLine = RhumbLine # PYCHOK see RhumbBase._RhumbLine 

322 

323 

324# Use I{Divided Differences} to determine (mu2 - mu1) / (psi2 - psi1) accurately. 

325# Definition: _Df(x,y,d) = (f(x) - f(y)) / (x - y), @see W. M. Kahan & R. J. 

326# Fateman, "Symbolic computation of Divided Differences", SIGSAM Bull. 33(3), 

327# 7-28 (1999). U{ACM<https://DL.ACM.org/doi/pdf/10.1145/334714.334716> and @see 

328# U{UCB<https://www.CS.Berkeley.edu/~fateman/papers/divdiff.pdf>}, Dec 8, 1999. 

329 

330def _Dasinh(x, y): 

331 hx = hypot1(x) 

332 d = x - y 

333 if d: 

334 hx *= y 

335 hy = x * hypot1(y) 

336 t = (d * (x + y) / (hy + hx)) if (x * y) > 0 else (hy - hx) 

337 r = asinh(t) / d 

338 else: 

339 r = _1_0 / hx 

340 return r 

341 

342 

343def _Datan(x, y): 

344 xy = x * y 

345 r = xy + _1_0 

346 d = x - y 

347 if d: # 2 * xy > -1 == 2 * xy + 1 > 0 == xy + r > 0 == xy > -r 

348 r = (atan1(d, r) if xy > -r else (atan1(x) - atan1(y))) / d 

349 else: 

350 r = _1_over(r) 

351 return r 

352 

353 

354def _Dcosh(x, y): 

355 return _Dsincos(x, y, sinh, sinh) 

356 

357 

358def _DeatanhE(x, y, E): # see .albers._Datanhee 

359 # Deatanhe(x, y) = eatanhe((x - y) / (1 - e^2 * x * y)) / (x - y) 

360 e = _1_0 - E.e2 * x * y 

361 if e: # assert not isnear0(e) 

362 d = x - y 

363 e = (E._es_atanh(d / e) / d) if d else (E.e2 / e) 

364 return e 

365 

366 

367def _DfEt(tx, ty, eF): # tangents 

368 # eF = Elliptic(-E.e12) # -E.e2 / (1 - E.e2) 

369 r, x, y, = _1_0, atan(tx), atan(ty) 

370 d = x - y 

371 if (x * y) > 0: 

372 # See U{DLMF<https://DLMF.NIST.gov/19.11>}: 19.11.2 and 19.11.4 

373 # letting theta -> x, phi -> -y, psi -> z 

374 # (E(x) - E(y)) / d = E(z)/d - k2 * sin(x) * sin(y) * sin(z)/d 

375 # tan(z/2) = (sin(x)*Delta(y) - sin(y)*Delta(x)) / (cos(x) + cos(y)) 

376 # = d * Dsin(x,y) * (sin(x) + sin(y))/(cos(x) + cos(y)) / 

377 # (sin(x)*Delta(y) + sin(y)*Delta(x)) 

378 # = t = d * Dt 

379 # sin(z) = 2*t/(1+t^2); cos(z) = (1-t^2)/(1+t^2) 

380 # Alt (this only works for |z| <= pi/2 -- however, this conditions 

381 # holds if x*y > 0): 

382 # sin(z) = d * Dsin(x,y) * (sin(x) + sin(y)) / 

383 # (sin(x)*cos(y)*Delta(y) + sin(y)*cos(x)*Delta(x)) 

384 # cos(z) = sqrt((1-sin(z))*(1+sin(z))) 

385 sx, cx, sy, cy = sincos2_(x, y) 

386 D = (cx + cy) * (eF.fDelta(sy, cy) * sx + 

387 eF.fDelta(sx, cx) * sy) 

388 D = (sx + sy) * _Dsin(x, y) / D 

389 t = D * d 

390 t2 = _1_0 + t**2 

391 D *= _2_0 / t2 

392 s = D * d 

393 if s: 

394 c = (t + _1_0) * (_1_0 - t) / t2 

395 r = eF.fE(s, c, eF.fDelta(s, c)) / s 

396 r = D * (r - eF.k2 * sx * sy) 

397 elif d: 

398 r = (eF.fE(x) - eF.fE(y)) / d 

399 return r 

400 

401 

402def _Dgd(x, y): 

403 return _Datan(sinh(x), sinh(y)) * _Dsinh(x, y) 

404 

405 

406def _Dgdinv(x, y): # x, y are tangents 

407 return _Dasinh(x, y) / _Datan(x, y) 

408 

409 

410def _Dlog(x, y): 

411 d = (x - y) * _0_5 

412 # Changed atanh(t / (x + y)) to asinh(t / (2 * sqrt(x*y))) to 

413 # avoid taking atanh(1) when x is large and y is 1. This also 

414 # fixes bogus results being returned for the area when an endpoint 

415 # is at a pole. N.B. this routine is invoked with positive x 

416 # and y, so the sqrt is always taken of a positive quantity. 

417 return (asinh(d / sqrt(x * y)) / d) if d else _1_over(x) 

418 

419 

420def _Dsin(x, y): 

421 return _Dsincos(x, y, sin, cos) 

422 

423 

424def _Dsincos(x, y, sin_, cos_): 

425 r = cos_((x + y) * _0_5) 

426 d = (x - y) * _0_5 

427 if d: 

428 r *= sin_(d) / d 

429 return r 

430 

431 

432def _Dsinh(x, y): 

433 return _Dsincos(x, y, sinh, cosh) 

434 

435 

436def _Dtan(x, y): # PYCHOK no cover 

437 return _Dtant(x - y, tan(x), tan(y)) 

438 

439 

440def _Dtant(dxy, tx, ty): 

441 txy = tx * ty 

442 r = txy + _1_0 

443 if dxy: # 2 * txy > -1 == 2 * txy + 1 > 0 == txy + r > 0 == txy > -r 

444 r = ((tan(dxy) * r) if txy > -r else (tx - ty)) / dxy 

445 return r 

446 

447 

448def _Eaux4(E_aux, mu_psi_x, mu_psi_y): # degrees 

449 # get inverse auxiliary lats in radians and tangents 

450 phix = radians(E_aux(mu_psi_x, inverse=True)) 

451 phiy = radians(E_aux(mu_psi_y, inverse=True)) 

452 return phix, phiy, tan(phix), tan(phiy) 

453 

454 

455def _gd(x): 

456 return atan(sinh(x)) 

457 

458 

459def _sincosSeries(sinp, x, y, C, n): 

460 # N.B. C[] has n+1 elements of which 

461 # C[0] is ignored and n >= 0 

462 # Use Clenshaw summation to evaluate 

463 # m = (g(x) + g(y)) / 2 -- mean value 

464 # s = (g(x) - g(y)) / (x - y) -- average slope 

465 # where 

466 # g(x) = sum(C[j] * SC(2 * j * x), j = 1..n) 

467 # SC = sinp ? sin : cos 

468 # CS = sinp ? cos : sin 

469 # ... 

470 d, _neg = (x - y), neg 

471 sp, cp, sd, cd = sincos2_(x + y, d) 

472 sd = (sd / d) if d else _1_0 

473 s = _neg(sp * sd) # negative 

474 # 2x2 matrices in row-major order 

475 a1 = s * d**2 

476 a2 = s * _4_0 

477 a0 = a3 = _2_0 * cp * cd # m 

478 b2 = b1 = _0_0s(4) 

479 if n > 0: 

480 b1 = C[n], _0_0, _0_0, C[n] 

481 

482 _fsum = _MODS.fsums.fsum1f_ 

483 for j in range(n - 1, 0, -1): # C[0] unused 

484 b1, b2, Cj = b2, b1, C[j] 

485 # b1 = a * b2 - b1 + C[j] * I 

486 m0, m1, m2, m3 = b2 

487 n0, n1, n2, n3 = map(_neg, b1) 

488 b1 = (_fsum(a0 * m0, a1 * m2, n0, Cj), 

489 _fsum(a0 * m1, a1 * m3, n1), 

490 _fsum(a2 * m0, a3 * m2, n2), 

491 _fsum(a2 * m1, a3 * m3, n3, Cj)) 

492 # Here are the full expressions for m and s 

493 # f01, f02, f11, f12 = (0, 0, cd * sp, 2 * sd * cp) if sinp else \ 

494 # (1, 0, cd * cp, -2 * sd * sp) 

495 # m = -b2[1] * f02 + (C[0] - b2[0]) * f01 + b1[0] * f11 + b1[1] * f12 

496 # s = -b2[2] * f01 + (C[0] - b2[3]) * f02 + b1[2] * f11 + b1[3] * f12 

497 cd *= b1[2] 

498 sd *= b1[3] * _2_0 

499 s = _fsum(cd * sp, sd * cp) if sinp else \ 

500 _fsum(cd * cp, _neg(sd * sp), _neg(b2[2])) 

501 return s 

502 

503 

504_RACoeffs = { # Generated by Maxima on 2015-05-15 08:24:04-04:00 

505 4: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 4 

506 691, 7860, -20160, 18900, 0, 56700, # R[0]/n^0, polynomial(n), order 4 

507 1772, -5340, 6930, -4725, 14175, # R[1]/n^1, polynomial(n), order 3 

508 -1747, 1590, -630, 4725, # PYCHOK R[2]/n^2, polynomial(n), order 2 

509 104, -31, 315, # R[3]/n^3, polynomial(n), order 1 

510 -41, 420), # PYCHOK R[4]/n^4, polynomial(n), order 0, count = 20 

511 5: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 5 

512 -79036, 22803, 259380, -665280, 623700, 0, 1871100, # PYCHOK R[0]/n^0, polynomial(n), order 5 

513 41662, 58476, -176220, 228690, -155925, 467775, # PYCHOK R[1]/n^1, polynomial(n), order 4 

514 18118, -57651, 52470, -20790, 155925, # PYCHOK R[2]/n^2, polynomial(n), order 3 

515 -23011, 17160, -5115, 51975, # PYCHOK R[3]/n^3, polynomial(n), order 2 

516 5480, -1353, 13860, # PYCHOK R[4]/n^4, polynomial(n), order 1 

517 -668, 5775), # PYCHOK R[5]/n^5, polynomial(n), order 0, count = 27 

518 6: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 6 

519 128346268, -107884140, 31126095, 354053700, -908107200, 851350500, 0, 2554051500, # R[0]/n^0, polynomial(n), order 6 

520 -114456994, 56868630, 79819740, -240540300, 312161850, -212837625, 638512875, # PYCHOK R[1]/n^1, polynomial(n), order 5 

521 51304574, 24731070, -78693615, 71621550, -28378350, 212837625, # R[2]/n^2, polynomial(n), order 4 

522 1554472, -6282003, 4684680, -1396395, 14189175, # R[3]/n^3, polynomial(n), order 3 

523 -4913956, 3205800, -791505, 8108100, # PYCHOK R[4]/n^4, polynomial(n), order 2 

524 1092376, -234468, 2027025, # R[5]/n^5, polynomial(n), order 1 

525 -313076, 2027025), # PYCHOK R[6]/n^6, polynomial(n), order 0, count = 35 

526 7: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 7 

527 -317195588, 385038804, -323652420, 93378285, 1062161100, -2724321600, 2554051500, 0, 7662154500, # PYCHOK R[0]/n^0, polynomial(n), order 7 

528 258618446, -343370982, 170605890, 239459220, -721620900, 936485550, -638512875, 1915538625, # PYCHOK R[1]/n^1, polynomial(n), order 6 

529 -248174686, 153913722, 74193210, -236080845, 214864650, -85135050, 638512875, # PYCHOK R[2]/n^2, polynomial(n), order 5 

530 114450437, 23317080, -94230045, 70270200, -20945925, 212837625, # PYCHOK R[3]/n^3, polynomial(n), order 4 

531 15445736, -103193076, 67321800, -16621605, 170270100, # PYCHOK R[4]/n^4, polynomial(n), order 3 

532 -27766753, 16385640, -3517020, 30405375, # PYCHOK R[4]/n^4, polynomial(n), order 3 

533 4892722, -939228, 6081075, # PYCHOK R[4]/n^4, polynomial(n), order 3 

534 -3189007, 14189175), # PYCHOK R[7]/n^7, polynomial(n), order 0, count = 44 

535 8: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 8 

536 71374704821, -161769749880, 196369790040, -165062734200, 47622925350, 541702161000, -1389404016000, 1302566265000, 0, 3907698795000, # R[0]/n^0, polynomial(n), order 8 

537 -13691187484, 65947703730, -87559600410, 43504501950, 61062101100, -184013329500, 238803815250, -162820783125, 488462349375, # PYCHOK R[1]/n^1, polynomial(n), order 7 

538 30802104839, -63284544930, 39247999110, 18919268550, -60200615475, 54790485750, -21709437750, 162820783125, # R[2]/n^2, polynomial(n), order 6 

539 -8934064508, 5836972287, 1189171080, -4805732295, 3583780200, -1068242175, 10854718875, # PYCHOK R[3]/n^3, polynomial(n), order 5 

540 50072287748, 3938662680, -26314234380, 17167059000, -4238509275, 43418875500, # R[4]/n^4, polynomial(n), order 4 

541 359094172, -9912730821, 5849673480, -1255576140, 10854718875, # R[5]/n^5, polynomial(n), order 3 

542 -16053944387, 8733508770, -1676521980, 10854718875, # PYCHOK R[6]/n^6, polynomial(n), order 2 

543 930092876, -162639357, 723647925, # R[7]/n^7, polynomial(n), order 1 

544 -673429061, 1929727800) # PYCHOK R[8]/n^8, polynomial(n), order 0, count = 54 

545} 

546 

547__all__ += _ALL_DOCS(Caps) 

548 

549# **) MIT License 

550# 

551# Copyright (C) 2022-2024 -- mrJean1 at Gmail -- All Rights Reserved. 

552# 

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

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

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

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

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

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

559# 

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

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

562# 

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

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

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

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

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

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

569# OTHER DEALINGS IN THE SOFTWARE.