Coverage for pygeodesy/sphericalTrigonometry.py: 94%

390 statements  

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

1 

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

3 

4u'''Spherical, C{trigonometry}-based geodesy. 

5 

6Trigonometric classes geodetic (lat-/longitude) L{LatLon} and 

7geocentric (ECEF) L{Cartesian} and functions L{areaOf}, L{intersection}, 

8L{intersections2}, L{isPoleEnclosedBy}, L{meanOf}, L{nearestOn3} and 

9L{perimeterOf}, I{all spherical}. 

10 

11Pure Python implementation of geodetic (lat-/longitude) methods using 

12spherical trigonometry, transcoded from JavaScript originals by 

13I{(C) Chris Veness 2011-2016} published under the same MIT Licence**, see 

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

15''' 

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

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

18 

19from pygeodesy.basics import copysign0, isscalar, map1, signOf 

20from pygeodesy.constants import EPS, EPS1, EPS4, PI, PI2, PI_2, PI_4, R_M, \ 

21 isnear0, isnear1, isnon0, _0_0, _0_5, \ 

22 _1_0, _2_0, _90_0 

23from pygeodesy.datums import _ellipsoidal_datum, _mean_radius 

24from pygeodesy.errors import _AssertionError, CrossError, crosserrors, \ 

25 _ValueError, IntersectionError, _xError, \ 

26 _xkwds, _xkwds_get, _xkwds_pop 

27from pygeodesy.fmath import favg, fdot, fmean, hypot 

28from pygeodesy.fsums import Fsum, fsum, fsumf_ 

29from pygeodesy.formy import antipode_, bearing_, _bearingTo2, excessAbc_, \ 

30 excessGirard_, excessLHuilier_, opposing_, _radical2, \ 

31 vincentys_ 

32from pygeodesy.interns import _1_, _2_, _coincident_, _composite_, _colinear_, \ 

33 _concentric_, _convex_, _end_, _infinite_, _invalid_,\ 

34 _line_, _near_, _not_, _null_, _parallel_, _point_, \ 

35 _SPACE_, _too_ 

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

37# from pygeodesy.named import notImplemented # from .points 

38# from pygeodesy.nvectorBase import NvectorBase, sumOf # _MODE 

39from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, NearestOn3Tuple, \ 

40 Triangle7Tuple, Triangle8Tuple 

41from pygeodesy.points import ispolar, nearestOn5 as _nearestOn5, \ 

42 Fmt as _Fmt, notImplemented # XXX shadowed 

43from pygeodesy.props import deprecated_function, deprecated_method 

44from pygeodesy.sphericalBase import _m2radians, CartesianSphericalBase, \ 

45 _intersecant2, LatLonSphericalBase, \ 

46 _rads3, _radians2m, _trilaterate5 

47# from pygeodesy.streprs import Fmt as _Fmt # from .points XXX shadowed 

48from pygeodesy.units import Bearing_, Height, Lam_, Phi_, Radius_, Scalar 

49from pygeodesy.utily import acos1, asin1, atan1d, atan2d, degrees90, degrees180, \ 

50 degrees2m, m2radians, radiansPI2, sincos2_, tan_2, \ 

51 unrollPI, _unrollon, _unrollon3, _Wrap, wrap180, wrapPI 

52from pygeodesy.vector3d import sumOf, Vector3d 

53 

54from math import asin, atan2, cos, degrees, fabs, radians, sin 

55 

56__all__ = _ALL_LAZY.sphericalTrigonometry 

57__version__ = '23.10.24' 

58 

59_PI_EPS4 = PI - EPS4 

60if _PI_EPS4 >= PI: 

61 raise _AssertionError(EPS4=EPS4, PI=PI, PI_EPS4=_PI_EPS4) 

62 

63 

64class Cartesian(CartesianSphericalBase): 

65 '''Extended to convert geocentric, L{Cartesian} points to 

66 spherical, geodetic L{LatLon}. 

67 ''' 

68 

69 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon 

70 '''Convert this cartesian point to a C{spherical} geodetic point. 

71 

72 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword 

73 arguments. Use C{B{LatLon}=...} to override 

74 this L{LatLon} class or specify C{B{LatLon}=None}. 

75 

76 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is C{None}, 

77 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} 

78 with C{C} and C{M} if available. 

79 

80 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument. 

81 ''' 

82 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum) 

83 return CartesianSphericalBase.toLatLon(self, **kwds) 

84 

85 

86class LatLon(LatLonSphericalBase): 

87 '''New point on spherical model earth model. 

88 

89 @example: 

90 

91 >>> p = LatLon(52.205, 0.119) # height=0 

92 ''' 

93 

94 def _ab1_ab2_db5(self, other, wrap): 

95 '''(INTERNAL) Helper for several methods. 

96 ''' 

97 a1, b1 = self.philam 

98 a2, b2 = self.others(other, up=2).philam 

99 if wrap: 

100 a2, b2 = _Wrap.philam(a2, b2) 

101 db, b2 = unrollPI(b1, b2, wrap=wrap) 

102 else: # unrollPI shortcut 

103 db = b2 - b1 

104 return a1, b1, a2, b2, db 

105 

106 def alongTrackDistanceTo(self, start, end, radius=R_M, wrap=False): 

107 '''Compute the (signed) distance from the start to the closest 

108 point on the great circle line defined by a start and an 

109 end point. 

110 

111 That is, if a perpendicular is drawn from this point to the 

112 great circle line, the along-track distance is the distance 

113 from the start point to the point where the perpendicular 

114 crosses the line. 

115 

116 @arg start: Start point of the great circle line (L{LatLon}). 

117 @arg end: End point of the great circle line (L{LatLon}). 

118 @kwarg radius: Mean earth radius (C{meter}) or C{None}. 

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

120 the B{C{start}} and B{C{end}} point (C{bool}). 

121 

122 @return: Distance along the great circle line (C{radians} 

123 if C{B{radius} is None} or C{meter}, same units 

124 as B{C{radius}}), positive if I{after} the 

125 B{C{start}} toward the B{C{end}} point of the 

126 line, I{negative} if before or C{0} if at the 

127 B{C{start}} point. 

128 

129 @raise TypeError: Invalid B{C{start}} or B{C{end}} point. 

130 

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

132 

133 @example: 

134 

135 >>> p = LatLon(53.2611, -0.7972) 

136 

137 >>> s = LatLon(53.3206, -1.7297) 

138 >>> e = LatLon(53.1887, 0.1334) 

139 >>> d = p.alongTrackDistanceTo(s, e) # 62331.58 

140 ''' 

141 r, x, b = self._a_x_b3(start, end, radius, wrap) 

142 cx = cos(x) 

143 return _0_0 if isnear0(cx) else \ 

144 _radians2m(copysign0(acos1(cos(r) / cx), cos(b)), radius) 

145 

146 def _a_x_b3(self, start, end, radius, wrap): 

147 '''(INTERNAL) Helper for .along-/crossTrackDistanceTo. 

148 ''' 

149 s = self.others(start=start) 

150 e = self.others(end=end) 

151 s, e, w = _unrollon3(self, s, e, wrap) 

152 

153 r = Radius_(radius) 

154 r = s.distanceTo(self, r, wrap=w) / r 

155 

156 b = radians(s.initialBearingTo(self, wrap=w) 

157 - s.initialBearingTo(e, wrap=w)) 

158 x = asin(sin(r) * sin(b)) 

159 return r, x, -b 

160 

161 @deprecated_method 

162 def bearingTo(self, other, wrap=False, raiser=False): # PYCHOK no cover 

163 '''DEPRECATED, use method L{initialBearingTo}. 

164 ''' 

165 return self.initialBearingTo(other, wrap=wrap, raiser=raiser) 

166 

167 def crossingParallels(self, other, lat, wrap=False): 

168 '''Return the pair of meridians at which a great circle defined 

169 by this and an other point crosses the given latitude. 

170 

171 @arg other: The other point defining great circle (L{LatLon}). 

172 @arg lat: Latitude at the crossing (C{degrees}). 

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

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

175 

176 @return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or 

177 C{None} if the great circle doesn't reach B{C{lat}}. 

178 ''' 

179 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap) 

180 sa, ca, sa1, ca1, \ 

181 sa2, ca2, sdb, cdb = sincos2_(radians(lat), a1, a2, db) 

182 sa1 *= ca2 * ca 

183 

184 x = sa1 * sdb 

185 y = sa1 * cdb - ca1 * sa2 * ca 

186 z = ca1 * sdb * ca2 * sa 

187 

188 h = hypot(x, y) 

189 if h < EPS or fabs(z) > h: # PYCHOK no cover 

190 return None # great circle doesn't reach latitude 

191 

192 m = atan2(-y, x) + b1 # longitude at max latitude 

193 d = acos1(z / h) # delta longitude to intersections 

194 return degrees180(m - d), degrees180(m + d) 

195 

196 def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False): 

197 '''Compute the (signed) distance from this point to 

198 the great circle defined by a start and an end point. 

199 

200 @arg start: Start point of the great circle line (L{LatLon}). 

201 @arg end: End point of the great circle line (L{LatLon}). 

202 @kwarg radius: Mean earth radius (C{meter}) or C{None}. 

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

204 the B{C{start}} and B{C{end}} point (C{bool}). 

205 

206 @return: Distance to the great circle (C{radians} if 

207 B{C{radius}} or C{meter}, same units as 

208 B{C{radius}}), I{negative} if to the left or 

209 I{positive} if to the right of the line. 

210 

211 @raise TypeError: If B{C{start}} or B{C{end}} is not L{LatLon}. 

212 

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

214 

215 @example: 

216 

217 >>> p = LatLon(53.2611, -0.7972) 

218 

219 >>> s = LatLon(53.3206, -1.7297) 

220 >>> e = LatLon(53.1887, 0.1334) 

221 >>> d = p.crossTrackDistanceTo(s, e) # -307.5 

222 ''' 

223 _, x, _ = self._a_x_b3(start, end, radius, wrap) 

224 return _radians2m(x, radius) 

225 

226 def destination(self, distance, bearing, radius=R_M, height=None): 

227 '''Locate the destination from this point after having 

228 travelled the given distance on the given initial bearing. 

229 

230 @arg distance: Distance travelled (C{meter}, same units as 

231 B{C{radius}}). 

232 @arg bearing: Bearing from this point (compass C{degrees360}). 

233 @kwarg radius: Mean earth radius (C{meter}). 

234 @kwarg height: Optional height at destination (C{meter}, same 

235 units a B{C{radius}}). 

236 

237 @return: Destination point (L{LatLon}). 

238 

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

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

241 

242 @example: 

243 

244 >>> p1 = LatLon(51.4778, -0.0015) 

245 >>> p2 = p1.destination(7794, 300.7) 

246 >>> p2.toStr() # '51.5135°N, 000.0983°W' 

247 ''' 

248 a, b = self.philam 

249 r, t = _m2radians(distance, radius, low=None), Bearing_(bearing) 

250 

251 a, b = _destination2(a, b, r, t) 

252 h = self._heigHt(height) 

253 return self.classof(degrees90(a), degrees180(b), height=h) 

254 

255 def distanceTo(self, other, radius=R_M, wrap=False): 

256 '''Compute the (angular) distance from this to an other point. 

257 

258 @arg other: The other point (L{LatLon}). 

259 @kwarg radius: Mean earth radius (C{meter}) or C{None}. 

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

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

262 

263 @return: Distance between this and the B{C{other}} point 

264 (C{meter}, same units as B{C{radius}} or 

265 C{radians} if B{C{radius}} is C{None}). 

266 

267 @raise TypeError: The B{C{other}} point is not L{LatLon}. 

268 

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

270 

271 @example: 

272 

273 >>> p1 = LatLon(52.205, 0.119) 

274 >>> p2 = LatLon(48.857, 2.351); 

275 >>> d = p1.distanceTo(p2) # 404300 

276 ''' 

277 a1, _, a2, _, db = self._ab1_ab2_db5(other, wrap) 

278 return _radians2m(vincentys_(a2, a1, db), radius) 

279 

280# @Property_RO 

281# def Ecef(self): 

282# '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}. 

283# ''' 

284# return _MODS.ecef.EcefKarney 

285 

286 def greatCircle(self, bearing, Vector=Vector3d, **Vector_kwds): 

287 '''Compute the vector normal to great circle obtained by heading 

288 on the given initial bearing from this point. 

289 

290 Direction of vector is such that initial bearing vector 

291 b = c × n, where n is an n-vector representing this point. 

292 

293 @arg bearing: Bearing from this point (compass C{degrees360}). 

294 @kwarg Vector: Vector class to return the great circle, 

295 overriding the default L{Vector3d}. 

296 @kwarg Vector_kwds: Optional, additional keyword argunents 

297 for B{C{Vector}}. 

298 

299 @return: Vector representing great circle (C{Vector}). 

300 

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

302 

303 @example: 

304 

305 >>> p = LatLon(53.3206, -1.7297) 

306 >>> g = p.greatCircle(96.0) 

307 >>> g.toStr() # (-0.794, 0.129, 0.594) 

308 ''' 

309 a, b = self.philam 

310 sa, ca, sb, cb, st, ct = sincos2_(a, b, Bearing_(bearing)) 

311 

312 return Vector(sb * ct - cb * sa * st, 

313 -cb * ct - sb * sa * st, 

314 ca * st, **Vector_kwds) # XXX .unit()? 

315 

316 def initialBearingTo(self, other, wrap=False, raiser=False): 

317 '''Compute the initial bearing (forward azimuth) from this 

318 to an other point. 

319 

320 @arg other: The other point (spherical L{LatLon}). 

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

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

323 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}), 

324 use C{B{raiser}=True} for behavior like 

325 C{sphericalNvector.LatLon.initialBearingTo}. 

326 

327 @return: Initial bearing (compass C{degrees360}). 

328 

329 @raise CrossError: If this and the B{C{other}} point coincide, 

330 provided both B{C{raiser}} is C{True} and 

331 L{pygeodesy.crosserrors} is C{True}. 

332 

333 @raise TypeError: The B{C{other}} point is not L{LatLon}. 

334 

335 @example: 

336 

337 >>> p1 = LatLon(52.205, 0.119) 

338 >>> p2 = LatLon(48.857, 2.351) 

339 >>> b = p1.initialBearingTo(p2) # 156.2 

340 ''' 

341 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap) 

342 # XXX behavior like sphericalNvector.LatLon.initialBearingTo 

343 if raiser and crosserrors() and max(fabs(a2 - a1), fabs(db)) < EPS: 

344 raise CrossError(_point_, self, other=other, wrap=wrap, txt=_coincident_) 

345 

346 return degrees(bearing_(a1, b1, a2, b2, final=False)) 

347 

348 def intermediateTo(self, other, fraction, height=None, wrap=False): 

349 '''Locate the point at given fraction between (or along) this 

350 and an other point. 

351 

352 @arg other: The other point (L{LatLon}). 

353 @arg fraction: Fraction between both points (C{scalar}, 

354 0.0 at this and 1.0 at the other point). 

355 @kwarg height: Optional height, overriding the intermediate 

356 height (C{meter}). 

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

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

359 

360 @return: Intermediate point (L{LatLon}). 

361 

362 @raise TypeError: The B{C{other}} point is not L{LatLon}. 

363 

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

365 

366 @see: Methods C{midpointTo} and C{rhumbMidpointTo}. 

367 

368 @example: 

369 

370 >>> p1 = LatLon(52.205, 0.119) 

371 >>> p2 = LatLon(48.857, 2.351) 

372 >>> p = p1.intermediateTo(p2, 0.25) # 51.3721°N, 000.7073°E 

373 ''' 

374 p = self 

375 f = Scalar(fraction=fraction) 

376 if not isnear0(f): 

377 p = p.others(other) 

378 if wrap: 

379 p = _Wrap.point(p) 

380 if not isnear1(f): # and not near0 

381 a1, b1 = self.philam 

382 a2, b2 = p.philam 

383 db, b2 = unrollPI(b1, b2, wrap=wrap) 

384 r = vincentys_(a2, a1, db) 

385 sr = sin(r) 

386 if isnon0(sr): 

387 sa1, ca1, sa2, ca2, \ 

388 sb1, cb1, sb2, cb2 = sincos2_(a1, a2, b1, b2) 

389 

390 t = f * r 

391 a = sin(r - t) # / sr superflous 

392 b = sin( t) # / sr superflous 

393 

394 x = a * ca1 * cb1 + b * ca2 * cb2 

395 y = a * ca1 * sb1 + b * ca2 * sb2 

396 z = a * sa1 + b * sa2 

397 

398 a = atan1d(z, hypot(x, y)) 

399 b = atan2d(y, x) 

400 

401 else: # PYCHOK no cover 

402 a = degrees90( favg(a1, a2, f=f)) # coincident 

403 b = degrees180(favg(b1, b2, f=f)) 

404 

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

406 p = self.classof(a, b, height=h) 

407 return p 

408 

409 def intersection(self, end1, other, end2, height=None, wrap=False): 

410 '''Compute the intersection point of two lines, each defined by 

411 two points or a start point and bearing from North. 

412 

413 @arg end1: End point of this line (L{LatLon}) or the initial 

414 bearing at this point (compass C{degrees360}). 

415 @arg other: Start point of the other line (L{LatLon}). 

416 @arg end2: End point of the other line (L{LatLon}) or the 

417 initial bearing at the B{C{other}} point (compass 

418 C{degrees360}). 

419 @kwarg height: Optional height for intersection point, 

420 overriding the mean height (C{meter}). 

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

422 B{C{start2}} and both B{C{end*}} points (C{bool}). 

423 

424 @return: The intersection point (L{LatLon}). An alternate 

425 intersection point might be the L{antipode} to 

426 the returned result. 

427 

428 @raise IntersectionError: Ambiguous or infinite intersection 

429 or colinear, parallel or otherwise 

430 non-intersecting lines. 

431 

432 @raise TypeError: If B{C{other}} is not L{LatLon} or B{C{end1}} 

433 or B{C{end2}} not C{scalar} nor L{LatLon}. 

434 

435 @raise ValueError: Invalid B{C{height}} or C{null} line. 

436 

437 @example: 

438 

439 >>> p = LatLon(51.8853, 0.2545) 

440 >>> s = LatLon(49.0034, 2.5735) 

441 >>> i = p.intersection(108.547, s, 32.435) # '50.9078°N, 004.5084°E' 

442 ''' 

443 try: 

444 s2 = self.others(other) 

445 return _intersect(self, end1, s2, end2, height=height, wrap=wrap, 

446 LatLon=self.classof) 

447 except (TypeError, ValueError) as x: 

448 raise _xError(x, start1=self, end1=end1, 

449 other=other, end2=end2, wrap=wrap) 

450 

451 def intersections2(self, rad1, other, rad2, radius=R_M, eps=_0_0, 

452 height=None, wrap=True): 

453 '''Compute the intersection points of two circles, each defined 

454 by a center point and radius. 

455 

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

457 see B{C{radius}}). 

458 @arg other: Center point of the other circle (L{LatLon}). 

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

460 see B{C{radius}}). 

461 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}}, 

462 B{C{rad2}} and B{C{eps}} are given in C{radians}). 

463 @kwarg eps: Required overlap (C{meter} or C{radians}, see 

464 B{C{radius}}). 

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

466 conventionally) or C{None} for the I{"radical height"} 

467 at the I{radical line} between both centers. 

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

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

470 

471 @return: 2-Tuple of the intersection points, each a L{LatLon} 

472 instance. For abutting circles, both intersection 

473 points are the same instance, aka the I{radical center}. 

474 

475 @raise IntersectionError: Concentric, antipodal, invalid or 

476 non-intersecting circles. 

477 

478 @raise TypeError: If B{C{other}} is not L{LatLon}. 

479 

480 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}}, 

481 B{C{eps}} or B{C{height}}. 

482 ''' 

483 try: 

484 c2 = self.others(other) 

485 return _intersects2(self, rad1, c2, rad2, radius=radius, eps=eps, 

486 height=height, wrap=wrap, 

487 LatLon=self.classof) 

488 except (TypeError, ValueError) as x: 

489 raise _xError(x, center=self, rad1=rad1, 

490 other=other, rad2=rad2, wrap=wrap) 

491 

492 @deprecated_method 

493 def isEnclosedBy(self, points): # PYCHOK no cover 

494 '''DEPRECATED, use method C{isenclosedBy}.''' 

495 return self.isenclosedBy(points) 

496 

497 def isenclosedBy(self, points, wrap=False): 

498 '''Check whether a (convex) polygon or composite encloses this point. 

499 

500 @arg points: The polygon points or composite (L{LatLon}[], 

501 L{BooleanFHP} or L{BooleanGH}). 

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

503 B{C{points}} (C{bool}). 

504 

505 @return: C{True} if this point is inside the polygon or 

506 composite, C{False} otherwise. 

507 

508 @raise PointsError: Insufficient number of B{C{points}}. 

509 

510 @raise TypeError: Some B{C{points}} are not L{LatLon}. 

511 

512 @raise ValueError: Invalid B{C{points}}, non-convex polygon. 

513 

514 @see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy} 

515 and L{pygeodesy.ispolar} especially if the B{C{points}} may 

516 enclose a pole or wrap around the earth I{longitudinally}. 

517 ''' 

518 if _MODS.booleans.isBoolean(points): 

519 return points._encloses(self.lat, self.lon, wrap=wrap) 

520 

521 Ps = self.PointsIter(points, loop=2, dedup=True, wrap=wrap) 

522 n0 = self._N_vector 

523 

524 v2 = Ps[0]._N_vector 

525 p1 = Ps[1] 

526 v1 = p1._N_vector 

527 # check whether this point on same side of all 

528 # polygon edges (to the left or right depending 

529 # on the anti-/clockwise polygon direction) 

530 gc1 = v2.cross(v1) 

531 t0 = gc1.angleTo(n0) > PI_2 

532 s0 = None 

533 # get great-circle vector for each edge 

534 for i, p2 in Ps.enumerate(closed=True): 

535 if wrap and not Ps.looped: 

536 p2 = _unrollon(p1, p2) 

537 p1 = p2 

538 v2 = p2._N_vector 

539 gc = v1.cross(v2) 

540 t = gc.angleTo(n0) > PI_2 

541 if t != t0: # different sides of edge i 

542 return False # outside 

543 

544 # check for convex polygon: angle between 

545 # gc vectors, signed by direction of n0 

546 # (otherwise the test above is not reliable) 

547 s = signOf(gc1.angleTo(gc, vSign=n0)) 

548 if s != s0: 

549 if s0 is None: 

550 s0 = s 

551 else: 

552 t = _Fmt.SQUARE(points=i) 

553 raise _ValueError(t, p2, wrap=wrap, txt=_not_(_convex_)) 

554 gc1, v1 = gc, v2 

555 

556 return True # inside 

557 

558 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False): 

559 '''Find the midpoint between this and an other point. 

560 

561 @arg other: The other point (L{LatLon}). 

562 @kwarg height: Optional height for midpoint, overriding 

563 the mean height (C{meter}). 

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

565 may be negative or greater than 1.0. 

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

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

568 

569 @return: Midpoint (L{LatLon}). 

570 

571 @raise TypeError: The B{C{other}} point is not L{LatLon}. 

572 

573 @raise ValueError: Invalid B{C{height}}. 

574 

575 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}. 

576 

577 @example: 

578 

579 >>> p1 = LatLon(52.205, 0.119) 

580 >>> p2 = LatLon(48.857, 2.351) 

581 >>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E' 

582 ''' 

583 if fraction is _0_5: 

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

585 a1, b, a2, _, db = self._ab1_ab2_db5(other, wrap) 

586 sa1, ca1, sa2, ca2, sdb, cdb = sincos2_(a1, a2, db) 

587 

588 x = ca2 * cdb + ca1 

589 y = ca2 * sdb 

590 

591 a = atan1d(sa1 + sa2, hypot(x, y)) 

592 b = degrees180(b + atan2(y, x)) 

593 

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

595 r = self.classof(a, b, height=h) 

596 else: 

597 r = self.intermediateTo(other, fraction, height=height, wrap=wrap) 

598 return r 

599 

600 def nearestOn(self, point1, point2, radius=R_M, **wrap_adjust_limit): 

601 '''Locate the point between two points closest to this point. 

602 

603 Distances are approximated by function L{pygeodesy.equirectangular_}, 

604 subject to the supplied B{C{options}}. 

605 

606 @arg point1: Start point (L{LatLon}). 

607 @arg point2: End point (L{LatLon}). 

608 @kwarg radius: Mean earth radius (C{meter}). 

609 @kwarg wrap_adjust_limit: Optional keyword arguments for functions 

610 L{sphericalTrigonometry.nearestOn3} and 

611 L{pygeodesy.equirectangular_}, 

612 

613 @return: Closest point on the great circle line (L{LatLon}). 

614 

615 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}}, 

616 see function L{pygeodesy.equirectangular_}. 

617 

618 @raise NotImplementedError: Keyword argument C{B{within}=False} 

619 is not (yet) supported. 

620 

621 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}. 

622 

623 @raise ValueError: Invalid B{C{radius}} or B{C{options}}. 

624 

625 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5} 

626 and method L{sphericalTrigonometry.LatLon.nearestOn3}. 

627 ''' 

628 # remove kwarg B{C{within}} if present 

629 w = _xkwds_pop(wrap_adjust_limit, within=True) 

630 if not w: 

631 notImplemented(self, within=w) 

632 

633# # UNTESTED - handle C{B{within}=False} and C{B{within}=True} 

634# wrap = _xkwds_get(options, wrap=False) 

635# a = self.alongTrackDistanceTo(point1, point2, radius=radius, wrap=wrap) 

636# if fabs(a) < EPS or (within and a < EPS): 

637# return point1 

638# d = point1.distanceTo(point2, radius=radius, wrap=wrap) 

639# if isnear0(d): 

640# return point1 # or point2 

641# elif fabs(d - a) < EPS or (a + EPS) > d: 

642# return point2 

643# f = a / d 

644# if within: 

645# if f > EPS1: 

646# return point2 

647# elif f < EPS: 

648# return point1 

649# return point1.intermediateTo(point2, f, wrap=wrap) 

650 

651 # without kwarg B{C{within}}, use backward compatible .nearestOn3 

652 return self.nearestOn3([point1, point2], closed=False, radius=radius, 

653 **wrap_adjust_limit)[0] 

654 

655 @deprecated_method 

656 def nearestOn2(self, points, closed=False, radius=R_M, **options): # PYCHOK no cover 

657 '''DEPRECATED, use method L{sphericalTrigonometry.LatLon.nearestOn3}. 

658 

659 @return: ... 2-Tuple C{(closest, distance)} of the closest 

660 point (L{LatLon}) on the polygon and the distance 

661 to that point from this point in C{meter}, same 

662 units of B{C{radius}}. 

663 ''' 

664 r = self.nearestOn3(points, closed=closed, radius=radius, **options) 

665 return r.closest, r.distance 

666 

667 def nearestOn3(self, points, closed=False, radius=R_M, **wrap_adjust_limit): 

668 '''Locate the point on a polygon closest to this point. 

669 

670 Distances are approximated by function L{pygeodesy.equirectangular_}, 

671 subject to the supplied B{C{options}}. 

672 

673 @arg points: The polygon points (L{LatLon}[]). 

674 @kwarg closed: Optionally, close the polygon (C{bool}). 

675 @kwarg radius: Mean earth radius (C{meter}). 

676 @kwarg wrap_adjust_limit: Optional keyword arguments for function 

677 L{sphericalTrigonometry.nearestOn3} and 

678 L{pygeodesy.equirectangular_}, 

679 

680 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of the 

681 C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular_} 

682 C{distance} between this and the C{closest} point converted to 

683 C{meter}, same units as B{C{radius}}. The C{angle} from this 

684 to the C{closest} point is in compass C{degrees360}, like 

685 function L{pygeodesy.compassAngle}. 

686 

687 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}}, 

688 see function L{pygeodesy.equirectangular_}. 

689 

690 @raise PointsError: Insufficient number of B{C{points}}. 

691 

692 @raise TypeError: Some B{C{points}} are not C{LatLon}. 

693 

694 @raise ValueError: Invalid B{C{radius}} or B{C{options}}. 

695 

696 @see: Functions L{pygeodesy.compassAngle}, L{pygeodesy.equirectangular_} 

697 and L{pygeodesy.nearestOn5}. 

698 ''' 

699 return nearestOn3(self, points, closed=closed, radius=radius, 

700 LatLon=self.classof, **wrap_adjust_limit) 

701 

702 def toCartesian(self, **Cartesian_datum_kwds): # PYCHOK Cartesian=Cartesian, datum=None 

703 '''Convert this point to C{Karney}-based cartesian (ECEF) 

704 coordinates. 

705 

706 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}} 

707 and other keyword arguments, ignored 

708 if C{B{Cartesian} is None}. Use 

709 C{B{Cartesian}=...} to override 

710 this L{Cartesian} class or specify 

711 C{B{Cartesian}=None}. 

712 

713 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}} 

714 is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, 

715 C, M, datum)} with C{C} and C{M} if available. 

716 

717 @raise TypeError: Invalid B{C{Cartesian_datum_kwds}} argument. 

718 ''' 

719 kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian, datum=self.datum) 

720 return LatLonSphericalBase.toCartesian(self, **kwds) 

721 

722 def triangle7(self, otherB, otherC, radius=R_M, wrap=False): 

723 '''Compute the angles, sides and area of a spherical triangle. 

724 

725 @arg otherB: Second triangle point (C{LatLon}). 

726 @arg otherC: Third triangle point (C{LatLon}). 

727 @kwarg radius: Mean earth radius, ellipsoid or datum 

728 (C{meter}, L{Ellipsoid}, L{Ellipsoid2}, 

729 L{Datum} or L{a_f2Tuple}) or C{None}. 

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

731 B{C{otherB}} and B{C{otherC}} points (C{bool}). 

732 

733 @return: L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} or if 

734 B{C{radius}} is C{None}, a L{Triangle8Tuple}C{(A, 

735 a, B, b, C, c, D, E)}. 

736 

737 @see: Function L{triangle7} and U{Spherical trigonometry 

738 <https://WikiPedia.org/wiki/Spherical_trigonometry>}. 

739 ''' 

740 B = self.others(otherB=otherB) 

741 C = self.others(otherC=otherC) 

742 B, C, _ = _unrollon3(self, B, C, wrap) 

743 

744 r = self.philam + B.philam + C.philam 

745 t = triangle8_(*r, wrap=wrap) 

746 return self._xnamed(_t7Tuple(t, radius)) 

747 

748 def trilaterate5(self, distance1, point2, distance2, point3, distance3, 

749 area=True, eps=EPS1, radius=R_M, wrap=False): 

750 '''Trilaterate three points by I{area overlap} or I{perimeter 

751 intersection} of three corresponding circles. 

752 

753 @arg distance1: Distance to this point (C{meter}, same units 

754 as B{C{radius}}). 

755 @arg point2: Second center point (C{LatLon}). 

756 @arg distance2: Distance to point2 (C{meter}, same units as 

757 B{C{radius}}). 

758 @arg point3: Third center point (C{LatLon}). 

759 @arg distance3: Distance to point3 (C{meter}, same units as 

760 B{C{radius}}). 

761 @kwarg area: If C{True} compute the area overlap, otherwise the 

762 perimeter intersection of the circles (C{bool}). 

763 @kwarg eps: The required I{minimal overlap} for C{B{area}=True} 

764 or the I{intersection margin} for C{B{area}=False} 

765 (C{meter}, same units as B{C{radius}}). 

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

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

768 B{C{point2}} and B{C{point3}} (C{bool}). 

769 

770 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)} 

771 with C{min} and C{max} in C{meter}, same units as B{C{eps}}, 

772 the corresponding trilaterated points C{minPoint} and 

773 C{maxPoint} as I{spherical} C{LatLon} and C{n}, the number 

774 of trilatered points found for the given B{C{eps}}. 

775 

776 If only a single trilaterated point is found, C{min I{is} 

777 max}, C{minPoint I{is} maxPoint} and C{n = 1}. 

778 

779 For C{B{area}=True}, C{min} and C{max} are the smallest 

780 respectively largest I{radial} overlap found. 

781 

782 For C{B{area}=False}, C{min} and C{max} represent the 

783 nearest respectively farthest intersection margin. 

784 

785 If C{B{area}=True} and all 3 circles are concentric, C{n = 

786 0} and C{minPoint} and C{maxPoint} are both the B{C{point#}} 

787 with the smallest B{C{distance#}} C{min} and C{max} the 

788 largest B{C{distance#}}. 

789 

790 @raise IntersectionError: Trilateration failed for the given B{C{eps}}, 

791 insufficient overlap for C{B{area}=True} or 

792 no intersection or all (near-)concentric for 

793 C{B{area}=False}. 

794 

795 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}. 

796 

797 @raise ValueError: Coincident B{C{point2}} or B{C{point3}} or invalid 

798 B{C{distance1}}, B{C{distance2}}, B{C{distance3}} 

799 or B{C{radius}}. 

800 ''' 

801 return _trilaterate5(self, distance1, 

802 self.others(point2=point2), distance2, 

803 self.others(point3=point3), distance3, 

804 area=area, radius=radius, eps=eps, wrap=wrap) 

805 

806 

807_T00 = LatLon(0, 0, name='T00') # reference instance (L{LatLon}) 

808 

809 

810def areaOf(points, radius=R_M, wrap=False): # was=True 

811 '''Calculate the area of a (spherical) polygon or composite 

812 (with the pointsjoined by great circle arcs). 

813 

814 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP} 

815 or L{BooleanGH}). 

816 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter}, 

817 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}) 

818 or C{None}. 

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

820 (C{bool}). 

821 

822 @return: Polygon area (C{meter} I{quared}, same units as B{C{radius}} 

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

824 

825 @raise PointsError: Insufficient number of B{C{points}}. 

826 

827 @raise TypeError: Some B{C{points}} are not L{LatLon}. 

828 

829 @raise ValueError: Invalid B{C{radius}} or semi-circular polygon edge. 

830 

831 @note: The area is based on I{Karney}'s U{'Area of a spherical 

832 polygon'<https://MathOverflow.net/questions/97711/ 

833 the-area-of-spherical-polygons>}, 3rd Answer. 

834 

835 @see: Functions L{pygeodesy.areaOf}, L{sphericalNvector.areaOf}, 

836 L{pygeodesy.excessKarney}, L{ellipsoidalExact.areaOf} and 

837 L{ellipsoidalKarney.areaOf}. 

838 

839 @example: 

840 

841 >>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1) 

842 >>> areaOf(b) # 8666058750.718977 

843 

844 >>> c = LatLon(0, 0), LatLon(1, 0), LatLon(0, 1) 

845 >>> areaOf(c) # 6.18e9 

846 ''' 

847 if _MODS.booleans.isBoolean(points): 

848 return points._sum2(LatLon, areaOf, radius=radius, wrap=wrap) 

849 

850 _at2, _t_2 = atan2, tan_2 

851 _un, _w180 = unrollPI, wrap180 

852 

853 Ps = _T00.PointsIter(points, loop=1, wrap=wrap) 

854 p1 = p2 = Ps[0] 

855 a1, b1 = p1.philam 

856 ta1, z1 = _t_2(a1), None 

857 

858 A = Fsum() # mean phi 

859 R = Fsum() # see L{pygeodesy.excessKarney_} 

860 # ispolar: Summation of course deltas around pole is 0° rather than normally ±360° 

861 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html> 

862 # XXX duplicate of function C{points.ispolar} to avoid copying all iterated points 

863 D = Fsum() 

864 for i, p2 in Ps.enumerate(closed=True): 

865 a2, b2 = p2.philam 

866 db, b2 = _un(b1, b2, wrap=wrap and not Ps.looped) 

867 A += a2 

868 ta2 = _t_2(a2) 

869 tdb = _t_2(db, points=i) 

870 R += _at2(tdb * (ta1 + ta2), 

871 _1_0 + ta1 * ta2) 

872 ta1, b1 = ta2, b2 

873 

874 if not p2.isequalTo(p1, eps=EPS): 

875 z, z2 = _bearingTo2(p1, p2, wrap=wrap) 

876 if z1 is not None: 

877 D += _w180(z - z1) # (z - z1 + 540) ... 

878 D += _w180(z2 - z) # (z2 - z + 540) % 360 - 180 

879 p1, z1 = p2, z2 

880 

881 R = abs(R * _2_0) 

882 if abs(D) < _90_0: # ispolar(points) 

883 R = abs(R - PI2) 

884 if radius: 

885 a = degrees(A.fover(len(A))) # mean lat 

886 R *= _mean_radius(radius, a)**2 

887 return float(R) 

888 

889 

890def _destination2(a, b, r, t): 

891 '''(INTERNAL) Destination lat- and longitude in C{radians}. 

892 

893 @arg a: Latitude (C{radians}). 

894 @arg b: Longitude (C{radians}). 

895 @arg r: Angular distance (C{radians}). 

896 @arg t: Bearing (compass C{radians}). 

897 

898 @return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}). 

899 ''' 

900 # see <https://www.EdWilliams.org/avform.htm#LL> 

901 sa, ca, sr, cr, st, ct = sincos2_(a, r, t) 

902 ca *= sr 

903 

904 a = asin1(ct * ca + cr * sa) 

905 d = atan2(st * ca, cr - sa * sin(a)) 

906 # note, in EdWilliams.org/avform.htm W is + and E is - 

907 return a, (b + d) # (mod(b + d + PI, PI2) - PI) 

908 

909 

910def _int3d2(s, end, wrap, _i_, Vector, hs): 

911 # see <https://www.EdWilliams.org/intersect.htm> (5) ff 

912 # and similar logic in .ellipsoidalBaseDI._intersect3 

913 a1, b1 = s.philam 

914 

915 if isscalar(end): # bearing, get pseudo-end point 

916 a2, b2 = _destination2(a1, b1, PI_4, radians(end)) 

917 else: # must be a point 

918 s.others(end, name=_end_ + _i_) 

919 hs.append(end.height) 

920 a2, b2 = end.philam 

921 if wrap: 

922 a2, b2 = _Wrap.philam(a2, b2) 

923 

924 db, b2 = unrollPI(b1, b2, wrap=wrap) 

925 if max(fabs(db), fabs(a2 - a1)) < EPS: 

926 raise _ValueError(_SPACE_(_line_ + _i_, _null_)) 

927 # note, in EdWilliams.org/avform.htm W is + and E is - 

928 sb21, cb21, sb12, cb12 = sincos2_(db * _0_5, 

929 -(b1 + b2) * _0_5) 

930 cb21 *= sin(a1 - a2) # sa21 

931 sb21 *= sin(a1 + a2) # sa12 

932 x = Vector(sb12 * cb21 - cb12 * sb21, 

933 cb12 * cb21 + sb12 * sb21, 

934 cos(a1) * cos(a2) * sin(db)) # ll=start 

935 return x.unit(), (db, (a2 - a1)) # negated d 

936 

937 

938def _intdot(ds, a1, b1, a, b, wrap): 

939 # compute dot product ds . (-b + b1, a - a1) 

940 db, _ = unrollPI(b1, b, wrap=wrap) 

941 return fdot(ds, db, a - a1) 

942 

943 

944def intersecant2(center, circle, point, other, **radius_exact_height_wrap): 

945 '''Compute the intersections of a circle and a (great circle) line given as 

946 two points or as a point and bearing. 

947 

948 @arg center: Center of the circle (L{LatLon}). 

949 @arg circle: Radius of the circle (C{meter}, same units as B{C{radius}}) 

950 or a point on the circle (L{LatLon}). 

951 @arg point: A point on the (great circle) line (L{LatLon}). 

952 @arg other: An other point on the (great circle) line (L{LatLon}) or 

953 the bearing at the B{C{point}} (compass C{degrees360}). 

954 @kwarg radius_exact_height_wrap: Optional keyword arguments, see 

955 method L{LatLon.intersecant2} for further details. 

956 

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

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

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

960 

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

962 

963 @raise TypeError: If B{C{center}} or B{C{point}} not L{LatLon} or 

964 B{C{circle}} or B{C{other}} invalid. 

965 

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

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

968 ''' 

969 c = _T00.others(center=center) 

970 p = _T00.others(point=point) 

971 try: 

972 return _intersecant2(c, circle, p, other, **radius_exact_height_wrap) 

973 except (TypeError, ValueError) as x: 

974 raise _xError(x, center=center, circle=circle, point=point, other=other, 

975 **radius_exact_height_wrap) 

976 

977 

978def _intersect(start1, end1, start2, end2, height=None, wrap=False, # in.ellipsoidalBaseDI._intersect3 

979 LatLon=None, **LatLon_kwds): 

980 # (INTERNAL) Intersect two (spherical) lines, see L{intersection} 

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

982 

983 s1, s2 = start1, start2 

984 if wrap: 

985 s2 = _Wrap.point(s2) 

986 hs = [s1.height, s2.height] 

987 

988 a1, b1 = s1.philam 

989 a2, b2 = s2.philam 

990 db, b2 = unrollPI(b1, b2, wrap=wrap) 

991 r12 = vincentys_(a2, a1, db) 

992 if fabs(r12) < EPS: # [nearly] coincident points 

993 a, b = favg(a1, a2), favg(b1, b2) 

994 

995 # see <https://www.EdWilliams.org/avform.htm#Intersection> 

996 elif isscalar(end1) and isscalar(end2): # both bearings 

997 sa1, ca1, sa2, ca2, sr12, cr12 = sincos2_(a1, a2, r12) 

998 

999 x1, x2 = (sr12 * ca1), (sr12 * ca2) 

1000 if isnear0(x1) or isnear0(x2): 

1001 raise IntersectionError(_parallel_) 

1002 # handle domain error for equivalent longitudes, 

1003 # see also functions asin_safe and acos_safe at 

1004 # <https://www.EdWilliams.org/avform.htm#Math> 

1005 t12, t13 = acos1((sa2 - sa1 * cr12) / x1), radiansPI2(end1) 

1006 t21, t23 = acos1((sa1 - sa2 * cr12) / x2), radiansPI2(end2) 

1007 if sin(db) > 0: 

1008 t21 = PI2 - t21 

1009 else: 

1010 t12 = PI2 - t12 

1011 sx1, cx1, sx2, cx2 = sincos2_(wrapPI(t13 - t12), # angle 2-1-3 

1012 wrapPI(t21 - t23)) # angle 1-2-3) 

1013 if isnear0(sx1) and isnear0(sx2): 

1014 raise IntersectionError(_infinite_) 

1015 sx3 = sx1 * sx2 

1016# XXX if sx3 < 0: 

1017# XXX raise ValueError(_ambiguous_) 

1018 x3 = acos1(cr12 * sx3 - cx2 * cx1) 

1019 r13 = atan2(sr12 * sx3, cx2 + cx1 * cos(x3)) 

1020 

1021 a, b = _destination2(a1, b1, r13, t13) 

1022 # like .ellipsoidalBaseDI,_intersect3, if this intersection 

1023 # is "before" the first point, use the antipodal intersection 

1024 if opposing_(t13, bearing_(a1, b1, a, b, wrap=wrap)): 

1025 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple 

1026 

1027 else: # end point(s) or bearing(s) 

1028 _N_vector_ = _MODS.nvectorBase._N_vector_ 

1029 

1030 x1, d1 = _int3d2(s1, end1, wrap, _1_, _N_vector_, hs) 

1031 x2, d2 = _int3d2(s2, end2, wrap, _2_, _N_vector_, hs) 

1032 x = x1.cross(x2) 

1033 if x.length < EPS: # [nearly] colinear or parallel lines 

1034 raise IntersectionError(_colinear_) 

1035 a, b = x.philam 

1036 # choose intersection similar to sphericalNvector 

1037 if not (_intdot(d1, a1, b1, a, b, wrap) * 

1038 _intdot(d2, a2, b2, a, b, wrap)) > 0: 

1039 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple 

1040 

1041 h = fmean(hs) if height is None else Height(height) 

1042 return _LL3Tuple(degrees90(a), degrees180(b), h, 

1043 intersection, LatLon, LatLon_kwds) 

1044 

1045 

1046def intersection(start1, end1, start2, end2, height=None, wrap=False, 

1047 LatLon=LatLon, **LatLon_kwds): 

1048 '''Compute the intersection point of two lines, each defined 

1049 by two points or a start point and bearing from North. 

1050 

1051 @arg start1: Start point of the first line (L{LatLon}). 

1052 @arg end1: End point of the first line (L{LatLon}) or 

1053 the initial bearing at the first start point 

1054 (compass C{degrees360}). 

1055 @arg start2: Start point of the second line (L{LatLon}). 

1056 @arg end2: End point of the second line (L{LatLon}) or 

1057 the initial bearing at the second start point 

1058 (compass C{degrees360}). 

1059 @kwarg height: Optional height for the intersection point, 

1060 overriding the mean height (C{meter}). 

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

1062 B{C{start2}} and both B{C{end*}} points (C{bool}). 

1063 @kwarg LatLon: Optional class to return the intersection 

1064 point (L{LatLon}) or C{None}. 

1065 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword 

1066 arguments, ignored if C{B{LatLon} is None}. 

1067 

1068 @return: The intersection point as a (B{C{LatLon}}) or if 

1069 C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, lon, 

1070 height)}. An alternate intersection point might 

1071 be the L{antipode} to the returned result. 

1072 

1073 @raise IntersectionError: Ambiguous or infinite intersection 

1074 or colinear, parallel or otherwise 

1075 non-intersecting lines. 

1076 

1077 @raise TypeError: A B{C{start1}}, B{C{end1}}, B{C{start2}} 

1078 or B{C{end2}} point not L{LatLon}. 

1079 

1080 @raise ValueError: Invalid B{C{height}} or C{null} line. 

1081 

1082 @example: 

1083 

1084 >>> p = LatLon(51.8853, 0.2545) 

1085 >>> s = LatLon(49.0034, 2.5735) 

1086 >>> i = intersection(p, 108.547, s, 32.435) # '50.9078°N, 004.5084°E' 

1087 ''' 

1088 s1 = _T00.others(start1=start1) 

1089 s2 = _T00.others(start2=start2) 

1090 try: 

1091 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap, 

1092 LatLon=LatLon, **LatLon_kwds) 

1093 except (TypeError, ValueError) as x: 

1094 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2) 

1095 

1096 

1097def intersections2(center1, rad1, center2, rad2, radius=R_M, eps=_0_0, 

1098 height=None, wrap=False, # was=True 

1099 LatLon=LatLon, **LatLon_kwds): 

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

1101 by a center point and a radius. 

1102 

1103 @arg center1: Center of the first circle (L{LatLon}). 

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

1105 see B{C{radius}}). 

1106 @arg center2: Center of the second circle (L{LatLon}). 

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

1108 see B{C{radius}}). 

1109 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}}, 

1110 B{C{rad2}} and B{C{eps}} are given in C{radians}). 

1111 @kwarg eps: Required overlap (C{meter} or C{radians}, see 

1112 B{C{radius}}). 

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

1114 conventionally) or C{None} for the I{"radical height"} 

1115 at the I{radical line} between both centers. 

1116 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{center2}} 

1117 (C{bool}). 

1118 @kwarg LatLon: Optional class to return the intersection 

1119 points (L{LatLon}) or C{None}. 

1120 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword 

1121 arguments, ignored if C{B{LatLon} is None}. 

1122 

1123 @return: 2-Tuple of the intersection points, each a B{C{LatLon}} 

1124 instance or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, 

1125 lon, height)}. For abutting circles, both intersection 

1126 points are the same instance, aka the I{radical center}. 

1127 

1128 @raise IntersectionError: Concentric, antipodal, invalid or 

1129 non-intersecting circles. 

1130 

1131 @raise TypeError: If B{C{center1}} or B{C{center2}} not L{LatLon}. 

1132 

1133 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}}, 

1134 B{C{eps}} or B{C{height}}. 

1135 

1136 @note: Courtesy of U{Samuel Čavoj<https://GitHub.com/mrJean1/PyGeodesy/issues/41>}. 

1137 

1138 @see: This U{Answer<https://StackOverflow.com/questions/53324667/ 

1139 find-intersection-coordinates-of-two-circles-on-earth/53331953>}. 

1140 ''' 

1141 c1 = _T00.others(center1=center1) 

1142 c2 = _T00.others(center2=center2) 

1143 try: 

1144 return _intersects2(c1, rad1, c2, rad2, radius=radius, eps=eps, 

1145 height=height, wrap=wrap, 

1146 LatLon=LatLon, **LatLon_kwds) 

1147 except (TypeError, ValueError) as x: 

1148 raise _xError(x, center1=center1, rad1=rad1, 

1149 center2=center2, rad2=rad2, wrap=wrap) 

1150 

1151 

1152def _intersects2(c1, rad1, c2, rad2, radius=R_M, eps=_0_0, # in .ellipsoidalBaseDI._intersects2 

1153 height=None, too_d=None, wrap=False, # was=True 

1154 LatLon=LatLon, **LatLon_kwds): 

1155 # (INTERNAL) Intersect two spherical circles, see L{intersections2} 

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

1157 

1158 def _dest3(bearing, h): 

1159 a, b = _destination2(a1, b1, r1, bearing) 

1160 return _LL3Tuple(degrees90(a), degrees180(b), h, 

1161 intersections2, LatLon, LatLon_kwds) 

1162 

1163 a1, b1 = c1.philam 

1164 a2, b2 = c2.philam 

1165 if wrap: 

1166 a2, b2 = _Wrap.philam(a2, b2) 

1167 

1168 r1, r2, f = _rads3(rad1, rad2, radius) 

1169 if f: # swapped radii, swap centers 

1170 a1, a2 = a2, a1 # PYCHOK swap! 

1171 b1, b2 = b2, b1 # PYCHOK swap! 

1172 

1173 db, b2 = unrollPI(b1, b2, wrap=wrap) 

1174 d = vincentys_(a2, a1, db) # radians 

1175 if d < max(r1 - r2, EPS): 

1176 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError? 

1177 

1178 r = eps if radius is None else (m2radians( 

1179 eps, radius=radius) if eps else _0_0) 

1180 if r < _0_0: 

1181 raise _ValueError(eps=r) 

1182 

1183 x = fsumf_(r1, r2, -d) # overlap 

1184 if x > max(r, EPS): 

1185 sd, cd, sr1, cr1, _, cr2 = sincos2_(d, r1, r2) 

1186 x = sd * sr1 

1187 if isnear0(x): 

1188 raise _ValueError(_invalid_) 

1189 x = acos1((cr2 - cd * cr1) / x) # 0 <= x <= PI 

1190 

1191 elif x < r: # PYCHOK no cover 

1192 t = (d * radius) if too_d is None else too_d 

1193 raise IntersectionError(_too_(_Fmt.distant(t))) 

1194 

1195 if height is None: # "radical height" 

1196 f = _radical2(d, r1, r2).ratio 

1197 h = Height(favg(c1.height, c2.height, f=f)) 

1198 else: 

1199 h = Height(height) 

1200 

1201 b = bearing_(a1, b1, a2, b2, final=False, wrap=wrap) 

1202 if x < EPS4: # externally ... 

1203 r = _dest3(b, h) 

1204 elif x > _PI_EPS4: # internally ... 

1205 r = _dest3(b + PI, h) 

1206 else: 

1207 return _dest3(b + x, h), _dest3(b - x, h) 

1208 return r, r # ... abutting circles 

1209 

1210 

1211@deprecated_function 

1212def isPoleEnclosedBy(points, wrap=False): # PYCHOK no cover 

1213 '''DEPRECATED, use function L{pygeodesy.ispolar}. 

1214 ''' 

1215 return ispolar(points, wrap=wrap) 

1216 

1217 

1218def _LL3Tuple(lat, lon, height, where, LatLon, LatLon_kwds): 

1219 '''(INTERNAL) Helper for L{intersection}, L{intersections2} and L{meanOf}. 

1220 ''' 

1221 n = where.__name__ 

1222 if LatLon is None: 

1223 r = LatLon3Tuple(lat, lon, height, name=n) 

1224 else: 

1225 kwds = _xkwds(LatLon_kwds, height=height, name=n) 

1226 r = LatLon(lat, lon, **kwds) 

1227 return r 

1228 

1229 

1230def meanOf(points, height=None, wrap=False, LatLon=LatLon, **LatLon_kwds): 

1231 '''Compute the I{geographic} mean of several points. 

1232 

1233 @arg points: Points to be averaged (L{LatLon}[]). 

1234 @kwarg height: Optional height at mean point, overriding the mean 

1235 height (C{meter}). 

1236 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{points}} 

1237 (C{bool}). 

1238 @kwarg LatLon: Optional class to return the mean point (L{LatLon}) 

1239 or C{None}. 

1240 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword 

1241 arguments, ignored if C{B{LatLon} is None}. 

1242 

1243 @return: The geographic mean and height (B{C{LatLon}}) or a 

1244 L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}} 

1245 is C{None}. 

1246 

1247 @raise TypeError: Some B{C{points}} are not L{LatLon}. 

1248 

1249 @raise ValueError: No B{C{points}} or invalid B{C{height}}. 

1250 ''' 

1251 def _N_vs(ps, w): 

1252 Ps = _T00.PointsIter(ps, wrap=w) 

1253 for p in Ps.iterate(closed=False): 

1254 yield p._N_vector 

1255 

1256 m = _MODS.nvectorBase 

1257 # geographic, vectorial mean 

1258 n = m.sumOf(_N_vs(points, wrap), h=height, Vector=m.NvectorBase) 

1259 lat, lon, h = n.latlonheight 

1260 return _LL3Tuple(lat, lon, h, meanOf, LatLon, LatLon_kwds) 

1261 

1262 

1263@deprecated_function 

1264def nearestOn2(point, points, **closed_radius_LatLon_options): # PYCHOK no cover 

1265 '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}. 

1266 

1267 @return: ... 2-tuple C{(closest, distance)} of the C{closest} 

1268 point (L{LatLon}) on the polygon and the C{distance} 

1269 between the C{closest} and the given B{C{point}}. The 

1270 C{closest} is a B{C{LatLon}} or a L{LatLon2Tuple}C{(lat, 

1271 lon)} if B{C{LatLon}} is C{None} ... 

1272 ''' 

1273 ll, d, _ = nearestOn3(point, points, **closed_radius_LatLon_options) # PYCHOK 3-tuple 

1274 if _xkwds_get(closed_radius_LatLon_options, LatLon=LatLon) is None: 

1275 ll = LatLon2Tuple(ll.lat, ll.lon) 

1276 return ll, d 

1277 

1278 

1279def nearestOn3(point, points, closed=False, radius=R_M, wrap=False, adjust=True, 

1280 limit=9, **LatLon_and_kwds): 

1281 '''Locate the point on a path or polygon closest to a reference point. 

1282 

1283 Distances are I{approximated} using function L{pygeodesy.equirectangular_}, 

1284 subject to the supplied B{C{options}}. 

1285 

1286 @arg point: The reference point (L{LatLon}). 

1287 @arg points: The path or polygon points (L{LatLon}[]). 

1288 @kwarg closed: Optionally, close the polygon (C{bool}). 

1289 @kwarg radius: Mean earth radius (C{meter}). 

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

1291 B{C{points}} (C{bool}). 

1292 @kwarg adjust: See function L{pygeodesy.equirectangular_} (C{bool}). 

1293 @kwarg limit: See function L{pygeodesy.equirectangular_} (C{degrees}), 

1294 default C{9 degrees} is about C{1,000 Kmeter} (for mean 

1295 spherical earth radius L{R_KM}). 

1296 @kwarg LatLon: Optional class to return the closest point (L{LatLon}) 

1297 or C{None}. 

1298 @kwarg options: Optional keyword arguments for function 

1299 L{pygeodesy.equirectangular_}. 

1300 

1301 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the 

1302 C{closest} point as B{C{LatLon}} or L{LatLon3Tuple}C{(lat, 

1303 lon, height)} if B{C{LatLon}} is C{None}. The C{distance} 

1304 is the L{pygeodesy.equirectangular_} distance between the 

1305 C{closest} and the given B{C{point}} converted to C{meter}, 

1306 same units as B{C{radius}}. The C{angle} from the given 

1307 B{C{point}} to the C{closest} is in compass C{degrees360}, 

1308 like function L{pygeodesy.compassAngle}. The C{height} is 

1309 the (interpolated) height at the C{closest} point. 

1310 

1311 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}}, 

1312 see function L{pygeodesy.equirectangular_}. 

1313 

1314 @raise PointsError: Insufficient number of B{C{points}}. 

1315 

1316 @raise TypeError: Some B{C{points}} are not C{LatLon}. 

1317 

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

1319 

1320 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}. 

1321 ''' 

1322 t = _nearestOn5(point, points, closed=closed, wrap=wrap, 

1323 adjust=adjust, limit=limit) 

1324 d = degrees2m(t.distance, radius=radius) 

1325 h = t.height 

1326 n = nearestOn3.__name__ 

1327 

1328 kwds = _xkwds(LatLon_and_kwds, height=h, name=n) 

1329 LL = _xkwds_pop(kwds, LatLon=LatLon) 

1330 r = LatLon3Tuple(t.lat, t.lon, h, name=n) if LL is None else \ 

1331 LL(t.lat, t.lon, **kwds) 

1332 return NearestOn3Tuple(r, d, t.angle, name=n) 

1333 

1334 

1335def perimeterOf(points, closed=False, radius=R_M, wrap=True): 

1336 '''Compute the perimeter of a (spherical) polygon or composite 

1337 (with great circle arcs joining the points). 

1338 

1339 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP} 

1340 or L{BooleanGH}). 

1341 @kwarg closed: Optionally, close the polygon (C{bool}). 

1342 @kwarg radius: Mean earth radius (C{meter}) or C{None}. 

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

1344 B{C{points}} (C{bool}). 

1345 

1346 @return: Polygon perimeter (C{meter}, same units as B{C{radius}} 

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

1348 

1349 @raise PointsError: Insufficient number of B{C{points}}. 

1350 

1351 @raise TypeError: Some B{C{points}} are not L{LatLon}. 

1352 

1353 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with 

1354 C{B{points}} a composite. 

1355 

1356 @note: Distances are based on function L{pygeodesy.vincentys_}. 

1357 

1358 @see: Functions L{pygeodesy.perimeterOf}, L{sphericalNvector.perimeterOf} 

1359 and L{ellipsoidalKarney.perimeterOf}. 

1360 ''' 

1361 def _rads(ps, c, w): # angular edge lengths in radians 

1362 Ps = _T00.PointsIter(ps, loop=1, wrap=w) 

1363 a1, b1 = Ps[0].philam 

1364 for p in Ps.iterate(closed=c): 

1365 a2, b2 = p.philam 

1366 db, b2 = unrollPI(b1, b2, wrap=w and not (c and Ps.looped)) 

1367 yield vincentys_(a2, a1, db) 

1368 a1, b1 = a2, b2 

1369 

1370 if _MODS.booleans.isBoolean(points): 

1371 if not closed: 

1372 raise _ValueError(closed=closed, points=_composite_) 

1373 r = points._sum2(LatLon, perimeterOf, closed=True, radius=radius, wrap=wrap) 

1374 else: 

1375 r = fsum(_rads(points, closed, wrap), floats=True) 

1376 return _radians2m(r, radius) 

1377 

1378 

1379def triangle7(latA, lonA, latB, lonB, latC, lonC, radius=R_M, 

1380 excess=excessAbc_, 

1381 wrap=False): 

1382 '''Compute the angles, sides, and area of a (spherical) triangle. 

1383 

1384 @arg latA: First corner latitude (C{degrees}). 

1385 @arg lonA: First corner longitude (C{degrees}). 

1386 @arg latB: Second corner latitude (C{degrees}). 

1387 @arg lonB: Second corner longitude (C{degrees}). 

1388 @arg latC: Third corner latitude (C{degrees}). 

1389 @arg lonC: Third corner longitude (C{degrees}). 

1390 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter}, 

1391 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}) 

1392 or C{None}. 

1393 @kwarg excess: I{Spherical excess} callable (L{excessAbc_}, 

1394 L{excessGirard_} or L{excessLHuilier_}). 

1395 @kwarg wrap: If C{True}, wrap and L{pygeodesy.unroll180} 

1396 longitudes (C{bool}). 

1397 

1398 @return: A L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} with 

1399 spherical angles C{A}, C{B} and C{C}, angular sides 

1400 C{a}, C{b} and C{c} all in C{degrees} and C{area} 

1401 in I{square} C{meter} or same units as B{C{radius}} 

1402 I{squared} or if C{B{radius}=0} or C{None}, a 

1403 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} all in 

1404 C{radians} with the I{spherical excess} C{E} as the 

1405 C{unit area} in C{radians}. 

1406 ''' 

1407 t = triangle8_(Phi_(latA=latA), Lam_(lonA=lonA), 

1408 Phi_(latB=latB), Lam_(lonB=lonB), 

1409 Phi_(latC=latC), Lam_(lonC=lonC), 

1410 excess=excess, wrap=wrap) 

1411 return _t7Tuple(t, radius) 

1412 

1413 

1414def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc_, 

1415 wrap=False): 

1416 '''Compute the angles, sides, I{spherical deficit} and I{spherical 

1417 excess} of a (spherical) triangle. 

1418 

1419 @arg phiA: First corner latitude (C{radians}). 

1420 @arg lamA: First corner longitude (C{radians}). 

1421 @arg phiB: Second corner latitude (C{radians}). 

1422 @arg lamB: Second corner longitude (C{radians}). 

1423 @arg phiC: Third corner latitude (C{radians}). 

1424 @arg lamC: Third corner longitude (C{radians}). 

1425 @kwarg excess: I{Spherical excess} callable (L{excessAbc_}, 

1426 L{excessGirard_} or L{excessLHuilier_}). 

1427 @kwarg wrap: If C{True}, L{pygeodesy.unrollPI} the 

1428 longitudinal deltas (C{bool}). 

1429 

1430 @return: A L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with 

1431 spherical angles C{A}, C{B} and C{C}, angular sides 

1432 C{a}, C{b} and C{c}, I{spherical deficit} C{D} and 

1433 I{spherical excess} C{E}, all in C{radians}. 

1434 ''' 

1435 def _a_r(w, phiA, lamA, phiB, lamB, phiC, lamC): 

1436 d, _ = unrollPI(lamB, lamC, wrap=w) 

1437 a = vincentys_(phiC, phiB, d) 

1438 return a, (phiB, lamB, phiC, lamC, phiA, lamA) # rotate A, B, C 

1439 

1440 def _A_r(a, sa, ca, sb, cb, sc, cc): 

1441 s = sb * sc 

1442 A = acos1((ca - cb * cc) / s) if isnon0(s) else a 

1443 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2_'s 

1444 

1445 # notation: side C{a} is oposite to corner C{A}, etc. 

1446 a, r = _a_r(wrap, phiA, lamA, phiB, lamB, phiC, lamC) 

1447 b, r = _a_r(wrap, *r) 

1448 c, _ = _a_r(wrap, *r) 

1449 

1450 A, r = _A_r(a, *sincos2_(a, b, c)) 

1451 B, r = _A_r(b, *r) 

1452 C, _ = _A_r(c, *r) 

1453 

1454 D = fsumf_(PI2, -a, -b, -c) # deficit aka defect 

1455 E = excessGirard_(A, B, C) if excess in (excessGirard_, True) else ( 

1456 excessLHuilier_(a, b, c) if excess in (excessLHuilier_, False) else 

1457 excessAbc_(*max((A, b, c), (B, c, a), (C, a, b)))) 

1458 

1459 return Triangle8Tuple(A, a, B, b, C, c, D, E) 

1460 

1461 

1462def _t7Tuple(t, radius): 

1463 '''(INTERNAL) Convert a L{Triangle8Tuple} to L{Triangle7Tuple}. 

1464 ''' 

1465 if radius: # not in (None, _0_0) 

1466 r = radius if isscalar(radius) else \ 

1467 _ellipsoidal_datum(radius).ellipsoid.Rmean 

1468 A, B, C = map1(degrees, t.A, t.B, t.C) 

1469 t = Triangle7Tuple(A, (r * t.a), 

1470 B, (r * t.b), 

1471 C, (r * t.c), t.E * r**2) 

1472 return t 

1473 

1474 

1475__all__ += _ALL_OTHER(Cartesian, LatLon, # classes 

1476 areaOf, # functions 

1477 intersecant2, intersection, intersections2, ispolar, 

1478 isPoleEnclosedBy, # DEPRECATED, use ispolar 

1479 meanOf, 

1480 nearestOn2, nearestOn3, 

1481 perimeterOf, 

1482 sumOf, # XXX == vector3d.sumOf 

1483 triangle7, triangle8_) 

1484 

1485# **) MIT License 

1486# 

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

1488# 

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

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

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

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

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

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

1495# 

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

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

1498# 

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

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

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

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

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

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

1505# OTHER DEALINGS IN THE SOFTWARE.