Coverage for pygeodesy/sphericalBase.py: 95%

217 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-10-04 12:08 -0400

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, map1 

16from pygeodesy.cartesianBase import CartesianBase, Bearing2Tuple 

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

18 _umod_360, isnear0, isnon0, _over, \ 

19 _0_0, _0_5, _1_0, _180_0 

20from pygeodesy.datums import Datums, _spherical_datum 

21from pygeodesy.errors import IntersectionError, _ValueError, _xError 

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

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

24 _distant_, _exceed_PI_radians_, _name_, \ 

25 _near_, _radius_, _too_ 

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

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

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

29from pygeodesy.nvectorBase import NvectorBase, Fmt, _xattrs 

30from pygeodesy.props import deprecated_method, property_doc_, \ 

31 property_RO, _update_all 

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

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

34 Radius_, Scalar_ 

35from pygeodesy.utily import acos1, atan2b, atan2d, degrees90, \ 

36 degrees180, sincos2, sincos2d, \ 

37 tanPI_2_2, _unrollon, wrap360, wrapPI 

38 

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

40 

41__all__ = _ALL_LAZY.sphericalBase 

42__version__ = '23.09.28' 

43 

44 

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

46 '''(INTERNAL) Return an angular distance in C{radians}. 

47 

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

49 ''' 

50 r = float(distance) 

51 if radius: 

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

53 if low is not None: 

54 # small near0 values from .rhumbDestination not exact OK 

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

56 return r 

57 

58 

59def _logPI_2_2(a2, a1): 

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

61 ''' 

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

63 

64 

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

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

67 ''' 

68 r1 = Radius_(rad1=rad1) 

69 r2 = Radius_(rad2=rad2) 

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

71 r1 = _angular(r1, radius) 

72 r2 = _angular(r2, radius) 

73 

74 x = r1 < r2 

75 if x: 

76 r1, r2 = r2, r1 

77 if r1 > PI: 

78 raise IntersectionError(rad1=rad1, rad2=rad2, 

79 txt=_exceed_PI_radians_) 

80 return r1, r2, x 

81 

82 

83class CartesianSphericalBase(CartesianBase): 

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

85 ''' 

86 _datum = Datums.Sphere # L{Datum} 

87 

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

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

90 by a center point and a radius. 

91 

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

93 see B{C{radius}}). 

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

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

96 see B{C{radius}}). 

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

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

99 

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

101 For abutting circles, the intersection points are the 

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

103 

104 @raise IntersectionError: Concentric, antipodal, invalid or 

105 non-intersecting circles. 

106 

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

108 

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

110 

111 @see: U{Calculating intersection of two Circles 

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

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

114 or function C{trilaterate3d2}. 

115 ''' 

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

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

118 if x: 

119 x1, x2 = x2, x1 

120 try: 

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

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

123 if n2 < EPS or isnear0(q1): 

124 raise ValueError(_near_(_concentric_)) 

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

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

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

128 n1 = _1_0 - x0.length2 

129 if n1 < EPS: 

130 raise ValueError(_too_(_distant_)) 

131 except ValueError as x: 

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

133 other=other, rad2=rad2, cause=x) 

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

135 if n.length > EPS: 

136 x1 = x0.plus(n) 

137 x2 = x0.minus(n) 

138 else: # abutting circles 

139 x1 = x2 = x0 

140 

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

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

143 

144 

145class LatLonSphericalBase(LatLonBase): 

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

147 ''' 

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

149 

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

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

152 height on the given datum. 

153 

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

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

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

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

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

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

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

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

162 conventionally). 

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

164 (C{bool}). 

165 @kwarg name: Optional name (string). 

166 

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

168 spherical. 

169 ''' 

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

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

172 self.datum = datum 

173 

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

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

176 azimuth) from this to an other point. 

177 

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

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

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

181 

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

183 

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

185 

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

187 ''' 

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

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

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

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

192 

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

194 def datum(self): 

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

196 ''' 

197 return self._datum 

198 

199 @datum.setter # PYCHOK setter! 

200 def datum(self, datum): 

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

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

203 

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

205 ''' 

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

207 if self._datum != d: 

208 _update_all(self) 

209 self._datum = d 

210 

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

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

213 an other point. 

214 

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

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

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

218 

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

220 

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

222 

223 @example: 

224 

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

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

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

228 ''' 

229 p = self.others(other) 

230 if wrap: 

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

232 # final bearing is the reverse of the other, initial one; 

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

234 b = p.initialBearingTo(self, wrap=False, raiser=raiser) 

235 return _umod_360(b + _180_0) 

236 

237 def intersecant2(self, circle, point, bearing, radius=R_M, exact=False, 

238 height=None, wrap=False): # was=True 

239 '''Compute the intersections of a circle and a line. 

240 

241 @arg circle: Radius of the circle centered at this location 

242 (C{meter}, same units as B{C{radius}}) or a point 

243 on the circle (this C{LatLon}). 

244 @arg point: An other point in- or outside the circle (this C{LatLon}). 

245 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360}) 

246 or a second point on the line (this C{LatLon}). 

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

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

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

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

251 circle} methods. 

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

253 conventionally) or C{None}. 

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

255 B{C{point}}, B{C{circle}} and/or B{C{bearing}} (C{bool}). 

256 

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

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

259 points are the same instance, the B{C{point}} or wrapped 

260 or I{normalized}. 

261 

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

263 

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

265 or B{C{bearing}} invalid. 

266 

267 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}}, 

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

269 ''' 

270 p = self.others(point=point) 

271 try: 

272 return _intersecant2(self, circle, p, bearing, radius=radius, exact=exact, 

273 height=height, wrap=wrap) 

274 except (TypeError, ValueError) as x: 

275 raise _xError(x, center=self, circle=circle, point=point, bearing=bearing, 

276 exact=exact, wrap=wrap) 

277 

278 def maxLat(self, bearing): 

279 '''Return the maximum latitude reached when travelling 

280 on a great circle on given bearing from this point 

281 based on Clairaut's formula. 

282 

283 The maximum latitude is independent of longitude 

284 and the same for all points on a given latitude. 

285 

286 Negate the result for the minimum latitude (on the 

287 Southern hemisphere). 

288 

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

290 

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

292 

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

294 ''' 

295 m = acos1(fabs(sin(Bearing_(bearing)) * cos(self.phi))) 

296 return degrees90(m) 

297 

298 def minLat(self, bearing): 

299 '''Return the minimum latitude reached when travelling 

300 on a great circle on given bearing from this point. 

301 

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

303 

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

305 

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

307 

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

309 ''' 

310 return -self.maxLat(bearing) 

311 

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

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

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

315 

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

317 see function L{pygeodesy.parse3llh}. 

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

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

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

321 overriding this name. 

322 

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

324 

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

326 ''' 

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

328 r = self.classof(*t) 

329 if name: 

330 r.rename(name) 

331 return r 

332 

333 @property_RO 

334 def _radius(self): 

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

336 ''' 

337 return self.datum.ellipsoid.equatoradius 

338 

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

340 '''(INTERNAL) Rhumb_ helper function. 

341 

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

343 ''' 

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

345 if wrap: 

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

347 a2, b2 = p.philam 

348 a1, b1 = self.philam 

349 # if |db| > 180 take shorter rhumb 

350 # line across the anti-meridian 

351 db = wrapPI(b2 - b1) 

352 dp = _logPI_2_2(a2, a1) 

353 da = a2 - a1 

354 if r: 

355 # on Mercator projection, longitude distances shrink 

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

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

358 # tolerance to avoid it 

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

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

361 return da, db, dp 

362 

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

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

365 between this and an other (spherical) point. 

366 

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

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

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

370 @kwarg exact: If C{True}, use I{Krüger} L{rhumbx} (C{bool}), 

371 default C{False} for backward compatibility. 

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

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

374 

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

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) 

389 else: 

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

391 z = atan2d(db, dp) # see .rhumbx.Rhumb.Inverse 

392 return z 

393 

394 @deprecated_method 

395 def rhumbBearingTo(self, other): # unwrapped 

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

397 return wrap360(self.rhumbAzimuthTo(other)) # [0..360) 

398 

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

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

401 from this point along a rhumb line (loxodrome) at the given bearing. 

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 bearing: Bearing (azimuth) at this point (compass C{degrees360}). 

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 

410 (C{meter}, same unit as B{C{radius}}). 

411 @kwarg exact: If C{True}, use I{Krüger} L{rhumbx} (C{bool}), 

412 default C{False} for backward compatibility. 

413 

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

415 

416 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}}, 

417 B{C{radius}} or B{C{height}}. 

418 

419 @example: 

420 

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

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

423 ''' 

424 if exact: # use series, always 

425 r = LatLonBase.rhumbDestination(self, distance, bearing, exact=False, # Krüger 

426 radius=radius, height=height) 

427 else: # radius=None from .rhumbMidpointTo 

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

429 d, r = self.datum, radius 

430 else: 

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

432 r = d.ellipsoid.equatoradius 

433 r = _angular(distance, r, low=-EPS) # distance=0 from .rhumbMidpointTo 

434 

435 a1, b1 = self.philam 

436 sb, cb = sincos2(Bearing_(bearing)) 

437 

438 da = r * cb 

439 a2 = a1 + da 

440 # normalize latitude if past pole 

441 if fabs(a2) > PI_2: 

442 a2 = _copysign(PI, a2) - a2 

443 

444 dp = _logPI_2_2(a2, a1) 

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

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

447 b2 = (b1 + r * sb / q) if fabs(q) > EPS else b1 

448 

449 h = self._heigHt(height) 

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

451 return r 

452 

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

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

455 a rhumb line (loxodrome). 

456 

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

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

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

460 C{B{exact}=True}. 

461 @kwarg exact: If C{True}, use I{Krüger} L{rhumbx} (C{bool}), 

462 default C{False} for backward compatibility. 

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

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

465 

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

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

468 

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

470 

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

472 

473 @example: 

474 

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

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

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

478 ''' 

479 if exact: # use series, always 

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

481 radius=radius, wrap=wrap) 

482 if radius is None: # angular distance in radians 

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

484 else: 

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

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

487 if radius is not None: 

488 r *= Radius(radius) 

489 return r 

490 

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

492 fraction=_0_5, wrap=False): 

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

494 this and an other point. 

495 

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

497 @kwarg height: Optional height, overriding the mean height 

498 (C{meter}). 

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

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

501 @kwarg exact: If C{True}, use I{Krüger} L{rhumbx} (C{bool}), 

502 default C{False} for backward compatibility. 

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

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

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

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

507 

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

509 the rhumb line (spherical C{LatLon}). 

510 

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

512 

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

514 

515 @example: 

516 

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

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

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

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

521 ''' 

522 if exact: # use series, always 

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

524 radius=radius, height=height, 

525 fraction=fraction, wrap=wrap) 

526 elif fraction is not _0_5: 

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

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

529 z = atan2b(db, dp) 

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

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

532 

533 else: # for backward compatibility, unwrapped 

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

535 a1, b1 = self.philam 

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

537 

538 if fabs(b2 - b1) > PI: 

539 b1 += PI2 # crossing anti-meridian 

540 

541 a3 = favg(a1, a2) 

542 b3 = favg(b1, b2) 

543 

544 f1 = tanPI_2_2(a1) 

545 if isnon0(f1): 

546 f2 = tanPI_2_2(a2) 

547 f = f2 / f1 

548 if isnon0(f): 

549 f = log(f) 

550 if isnon0(f): 

551 f3 = tanPI_2_2(a3) 

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

553 -b2, b1, b2 - b1) / f 

554 

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

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

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

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

559 return r 

560 

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

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

563 height}. 

564 

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

566 keyword arguments, ignored if 

567 C{B{Nvector} is None}. 

568 

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

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

571 

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

573 ''' 

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

575 

576 

577def _intersecant2(c, r, p, b, radius=R_M, exact=False, 

578 height=None, wrap=False): 

579 # (INTERNAL) Intersect a circle and bearing, see L{intersecant2} 

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

581 

582 if wrap: 

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

584 nonexact = exact is None 

585 

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

587 r = Radius_(circle=r) 

588 elif nonexact: 

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

590 elif isbool(exact): 

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

592 else: 

593 raise _ValueError(exact=exact) 

594 

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

596 b = Bearing(b) 

597 elif nonexact: 

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

599 else: 

600 b = p.rhumbAzimuthTo(b, radius=radius, exact=exact, wrap=wrap) 

601 

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

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

604 if d > EPS: 

605 a = p.initialBearingTo(c) if nonexact else wrap360( 

606 p.rhumbAzimuthTo(c, radius=radius, exact=exact)) 

607 s, c = sincos2d(b - a) 

608 s = sqrt_a(r, fabs(s * d)) 

609 if s > r: 

610 raise IntersectionError(_too_(Fmt.distant(s))) 

611 elif (r - s) < EPS: 

612 return p, p # tangent 

613 c *= d 

614 else: # p and c coincide 

615 s, c = r, 0 

616 t = () 

617 for d, b in ((s + c, b), (s - c, b + _180_0)): # bearing direction first 

618 t += (p.destination(d, b, radius=radius, height=height) if nonexact else 

619 p.rhumbDestination(d, b, radius=radius, height=height, exact=exact)), 

620 return t 

621 

622 

623def _r2m(r, radius): 

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

625 ''' 

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

627 r *= R_M if radius is R_M else Radius(radius) 

628 return r 

629 

630 

631__all__ += _ALL_DOCS(CartesianSphericalBase, LatLonSphericalBase) 

632 

633# **) MIT License 

634# 

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

636# 

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

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

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

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

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

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

643# 

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

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

646# 

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

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

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

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

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

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

653# OTHER DEALINGS IN THE SOFTWARE.