Coverage for pygeodesy/sphericalTrigonometry.py: 94%

390 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-10-04 14:05 -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_, _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 _angular, CartesianSphericalBase, _intersecant2, \ 

45 LatLonSphericalBase, _rads3, _r2m, _trilaterate5 

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

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

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

49 degrees2m, m2radians, radiansPI2, sincos2_, tan_2, \ 

50 unrollPI, _unrollon, _unrollon3, _Wrap, wrap180, wrapPI 

51from pygeodesy.vector3d import sumOf, Vector3d 

52 

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

54 

55__all__ = _ALL_LAZY.sphericalTrigonometry 

56__version__ = '23.09.29' 

57 

58_PI_EPS4 = PI - EPS4 

59if _PI_EPS4 >= PI: 

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

61 

62 

63class Cartesian(CartesianSphericalBase): 

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

65 spherical, geodetic L{LatLon}. 

66 ''' 

67 

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

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

70 

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

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

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

74 

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

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

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

78 

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

80 ''' 

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

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

83 

84 

85class LatLon(LatLonSphericalBase): 

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

87 

88 @example: 

89 

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

91 ''' 

92 

93 def _ab1_ab2_db5(self, other, wrap): 

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

95 ''' 

96 a1, b1 = self.philam 

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

98 if wrap: 

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

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

101 else: # unrollPI shortcut 

102 db = b2 - b1 

103 return a1, b1, a2, b2, db 

104 

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

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

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

108 end point. 

109 

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

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

112 from the start point to the point where the perpendicular 

113 crosses the line. 

114 

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

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

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

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

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

120 

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

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

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

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

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

126 B{C{start}} point. 

127 

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

129 

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

131 

132 @example: 

133 

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

135 

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

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

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

139 ''' 

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

141 cx = cos(x) 

142 return _0_0 if isnear0(cx) else \ 

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

144 

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

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

147 ''' 

148 s = self.others(start=start) 

149 e = self.others(end=end) 

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

151 

152 r = Radius_(radius) 

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

154 

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

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

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

158 return r, x, -b 

159 

160 @deprecated_method 

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

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

163 ''' 

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

165 

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

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

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

169 

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

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

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

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

174 

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

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

177 ''' 

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

179 sa, ca, sa1, ca1, \ 

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

181 sa1 *= ca2 * ca 

182 

183 x = sa1 * sdb 

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

185 z = ca1 * sdb * ca2 * sa 

186 

187 h = hypot(x, y) 

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

189 return None # great circle doesn't reach latitude 

190 

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

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

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

194 

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

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

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

198 

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

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

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

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

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

204 

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

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

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

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

209 

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

211 

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

213 

214 @example: 

215 

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

217 

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

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

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

221 ''' 

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

223 return _r2m(x, radius) 

224 

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

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

227 travelled the given distance on the given initial bearing. 

228 

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

230 B{C{radius}}). 

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

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

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

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

235 

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

237 

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

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

240 

241 @example: 

242 

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

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

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

246 ''' 

247 a, b = self.philam 

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

249 

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

251 h = self._heigHt(height) 

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

253 

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

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

256 

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

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

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

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

261 

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

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

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

265 

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

267 

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

269 

270 @example: 

271 

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

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

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

275 ''' 

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

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

278 

279# @Property_RO 

280# def Ecef(self): 

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

282# ''' 

283# return _MODS.ecef.EcefKarney 

284 

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

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

287 on the given initial bearing from this point. 

288 

289 Direction of vector is such that initial bearing vector 

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

291 

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

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

294 overriding the default L{Vector3d}. 

295 @kwarg Vector_kwds: Optional, additional keyword argunents 

296 for B{C{Vector}}. 

297 

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

299 

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

301 

302 @example: 

303 

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

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

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

307 ''' 

308 a, b = self.philam 

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

310 

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

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

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

314 

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

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

317 to an other point. 

318 

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

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

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

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

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

324 C{sphericalNvector.LatLon.initialBearingTo}. 

325 

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

327 

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

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

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

331 

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

333 

334 @example: 

335 

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

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

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

339 ''' 

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

341 # XXX behavior like sphericalNvector.LatLon.initialBearingTo 

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

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

344 

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

346 

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

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

349 and an other point. 

350 

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

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

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

354 @kwarg height: Optional height, overriding the intermediate 

355 height (C{meter}). 

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

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

358 

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

360 

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

362 

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

364 

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

366 

367 @example: 

368 

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

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

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

372 ''' 

373 p = self 

374 f = Scalar(fraction=fraction) 

375 if not isnear0(f): 

376 p = p.others(other) 

377 if wrap: 

378 p = _Wrap.point(p) 

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

380 a1, b1 = self.philam 

381 a2, b2 = p.philam 

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

383 r = vincentys_(a2, a1, db) 

384 sr = sin(r) 

385 if isnon0(sr): 

386 sa1, ca1, sa2, ca2, \ 

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

388 

389 t = f * r 

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

391 b = sin( t) # / sr superflous 

392 

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

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

395 z = a * sa1 + b * sa2 

396 

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

398 b = atan2d(y, x) 

399 

400 else: # PYCHOK no cover 

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

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

403 

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

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

406 return p 

407 

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

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

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

411 

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

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

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

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

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

417 C{degrees360}). 

418 @kwarg height: Optional height for intersection point, 

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

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

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

422 

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

424 intersection point might be the L{antipode} to 

425 the returned result. 

426 

427 @raise IntersectionError: Ambiguous or infinite intersection 

428 or colinear, parallel or otherwise 

429 non-intersecting lines. 

430 

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

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

433 

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

435 

436 @example: 

437 

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

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

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

441 ''' 

442 try: 

443 s2 = self.others(other) 

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

445 LatLon=self.classof) 

446 except (TypeError, ValueError) as x: 

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

448 other=other, end2=end2, wrap=wrap) 

449 

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

451 height=None, wrap=True): 

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

453 by a center point and radius. 

454 

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

456 see B{C{radius}}). 

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

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

459 see B{C{radius}}). 

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

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

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

463 B{C{radius}}). 

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

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

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

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

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

469 

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

471 instance. For abutting circles, both intersection 

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

473 

474 @raise IntersectionError: Concentric, antipodal, invalid or 

475 non-intersecting circles. 

476 

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

478 

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

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

481 ''' 

482 try: 

483 c2 = self.others(other) 

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

485 height=height, wrap=wrap, 

486 LatLon=self.classof) 

487 except (TypeError, ValueError) as x: 

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

489 other=other, rad2=rad2, wrap=wrap) 

490 

491 @deprecated_method 

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

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

494 return self.isenclosedBy(points) 

495 

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

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

498 

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

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

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

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

503 

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

505 composite, C{False} otherwise. 

506 

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

508 

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

510 

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

512 

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

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

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

516 ''' 

517 if _MODS.booleans.isBoolean(points): 

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

519 

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

521 n0 = self._N_vector 

522 

523 v2 = Ps[0]._N_vector 

524 p1 = Ps[1] 

525 v1 = p1._N_vector 

526 # check whether this point on same side of all 

527 # polygon edges (to the left or right depending 

528 # on the anti-/clockwise polygon direction) 

529 gc1 = v2.cross(v1) 

530 t0 = gc1.angleTo(n0) > PI_2 

531 s0 = None 

532 # get great-circle vector for each edge 

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

534 if wrap and not Ps.looped: 

535 p2 = _unrollon(p1, p2) 

536 p1 = p2 

537 v2 = p2._N_vector 

538 gc = v1.cross(v2) 

539 t = gc.angleTo(n0) > PI_2 

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

541 return False # outside 

542 

543 # check for convex polygon: angle between 

544 # gc vectors, signed by direction of n0 

545 # (otherwise the test above is not reliable) 

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

547 if s != s0: 

548 if s0 is None: 

549 s0 = s 

550 else: 

551 t = _Fmt.SQUARE(points=i) 

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

553 gc1, v1 = gc, v2 

554 

555 return True # inside 

556 

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

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

559 

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

561 @kwarg height: Optional height for midpoint, overriding 

562 the mean height (C{meter}). 

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

564 may be negative or greater than 1.0. 

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

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

567 

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

569 

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

571 

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

573 

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

575 

576 @example: 

577 

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

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

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

581 ''' 

582 if fraction is _0_5: 

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

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

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

586 

587 x = ca2 * cdb + ca1 

588 y = ca2 * sdb 

589 

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

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

592 

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

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

595 else: 

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

597 return r 

598 

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

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

601 

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

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

604 

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

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

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

608 @kwarg wrap_adjust_limit: Optional keyword arguments for functions 

609 L{sphericalTrigonometry.nearestOn3} and 

610 L{pygeodesy.equirectangular_}, 

611 

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

613 

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

615 see function L{pygeodesy.equirectangular_}. 

616 

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

618 is not (yet) supported. 

619 

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

621 

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

623 

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

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

626 ''' 

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

628 w = _xkwds_pop(wrap_adjust_limit, within=True) 

629 if not w: 

630 notImplemented(self, within=w) 

631 

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

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

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

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

636# return point1 

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

638# if isnear0(d): 

639# return point1 # or point2 

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

641# return point2 

642# f = a / d 

643# if within: 

644# if f > EPS1: 

645# return point2 

646# elif f < EPS: 

647# return point1 

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

649 

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

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

652 **wrap_adjust_limit)[0] 

653 

654 @deprecated_method 

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

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

657 

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

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

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

661 units of B{C{radius}}. 

662 ''' 

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

664 return r.closest, r.distance 

665 

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

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

668 

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

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

671 

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

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

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

675 @kwarg wrap_adjust_limit: Optional keyword arguments for function 

676 L{sphericalTrigonometry.nearestOn3} and 

677 L{pygeodesy.equirectangular_}, 

678 

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

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

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

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

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

684 function L{pygeodesy.compassAngle}. 

685 

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

687 see function L{pygeodesy.equirectangular_}. 

688 

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

690 

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

692 

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

694 

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

696 and L{pygeodesy.nearestOn5}. 

697 ''' 

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

699 LatLon=self.classof, **wrap_adjust_limit) 

700 

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

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

703 coordinates. 

704 

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

706 and other keyword arguments, ignored 

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

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

709 this L{Cartesian} class or specify 

710 C{B{Cartesian}=None}. 

711 

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

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

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

715 

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

717 ''' 

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

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

720 

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

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

723 

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

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

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

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

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

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

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

731 

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

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

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

735 

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

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

738 ''' 

739 B = self.others(otherB=otherB) 

740 C = self.others(otherC=otherC) 

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

742 

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

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

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

746 

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

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

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

750 intersection} of three corresponding circles. 

751 

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

753 as B{C{radius}}). 

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

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

756 B{C{radius}}). 

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

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

759 B{C{radius}}). 

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

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

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

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

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

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

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

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

768 

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

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

771 the corresponding trilaterated points C{minPoint} and 

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

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

774 

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

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

777 

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

779 respectively largest I{radial} overlap found. 

780 

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

782 nearest respectively farthest intersection margin. 

783 

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

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

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

787 largest B{C{distance#}}. 

788 

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

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

791 no intersection or all (near-)concentric for 

792 C{B{area}=False}. 

793 

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

795 

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

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

798 or B{C{radius}}. 

799 ''' 

800 return _trilaterate5(self, distance1, 

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

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

803 area=area, radius=radius, eps=eps, wrap=wrap) 

804 

805 

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

807 

808 

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

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

811 (with the pointsjoined by great circle arcs). 

812 

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

814 or L{BooleanGH}). 

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

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

817 or C{None}. 

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

819 (C{bool}). 

820 

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

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

823 

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

825 

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

827 

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

829 

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

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

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

833 

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

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

836 L{ellipsoidalKarney.areaOf}. 

837 

838 @example: 

839 

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

841 >>> areaOf(b) # 8666058750.718977 

842 

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

844 >>> areaOf(c) # 6.18e9 

845 ''' 

846 if _MODS.booleans.isBoolean(points): 

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

848 

849 _at2, _t_2 = atan2, tan_2 

850 _un, _w180 = unrollPI, wrap180 

851 

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

853 p1 = p2 = Ps[0] 

854 a1, b1 = p1.philam 

855 ta1, z1 = _t_2(a1), None 

856 

857 A = Fsum() # mean phi 

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

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

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

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

862 D = Fsum() 

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

864 a2, b2 = p2.philam 

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

866 A += a2 

867 ta2 = _t_2(a2) 

868 tdb = _t_2(db, points=i) 

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

870 _1_0 + ta1 * ta2) 

871 ta1, b1 = ta2, b2 

872 

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

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

875 if z1 is not None: 

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

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

878 p1, z1 = p2, z2 

879 

880 R = abs(R * _2_0) 

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

882 R = abs(R - PI2) 

883 if radius: 

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

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

886 return float(R) 

887 

888 

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

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

891 

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

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

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

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

896 

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

898 ''' 

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

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

901 ca *= sr 

902 

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

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

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

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

907 

908 

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

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

911 # and similar logic in .ellipsoidalBaseDI._intersect3 

912 a1, b1 = s.philam 

913 

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

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

916 else: # must be a point 

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

918 hs.append(end.height) 

919 a2, b2 = end.philam 

920 if wrap: 

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

922 

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

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

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

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

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

928 -(b1 + b2) * _0_5) 

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

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

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

932 cb12 * cb21 + sb12 * sb21, 

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

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

935 

936 

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

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

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

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

941 

942 

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

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

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

946 

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

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

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

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

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

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

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

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

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

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

957 circle} methods. 

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

959 conventionally) or C{None}. 

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

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

962 

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

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

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

966 or I{normalized}. 

967 

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

969 

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

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

972 

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

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

975 ''' 

976 c = _T00.others(center=center) 

977 p = _T00.others(point=point) 

978 try: 

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

980 height=height, wrap=wrap) 

981 except (TypeError, ValueError) as x: 

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

983 

984 

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

986 LatLon=None, **LatLon_kwds): 

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

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

989 

990 s1, s2 = start1, start2 

991 if wrap: 

992 s2 = _Wrap.point(s2) 

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

994 

995 a1, b1 = s1.philam 

996 a2, b2 = s2.philam 

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

998 r12 = vincentys_(a2, a1, db) 

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

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

1001 

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

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

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

1005 

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

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

1008 raise IntersectionError(_parallel_) 

1009 # handle domain error for equivalent longitudes, 

1010 # see also functions asin_safe and acos_safe at 

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

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

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

1014 if sin(db) > 0: 

1015 t21 = PI2 - t21 

1016 else: 

1017 t12 = PI2 - t12 

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

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

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

1021 raise IntersectionError(_infinite_) 

1022 sx3 = sx1 * sx2 

1023# XXX if sx3 < 0: 

1024# XXX raise ValueError(_ambiguous_) 

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

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

1027 

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

1029 # like .ellipsoidalBaseDI,_intersect3, if this intersection 

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

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

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

1033 

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

1035 _N_vector_ = _MODS.nvectorBase._N_vector_ 

1036 

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

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

1039 x = x1.cross(x2) 

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

1041 raise IntersectionError(_colinear_) 

1042 a, b = x.philam 

1043 # choose intersection similar to sphericalNvector 

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

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

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

1047 

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

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

1050 intersection, LatLon, LatLon_kwds) 

1051 

1052 

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

1054 LatLon=LatLon, **LatLon_kwds): 

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

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

1057 

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

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

1060 the initial bearing at the first start point 

1061 (compass C{degrees360}). 

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

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

1064 the initial bearing at the second start point 

1065 (compass C{degrees360}). 

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

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

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

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

1070 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1074 

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

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

1077 height)}. An alternate intersection point might 

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

1079 

1080 @raise IntersectionError: Ambiguous or infinite intersection 

1081 or colinear, parallel or otherwise 

1082 non-intersecting lines. 

1083 

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

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

1086 

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

1088 

1089 @example: 

1090 

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

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

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

1094 ''' 

1095 s1 = _T00.others(start1=start1) 

1096 s2 = _T00.others(start2=start2) 

1097 try: 

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

1099 LatLon=LatLon, **LatLon_kwds) 

1100 except (TypeError, ValueError) as x: 

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

1102 

1103 

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

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

1106 LatLon=LatLon, **LatLon_kwds): 

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

1108 by a center point and a radius. 

1109 

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

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

1112 see B{C{radius}}). 

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

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

1115 see B{C{radius}}). 

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

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

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

1119 B{C{radius}}). 

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

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

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

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

1124 (C{bool}). 

1125 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1129 

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

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

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

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

1134 

1135 @raise IntersectionError: Concentric, antipodal, invalid or 

1136 non-intersecting circles. 

1137 

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

1139 

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

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

1142 

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

1144 

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

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

1147 ''' 

1148 c1 = _T00.others(center1=center1) 

1149 c2 = _T00.others(center2=center2) 

1150 try: 

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

1152 height=height, wrap=wrap, 

1153 LatLon=LatLon, **LatLon_kwds) 

1154 except (TypeError, ValueError) as x: 

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

1156 center2=center2, rad2=rad2, wrap=wrap) 

1157 

1158 

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

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

1161 LatLon=LatLon, **LatLon_kwds): 

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

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

1164 

1165 def _dest3(bearing, h): 

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

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

1168 intersections2, LatLon, LatLon_kwds) 

1169 

1170 a1, b1 = c1.philam 

1171 a2, b2 = c2.philam 

1172 if wrap: 

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

1174 

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

1176 if f: # swapped radii, swap centers 

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

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

1179 

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

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

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

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

1184 

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

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

1187 if r < _0_0: 

1188 raise _ValueError(eps=r) 

1189 

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

1191 if x > max(r, EPS): 

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

1193 x = sd * sr1 

1194 if isnear0(x): 

1195 raise _ValueError(_invalid_) 

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

1197 

1198 elif x < r: # PYCHOK no cover 

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

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

1201 

1202 if height is None: # "radical height" 

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

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

1205 else: 

1206 h = Height(height) 

1207 

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

1209 if x < EPS4: # externally ... 

1210 r = _dest3(b, h) 

1211 elif x > _PI_EPS4: # internally ... 

1212 r = _dest3(b + PI, h) 

1213 else: 

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

1215 return r, r # ... abutting circles 

1216 

1217 

1218@deprecated_function 

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

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

1221 ''' 

1222 return ispolar(points, wrap=wrap) 

1223 

1224 

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

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

1227 ''' 

1228 n = where.__name__ 

1229 if LatLon is None: 

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

1231 else: 

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

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

1234 return r 

1235 

1236 

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

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

1239 

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

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

1242 height (C{meter}). 

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

1244 (C{bool}). 

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

1246 or C{None}. 

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

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

1249 

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

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

1252 is C{None}. 

1253 

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

1255 

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

1257 ''' 

1258 def _N_vs(ps, w): 

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

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

1261 yield p._N_vector 

1262 

1263 m = _MODS.nvectorBase 

1264 # geographic, vectorial mean 

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

1266 lat, lon, h = n.latlonheight 

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

1268 

1269 

1270@deprecated_function 

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

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

1273 

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

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

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

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

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

1279 ''' 

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

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

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

1283 return ll, d 

1284 

1285 

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

1287 limit=9, **LatLon_and_kwds): 

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

1289 

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

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

1292 

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

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

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

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

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

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

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

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

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

1302 spherical earth radius L{R_KM}). 

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

1304 or C{None}. 

1305 @kwarg options: Optional keyword arguments for function 

1306 L{pygeodesy.equirectangular_}. 

1307 

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

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

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

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

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

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

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

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

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

1317 

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

1319 see function L{pygeodesy.equirectangular_}. 

1320 

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

1322 

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

1324 

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

1326 

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

1328 ''' 

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

1330 adjust=adjust, limit=limit) 

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

1332 h = t.height 

1333 n = nearestOn3.__name__ 

1334 

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

1336 LL = _xkwds_pop(kwds, LatLon=LatLon) 

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

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

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

1340 

1341 

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

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

1344 (with great circle arcs joining the points). 

1345 

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

1347 or L{BooleanGH}). 

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

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

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

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

1352 

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

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

1355 

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

1357 

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

1359 

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

1361 C{B{points}} a composite. 

1362 

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

1364 

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

1366 and L{ellipsoidalKarney.perimeterOf}. 

1367 ''' 

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

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

1370 a1, b1 = Ps[0].philam 

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

1372 a2, b2 = p.philam 

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

1374 yield vincentys_(a2, a1, db) 

1375 a1, b1 = a2, b2 

1376 

1377 if _MODS.booleans.isBoolean(points): 

1378 if not closed: 

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

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

1381 else: 

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

1383 return _r2m(r, radius) 

1384 

1385 

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

1387 excess=excessAbc_, 

1388 wrap=False): 

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

1390 

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

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

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

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

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

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

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

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

1399 or C{None}. 

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

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

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

1403 longitudes (C{bool}). 

1404 

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

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

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

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

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

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

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

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

1413 ''' 

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

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

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

1417 excess=excess, wrap=wrap) 

1418 return _t7Tuple(t, radius) 

1419 

1420 

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

1422 wrap=False): 

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

1424 excess} of a (spherical) triangle. 

1425 

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

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

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

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

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

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

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

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

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

1435 longitudinal deltas (C{bool}). 

1436 

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

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

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

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

1441 ''' 

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

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

1444 a = vincentys_(phiC, phiB, d) 

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

1446 

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

1448 s = sb * sc 

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

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

1451 

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

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

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

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

1456 

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

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

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

1460 

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

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

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

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

1465 

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

1467 

1468 

1469def _t7Tuple(t, radius): 

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

1471 ''' 

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

1473 r = radius if isscalar(radius) else \ 

1474 _ellipsoidal_datum(radius).ellipsoid.Rmean 

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

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

1477 B, (r * t.b), 

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

1479 return t 

1480 

1481 

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

1483 areaOf, # functions 

1484 intersecant2, intersection, intersections2, ispolar, 

1485 isPoleEnclosedBy, # DEPRECATED, use ispolar 

1486 meanOf, 

1487 nearestOn2, nearestOn3, 

1488 perimeterOf, 

1489 sumOf, # XXX == vector3d.sumOf 

1490 triangle7, triangle8_) 

1491 

1492# **) MIT License 

1493# 

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

1495# 

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

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

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

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

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

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

1502# 

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

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

1505# 

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

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

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

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

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

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

1512# OTHER DEALINGS IN THE SOFTWARE.