Coverage for pygeodesy/sphericalBase.py: 94%

251 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-12-02 13:46 -0500

1 

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

3 

4u'''(INTERNAL) Private spherical base classes C{CartesianSphericalBase} and 

5C{LatLonSphericalBase} for L{sphericalNvector} and L{sphericalTrigonometry}. 

6 

7A pure Python implementation of geodetic (lat-/longitude) functions, 

8transcoded in part from JavaScript originals by I{(C) Chris Veness 2011-2016} 

9and published under the same MIT Licence**, see 

10U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>}. 

11''' 

12# make sure int/int division yields float quotient, see .basics 

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

14 

15from pygeodesy.basics import _copysign, isbool, isinstanceof, isscalar, map1 

16from pygeodesy.cartesianBase import CartesianBase, Bearing2Tuple 

17from pygeodesy.constants import EPS, EPS0, PI, PI2, PI_2, R_M, \ 

18 _0_0, _0_5, _1_0, _180_0, _360_0, \ 

19 _over, isnear0, isnon0 

20from pygeodesy.datums import Datums, _earth_ellipsoid, _spherical_datum 

21from pygeodesy.errors import IntersectionError, _ValueError, \ 

22 _xattr, _xError 

23from pygeodesy.fmath import favg, fdot, hypot, sqrt_a 

24from pygeodesy.interns import NN, _COMMA_, _concentric_, _datum_, \ 

25 _distant_, _exceed_PI_radians_, _name_, \ 

26 _near_, _radius_, _too_ 

27from pygeodesy.latlonBase import LatLonBase, _trilaterate5 # PYCHOK passed 

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

29# from pygeodesy.namedTuples import Bearing2Tuple # from .cartesianBase 

30from pygeodesy.nvectorBase import NvectorBase, Fmt, _xattrs 

31from pygeodesy.props import deprecated_method, property_doc_, \ 

32 property_RO, _update_all 

33# from pygeodesy.streprs import Fmt, _xattrs # from .nvectorBase 

34from pygeodesy.units import Bearing, Bearing_, Radians_, Radius, \ 

35 Radius_, Scalar_, _100km 

36from pygeodesy.utily import acos1, asin1, atan2b, atan2d, degrees90, \ 

37 degrees180, sincos2, sincos2d, _unrollon, \ 

38 tanPI_2_2, wrapPI 

39 

40from math import cos, fabs, log, sin, sqrt 

41 

42__all__ = _ALL_LAZY.sphericalBase 

43__version__ = '23.12.01' 

44 

45 

46class CartesianSphericalBase(CartesianBase): 

47 '''(INTERNAL) Base class for spherical C{Cartesian}s. 

48 ''' 

49 _datum = Datums.Sphere # L{Datum} 

50 

51 def intersections2(self, rad1, other, rad2, radius=R_M): 

52 '''Compute the intersection points of two circles each defined 

53 by a center point and a radius. 

54 

55 @arg rad1: Radius of the this circle (C{meter} or C{radians}, 

56 see B{C{radius}}). 

57 @arg other: Center of the other circle (C{Cartesian}). 

58 @arg rad2: Radius of the other circle (C{meter} or C{radians}, 

59 see B{C{radius}}). 

60 @kwarg radius: Mean earth radius (C{meter} or C{None} if both 

61 B{C{rad1}} and B{C{rad2}} are given in C{radians}). 

62 

63 @return: 2-Tuple of the intersection points, each C{Cartesian}. 

64 For abutting circles, the intersection points are the 

65 same C{Cartesian} instance, aka the I{radical center}. 

66 

67 @raise IntersectionError: Concentric, antipodal, invalid or 

68 non-intersecting circles. 

69 

70 @raise TypeError: If B{C{other}} is not C{Cartesian}. 

71 

72 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}} or B{C{radius}}. 

73 

74 @see: U{Calculating intersection of two Circles 

75 <https://GIS.StackExchange.com/questions/48937/ 

76 calculating-intersection-of-two-circles>} and method 

77 or function C{trilaterate3d2}. 

78 ''' 

79 x1, x2 = self, self.others(other) 

80 r1, r2, x = _rads3(rad1, rad2, radius) 

81 if x: 

82 x1, x2 = x2, x1 

83 try: 

84 n, q = x1.cross(x2), x1.dot(x2) 

85 n2, q1 = n.length2, (_1_0 - q**2) 

86 if n2 < EPS or isnear0(q1): 

87 raise ValueError(_near_(_concentric_)) 

88 c1, c2 = cos(r1), cos(r2) 

89 x0 = x1.times((c1 - q * c2) / q1).plus( 

90 x2.times((c2 - q * c1) / q1)) 

91 n1 = _1_0 - x0.length2 

92 if n1 < EPS: 

93 raise ValueError(_too_(_distant_)) 

94 except ValueError as x: 

95 raise IntersectionError(center=self, rad1=rad1, 

96 other=other, rad2=rad2, cause=x) 

97 n = n.times(sqrt(n1 / n2)) 

98 if n.length > EPS: 

99 x1 = x0.plus(n) 

100 x2 = x0.minus(n) 

101 else: # abutting circles 

102 x1 = x2 = x0 

103 

104 return (_xattrs(x1, self, _datum_, _name_), 

105 _xattrs(x2, self, _datum_, _name_)) 

106 

107 @property_RO 

108 def sphericalCartesian(self): 

109 '''Get this C{Cartesian}'s spherical class. 

110 ''' 

111 return type(self) 

112 

113 

114class LatLonSphericalBase(LatLonBase): 

115 '''(INTERNAL) Base class for spherical C{LatLon}s. 

116 ''' 

117 _datum = Datums.Sphere # spherical L{Datum} 

118 _napieradius = _100km 

119 

120 def __init__(self, latlonh, lon=None, height=0, datum=None, wrap=False, name=NN): 

121 '''Create a spherical C{LatLon} point frome the given lat-, longitude and 

122 height on the given datum. 

123 

124 @arg latlonh: Latitude (C{degrees} or DMS C{str} with N or S suffix) or 

125 a previous C{LatLon} instance provided C{B{lon}=None}. 

126 @kwarg lon: Longitude (C{degrees} or DMS C{str} with E or W suffix) or 

127 C(None), indicating B{C{latlonh}} is a C{LatLon}. 

128 @kwarg height: Optional height above (or below) the earth surface (C{meter}, 

129 same units as the datum's ellipsoid axes or radius). 

130 @kwarg datum: Optional, spherical datum to use (L{Datum}, L{Ellipsoid}, 

131 L{Ellipsoid2}, L{a_f2Tuple}) or earth radius in C{meter}, 

132 conventionally). 

133 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{lat}} and B{C{lon}} 

134 (C{bool}). 

135 @kwarg name: Optional name (string). 

136 

137 @raise TypeError: If B{C{latlonh}} is not a C{LatLon} or B{C{datum}} not 

138 spherical. 

139 ''' 

140 LatLonBase.__init__(self, latlonh, lon=lon, height=height, wrap=wrap, name=name) 

141 if datum not in (None, self.datum): 

142 self.datum = datum 

143 

144 def bearingTo2(self, other, wrap=False, raiser=False): 

145 '''Return the initial and final bearing (forward and reverse 

146 azimuth) from this to an other point. 

147 

148 @arg other: The other point (C{LatLon}). 

149 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the 

150 B{C{other}} point (C{bool}). 

151 

152 @return: A L{Bearing2Tuple}C{(initial, final)}. 

153 

154 @raise TypeError: The B{C{other}} point is not spherical. 

155 

156 @see: Methods C{initialBearingTo} and C{finalBearingTo}. 

157 ''' 

158 # .initialBearingTo is inside .-Nvector and .-Trigonometry 

159 i = self.initialBearingTo(other, wrap=wrap, raiser=raiser) # PYCHOK .initialBearingTo 

160 f = self.finalBearingTo( other, wrap=wrap, raiser=raiser) 

161 return Bearing2Tuple(i, f, name=self.name) 

162 

163 @property_doc_(''' this point's datum (L{Datum}).''') 

164 def datum(self): 

165 '''Get this point's datum (L{Datum}). 

166 ''' 

167 return self._datum 

168 

169 @datum.setter # PYCHOK setter! 

170 def datum(self, datum): 

171 '''Set this point's datum I{without conversion} (L{Datum}, L{Ellipsoid}, 

172 L{Ellipsoid2}, L{a_f2Tuple}) or C{scalar} spherical earth radius). 

173 

174 @raise TypeError: If B{C{datum}} invalid or not not spherical. 

175 ''' 

176 d = _spherical_datum(datum, name=self.name, raiser=_datum_) 

177 if self._datum != d: 

178 _update_all(self) 

179 self._datum = d 

180 

181 def finalBearingTo(self, other, wrap=False, raiser=False): 

182 '''Return the final bearing (reverse azimuth) from this to 

183 an other point. 

184 

185 @arg other: The other point (spherical C{LatLon}). 

186 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll 

187 the B{C{other}} point (C{bool}). 

188 

189 @return: Final bearing (compass C{degrees360}). 

190 

191 @raise TypeError: The B{C{other}} point is not spherical. 

192 

193 @example: 

194 

195 >>> p = LatLon(52.205, 0.119) 

196 >>> q = LatLon(48.857, 2.351) 

197 >>> b = p.finalBearingTo(q) # 157.9 

198 ''' 

199 p = self.others(other) 

200 if wrap: 

201 p = _unrollon(self, p, wrap=wrap) 

202 # final bearing is the reverse of the other, initial one 

203 b = p.initialBearingTo(self, wrap=False, raiser=raiser) + _180_0 

204 return b if b < 360 else (b - _360_0) 

205 

206 def intersecant2(self, circle, point, other, radius=R_M, exact=False, # PYCHOK signature 

207 height=None, wrap=False): 

208 '''Compute the intersections of a circle and a (great circle) line 

209 given as two points or as a point and bearing. 

210 

211 @arg circle: Radius of the circle centered at this location (C{meter}, 

212 same units as B{C{radius}}) or a point on the circle 

213 (this C{LatLon}). 

214 @arg point: A point on the (great circle) line (this C{LatLon}). 

215 @arg other: An other point I{on} (this {LatLon}) or the bearing at 

216 B{C{point}} I{of} the (great circle) line (compass 

217 C{degrees}). 

218 @kwarg radius: Mean earth radius (C{meter}, conventionally). 

219 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth, 

220 destination and distance, if C{False} use the basic 

221 rhumb methods (C{bool}) or if C{None} use the I{great 

222 circle} methods. 

223 @kwarg height: Optional height for the intersection points (C{meter}, 

224 conventionally) or C{None} for interpolated heights. 

225 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the points 

226 B{C{circle}}, B{C{point}} and/or B{C{other}} (C{bool}). 

227 

228 @return: 2-Tuple of the intersection points (representing a chord), each 

229 an instance of the B{C{point}} class. Both points are the same 

230 instance if the (great circle) line is tangent to the circle. 

231 

232 @raise IntersectionError: The circle and line do not intersect. 

233 

234 @raise TypeError: If B{C{point}} is not this C{LatLon} or B{C{circle}} 

235 or B{C{other}} invalid. 

236 

237 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{radius}}, 

238 B{C{exact}}, B{C{height}} or B{C{napieradius}}. 

239 ''' 

240 p = self.others(point=point) 

241 try: 

242 return _intersecant2(self, circle, p, other, radius=radius, exact=exact, 

243 height=height, wrap=wrap) 

244 except (TypeError, ValueError) as x: 

245 raise _xError(x, center=self, circle=circle, point=point, other=other, 

246 radius=radius, exact=exact, height=height, wrap=wrap) 

247 

248 def maxLat(self, bearing): 

249 '''Return the maximum latitude reached when travelling on a great circle 

250 on given bearing from this point based on Clairaut's formula. 

251 

252 The maximum latitude is independent of longitude and the same for all 

253 points on a given latitude. 

254 

255 Negate the result for the minimum latitude (on the Southern hemisphere). 

256 

257 @arg bearing: Initial bearing (compass C{degrees360}). 

258 

259 @return: Maximum latitude (C{degrees90}). 

260 

261 @raise ValueError: Invalid B{C{bearing}}. 

262 ''' 

263 r = acos1(fabs(sin(Bearing_(bearing)) * cos(self.phi))) 

264 return degrees90(r) 

265 

266 def minLat(self, bearing): 

267 '''Return the minimum latitude reached when travelling on a great circle 

268 on given bearing from this point. 

269 

270 @arg bearing: Initial bearing (compass C{degrees360}). 

271 

272 @return: Minimum latitude (C{degrees90}). 

273 

274 @see: Method L{maxLat} for more details. 

275 

276 @raise ValueError: Invalid B{C{bearing}}. 

277 ''' 

278 return -self.maxLat(bearing) 

279 

280 def _mpr(self, radius=R_M, exact=None): # meter per radian 

281 if exact and not isscalar(radius): # see .rhumb.ekx.Rhumb._mpr 

282 radius = _earth_ellipsoid(radius)._Lpr 

283 return radius 

284 

285 @property_doc_(''' the I{Napier} radius to apply spherical trigonometry.''') 

286 def napieradius(self): 

287 '''Get the I{Napier} radius (C{meter}, conventionally). 

288 ''' 

289 return self._napieradius 

290 

291 @napieradius.setter # PYCHOK setter! 

292 def napieradius(self, radius): 

293 '''Set this I{Napier} radius (C{meter}, conventionally) or C{0}. 

294 

295 In methods L{intersecant2} and L{rhumbIntersecant2}, I{Napier}'s 

296 spherical trigonometry is applied if the circle radius exceeds 

297 the I{Napier} radius, otherwise planar trigonometry is used. 

298 

299 @raise UnitError: Invalid B{C{radius}}. 

300 ''' 

301 self._napieradius = Radius(napieradius=radius or 0) 

302 

303# def nearestTo(self, point, other, **radius_exact_height_wrap): # PYCHOK signature 

304# p = self.others(point=point) 

305# try: 

306# p, q = _intersecant2(self, p, p, other, **radius_exact_height_wrap) 

307# except (TypeError, ValueError) as x: 

308# raise _xError(x, this=self, point=point, other=other, **radius_exact_height_wrap) 

309# return p.midpointTo(q) 

310 

311 def parse(self, strllh, height=0, sep=_COMMA_, name=NN): 

312 '''Parse a string representing a similar, spherical C{LatLon} 

313 point, consisting of C{"lat, lon[, height]"}. 

314 

315 @arg strllh: Lat, lon and optional height (C{str}), 

316 see function L{pygeodesy.parse3llh}. 

317 @kwarg height: Optional, default height (C{meter}). 

318 @kwarg sep: Optional separator (C{str}). 

319 @kwarg name: Optional instance name (C{str}), 

320 overriding this name. 

321 

322 @return: The similar point (spherical C{LatLon}). 

323 

324 @raise ParseError: Invalid B{C{strllh}}. 

325 ''' 

326 t = _MODS.dms.parse3llh(strllh, height=height, sep=sep) 

327 r = self.classof(*t) 

328 if name: 

329 r.rename(name) 

330 return r 

331 

332 @property_RO 

333 def _radius(self): 

334 '''(INTERNAL) Get this sphere's radius. 

335 ''' 

336 return self.datum.ellipsoid.equatoradius 

337 

338 def _rhumbs3(self, other, wrap, r=False): # != .latlonBase._rhumbx3 

339 '''(INTERNAL) Rhumb_ helper function. 

340 

341 @arg other: The other point (spherical C{LatLon}). 

342 ''' 

343 p = self.others(other, up=2) 

344 if wrap: 

345 p = _unrollon(self, p, wrap=wrap) 

346 a2, b2 = p.philam 

347 a1, b1 = self.philam 

348 # if |db| > 180 take shorter rhumb 

349 # line across the anti-meridian 

350 db = wrapPI(b2 - b1) 

351 dp = _logPI_2_2(a2, a1) 

352 da = a2 - a1 

353 if r: 

354 # on Mercator projection, longitude distances shrink 

355 # by latitude; the 'stretch factor' q becomes ill- 

356 # conditioned along E-W line (0/0); use an empirical 

357 # tolerance to avoid it 

358 q = (da / dp) if fabs(dp) > EPS else cos(a1) 

359 da = hypot(da, q * db) # angular distance radians 

360 return da, db, dp 

361 

362 def rhumbAzimuthTo(self, other, radius=R_M, exact=False, wrap=False, b360=False): 

363 '''Return the azimuth (bearing) of a rhumb line (loxodrome) between 

364 this and an other (spherical) point. 

365 

366 @arg other: The other point (spherical C{LatLon}). 

367 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum}, 

368 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}). 

369 @kwarg exact: If C{True}, use I{Elliptic, Krüger} L{Rhumb} (C{bool}), 

370 default C{False} for backward compatibility. 

371 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the 

372 B{C{other}} point (C{bool}). 

373 @kwarg b360: If C{True}, return the azimuth in the bearing range. 

374 

375 @return: Rhumb azimuth (compass C{degrees180} or C{degrees360}). 

376 

377 @raise TypeError: The B{C{other}} point is incompatible or 

378 B{C{radius}} is invalid. 

379 

380 @example: 

381 

382 >>> p = LatLon(51.127, 1.338) 

383 >>> q = LatLon(50.964, 1.853) 

384 >>> b = p.rhumbBearingTo(q) # 116.7 

385 ''' 

386 if exact: # use series, always 

387 z = LatLonBase.rhumbAzimuthTo(self, other, exact=False, # Krüger 

388 radius=radius, wrap=wrap, b360=b360) 

389 else: 

390 _, db, dp = self._rhumbs3(other, wrap) 

391 z = (atan2b if b360 else atan2d)(db, dp) # see .rhumbBase.RhumbBase.Inverse 

392 return z 

393 

394 @deprecated_method 

395 def rhumbBearingTo(self, other): # unwrapped 

396 '''DEPRECATED, use method C{.rhumbAzimuthTo}.''' 

397 return self.rhumbAzimuthTo(other, b360=True) # [0..360) 

398 

399 def rhumbDestination(self, distance, azimuth, radius=R_M, height=None, exact=False): 

400 '''Return the destination point having travelled the given distance from 

401 this point along a rhumb line (loxodrome) of the given azimuth. 

402 

403 @arg distance: Distance travelled (C{meter}, same units as B{C{radius}}), 

404 may be negative if C{B{exact}=True}. 

405 @arg azimuth: Azimuth (bearing) of the rhumb line (compass C{degrees}). 

406 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum}, 

407 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}) if 

408 C{B{exact}=True}. 

409 @kwarg height: Optional height, overriding the default height (C{meter}. 

410 @kwarg exact: If C{True}, use I{Elliptic, Krüger} L{Rhumb} (C{bool}), 

411 default C{False} for backward compatibility. 

412 

413 @return: The destination point (spherical C{LatLon}). 

414 

415 @raise ValueError: Invalid B{C{distance}}, B{C{azimuth}}, B{C{radius}} 

416 or B{C{height}}. 

417 

418 @example: 

419 

420 >>> p = LatLon(51.127, 1.338) 

421 >>> q = p.rhumbDestination(40300, 116.7) # 50.9642°N, 001.8530°E 

422 ''' 

423 if exact: # use series, always 

424 r = LatLonBase.rhumbDestination(self, distance, azimuth, exact=False, # Krüger 

425 radius=radius, height=height) 

426 else: # radius=None from .rhumbMidpointTo 

427 if radius in (None, self._radius): 

428 d, r = self.datum, radius 

429 else: 

430 d = _spherical_datum(radius, raiser=_radius_) # spherical only 

431 r = d.ellipsoid.equatoradius 

432 r = _m2radians(distance, r, low=-EPS) # distance=0 from .rhumbMidpointTo 

433 

434 a1, b1 = self.philam 

435 sb, cb = sincos2(Bearing_(azimuth)) # radians 

436 

437 da = r * cb 

438 a2 = a1 + da 

439 # normalize latitude if past pole 

440 if fabs(a2) > PI_2: 

441 a2 = _copysign(PI, a2) - a2 

442 

443 dp = _logPI_2_2(a2, a1) 

444 # q becomes ill-conditioned on E-W course 0/0 

445 q = cos(a1) if isnear0(dp) else (da / dp) 

446 b2 = b1 if isnear0(q) else (b1 + r * sb / q) 

447 

448 h = self._heigHt(height) 

449 r = self.classof(degrees90(a2), degrees180(b2), datum=d, height=h) 

450 return r 

451 

452 def rhumbDistanceTo(self, other, radius=R_M, exact=False, wrap=False): 

453 '''Return the distance from this to an other point along 

454 a rhumb line (loxodrome). 

455 

456 @arg other: The other point (spherical C{LatLon}). 

457 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum}, 

458 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}) if 

459 C{B{exact}=True}. 

460 @kwarg exact: If C{True}, use I{Elliptic, Krüger} L{Rhumb} (C{bool}), 

461 default C{False} for backward compatibility. 

462 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the 

463 B{C{other}} point (C{bool}). 

464 

465 @return: Distance (C{meter}, the same units as B{C{radius}} 

466 or C{radians} if B{C{radius}} is C{None}). 

467 

468 @raise TypeError: The B{C{other}} point is incompatible. 

469 

470 @raise ValueError: Invalid B{C{radius}}. 

471 

472 @example: 

473 

474 >>> p = LatLon(51.127, 1.338) 

475 >>> q = LatLon(50.964, 1.853) 

476 >>> d = p.rhumbDistanceTo(q) # 403100 

477 ''' 

478 if exact: # use series, always 

479 r = LatLonBase.rhumbDistanceTo(self, other, exact=False, # Krüger 

480 radius=radius, wrap=wrap) 

481 if radius is None: # angular distance in radians 

482 r = r / self._radius # /= chokes PyChecker 

483 else: 

484 # see <https://www.EdWilliams.org/avform.htm#Rhumb> 

485 r, _, _ = self._rhumbs3(other, wrap, r=True) 

486 if radius is not None: 

487 r *= Radius(radius) 

488 return r 

489 

490 def rhumbIntersecant2(self, circle, point, other, radius=R_M, exact=True, # PYCHOK signature 

491 height=None, wrap=False): 

492 '''Compute the intersections of a circle and a rhumb line given as two 

493 points and as a point and azimuth. 

494 

495 @arg circle: Radius of the circle centered at this location (C{meter}, 

496 same units as B{C{radius}}) or a point on the circle 

497 (this C{LatLon}). 

498 @arg point: The rhumb line's start point (this C{LatLon}). 

499 @arg other: An other point (this I{on} C{LatLon}) or the azimuth I{of} 

500 (compass C{degrees}) the rhumb line. 

501 @kwarg radius: Mean earth radius (C{meter}, conventionally). 

502 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth, 

503 destination and distance, if C{False} use the basic 

504 rhumb methods (C{bool}) or if C{None} use the I{great 

505 circle} methods. 

506 @kwarg height: Optional height for the intersection points (C{meter}, 

507 conventionally) or C{None}. 

508 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the points 

509 B{C{circle}}, B{C{point}} and/or B{C{other}} (C{bool}). 

510 

511 @return: 2-Tuple of the intersection points (representing a chord), 

512 each an instance of this class. For a tangent line, both 

513 points are the same instance, wrapped or I{normalized}. 

514 

515 @raise IntersectionError: The circle and line do not intersect. 

516 

517 @raise TypeError: If B{C{point}} is not this C{LatLon} or B{C{circle}} 

518 or B{C{other}} invalid. 

519 

520 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{radius}}, 

521 B{C{exact}} or B{C{height}}. 

522 ''' 

523 m = LatLonBase.rhumbIntersecant2 if exact else \ 

524 LatLonSphericalBase.intersecant2 

525 return m(self, circle, point, other, radius=radius, exact=exact, 

526 height=height, wrap=wrap) 

527 

528 def rhumbMidpointTo(self, other, height=None, radius=R_M, exact=False, 

529 fraction=_0_5, wrap=False): 

530 '''Return the (loxodromic) midpoint on the rhumb line between 

531 this and an other point. 

532 

533 @arg other: The other point (spherical LatLon). 

534 @kwarg height: Optional height, overriding the mean height (C{meter}). 

535 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum}, 

536 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}). 

537 @kwarg exact: If C{True}, use I{Elliptic, Krüger} L{Rhumb} (C{bool}), 

538 default C{False} for backward compatibility. 

539 @kwarg fraction: Midpoint location from this point (C{scalar}), may 

540 be negative if C{B{exact}=True}. 

541 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{other}} 

542 point (C{bool}). 

543 

544 @return: The (mid)point at the given B{C{fraction}} along the rhumb 

545 line (spherical C{LatLon}). 

546 

547 @raise TypeError: The B{C{other}} point is incompatible. 

548 

549 @raise ValueError: Invalid B{C{height}} or B{C{fraction}} 

550 

551 @example: 

552 

553 >>> p = LatLon(51.127, 1.338) 

554 >>> q = LatLon(50.964, 1.853) 

555 >>> m = p.rhumb_midpointTo(q) 

556 >>> m.toStr() # '51.0455°N, 001.5957°E' 

557 ''' 

558 if exact: # use series, always 

559 r = LatLonBase.rhumbMidpointTo(self, other, exact=False, # Krüger 

560 radius=radius, height=height, 

561 fraction=fraction, wrap=wrap) 

562 elif fraction is not _0_5: 

563 f = Scalar_(fraction=fraction) # low=_0_0 

564 r, db, dp = self._rhumbs3(other, wrap, r=True) # radians 

565 z = atan2b(db, dp) 

566 h = self._havg(other, f=f, h=height) 

567 r = self.rhumbDestination(r * f, z, radius=None, height=h) 

568 

569 else: # for backward compatibility, unwrapped 

570 # see <https://MathForum.org/library/drmath/view/51822.html> 

571 a1, b1 = self.philam 

572 a2, b2 = self.others(other).philam 

573 

574 if fabs(b2 - b1) > PI: 

575 b1 += PI2 # crossing anti-meridian 

576 

577 a3 = favg(a1, a2) 

578 b3 = favg(b1, b2) 

579 

580 f1 = tanPI_2_2(a1) 

581 if isnon0(f1): 

582 f2 = tanPI_2_2(a2) 

583 f = f2 / f1 

584 if isnon0(f): 

585 f = log(f) 

586 if isnon0(f): 

587 f3 = tanPI_2_2(a3) 

588 b3 = fdot(map1(log, f1, f2, f3), 

589 -b2, b1, b2 - b1) / f 

590 

591 d = self.datum if radius in (None, self._radius) else \ 

592 _spherical_datum(radius, name=self.name, raiser=_radius_) 

593 h = self._havg(other, h=height) 

594 r = self.classof(degrees90(a3), degrees180(b3), datum=d, height=h) 

595 return r 

596 

597 @property_RO 

598 def sphericalLatLon(self): 

599 '''Get this C{LatLon}'s spherical class. 

600 ''' 

601 return type(self) 

602 

603 def toNvector(self, Nvector=NvectorBase, **Nvector_kwds): # PYCHOK signature 

604 '''Convert this point to C{Nvector} components, I{including 

605 height}. 

606 

607 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}} 

608 keyword arguments, ignored if 

609 C{B{Nvector} is None}. 

610 

611 @return: An B{C{Nvector}} or a L{Vector4Tuple}C{(x, y, z, h)} 

612 if B{C{Nvector}} is C{None}. 

613 

614 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}. 

615 ''' 

616 return LatLonBase.toNvector(self, Nvector=Nvector, **Nvector_kwds) 

617 

618 

619def _intersecant2(c, r, p, b, radius=R_M, exact=False, height=None, wrap=False): 

620 # (INTERNAL) Intersect a circle and line, see L{intersecant2} 

621 # above, separated to allow callers to embellish any exceptions 

622 

623 if wrap: 

624 p = _unrollon(c, p, wrap=wrap) 

625 nonexact = exact is None 

626 

627 if not isinstanceof(r, c.__class__, p.__class__): 

628 r = Radius_(circle=r) 

629 elif nonexact: 

630 r = c.distanceTo(r, radius=radius, wrap=wrap) 

631 elif isbool(exact): 

632 r = c.rhumbDistanceTo(r, radius=radius, exact=exact, wrap=wrap) 

633 else: 

634 raise _ValueError(exact=exact) 

635 

636 if not isinstanceof(b, c.__class__, p.__class__): 

637 b = Bearing(b) 

638 elif nonexact: 

639 b = p.initialBearingTo(b, wrap=wrap) 

640 else: 

641 b = p.rhumbAzimuthTo(b, radius=radius, exact=exact, wrap=wrap, 

642 b360=True) 

643 

644 d = p.distanceTo(c, radius=radius) if nonexact else \ 

645 p.rhumbDistanceTo(c, radius=radius, exact=exact) 

646 if d > EPS0: 

647 n = _xattr(c, napieradius=0) 

648 a = p.initialBearingTo(c) if nonexact else \ 

649 p.rhumbAzimuthTo(c, radius=radius, exact=exact, b360=True) 

650 s, c = sincos2d(b - a) # Napier's sin(A), cos(A) 

651 if r > n: 

652 # Napier's right spherical triangle rules (R2) and (R1) 

653 # <https://WikiPedia.org/wiki/Spherical_trigonometry> 

654 m = p._mpr(radius=radius, exact=exact) # meter per radian 

655 if fabs(c) > EPS0: 

656 d = d / m # /= chokes PyChecker 

657 a = asin1(sin(d) * fabs(s)) # Napier's a 

658 c = _copysign(cos(a), c) 

659 d = acos1(cos(d) / c) * m 

660 a *= m # meter 

661 else: # point and chord center coincident 

662 a, d = d, 0 

663 c = cos(a / m) 

664 h = (acos1(cos(r / m) / c) * m) if a < r else 0 

665 else: # distance from the chord center to ... 

666 a = fabs(d * s) # ... the cicle center ... 

667 d *= c # ... and to the point 

668 h = sqrt_a(r, a) if a < r else 0 # half chord length 

669 if a > r: 

670 raise IntersectionError(_too_(Fmt.distant(a))) 

671 else: 

672 d, h = 0, r # point and circle center coincident 

673 

674 _intersecant1, kwds = (p.destination, {}) if nonexact else \ 

675 (p.rhumbDestination, dict(exact=exact)) 

676 kwds.update(radius=radius, height=height) 

677 t = (_intersecant1(d + h, b, **kwds),) 

678 if h: 

679 t += (_intersecant1(d - h, b, **kwds),) 

680 else: # same instance twice 

681 t *= 2 

682 return t 

683 

684 

685def _logPI_2_2(a2, a1): 

686 '''(INTERNAL) C{log} of C{tanPI_2_2}'s quotient. 

687 ''' 

688 return log(_over(tanPI_2_2(a2), tanPI_2_2(a1))) 

689 

690 

691def _m2radians(distance, radius, low=EPS): # PYCHOK in .spherical* 

692 '''(INTERNAL) Distance in C{meter} to angular distance in C{radians}. 

693 

694 @raise UnitError: Invalid B{C{distance}} or B{C{radius}}. 

695 ''' 

696 r = float(distance) 

697 if radius: 

698 r = r / Radius_(radius=radius) # /= chokes PyChecker 

699 if low is not None: 

700 # small near0 values from .rhumbDestination not exact OK 

701 r = _0_0 if low < 0 and r < 0 else Radians_(r, low=low) 

702 # _0_0 if low < 0 and low < r < 0 else Radians_(r, low=low) 

703 return r 

704 

705 

706def _radians2m(rad, radius): 

707 '''(INTERNAL) Angular distance in C{radians} to distance in C{meter}. 

708 ''' 

709 if radius is not None: # not in (None, _0_0) 

710 rad *= R_M if radius is R_M else Radius(radius) 

711 return rad 

712 

713 

714def _rads3(rad1, rad2, radius): # in .sphericalTrigonometry 

715 '''(INTERNAL) Convert radii to radians. 

716 ''' 

717 r1 = Radius_(rad1=rad1) 

718 r2 = Radius_(rad2=rad2) 

719 if radius is not None: # convert radii to radians 

720 r1 = _m2radians(r1, radius) 

721 r2 = _m2radians(r2, radius) 

722 

723 x = r1 < r2 

724 if x: 

725 r1, r2 = r2, r1 

726 if r1 > PI: 

727 raise IntersectionError(rad1=rad1, rad2=rad2, 

728 txt=_exceed_PI_radians_) 

729 return r1, r2, x 

730 

731 

732__all__ += _ALL_DOCS(CartesianSphericalBase, LatLonSphericalBase) 

733 

734# **) MIT License 

735# 

736# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

737# 

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

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

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

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

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

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

744# 

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

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

747# 

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

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

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

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

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

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

754# OTHER DEALINGS IN THE SOFTWARE.