Coverage for pygeodesy/sphericalTrigonometry.py: 94%

395 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-07-12 13:40 -0400

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_, \ 

34 _invalid_, _line_, _near_, _not_, _null_, \ 

35 _point_, _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, \ 

40 NearestOn3Tuple, Triangle7Tuple, \ 

41 Triangle8Tuple 

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

43 notImplemented, Fmt as _Fmt # XXX shadowed 

44from pygeodesy.props import deprecated_function, deprecated_method 

45from pygeodesy.sphericalBase import _angular, CartesianSphericalBase, _intersecant2, \ 

46 LatLonSphericalBase, _rads3, _trilaterate5 

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

48from pygeodesy.units import Bearing_, Height, Lam_, Phi_, Radius, \ 

49 Radius_, Scalar 

50from pygeodesy.utily import acos1, asin1, degrees90, degrees180, degrees2m, \ 

51 m2radians, radiansPI2, sincos2_, tan_2, _unrollon, \ 

52 unrollPI, _unrollon3, _Wrap, wrap180, wrapPI 

53from pygeodesy.vector3d import sumOf, Vector3d 

54 

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

56 

57__all__ = _ALL_LAZY.sphericalTrigonometry 

58__version__ = '23.06.12' 

59 

60_parallel_ = 'parallel' 

61 

62_PI_EPS4 = PI - EPS4 

63if _PI_EPS4 >= PI: 

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

65 

66 

67class Cartesian(CartesianSphericalBase): 

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

69 spherical, geodetic L{LatLon}. 

70 ''' 

71 

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

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

74 

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

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

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

78 

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

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

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

82 

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

84 ''' 

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

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

87 

88 

89class LatLon(LatLonSphericalBase): 

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

91 

92 @example: 

93 

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

95 ''' 

96 

97 def _ab1_ab2_db5(self, other, wrap): 

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

99 ''' 

100 a1, b1 = self.philam 

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

102 if wrap: 

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

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

105 else: # unrollPI shortcut 

106 db = b2 - b1 

107 return a1, b1, a2, b2, db 

108 

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

110 '''Compute the (angular) distance (signed) from the start to 

111 the closest point on the great circle line defined by a 

112 start and an end point. 

113 

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

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

116 from the start point to the point where the perpendicular 

117 crosses the line. 

118 

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

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

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

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

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

124 

125 @return: Distance along the great circle line (C{meter}, 

126 same units as B{C{radius}}) or C{radians} if 

127 C{B{radius} is None}, positive if I{after} the 

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

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

130 B{C{start}} point. 

131 

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

133 

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

135 

136 @example: 

137 

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

139 

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

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

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

143 ''' 

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

145 cx = cos(x) 

146 return _0_0 if isnear0(cx) else \ 

147 _r2m(copysign0(acos1(cos(r) / cx), cos(b)), radius) 

148 

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

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

151 ''' 

152 s = self.others(start=start) 

153 e = self.others(end=end) 

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

155 

156 r = Radius_(radius) 

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

158 

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

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

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

162 return r, x, -b 

163 

164 @deprecated_method 

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

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

167 ''' 

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

169 

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

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

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

173 

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

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

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

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

178 

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

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

181 ''' 

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

183 sa, ca, sa1, ca1, \ 

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

185 sa1 *= ca2 * ca 

186 

187 x = sa1 * sdb 

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

189 z = ca1 * sdb * ca2 * sa 

190 

191 h = hypot(x, y) 

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

193 return None # great circle doesn't reach latitude 

194 

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

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

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

198 

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

200 '''Compute the (signed, angular) distance from this point to 

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

202 

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

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

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

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

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

208 

209 @return: Distance to the great circle (I{negative} if to 

210 the left or I{positive} if to the right of the 

211 line) (C{meter}, same units as B{C{radius}} or 

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

213 

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

215 

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

217 

218 @example: 

219 

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

221 

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

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

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

225 ''' 

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

227 return _r2m(x, radius) 

228 

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

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

231 travelled the given distance on the given initial bearing. 

232 

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

234 B{C{radius}}). 

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

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

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

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

239 

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

241 

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

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

244 

245 @example: 

246 

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

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

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

250 ''' 

251 a, b = self.philam 

252 r, t = _angular(distance, radius, low=None), Bearing_(bearing) 

253 

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

255 h = self._heigHt(height) 

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

257 

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

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

260 

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

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

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

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

265 

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

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

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

269 

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

271 

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

273 

274 @example: 

275 

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

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

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

279 ''' 

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

281 return _r2m(vincentys_(a2, a1, db), radius) 

282 

283# @Property_RO 

284# def Ecef(self): 

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

286# ''' 

287# return _MODS.ecef.EcefKarney 

288 

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

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

291 on the given initial bearing from this point. 

292 

293 Direction of vector is such that initial bearing vector 

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

295 

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

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

298 overriding the default L{Vector3d}. 

299 @kwarg Vector_kwds: Optional, additional keyword argunents 

300 for B{C{Vector}}. 

301 

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

303 

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

305 

306 @example: 

307 

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

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

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

311 ''' 

312 a, b = self.philam 

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

314 

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

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

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

318 

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

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

321 to an other point. 

322 

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

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

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

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

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

328 C{sphericalNvector.LatLon.initialBearingTo}. 

329 

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

331 

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

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

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

335 

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

337 

338 @example: 

339 

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

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

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

343 ''' 

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

345 # XXX behavior like sphericalNvector.LatLon.initialBearingTo 

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

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

348 

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

350 

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

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

353 and an other point. 

354 

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

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

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

358 @kwarg height: Optional height, overriding the intermediate 

359 height (C{meter}). 

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

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

362 

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

364 

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

366 

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

368 

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

370 

371 @example: 

372 

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

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

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

376 ''' 

377 p = self 

378 f = Scalar(fraction=fraction) 

379 if not isnear0(f): 

380 p = p.others(other) 

381 if wrap: 

382 p = _Wrap.point(p) 

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

384 a1, b1 = self.philam 

385 a2, b2 = p.philam 

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

387 r = vincentys_(a2, a1, db) 

388 sr = sin(r) 

389 if isnon0(sr): 

390 sa1, ca1, sa2, ca2, \ 

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

392 

393 t = f * r 

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

395 b = sin( t) # / sr superflous 

396 

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

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

399 z = a * sa1 + b * sa2 

400 

401 a = atan2(z, hypot(x, y)) 

402 b = atan2(y, x) 

403 

404 else: # PYCHOK no cover 

405 a = favg(a1, a2, f=f) # coincident 

406 b = favg(b1, b2, f=f) 

407 

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

409 p = self.classof(degrees90(a), degrees180(b), height=h) 

410 return p 

411 

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

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

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

415 

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

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

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

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

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

421 C{degrees360}). 

422 @kwarg height: Optional height for intersection point, 

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

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

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

426 

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

428 intersection point might be the L{antipode} to 

429 the returned result. 

430 

431 @raise IntersectionError: Ambiguous or infinite intersection 

432 or colinear, parallel or otherwise 

433 non-intersecting lines. 

434 

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

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

437 

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

439 

440 @example: 

441 

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

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

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

445 ''' 

446 try: 

447 s2 = self.others(other) 

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

449 LatLon=self.classof) 

450 except (TypeError, ValueError) as x: 

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

452 other=other, end2=end2, wrap=wrap) 

453 

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

455 height=None, wrap=True): 

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

457 by a center point and radius. 

458 

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

460 see B{C{radius}}). 

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

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

463 see B{C{radius}}). 

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

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

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

467 B{C{radius}}). 

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

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

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

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

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

473 

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

475 instance. For abutting circles, both intersection 

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

477 

478 @raise IntersectionError: Concentric, antipodal, invalid or 

479 non-intersecting circles. 

480 

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

482 

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

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

485 ''' 

486 try: 

487 c2 = self.others(other) 

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

489 height=height, wrap=wrap, 

490 LatLon=self.classof) 

491 except (TypeError, ValueError) as x: 

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

493 other=other, rad2=rad2, wrap=wrap) 

494 

495 @deprecated_method 

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

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

498 return self.isenclosedBy(points) 

499 

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

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

502 

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

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

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

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

507 

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

509 composite, C{False} otherwise. 

510 

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

512 

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

514 

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

516 

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

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

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

520 ''' 

521 if _MODS.booleans.isBoolean(points): 

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

523 

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

525 n0 = self._N_vector 

526 

527 v2 = Ps[0]._N_vector 

528 p1 = Ps[1] 

529 v1 = p1._N_vector 

530 # check whether this point on same side of all 

531 # polygon edges (to the left or right depending 

532 # on the anti-/clockwise polygon direction) 

533 gc1 = v2.cross(v1) 

534 t0 = gc1.angleTo(n0) > PI_2 

535 s0 = None 

536 # get great-circle vector for each edge 

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

538 if wrap and not Ps.looped: 

539 p2 = _unrollon(p1, p2) 

540 p1 = p2 

541 v2 = p2._N_vector 

542 gc = v1.cross(v2) 

543 t = gc.angleTo(n0) > PI_2 

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

545 return False # outside 

546 

547 # check for convex polygon: angle between 

548 # gc vectors, signed by direction of n0 

549 # (otherwise the test above is not reliable) 

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

551 if s != s0: 

552 if s0 is None: 

553 s0 = s 

554 else: 

555 t = _Fmt.SQUARE(points=i) 

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

557 gc1, v1 = gc, v2 

558 

559 return True # inside 

560 

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

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

563 

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

565 @kwarg height: Optional height for midpoint, overriding 

566 the mean height (C{meter}). 

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

568 may be negative or greater than 1.0. 

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

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

571 

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

573 

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

575 

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

577 

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

579 

580 @example: 

581 

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

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

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

585 ''' 

586 if fraction is _0_5: 

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

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

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

590 

591 x = ca2 * cdb + ca1 

592 y = ca2 * sdb 

593 

594 a = atan2(sa1 + sa2, hypot(x, y)) 

595 b += atan2(y, x) 

596 

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

598 r = self.classof(degrees90(a), degrees180(b), height=h) 

599 else: 

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

601 return r 

602 

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

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

605 

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

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

608 

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

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

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

612 @kwarg wrap_adjust_limit: Optional keyword arguments for functions 

613 L{sphericalTrigonometry.nearestOn3} and 

614 L{pygeodesy.equirectangular_}, 

615 

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

617 

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

619 see function L{pygeodesy.equirectangular_}. 

620 

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

622 is not (yet) supported. 

623 

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

625 

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

627 

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

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

630 ''' 

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

632 w = _xkwds_pop(wrap_adjust_limit, within=True) 

633 if not w: 

634 notImplemented(self, within=w) 

635 

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

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

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

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

640# return point1 

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

642# if isnear0(d): 

643# return point1 # or point2 

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

645# return point2 

646# f = a / d 

647# if within: 

648# if f > EPS1: 

649# return point2 

650# elif f < EPS: 

651# return point1 

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

653 

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

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

656 **wrap_adjust_limit)[0] 

657 

658 @deprecated_method 

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

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

661 

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

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

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

665 units of B{C{radius}}. 

666 ''' 

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

668 return r.closest, r.distance 

669 

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

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

672 

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

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

675 

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

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

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

679 @kwarg wrap_adjust_limit: Optional keyword arguments for function 

680 L{sphericalTrigonometry.nearestOn3} and 

681 L{pygeodesy.equirectangular_}, 

682 

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

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

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

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

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

688 function L{pygeodesy.compassAngle}. 

689 

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

691 see function L{pygeodesy.equirectangular_}. 

692 

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

694 

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

696 

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

698 

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

700 and L{pygeodesy.nearestOn5}. 

701 ''' 

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

703 LatLon=self.classof, **wrap_adjust_limit) 

704 

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

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

707 coordinates. 

708 

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

710 and other keyword arguments, ignored 

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

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

713 this L{Cartesian} class or specify 

714 C{B{Cartesian}=None}. 

715 

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

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

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

719 

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

721 ''' 

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

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

724 

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

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

727 

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

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

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

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

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

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

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

735 

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

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

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

739 

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

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

742 ''' 

743 B = self.others(otherB=otherB) 

744 C = self.others(otherC=otherC) 

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

746 

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

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

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

750 

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

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

753 '''Trilaterate three points by area overlap or perimeter intersection 

754 of three corresponding circles. 

755 

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

757 as B{C{radius}}). 

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

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

760 B{C{radius}}). 

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

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

763 B{C{radius}}). 

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

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

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

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

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

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

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

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

772 

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

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

775 the corresponding trilaterated points C{minPoint} and 

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

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

778 

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

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

781 

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

783 respectively largest I{radial} overlap found. 

784 

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

786 nearest respectively farthest intersection margin. 

787 

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

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

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

791 largest B{C{distance#}}. 

792 

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

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

795 no intersection or all (near-)concentric for 

796 C{B{area}=False}. 

797 

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

799 

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

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

802 or B{C{radius}}. 

803 ''' 

804 return _trilaterate5(self, distance1, 

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

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

807 area=area, radius=radius, eps=eps, wrap=wrap) 

808 

809 

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

811 

812 

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

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

815 (with the pointsjoined by great circle arcs). 

816 

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

818 or L{BooleanGH}). 

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

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

821 or C{None}. 

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

823 (C{bool}). 

824 

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

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

827 

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

829 

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

831 

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

833 

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

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

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

837 

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

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

840 L{ellipsoidalKarney.areaOf}. 

841 

842 @example: 

843 

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

845 >>> areaOf(b) # 8666058750.718977 

846 

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

848 >>> areaOf(c) # 6.18e9 

849 ''' 

850 if _MODS.booleans.isBoolean(points): 

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

852 

853 _at2, _t_2 = atan2, tan_2 

854 _un, _w180 = unrollPI, wrap180 

855 

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

857 p1 = p2 = Ps[0] 

858 a1, b1 = p1.philam 

859 ta1, z1 = _t_2(a1), None 

860 

861 A = Fsum() # mean phi 

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

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

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

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

866 D = Fsum() 

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

868 a2, b2 = p2.philam 

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

870 A += a2 

871 ta2 = _t_2(a2) 

872 tdb = _t_2(db, points=i) 

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

874 _1_0 + ta1 * ta2) 

875 ta1, b1 = ta2, b2 

876 

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

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

879 if z1 is not None: 

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

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

882 p1, z1 = p2, z2 

883 

884 R = abs(R * _2_0) 

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

886 R = abs(R - PI2) 

887 if radius: 

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

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

890 return float(R) 

891 

892 

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

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

895 

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

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

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

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

900 

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

902 ''' 

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

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

905 ca *= sr 

906 

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

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

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

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

911 

912 

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

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

915 # and similar logic in .ellipsoidalBaseDI._intersect3 

916 a1, b1 = s.philam 

917 

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

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

920 else: # must be a point 

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

922 hs.append(end.height) 

923 a2, b2 = end.philam 

924 if wrap: 

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

926 

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

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

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

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

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

932 -(b1 + b2) * _0_5) 

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

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

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

936 cb12 * cb21 + sb12 * sb21, 

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

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

939 

940 

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

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

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

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

945 

946 

947def intersecant2(center, circle, point, bearing, radius=R_M, exact=False, 

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

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

950 

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

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

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

954 @arg point: A point in- or outside the circle (L{LatLon}). 

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

956 a second point on the line (L{LatLon}). 

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

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

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

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

961 circle} methods. 

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

963 conventionally) or C{None}. 

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

965 and the B{C{circle}} and B{C{bearing}} points (C{bool}). 

966 

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

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

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

970 or I{normalized}. 

971 

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

973 

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

975 B{C{circle}} or B{C{bearing}} invalid. 

976 

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

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

979 ''' 

980 c = _T00.others(center=center) 

981 p = _T00.others(point=point) 

982 try: 

983 return _intersecant2(c, circle, p, bearing, radius=radius, exact=exact, 

984 height=height, wrap=wrap) 

985 except (TypeError, ValueError) as x: 

986 raise _xError(x, center=center, circle=circle, point=point, bearing=bearing, exact=exact) 

987 

988 

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

990 LatLon=None, **LatLon_kwds): 

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

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

993 

994 s1, s2 = start1, start2 

995 if wrap: 

996 s2 = _Wrap.point(s2) 

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

998 

999 a1, b1 = s1.philam 

1000 a2, b2 = s2.philam 

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

1002 r12 = vincentys_(a2, a1, db) 

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

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

1005 

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

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

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

1009 

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

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

1012 raise IntersectionError(_parallel_) 

1013 # handle domain error for equivalent longitudes, 

1014 # see also functions asin_safe and acos_safe at 

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

1016 t1, t2 = acos1((sa2 - sa1 * cr12) / x1), \ 

1017 acos1((sa1 - sa2 * cr12) / x2) 

1018 if sin(db) > 0: 

1019 t12, t21 = t1, PI2 - t2 

1020 else: 

1021 t12, t21 = PI2 - t1, t2 

1022 t13, t23 = radiansPI2(end1), radiansPI2(end2) 

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

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

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

1026 raise IntersectionError(_infinite_) 

1027 sx3 = sx1 * sx2 

1028# XXX if sx3 < 0: 

1029# XXX raise ValueError(_ambiguous_) 

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

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

1032 

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

1034 # like .ellipsoidalBaseDI,_intersect3, if this intersection 

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

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

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

1038 

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

1040 _N_vector_ = _MODS.nvectorBase._N_vector_ 

1041 

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

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

1044 x = x1.cross(x2) 

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

1046 raise IntersectionError(_colinear_) 

1047 a, b = x.philam 

1048 # choose intersection similar to sphericalNvector 

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

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

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

1052 

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

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

1055 intersection, LatLon, LatLon_kwds) 

1056 

1057 

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

1059 LatLon=LatLon, **LatLon_kwds): 

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

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

1062 

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

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

1065 the initial bearing at the first start point 

1066 (compass C{degrees360}). 

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

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

1069 the initial bearing at the second start point 

1070 (compass C{degrees360}). 

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

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

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

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

1075 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1079 

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

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

1082 height)}. An alternate intersection point might 

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

1084 

1085 @raise IntersectionError: Ambiguous or infinite intersection 

1086 or colinear, parallel or otherwise 

1087 non-intersecting lines. 

1088 

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

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

1091 

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

1093 

1094 @example: 

1095 

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

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

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

1099 ''' 

1100 s1 = _T00.others(start1=start1) 

1101 s2 = _T00.others(start2=start2) 

1102 try: 

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

1104 LatLon=LatLon, **LatLon_kwds) 

1105 except (TypeError, ValueError) as x: 

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

1107 

1108 

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

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

1111 LatLon=LatLon, **LatLon_kwds): 

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

1113 by a center point and a radius. 

1114 

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

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

1117 see B{C{radius}}). 

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

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

1120 see B{C{radius}}). 

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

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

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

1124 B{C{radius}}). 

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

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

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

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

1129 (C{bool}). 

1130 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1134 

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

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

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

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

1139 

1140 @raise IntersectionError: Concentric, antipodal, invalid or 

1141 non-intersecting circles. 

1142 

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

1144 

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

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

1147 

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

1149 

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

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

1152 ''' 

1153 c1 = _T00.others(center1=center1) 

1154 c2 = _T00.others(center2=center2) 

1155 try: 

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

1157 height=height, wrap=wrap, 

1158 LatLon=LatLon, **LatLon_kwds) 

1159 except (TypeError, ValueError) as x: 

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

1161 center2=center2, rad2=rad2, wrap=wrap) 

1162 

1163 

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

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

1166 LatLon=LatLon, **LatLon_kwds): 

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

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

1169 

1170 def _dest3(bearing, h): 

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

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

1173 intersections2, LatLon, LatLon_kwds) 

1174 

1175 a1, b1 = c1.philam 

1176 a2, b2 = c2.philam 

1177 if wrap: 

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

1179 

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

1181 if f: # swapped radii, swap centers 

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

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

1184 

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

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

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

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

1189 

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

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

1192 if r < _0_0: 

1193 raise _ValueError(eps=r) 

1194 

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

1196 if x > max(r, EPS): 

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

1198 x = sd * sr1 

1199 if isnear0(x): 

1200 raise _ValueError(_invalid_) 

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

1202 

1203 elif x < r: # PYCHOK no cover 

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

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

1206 

1207 if height is None: # "radical height" 

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

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

1210 else: 

1211 h = Height(height) 

1212 

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

1214 if x < EPS4: # externally ... 

1215 r = _dest3(b, h) 

1216 elif x > _PI_EPS4: # internally ... 

1217 r = _dest3(b + PI, h) 

1218 else: 

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

1220 return r, r # ... abutting circles 

1221 

1222 

1223@deprecated_function 

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

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

1226 ''' 

1227 return ispolar(points, wrap=wrap) 

1228 

1229 

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

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

1232 ''' 

1233 n = where.__name__ 

1234 if LatLon is None: 

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

1236 else: 

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

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

1239 return r 

1240 

1241 

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

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

1244 

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

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

1247 height (C{meter}). 

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

1249 (C{bool}). 

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

1251 or C{None}. 

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

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

1254 

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

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

1257 is C{None}. 

1258 

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

1260 

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

1262 ''' 

1263 def _N_vs(ps, w): 

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

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

1266 yield p._N_vector 

1267 

1268 m = _MODS.nvectorBase 

1269 # geographic, vectorial mean 

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

1271 lat, lon, h = n.latlonheight 

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

1273 

1274 

1275@deprecated_function 

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

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

1278 

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

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

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

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

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

1284 ''' 

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

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

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

1288 return ll, d 

1289 

1290 

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

1292 limit=9, **LatLon_and_kwds): 

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

1294 

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

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

1297 

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

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

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

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

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

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

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

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

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

1307 spherical earth radius L{R_KM}). 

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

1309 or C{None}. 

1310 @kwarg options: Optional keyword arguments for function 

1311 L{pygeodesy.equirectangular_}. 

1312 

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

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

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

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

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

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

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

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

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

1322 

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

1324 see function L{pygeodesy.equirectangular_}. 

1325 

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

1327 

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

1329 

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

1331 

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

1333 ''' 

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

1335 adjust=adjust, limit=limit) 

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

1337 h = t.height 

1338 n = nearestOn3.__name__ 

1339 

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

1341 LL = _xkwds_pop(kwds, LatLon=LatLon) 

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

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

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

1345 

1346 

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

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

1349 (with great circle arcs joining the points). 

1350 

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

1352 or L{BooleanGH}). 

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

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

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

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

1357 

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

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

1360 

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

1362 

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

1364 

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

1366 C{B{points}} a composite. 

1367 

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

1369 

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

1371 and L{ellipsoidalKarney.perimeterOf}. 

1372 ''' 

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

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

1375 a1, b1 = Ps[0].philam 

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

1377 a2, b2 = p.philam 

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

1379 yield vincentys_(a2, a1, db) 

1380 a1, b1 = a2, b2 

1381 

1382 if _MODS.booleans.isBoolean(points): 

1383 if not closed: 

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

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

1386 else: 

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

1388 return _r2m(r, radius) 

1389 

1390 

1391def _r2m(r, radius): 

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

1393 ''' 

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

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

1396 return r 

1397 

1398 

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

1400 excess=excessAbc_, 

1401 wrap=False): 

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

1403 

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

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

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

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

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

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

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

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

1412 or C{None}. 

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

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

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

1416 longitudes (C{bool}). 

1417 

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

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

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

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

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

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

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

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

1426 ''' 

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

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

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

1430 excess=excess, wrap=wrap) 

1431 return _t7Tuple(t, radius) 

1432 

1433 

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

1435 wrap=False): 

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

1437 excess} of a (spherical) triangle. 

1438 

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

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

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

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

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

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

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

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

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

1448 longitudinal deltas (C{bool}). 

1449 

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

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

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

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

1454 ''' 

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

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

1457 a = vincentys_(phiC, phiB, d) 

1458 return a, (phiB, lamB, phiC, lamC, phiA, lamA) 

1459 

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

1461 s = sb * sc 

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

1463 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2's 

1464 

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

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

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

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

1469 

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

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

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

1473 

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

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

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

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

1478 

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

1480 

1481 

1482def _t7Tuple(t, radius): 

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

1484 ''' 

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

1486 r = radius if isscalar(radius) else \ 

1487 _ellipsoidal_datum(radius).ellipsoid.Rmean 

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

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

1490 B, (r * t.b), 

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

1492 return t 

1493 

1494 

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

1496 areaOf, # functions 

1497 intersecant2, intersection, intersections2, ispolar, 

1498 isPoleEnclosedBy, # DEPRECATED, use ispolar 

1499 meanOf, 

1500 nearestOn2, nearestOn3, 

1501 perimeterOf, 

1502 sumOf, # XXX == vector3d.sumOf 

1503 triangle7, triangle8_) 

1504 

1505# **) MIT License 

1506# 

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

1508# 

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

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

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

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

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

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

1515# 

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

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

1518# 

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

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

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

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

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

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

1525# OTHER DEALINGS IN THE SOFTWARE.