Coverage for pygeodesy/sphericalTrigonometry.py: 94%

390 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-28 15:52 -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 _parallel_, _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, _r2m, _trilaterate5 

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

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

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

50 m2radians, radiansPI2, sincos2_, tan_2, _unrollon, \ 

51 unrollPI, _unrollon3, _Wrap, wrap180, wrapPI 

52from pygeodesy.vector3d import sumOf, Vector3d 

53 

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

55 

56__all__ = _ALL_LAZY.sphericalTrigonometry 

57__version__ = '23.08.25' 

58 

59_PI_EPS4 = PI - EPS4 

60if _PI_EPS4 >= PI: 

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

62 

63 

64class Cartesian(CartesianSphericalBase): 

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

66 spherical, geodetic L{LatLon}. 

67 ''' 

68 

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

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

71 

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

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

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

75 

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

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

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

79 

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

81 ''' 

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

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

84 

85 

86class LatLon(LatLonSphericalBase): 

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

88 

89 @example: 

90 

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

92 ''' 

93 

94 def _ab1_ab2_db5(self, other, wrap): 

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

96 ''' 

97 a1, b1 = self.philam 

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

99 if wrap: 

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

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

102 else: # unrollPI shortcut 

103 db = b2 - b1 

104 return a1, b1, a2, b2, db 

105 

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

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

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

109 end point. 

110 

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

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

113 from the start point to the point where the perpendicular 

114 crosses the line. 

115 

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

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

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

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

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

121 

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

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

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

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

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

127 B{C{start}} point. 

128 

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

130 

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

132 

133 @example: 

134 

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

136 

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

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

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

140 ''' 

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

142 cx = cos(x) 

143 return _0_0 if isnear0(cx) else \ 

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

145 

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

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

148 ''' 

149 s = self.others(start=start) 

150 e = self.others(end=end) 

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

152 

153 r = Radius_(radius) 

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

155 

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

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

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

159 return r, x, -b 

160 

161 @deprecated_method 

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

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

164 ''' 

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

166 

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

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

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

170 

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

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

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

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

175 

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

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

178 ''' 

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

180 sa, ca, sa1, ca1, \ 

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

182 sa1 *= ca2 * ca 

183 

184 x = sa1 * sdb 

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

186 z = ca1 * sdb * ca2 * sa 

187 

188 h = hypot(x, y) 

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

190 return None # great circle doesn't reach latitude 

191 

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

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

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

195 

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

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

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

199 

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

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

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

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

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

205 

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

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

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

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

210 

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

212 

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

214 

215 @example: 

216 

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

218 

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

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

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

222 ''' 

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

224 return _r2m(x, radius) 

225 

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

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

228 travelled the given distance on the given initial bearing. 

229 

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

231 B{C{radius}}). 

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

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

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

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

236 

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

238 

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

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

241 

242 @example: 

243 

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

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

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

247 ''' 

248 a, b = self.philam 

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

250 

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

252 h = self._heigHt(height) 

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

254 

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

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

257 

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

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

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

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

262 

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

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

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

266 

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

268 

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

270 

271 @example: 

272 

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

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

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

276 ''' 

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

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

279 

280# @Property_RO 

281# def Ecef(self): 

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

283# ''' 

284# return _MODS.ecef.EcefKarney 

285 

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

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

288 on the given initial bearing from this point. 

289 

290 Direction of vector is such that initial bearing vector 

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

292 

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

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

295 overriding the default L{Vector3d}. 

296 @kwarg Vector_kwds: Optional, additional keyword argunents 

297 for B{C{Vector}}. 

298 

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

300 

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

302 

303 @example: 

304 

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

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

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

308 ''' 

309 a, b = self.philam 

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

311 

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

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

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

315 

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

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

318 to an other point. 

319 

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

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

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

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

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

325 C{sphericalNvector.LatLon.initialBearingTo}. 

326 

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

328 

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

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

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

332 

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

334 

335 @example: 

336 

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

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

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

340 ''' 

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

342 # XXX behavior like sphericalNvector.LatLon.initialBearingTo 

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

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

345 

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

347 

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

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

350 and an other point. 

351 

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

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

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

355 @kwarg height: Optional height, overriding the intermediate 

356 height (C{meter}). 

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

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

359 

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

361 

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

363 

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

365 

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

367 

368 @example: 

369 

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

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

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

373 ''' 

374 p = self 

375 f = Scalar(fraction=fraction) 

376 if not isnear0(f): 

377 p = p.others(other) 

378 if wrap: 

379 p = _Wrap.point(p) 

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

381 a1, b1 = self.philam 

382 a2, b2 = p.philam 

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

384 r = vincentys_(a2, a1, db) 

385 sr = sin(r) 

386 if isnon0(sr): 

387 sa1, ca1, sa2, ca2, \ 

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

389 

390 t = f * r 

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

392 b = sin( t) # / sr superflous 

393 

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

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

396 z = a * sa1 + b * sa2 

397 

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

399 b = atan2(y, x) 

400 

401 else: # PYCHOK no cover 

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

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

404 

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

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

407 return p 

408 

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

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

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

412 

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

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

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

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

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

418 C{degrees360}). 

419 @kwarg height: Optional height for intersection point, 

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

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

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

423 

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

425 intersection point might be the L{antipode} to 

426 the returned result. 

427 

428 @raise IntersectionError: Ambiguous or infinite intersection 

429 or colinear, parallel or otherwise 

430 non-intersecting lines. 

431 

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

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

434 

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

436 

437 @example: 

438 

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

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

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

442 ''' 

443 try: 

444 s2 = self.others(other) 

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

446 LatLon=self.classof) 

447 except (TypeError, ValueError) as x: 

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

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

450 

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

452 height=None, wrap=True): 

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

454 by a center point and radius. 

455 

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

457 see B{C{radius}}). 

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

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

460 see B{C{radius}}). 

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

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

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

464 B{C{radius}}). 

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

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

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

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

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

470 

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

472 instance. For abutting circles, both intersection 

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

474 

475 @raise IntersectionError: Concentric, antipodal, invalid or 

476 non-intersecting circles. 

477 

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

479 

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

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

482 ''' 

483 try: 

484 c2 = self.others(other) 

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

486 height=height, wrap=wrap, 

487 LatLon=self.classof) 

488 except (TypeError, ValueError) as x: 

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

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

491 

492 @deprecated_method 

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

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

495 return self.isenclosedBy(points) 

496 

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

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

499 

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

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

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

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

504 

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

506 composite, C{False} otherwise. 

507 

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

509 

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

511 

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

513 

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

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

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

517 ''' 

518 if _MODS.booleans.isBoolean(points): 

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

520 

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

522 n0 = self._N_vector 

523 

524 v2 = Ps[0]._N_vector 

525 p1 = Ps[1] 

526 v1 = p1._N_vector 

527 # check whether this point on same side of all 

528 # polygon edges (to the left or right depending 

529 # on the anti-/clockwise polygon direction) 

530 gc1 = v2.cross(v1) 

531 t0 = gc1.angleTo(n0) > PI_2 

532 s0 = None 

533 # get great-circle vector for each edge 

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

535 if wrap and not Ps.looped: 

536 p2 = _unrollon(p1, p2) 

537 p1 = p2 

538 v2 = p2._N_vector 

539 gc = v1.cross(v2) 

540 t = gc.angleTo(n0) > PI_2 

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

542 return False # outside 

543 

544 # check for convex polygon: angle between 

545 # gc vectors, signed by direction of n0 

546 # (otherwise the test above is not reliable) 

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

548 if s != s0: 

549 if s0 is None: 

550 s0 = s 

551 else: 

552 t = _Fmt.SQUARE(points=i) 

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

554 gc1, v1 = gc, v2 

555 

556 return True # inside 

557 

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

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

560 

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

562 @kwarg height: Optional height for midpoint, overriding 

563 the mean height (C{meter}). 

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

565 may be negative or greater than 1.0. 

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

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

568 

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

570 

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

572 

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

574 

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

576 

577 @example: 

578 

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

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

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

582 ''' 

583 if fraction is _0_5: 

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

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

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

587 

588 x = ca2 * cdb + ca1 

589 y = ca2 * sdb 

590 

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

592 b += atan2(y, x) 

593 

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

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

596 else: 

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

598 return r 

599 

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

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

602 

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

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

605 

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

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

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

609 @kwarg wrap_adjust_limit: Optional keyword arguments for functions 

610 L{sphericalTrigonometry.nearestOn3} and 

611 L{pygeodesy.equirectangular_}, 

612 

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

614 

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

616 see function L{pygeodesy.equirectangular_}. 

617 

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

619 is not (yet) supported. 

620 

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

622 

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

624 

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

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

627 ''' 

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

629 w = _xkwds_pop(wrap_adjust_limit, within=True) 

630 if not w: 

631 notImplemented(self, within=w) 

632 

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

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

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

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

637# return point1 

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

639# if isnear0(d): 

640# return point1 # or point2 

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

642# return point2 

643# f = a / d 

644# if within: 

645# if f > EPS1: 

646# return point2 

647# elif f < EPS: 

648# return point1 

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

650 

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

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

653 **wrap_adjust_limit)[0] 

654 

655 @deprecated_method 

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

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

658 

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

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

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

662 units of B{C{radius}}. 

663 ''' 

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

665 return r.closest, r.distance 

666 

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

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

669 

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

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

672 

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

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

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

676 @kwarg wrap_adjust_limit: Optional keyword arguments for function 

677 L{sphericalTrigonometry.nearestOn3} and 

678 L{pygeodesy.equirectangular_}, 

679 

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

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

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

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

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

685 function L{pygeodesy.compassAngle}. 

686 

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

688 see function L{pygeodesy.equirectangular_}. 

689 

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

691 

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

693 

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

695 

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

697 and L{pygeodesy.nearestOn5}. 

698 ''' 

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

700 LatLon=self.classof, **wrap_adjust_limit) 

701 

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

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

704 coordinates. 

705 

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

707 and other keyword arguments, ignored 

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

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

710 this L{Cartesian} class or specify 

711 C{B{Cartesian}=None}. 

712 

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

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

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

716 

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

718 ''' 

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

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

721 

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

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

724 

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

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

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

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

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

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

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

732 

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

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

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

736 

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

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

739 ''' 

740 B = self.others(otherB=otherB) 

741 C = self.others(otherC=otherC) 

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

743 

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

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

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

747 

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

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

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

751 of three corresponding circles. 

752 

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

754 as B{C{radius}}). 

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

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

757 B{C{radius}}). 

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

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

760 B{C{radius}}). 

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

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

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

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

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

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

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

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

769 

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

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

772 the corresponding trilaterated points C{minPoint} and 

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

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

775 

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

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

778 

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

780 respectively largest I{radial} overlap found. 

781 

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

783 nearest respectively farthest intersection margin. 

784 

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

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

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

788 largest B{C{distance#}}. 

789 

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

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

792 no intersection or all (near-)concentric for 

793 C{B{area}=False}. 

794 

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

796 

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

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

799 or B{C{radius}}. 

800 ''' 

801 return _trilaterate5(self, distance1, 

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

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

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

805 

806 

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

808 

809 

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

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

812 (with the pointsjoined by great circle arcs). 

813 

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

815 or L{BooleanGH}). 

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

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

818 or C{None}. 

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

820 (C{bool}). 

821 

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

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

824 

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

826 

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

828 

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

830 

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

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

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

834 

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

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

837 L{ellipsoidalKarney.areaOf}. 

838 

839 @example: 

840 

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

842 >>> areaOf(b) # 8666058750.718977 

843 

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

845 >>> areaOf(c) # 6.18e9 

846 ''' 

847 if _MODS.booleans.isBoolean(points): 

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

849 

850 _at2, _t_2 = atan2, tan_2 

851 _un, _w180 = unrollPI, wrap180 

852 

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

854 p1 = p2 = Ps[0] 

855 a1, b1 = p1.philam 

856 ta1, z1 = _t_2(a1), None 

857 

858 A = Fsum() # mean phi 

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

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

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

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

863 D = Fsum() 

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

865 a2, b2 = p2.philam 

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

867 A += a2 

868 ta2 = _t_2(a2) 

869 tdb = _t_2(db, points=i) 

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

871 _1_0 + ta1 * ta2) 

872 ta1, b1 = ta2, b2 

873 

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

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

876 if z1 is not None: 

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

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

879 p1, z1 = p2, z2 

880 

881 R = abs(R * _2_0) 

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

883 R = abs(R - PI2) 

884 if radius: 

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

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

887 return float(R) 

888 

889 

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

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

892 

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

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

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

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

897 

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

899 ''' 

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

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

902 ca *= sr 

903 

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

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

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

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

908 

909 

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

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

912 # and similar logic in .ellipsoidalBaseDI._intersect3 

913 a1, b1 = s.philam 

914 

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

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

917 else: # must be a point 

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

919 hs.append(end.height) 

920 a2, b2 = end.philam 

921 if wrap: 

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

923 

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

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

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

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

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

929 -(b1 + b2) * _0_5) 

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

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

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

933 cb12 * cb21 + sb12 * sb21, 

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

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

936 

937 

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

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

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

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

942 

943 

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

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

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

947 

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

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

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

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

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

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

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

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

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

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

958 circle} methods. 

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

960 conventionally) or C{None}. 

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

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

963 

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

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

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

967 or I{normalized}. 

968 

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

970 

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

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

973 

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

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

976 ''' 

977 c = _T00.others(center=center) 

978 p = _T00.others(point=point) 

979 try: 

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

981 height=height, wrap=wrap) 

982 except (TypeError, ValueError) as x: 

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

984 

985 

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

987 LatLon=None, **LatLon_kwds): 

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

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

990 

991 s1, s2 = start1, start2 

992 if wrap: 

993 s2 = _Wrap.point(s2) 

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

995 

996 a1, b1 = s1.philam 

997 a2, b2 = s2.philam 

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

999 r12 = vincentys_(a2, a1, db) 

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

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

1002 

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

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

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

1006 

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

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

1009 raise IntersectionError(_parallel_) 

1010 # handle domain error for equivalent longitudes, 

1011 # see also functions asin_safe and acos_safe at 

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

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

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

1015 if sin(db) > 0: 

1016 t21 = PI2 - t21 

1017 else: 

1018 t12 = PI2 - t12 

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

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

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

1022 raise IntersectionError(_infinite_) 

1023 sx3 = sx1 * sx2 

1024# XXX if sx3 < 0: 

1025# XXX raise ValueError(_ambiguous_) 

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

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

1028 

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

1030 # like .ellipsoidalBaseDI,_intersect3, if this intersection 

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

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

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

1034 

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

1036 _N_vector_ = _MODS.nvectorBase._N_vector_ 

1037 

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

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

1040 x = x1.cross(x2) 

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

1042 raise IntersectionError(_colinear_) 

1043 a, b = x.philam 

1044 # choose intersection similar to sphericalNvector 

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

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

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

1048 

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

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

1051 intersection, LatLon, LatLon_kwds) 

1052 

1053 

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

1055 LatLon=LatLon, **LatLon_kwds): 

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

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

1058 

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

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

1061 the initial bearing at the first start point 

1062 (compass C{degrees360}). 

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

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

1065 the initial bearing at the second start point 

1066 (compass C{degrees360}). 

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

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

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

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

1071 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1075 

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

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

1078 height)}. An alternate intersection point might 

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

1080 

1081 @raise IntersectionError: Ambiguous or infinite intersection 

1082 or colinear, parallel or otherwise 

1083 non-intersecting lines. 

1084 

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

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

1087 

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

1089 

1090 @example: 

1091 

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

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

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

1095 ''' 

1096 s1 = _T00.others(start1=start1) 

1097 s2 = _T00.others(start2=start2) 

1098 try: 

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

1100 LatLon=LatLon, **LatLon_kwds) 

1101 except (TypeError, ValueError) as x: 

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

1103 

1104 

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

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

1107 LatLon=LatLon, **LatLon_kwds): 

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

1109 by a center point and a radius. 

1110 

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

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

1113 see B{C{radius}}). 

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

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

1116 see B{C{radius}}). 

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

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

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

1120 B{C{radius}}). 

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

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

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

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

1125 (C{bool}). 

1126 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1130 

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

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

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

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

1135 

1136 @raise IntersectionError: Concentric, antipodal, invalid or 

1137 non-intersecting circles. 

1138 

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

1140 

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

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

1143 

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

1145 

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

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

1148 ''' 

1149 c1 = _T00.others(center1=center1) 

1150 c2 = _T00.others(center2=center2) 

1151 try: 

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

1153 height=height, wrap=wrap, 

1154 LatLon=LatLon, **LatLon_kwds) 

1155 except (TypeError, ValueError) as x: 

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

1157 center2=center2, rad2=rad2, wrap=wrap) 

1158 

1159 

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

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

1162 LatLon=LatLon, **LatLon_kwds): 

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

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

1165 

1166 def _dest3(bearing, h): 

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

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

1169 intersections2, LatLon, LatLon_kwds) 

1170 

1171 a1, b1 = c1.philam 

1172 a2, b2 = c2.philam 

1173 if wrap: 

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

1175 

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

1177 if f: # swapped radii, swap centers 

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

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

1180 

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

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

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

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

1185 

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

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

1188 if r < _0_0: 

1189 raise _ValueError(eps=r) 

1190 

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

1192 if x > max(r, EPS): 

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

1194 x = sd * sr1 

1195 if isnear0(x): 

1196 raise _ValueError(_invalid_) 

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

1198 

1199 elif x < r: # PYCHOK no cover 

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

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

1202 

1203 if height is None: # "radical height" 

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

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

1206 else: 

1207 h = Height(height) 

1208 

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

1210 if x < EPS4: # externally ... 

1211 r = _dest3(b, h) 

1212 elif x > _PI_EPS4: # internally ... 

1213 r = _dest3(b + PI, h) 

1214 else: 

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

1216 return r, r # ... abutting circles 

1217 

1218 

1219@deprecated_function 

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

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

1222 ''' 

1223 return ispolar(points, wrap=wrap) 

1224 

1225 

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

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

1228 ''' 

1229 n = where.__name__ 

1230 if LatLon is None: 

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

1232 else: 

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

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

1235 return r 

1236 

1237 

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

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

1240 

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

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

1243 height (C{meter}). 

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

1245 (C{bool}). 

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

1247 or C{None}. 

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

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

1250 

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

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

1253 is C{None}. 

1254 

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

1256 

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

1258 ''' 

1259 def _N_vs(ps, w): 

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

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

1262 yield p._N_vector 

1263 

1264 m = _MODS.nvectorBase 

1265 # geographic, vectorial mean 

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

1267 lat, lon, h = n.latlonheight 

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

1269 

1270 

1271@deprecated_function 

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

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

1274 

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

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

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

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

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

1280 ''' 

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

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

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

1284 return ll, d 

1285 

1286 

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

1288 limit=9, **LatLon_and_kwds): 

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

1290 

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

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

1293 

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

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

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

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

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

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

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

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

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

1303 spherical earth radius L{R_KM}). 

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

1305 or C{None}. 

1306 @kwarg options: Optional keyword arguments for function 

1307 L{pygeodesy.equirectangular_}. 

1308 

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

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

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

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

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

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

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

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

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

1318 

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

1320 see function L{pygeodesy.equirectangular_}. 

1321 

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

1323 

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

1325 

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

1327 

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

1329 ''' 

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

1331 adjust=adjust, limit=limit) 

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

1333 h = t.height 

1334 n = nearestOn3.__name__ 

1335 

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

1337 LL = _xkwds_pop(kwds, LatLon=LatLon) 

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

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

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

1341 

1342 

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

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

1345 (with great circle arcs joining the points). 

1346 

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

1348 or L{BooleanGH}). 

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

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

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

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

1353 

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

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

1356 

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

1358 

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

1360 

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

1362 C{B{points}} a composite. 

1363 

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

1365 

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

1367 and L{ellipsoidalKarney.perimeterOf}. 

1368 ''' 

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

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

1371 a1, b1 = Ps[0].philam 

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

1373 a2, b2 = p.philam 

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

1375 yield vincentys_(a2, a1, db) 

1376 a1, b1 = a2, b2 

1377 

1378 if _MODS.booleans.isBoolean(points): 

1379 if not closed: 

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

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

1382 else: 

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

1384 return _r2m(r, radius) 

1385 

1386 

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

1388 excess=excessAbc_, 

1389 wrap=False): 

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

1391 

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

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

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

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

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

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

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

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

1400 or C{None}. 

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

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

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

1404 longitudes (C{bool}). 

1405 

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

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

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

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

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

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

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

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

1414 ''' 

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

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

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

1418 excess=excess, wrap=wrap) 

1419 return _t7Tuple(t, radius) 

1420 

1421 

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

1423 wrap=False): 

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

1425 excess} of a (spherical) triangle. 

1426 

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

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

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

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

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

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

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

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

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

1436 longitudinal deltas (C{bool}). 

1437 

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

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

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

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

1442 ''' 

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

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

1445 a = vincentys_(phiC, phiB, d) 

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

1447 

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

1449 s = sb * sc 

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

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

1452 

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

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

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

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

1457 

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

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

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

1461 

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

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

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

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

1466 

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

1468 

1469 

1470def _t7Tuple(t, radius): 

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

1472 ''' 

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

1474 r = radius if isscalar(radius) else \ 

1475 _ellipsoidal_datum(radius).ellipsoid.Rmean 

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

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

1478 B, (r * t.b), 

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

1480 return t 

1481 

1482 

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

1484 areaOf, # functions 

1485 intersecant2, intersection, intersections2, ispolar, 

1486 isPoleEnclosedBy, # DEPRECATED, use ispolar 

1487 meanOf, 

1488 nearestOn2, nearestOn3, 

1489 perimeterOf, 

1490 sumOf, # XXX == vector3d.sumOf 

1491 triangle7, triangle8_) 

1492 

1493# **) MIT License 

1494# 

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

1496# 

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

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

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

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

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

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

1503# 

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

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

1506# 

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

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

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

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

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

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

1513# OTHER DEALINGS IN THE SOFTWARE.