Coverage for pygeodesy/sphericalTrigonometry.py: 94%

391 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-07 07:28 -0400

1 

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

3 

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

5 

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

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

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

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

10 

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

12spherical trigonometry, transcoded from JavaScript originals by 

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

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

15''' 

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

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

18 

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

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

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

22 _1_0, _2_0, _90_0 

23from pygeodesy.datums import _ellipsoidal_datum, _mean_radius 

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

25 _ValueError, IntersectionError, _xError, \ 

26 _xkwds, _xkwds_get, _xkwds_pop 

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

28from pygeodesy.fsums import Fsum, fsum, fsumf_ 

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

30 excessGirard_, excessLHuilier_, opposing_, _radical2, \ 

31 vincentys_ 

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

33 _concentric_, _convex_, _end_, _infinite_, \ 

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

35 _point_, _SPACE_, _too_ 

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

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

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

39from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, \ 

40 NearestOn3Tuple, Triangle7Tuple, \ 

41 Triangle8Tuple 

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

43 notImplemented, Fmt as _Fmt # XXX shadowed 

44from pygeodesy.props import deprecated_function, deprecated_method 

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

46 LatLonSphericalBase, _rads3, _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.07.01' 

58 

59_parallel_ = 'parallel' 

60 

61_PI_EPS4 = PI - EPS4 

62if _PI_EPS4 >= PI: 

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

64 

65 

66class Cartesian(CartesianSphericalBase): 

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

68 spherical, geodetic L{LatLon}. 

69 ''' 

70 

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

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

73 

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

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

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

77 

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

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

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

81 

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

83 ''' 

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

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

86 

87 

88class LatLon(LatLonSphericalBase): 

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

90 

91 @example: 

92 

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

94 ''' 

95 

96 def _ab1_ab2_db5(self, other, wrap): 

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

98 ''' 

99 a1, b1 = self.philam 

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

101 if wrap: 

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

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

104 else: # unrollPI shortcut 

105 db = b2 - b1 

106 return a1, b1, a2, b2, db 

107 

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

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

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

111 end point. 

112 

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

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

115 from the start point to the point where the perpendicular 

116 crosses the line. 

117 

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

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

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

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

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

123 

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

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

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

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

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

129 B{C{start}} point. 

130 

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

132 

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

134 

135 @example: 

136 

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

138 

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

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

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

142 ''' 

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

144 cx = cos(x) 

145 return _0_0 if isnear0(cx) else \ 

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

147 

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

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

150 ''' 

151 s = self.others(start=start) 

152 e = self.others(end=end) 

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

154 

155 r = Radius_(radius) 

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

157 

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

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

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

161 return r, x, -b 

162 

163 @deprecated_method 

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

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

166 ''' 

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

168 

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

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

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

172 

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

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

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

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

177 

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

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

180 ''' 

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

182 sa, ca, sa1, ca1, \ 

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

184 sa1 *= ca2 * ca 

185 

186 x = sa1 * sdb 

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

188 z = ca1 * sdb * ca2 * sa 

189 

190 h = hypot(x, y) 

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

192 return None # great circle doesn't reach latitude 

193 

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

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

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

197 

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

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

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

201 

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

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

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

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

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

207 

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

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

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

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

212 

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

214 

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

216 

217 @example: 

218 

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

220 

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

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

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

224 ''' 

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

226 return _r2m(x, radius) 

227 

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

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

230 travelled the given distance on the given initial bearing. 

231 

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

233 B{C{radius}}). 

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

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

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

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

238 

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

240 

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

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

243 

244 @example: 

245 

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

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

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

249 ''' 

250 a, b = self.philam 

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

252 

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

254 h = self._heigHt(height) 

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

256 

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

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

259 

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

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

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

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

264 

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

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

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

268 

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

270 

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

272 

273 @example: 

274 

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

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

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

278 ''' 

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

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

281 

282# @Property_RO 

283# def Ecef(self): 

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

285# ''' 

286# return _MODS.ecef.EcefKarney 

287 

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

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

290 on the given initial bearing from this point. 

291 

292 Direction of vector is such that initial bearing vector 

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

294 

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

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

297 overriding the default L{Vector3d}. 

298 @kwarg Vector_kwds: Optional, additional keyword argunents 

299 for B{C{Vector}}. 

300 

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

302 

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

304 

305 @example: 

306 

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

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

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

310 ''' 

311 a, b = self.philam 

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

313 

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

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

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

317 

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

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

320 to an other point. 

321 

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

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

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

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

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

327 C{sphericalNvector.LatLon.initialBearingTo}. 

328 

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

330 

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

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

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

334 

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

336 

337 @example: 

338 

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

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

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

342 ''' 

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

344 # XXX behavior like sphericalNvector.LatLon.initialBearingTo 

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

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

347 

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

349 

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

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

352 and an other point. 

353 

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

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

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

357 @kwarg height: Optional height, overriding the intermediate 

358 height (C{meter}). 

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

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

361 

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

363 

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

365 

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

367 

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

369 

370 @example: 

371 

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

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

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

375 ''' 

376 p = self 

377 f = Scalar(fraction=fraction) 

378 if not isnear0(f): 

379 p = p.others(other) 

380 if wrap: 

381 p = _Wrap.point(p) 

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

383 a1, b1 = self.philam 

384 a2, b2 = p.philam 

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

386 r = vincentys_(a2, a1, db) 

387 sr = sin(r) 

388 if isnon0(sr): 

389 sa1, ca1, sa2, ca2, \ 

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

391 

392 t = f * r 

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

394 b = sin( t) # / sr superflous 

395 

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

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

398 z = a * sa1 + b * sa2 

399 

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

401 b = atan2(y, x) 

402 

403 else: # PYCHOK no cover 

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

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

406 

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

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

409 return p 

410 

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

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

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

414 

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

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

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

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

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

420 C{degrees360}). 

421 @kwarg height: Optional height for intersection point, 

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

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

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

425 

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

427 intersection point might be the L{antipode} to 

428 the returned result. 

429 

430 @raise IntersectionError: Ambiguous or infinite intersection 

431 or colinear, parallel or otherwise 

432 non-intersecting lines. 

433 

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

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

436 

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

438 

439 @example: 

440 

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

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

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

444 ''' 

445 try: 

446 s2 = self.others(other) 

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

448 LatLon=self.classof) 

449 except (TypeError, ValueError) as x: 

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

451 other=other, end2=end2, wrap=wrap) 

452 

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

454 height=None, wrap=True): 

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

456 by a center point and radius. 

457 

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

459 see B{C{radius}}). 

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

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

462 see B{C{radius}}). 

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

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

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

466 B{C{radius}}). 

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

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

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

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

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

472 

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

474 instance. For abutting circles, both intersection 

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

476 

477 @raise IntersectionError: Concentric, antipodal, invalid or 

478 non-intersecting circles. 

479 

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

481 

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

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

484 ''' 

485 try: 

486 c2 = self.others(other) 

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

488 height=height, wrap=wrap, 

489 LatLon=self.classof) 

490 except (TypeError, ValueError) as x: 

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

492 other=other, rad2=rad2, wrap=wrap) 

493 

494 @deprecated_method 

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

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

497 return self.isenclosedBy(points) 

498 

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

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

501 

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

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

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

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

506 

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

508 composite, C{False} otherwise. 

509 

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

511 

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

513 

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

515 

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

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

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

519 ''' 

520 if _MODS.booleans.isBoolean(points): 

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

522 

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

524 n0 = self._N_vector 

525 

526 v2 = Ps[0]._N_vector 

527 p1 = Ps[1] 

528 v1 = p1._N_vector 

529 # check whether this point on same side of all 

530 # polygon edges (to the left or right depending 

531 # on the anti-/clockwise polygon direction) 

532 gc1 = v2.cross(v1) 

533 t0 = gc1.angleTo(n0) > PI_2 

534 s0 = None 

535 # get great-circle vector for each edge 

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

537 if wrap and not Ps.looped: 

538 p2 = _unrollon(p1, p2) 

539 p1 = p2 

540 v2 = p2._N_vector 

541 gc = v1.cross(v2) 

542 t = gc.angleTo(n0) > PI_2 

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

544 return False # outside 

545 

546 # check for convex polygon: angle between 

547 # gc vectors, signed by direction of n0 

548 # (otherwise the test above is not reliable) 

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

550 if s != s0: 

551 if s0 is None: 

552 s0 = s 

553 else: 

554 t = _Fmt.SQUARE(points=i) 

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

556 gc1, v1 = gc, v2 

557 

558 return True # inside 

559 

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

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

562 

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

564 @kwarg height: Optional height for midpoint, overriding 

565 the mean height (C{meter}). 

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

567 may be negative or greater than 1.0. 

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

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

570 

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

572 

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

574 

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

576 

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

578 

579 @example: 

580 

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

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

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

584 ''' 

585 if fraction is _0_5: 

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

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

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

589 

590 x = ca2 * cdb + ca1 

591 y = ca2 * sdb 

592 

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

594 b += atan2(y, x) 

595 

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

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

598 else: 

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

600 return r 

601 

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

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

604 

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

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

607 

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

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

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

611 @kwarg wrap_adjust_limit: Optional keyword arguments for functions 

612 L{sphericalTrigonometry.nearestOn3} and 

613 L{pygeodesy.equirectangular_}, 

614 

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

616 

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

618 see function L{pygeodesy.equirectangular_}. 

619 

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

621 is not (yet) supported. 

622 

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

624 

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

626 

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

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

629 ''' 

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

631 w = _xkwds_pop(wrap_adjust_limit, within=True) 

632 if not w: 

633 notImplemented(self, within=w) 

634 

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

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

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

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

639# return point1 

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

641# if isnear0(d): 

642# return point1 # or point2 

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

644# return point2 

645# f = a / d 

646# if within: 

647# if f > EPS1: 

648# return point2 

649# elif f < EPS: 

650# return point1 

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

652 

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

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

655 **wrap_adjust_limit)[0] 

656 

657 @deprecated_method 

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

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

660 

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

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

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

664 units of B{C{radius}}. 

665 ''' 

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

667 return r.closest, r.distance 

668 

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

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

671 

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

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

674 

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

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

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

678 @kwarg wrap_adjust_limit: Optional keyword arguments for function 

679 L{sphericalTrigonometry.nearestOn3} and 

680 L{pygeodesy.equirectangular_}, 

681 

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

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

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

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

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

687 function L{pygeodesy.compassAngle}. 

688 

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

690 see function L{pygeodesy.equirectangular_}. 

691 

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

693 

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

695 

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

697 

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

699 and L{pygeodesy.nearestOn5}. 

700 ''' 

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

702 LatLon=self.classof, **wrap_adjust_limit) 

703 

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

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

706 coordinates. 

707 

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

709 and other keyword arguments, ignored 

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

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

712 this L{Cartesian} class or specify 

713 C{B{Cartesian}=None}. 

714 

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

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

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

718 

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

720 ''' 

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

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

723 

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

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

726 

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

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

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

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

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

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

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

734 

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

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

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

738 

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

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

741 ''' 

742 B = self.others(otherB=otherB) 

743 C = self.others(otherC=otherC) 

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

745 

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

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

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

749 

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

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

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

753 of three corresponding circles. 

754 

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

756 as B{C{radius}}). 

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

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

759 B{C{radius}}). 

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

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

762 B{C{radius}}). 

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

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

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

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

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

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

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

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

771 

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

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

774 the corresponding trilaterated points C{minPoint} and 

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

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

777 

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

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

780 

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

782 respectively largest I{radial} overlap found. 

783 

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

785 nearest respectively farthest intersection margin. 

786 

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

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

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

790 largest B{C{distance#}}. 

791 

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

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

794 no intersection or all (near-)concentric for 

795 C{B{area}=False}. 

796 

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

798 

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

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

801 or B{C{radius}}. 

802 ''' 

803 return _trilaterate5(self, distance1, 

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

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

806 area=area, radius=radius, eps=eps, wrap=wrap) 

807 

808 

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

810 

811 

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

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

814 (with the pointsjoined by great circle arcs). 

815 

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

817 or L{BooleanGH}). 

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

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

820 or C{None}. 

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

822 (C{bool}). 

823 

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

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

826 

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

828 

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

830 

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

832 

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

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

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

836 

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

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

839 L{ellipsoidalKarney.areaOf}. 

840 

841 @example: 

842 

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

844 >>> areaOf(b) # 8666058750.718977 

845 

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

847 >>> areaOf(c) # 6.18e9 

848 ''' 

849 if _MODS.booleans.isBoolean(points): 

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

851 

852 _at2, _t_2 = atan2, tan_2 

853 _un, _w180 = unrollPI, wrap180 

854 

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

856 p1 = p2 = Ps[0] 

857 a1, b1 = p1.philam 

858 ta1, z1 = _t_2(a1), None 

859 

860 A = Fsum() # mean phi 

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

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

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

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

865 D = Fsum() 

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

867 a2, b2 = p2.philam 

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

869 A += a2 

870 ta2 = _t_2(a2) 

871 tdb = _t_2(db, points=i) 

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

873 _1_0 + ta1 * ta2) 

874 ta1, b1 = ta2, b2 

875 

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

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

878 if z1 is not None: 

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

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

881 p1, z1 = p2, z2 

882 

883 R = abs(R * _2_0) 

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

885 R = abs(R - PI2) 

886 if radius: 

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

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

889 return float(R) 

890 

891 

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

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

894 

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

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

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

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

899 

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

901 ''' 

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

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

904 ca *= sr 

905 

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

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

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

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

910 

911 

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

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

914 # and similar logic in .ellipsoidalBaseDI._intersect3 

915 a1, b1 = s.philam 

916 

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

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

919 else: # must be a point 

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

921 hs.append(end.height) 

922 a2, b2 = end.philam 

923 if wrap: 

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

925 

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

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

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

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

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

931 -(b1 + b2) * _0_5) 

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

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

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

935 cb12 * cb21 + sb12 * sb21, 

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

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

938 

939 

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

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

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

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

944 

945 

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

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

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

949 

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

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

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

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

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

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

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

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

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

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

960 circle} methods. 

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

962 conventionally) or C{None}. 

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

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

965 

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

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

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

969 or I{normalized}. 

970 

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

972 

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

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

975 

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

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

978 ''' 

979 c = _T00.others(center=center) 

980 p = _T00.others(point=point) 

981 try: 

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

983 height=height, wrap=wrap) 

984 except (TypeError, ValueError) as x: 

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

986 

987 

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

989 LatLon=None, **LatLon_kwds): 

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

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

992 

993 s1, s2 = start1, start2 

994 if wrap: 

995 s2 = _Wrap.point(s2) 

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

997 

998 a1, b1 = s1.philam 

999 a2, b2 = s2.philam 

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

1001 r12 = vincentys_(a2, a1, db) 

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

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

1004 

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

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

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

1008 

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

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

1011 raise IntersectionError(_parallel_) 

1012 # handle domain error for equivalent longitudes, 

1013 # see also functions asin_safe and acos_safe at 

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

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

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

1017 if sin(db) > 0: 

1018 t12, t21 = t1, PI2 - t2 

1019 else: 

1020 t12, t21 = PI2 - t1, t2 

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

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

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

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

1025 raise IntersectionError(_infinite_) 

1026 sx3 = sx1 * sx2 

1027# XXX if sx3 < 0: 

1028# XXX raise ValueError(_ambiguous_) 

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

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

1031 

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

1033 # like .ellipsoidalBaseDI,_intersect3, if this intersection 

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

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

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

1037 

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

1039 _N_vector_ = _MODS.nvectorBase._N_vector_ 

1040 

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

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

1043 x = x1.cross(x2) 

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

1045 raise IntersectionError(_colinear_) 

1046 a, b = x.philam 

1047 # choose intersection similar to sphericalNvector 

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

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

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

1051 

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

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

1054 intersection, LatLon, LatLon_kwds) 

1055 

1056 

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

1058 LatLon=LatLon, **LatLon_kwds): 

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

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

1061 

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

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

1064 the initial bearing at the first start point 

1065 (compass C{degrees360}). 

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

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

1068 the initial bearing at the second start point 

1069 (compass C{degrees360}). 

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

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

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

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

1074 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1078 

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

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

1081 height)}. An alternate intersection point might 

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

1083 

1084 @raise IntersectionError: Ambiguous or infinite intersection 

1085 or colinear, parallel or otherwise 

1086 non-intersecting lines. 

1087 

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

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

1090 

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

1092 

1093 @example: 

1094 

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

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

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

1098 ''' 

1099 s1 = _T00.others(start1=start1) 

1100 s2 = _T00.others(start2=start2) 

1101 try: 

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

1103 LatLon=LatLon, **LatLon_kwds) 

1104 except (TypeError, ValueError) as x: 

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

1106 

1107 

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

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

1110 LatLon=LatLon, **LatLon_kwds): 

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

1112 by a center point and a radius. 

1113 

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

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

1116 see B{C{radius}}). 

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

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

1119 see B{C{radius}}). 

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

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

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

1123 B{C{radius}}). 

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

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

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

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

1128 (C{bool}). 

1129 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1133 

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

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

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

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

1138 

1139 @raise IntersectionError: Concentric, antipodal, invalid or 

1140 non-intersecting circles. 

1141 

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

1143 

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

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

1146 

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

1148 

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

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

1151 ''' 

1152 c1 = _T00.others(center1=center1) 

1153 c2 = _T00.others(center2=center2) 

1154 try: 

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

1156 height=height, wrap=wrap, 

1157 LatLon=LatLon, **LatLon_kwds) 

1158 except (TypeError, ValueError) as x: 

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

1160 center2=center2, rad2=rad2, wrap=wrap) 

1161 

1162 

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

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

1165 LatLon=LatLon, **LatLon_kwds): 

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

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

1168 

1169 def _dest3(bearing, h): 

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

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

1172 intersections2, LatLon, LatLon_kwds) 

1173 

1174 a1, b1 = c1.philam 

1175 a2, b2 = c2.philam 

1176 if wrap: 

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

1178 

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

1180 if f: # swapped radii, swap centers 

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

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

1183 

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

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

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

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

1188 

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

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

1191 if r < _0_0: 

1192 raise _ValueError(eps=r) 

1193 

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

1195 if x > max(r, EPS): 

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

1197 x = sd * sr1 

1198 if isnear0(x): 

1199 raise _ValueError(_invalid_) 

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

1201 

1202 elif x < r: # PYCHOK no cover 

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

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

1205 

1206 if height is None: # "radical height" 

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

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

1209 else: 

1210 h = Height(height) 

1211 

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

1213 if x < EPS4: # externally ... 

1214 r = _dest3(b, h) 

1215 elif x > _PI_EPS4: # internally ... 

1216 r = _dest3(b + PI, h) 

1217 else: 

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

1219 return r, r # ... abutting circles 

1220 

1221 

1222@deprecated_function 

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

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

1225 ''' 

1226 return ispolar(points, wrap=wrap) 

1227 

1228 

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

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

1231 ''' 

1232 n = where.__name__ 

1233 if LatLon is None: 

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

1235 else: 

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

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

1238 return r 

1239 

1240 

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

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

1243 

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

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

1246 height (C{meter}). 

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

1248 (C{bool}). 

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

1250 or C{None}. 

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

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

1253 

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

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

1256 is C{None}. 

1257 

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

1259 

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

1261 ''' 

1262 def _N_vs(ps, w): 

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

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

1265 yield p._N_vector 

1266 

1267 m = _MODS.nvectorBase 

1268 # geographic, vectorial mean 

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

1270 lat, lon, h = n.latlonheight 

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

1272 

1273 

1274@deprecated_function 

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

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

1277 

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

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

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

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

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

1283 ''' 

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

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

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

1287 return ll, d 

1288 

1289 

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

1291 limit=9, **LatLon_and_kwds): 

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

1293 

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

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

1296 

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

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

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

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

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

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

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

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

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

1306 spherical earth radius L{R_KM}). 

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

1308 or C{None}. 

1309 @kwarg options: Optional keyword arguments for function 

1310 L{pygeodesy.equirectangular_}. 

1311 

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

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

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

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

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

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

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

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

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

1321 

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

1323 see function L{pygeodesy.equirectangular_}. 

1324 

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

1326 

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

1328 

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

1330 

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

1332 ''' 

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

1334 adjust=adjust, limit=limit) 

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

1336 h = t.height 

1337 n = nearestOn3.__name__ 

1338 

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

1340 LL = _xkwds_pop(kwds, LatLon=LatLon) 

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

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

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

1344 

1345 

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

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

1348 (with great circle arcs joining the points). 

1349 

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

1351 or L{BooleanGH}). 

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

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

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

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

1356 

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

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

1359 

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

1361 

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

1363 

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

1365 C{B{points}} a composite. 

1366 

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

1368 

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

1370 and L{ellipsoidalKarney.perimeterOf}. 

1371 ''' 

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

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

1374 a1, b1 = Ps[0].philam 

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

1376 a2, b2 = p.philam 

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

1378 yield vincentys_(a2, a1, db) 

1379 a1, b1 = a2, b2 

1380 

1381 if _MODS.booleans.isBoolean(points): 

1382 if not closed: 

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

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

1385 else: 

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

1387 return _r2m(r, radius) 

1388 

1389 

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

1391 excess=excessAbc_, 

1392 wrap=False): 

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

1394 

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

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

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

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

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

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

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

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

1403 or C{None}. 

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

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

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

1407 longitudes (C{bool}). 

1408 

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

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

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

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

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

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

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

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

1417 ''' 

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

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

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

1421 excess=excess, wrap=wrap) 

1422 return _t7Tuple(t, radius) 

1423 

1424 

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

1426 wrap=False): 

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

1428 excess} of a (spherical) triangle. 

1429 

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

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

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

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

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

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

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

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

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

1439 longitudinal deltas (C{bool}). 

1440 

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

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

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

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

1445 ''' 

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

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

1448 a = vincentys_(phiC, phiB, d) 

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

1450 

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

1452 s = sb * sc 

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

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

1455 

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

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

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

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

1460 

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

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

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

1464 

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

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

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

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

1469 

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

1471 

1472 

1473def _t7Tuple(t, radius): 

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

1475 ''' 

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

1477 r = radius if isscalar(radius) else \ 

1478 _ellipsoidal_datum(radius).ellipsoid.Rmean 

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

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

1481 B, (r * t.b), 

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

1483 return t 

1484 

1485 

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

1487 areaOf, # functions 

1488 intersecant2, intersection, intersections2, ispolar, 

1489 isPoleEnclosedBy, # DEPRECATED, use ispolar 

1490 meanOf, 

1491 nearestOn2, nearestOn3, 

1492 perimeterOf, 

1493 sumOf, # XXX == vector3d.sumOf 

1494 triangle7, triangle8_) 

1495 

1496# **) MIT License 

1497# 

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

1499# 

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

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

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

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

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

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

1506# 

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

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

1509# 

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

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

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

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

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

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

1516# OTHER DEALINGS IN THE SOFTWARE.