Coverage for pygeodesy/sphericalTrigonometry.py: 94%

390 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-10-11 16:04 -0400

1 

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

3 

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

5 

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

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

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

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

10 

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

12spherical trigonometry, transcoded from JavaScript originals by 

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

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

15''' 

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

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

18 

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

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

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

22 _1_0, _2_0, _90_0 

23from pygeodesy.datums import _ellipsoidal_datum, _mean_radius 

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

25 _ValueError, IntersectionError, _xError, \ 

26 _xkwds, _xkwds_get, _xkwds_pop 

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

28from pygeodesy.fsums import Fsum, fsum, fsumf_ 

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

30 excessGirard_, excessLHuilier_, opposing_, _radical2, \ 

31 vincentys_ 

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

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

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

35 _SPACE_, _too_ 

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

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

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

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

40 Triangle7Tuple, Triangle8Tuple 

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

42 Fmt as _Fmt, notImplemented # XXX shadowed 

43from pygeodesy.props import deprecated_function, deprecated_method 

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

45 LatLonSphericalBase, _rads3, _r2m, _trilaterate5 

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

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

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

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

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

51from pygeodesy.vector3d import sumOf, Vector3d 

52 

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

54 

55__all__ = _ALL_LAZY.sphericalTrigonometry 

56__version__ = '23.09.29' 

57 

58_PI_EPS4 = PI - EPS4 

59if _PI_EPS4 >= PI: 

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

61 

62 

63class Cartesian(CartesianSphericalBase): 

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

65 spherical, geodetic L{LatLon}. 

66 ''' 

67 

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

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

70 

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

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

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

74 

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

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

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

78 

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

80 ''' 

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

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

83 

84 

85class LatLon(LatLonSphericalBase): 

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

87 

88 @example: 

89 

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

91 ''' 

92 

93 def _ab1_ab2_db5(self, other, wrap): 

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

95 ''' 

96 a1, b1 = self.philam 

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

98 if wrap: 

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

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

101 else: # unrollPI shortcut 

102 db = b2 - b1 

103 return a1, b1, a2, b2, db 

104 

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

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

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

108 end point. 

109 

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

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

112 from the start point to the point where the perpendicular 

113 crosses the line. 

114 

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

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

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

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

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

120 

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

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

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

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

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

126 B{C{start}} point. 

127 

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

129 

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

131 

132 @example: 

133 

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

135 

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

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

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

139 ''' 

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

141 cx = cos(x) 

142 return _0_0 if isnear0(cx) else \ 

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

144 

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

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

147 ''' 

148 s = self.others(start=start) 

149 e = self.others(end=end) 

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

151 

152 r = Radius_(radius) 

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

154 

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

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

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

158 return r, x, -b 

159 

160 @deprecated_method 

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

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

163 ''' 

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

165 

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

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

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

169 

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

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

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

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

174 

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

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

177 ''' 

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

179 sa, ca, sa1, ca1, \ 

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

181 sa1 *= ca2 * ca 

182 

183 x = sa1 * sdb 

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

185 z = ca1 * sdb * ca2 * sa 

186 

187 h = hypot(x, y) 

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

189 return None # great circle doesn't reach latitude 

190 

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

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

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

194 

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

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

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

198 

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

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

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

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

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

204 

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

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

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

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

209 

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

211 

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

213 

214 @example: 

215 

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

217 

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

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

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

221 ''' 

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

223 return _r2m(x, radius) 

224 

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

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

227 travelled the given distance on the given initial bearing. 

228 

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

230 B{C{radius}}). 

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

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

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

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

235 

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

237 

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

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

240 

241 @example: 

242 

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

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

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

246 ''' 

247 a, b = self.philam 

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

249 

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

251 h = self._heigHt(height) 

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

253 

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

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

256 

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

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

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

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

261 

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

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

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

265 

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

267 

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

269 

270 @example: 

271 

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

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

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

275 ''' 

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

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

278 

279# @Property_RO 

280# def Ecef(self): 

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

282# ''' 

283# return _MODS.ecef.EcefKarney 

284 

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

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

287 on the given initial bearing from this point. 

288 

289 Direction of vector is such that initial bearing vector 

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

291 

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

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

294 overriding the default L{Vector3d}. 

295 @kwarg Vector_kwds: Optional, additional keyword argunents 

296 for B{C{Vector}}. 

297 

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

299 

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

301 

302 @example: 

303 

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

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

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

307 ''' 

308 a, b = self.philam 

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

310 

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

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

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

314 

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

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

317 to an other point. 

318 

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

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

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

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

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

324 C{sphericalNvector.LatLon.initialBearingTo}. 

325 

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

327 

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

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

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

331 

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

333 

334 @example: 

335 

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

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

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

339 ''' 

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

341 # XXX behavior like sphericalNvector.LatLon.initialBearingTo 

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

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

344 

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

346 

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

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

349 and an other point. 

350 

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

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

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

354 @kwarg height: Optional height, overriding the intermediate 

355 height (C{meter}). 

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

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

358 

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

360 

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

362 

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

364 

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

366 

367 @example: 

368 

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

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

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

372 ''' 

373 p = self 

374 f = Scalar(fraction=fraction) 

375 if not isnear0(f): 

376 p = p.others(other) 

377 if wrap: 

378 p = _Wrap.point(p) 

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

380 a1, b1 = self.philam 

381 a2, b2 = p.philam 

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

383 r = vincentys_(a2, a1, db) 

384 sr = sin(r) 

385 if isnon0(sr): 

386 sa1, ca1, sa2, ca2, \ 

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

388 

389 t = f * r 

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

391 b = sin( t) # / sr superflous 

392 

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

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

395 z = a * sa1 + b * sa2 

396 

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

398 b = atan2d(y, x) 

399 

400 else: # PYCHOK no cover 

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

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

403 

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

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

406 return p 

407 

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

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

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

411 

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

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

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

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

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

417 C{degrees360}). 

418 @kwarg height: Optional height for intersection point, 

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

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

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

422 

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

424 intersection point might be the L{antipode} to 

425 the returned result. 

426 

427 @raise IntersectionError: Ambiguous or infinite intersection 

428 or colinear, parallel or otherwise 

429 non-intersecting lines. 

430 

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

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

433 

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

435 

436 @example: 

437 

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

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

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

441 ''' 

442 try: 

443 s2 = self.others(other) 

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

445 LatLon=self.classof) 

446 except (TypeError, ValueError) as x: 

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

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

449 

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

451 height=None, wrap=True): 

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

453 by a center point and radius. 

454 

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

456 see B{C{radius}}). 

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

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

459 see B{C{radius}}). 

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

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

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

463 B{C{radius}}). 

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

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

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

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

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

469 

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

471 instance. For abutting circles, both intersection 

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

473 

474 @raise IntersectionError: Concentric, antipodal, invalid or 

475 non-intersecting circles. 

476 

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

478 

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

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

481 ''' 

482 try: 

483 c2 = self.others(other) 

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

485 height=height, wrap=wrap, 

486 LatLon=self.classof) 

487 except (TypeError, ValueError) as x: 

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

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

490 

491 @deprecated_method 

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

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

494 return self.isenclosedBy(points) 

495 

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

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

498 

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

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

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

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

503 

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

505 composite, C{False} otherwise. 

506 

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

508 

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

510 

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

512 

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

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

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

516 ''' 

517 if _MODS.booleans.isBoolean(points): 

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

519 

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

521 n0 = self._N_vector 

522 

523 v2 = Ps[0]._N_vector 

524 p1 = Ps[1] 

525 v1 = p1._N_vector 

526 # check whether this point on same side of all 

527 # polygon edges (to the left or right depending 

528 # on the anti-/clockwise polygon direction) 

529 gc1 = v2.cross(v1) 

530 t0 = gc1.angleTo(n0) > PI_2 

531 s0 = None 

532 # get great-circle vector for each edge 

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

534 if wrap and not Ps.looped: 

535 p2 = _unrollon(p1, p2) 

536 p1 = p2 

537 v2 = p2._N_vector 

538 gc = v1.cross(v2) 

539 t = gc.angleTo(n0) > PI_2 

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

541 return False # outside 

542 

543 # check for convex polygon: angle between 

544 # gc vectors, signed by direction of n0 

545 # (otherwise the test above is not reliable) 

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

547 if s != s0: 

548 if s0 is None: 

549 s0 = s 

550 else: 

551 t = _Fmt.SQUARE(points=i) 

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

553 gc1, v1 = gc, v2 

554 

555 return True # inside 

556 

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

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

559 

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

561 @kwarg height: Optional height for midpoint, overriding 

562 the mean height (C{meter}). 

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

564 may be negative or greater than 1.0. 

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

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

567 

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

569 

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

571 

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

573 

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

575 

576 @example: 

577 

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

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

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

581 ''' 

582 if fraction is _0_5: 

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

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

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

586 

587 x = ca2 * cdb + ca1 

588 y = ca2 * sdb 

589 

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

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

592 

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

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

595 else: 

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

597 return r 

598 

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

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

601 

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

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

604 

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

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

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

608 @kwarg wrap_adjust_limit: Optional keyword arguments for functions 

609 L{sphericalTrigonometry.nearestOn3} and 

610 L{pygeodesy.equirectangular_}, 

611 

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

613 

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

615 see function L{pygeodesy.equirectangular_}. 

616 

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

618 is not (yet) supported. 

619 

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

621 

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

623 

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

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

626 ''' 

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

628 w = _xkwds_pop(wrap_adjust_limit, within=True) 

629 if not w: 

630 notImplemented(self, within=w) 

631 

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

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

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

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

636# return point1 

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

638# if isnear0(d): 

639# return point1 # or point2 

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

641# return point2 

642# f = a / d 

643# if within: 

644# if f > EPS1: 

645# return point2 

646# elif f < EPS: 

647# return point1 

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

649 

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

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

652 **wrap_adjust_limit)[0] 

653 

654 @deprecated_method 

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

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

657 

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

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

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

661 units of B{C{radius}}. 

662 ''' 

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

664 return r.closest, r.distance 

665 

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

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

668 

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

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

671 

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

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

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

675 @kwarg wrap_adjust_limit: Optional keyword arguments for function 

676 L{sphericalTrigonometry.nearestOn3} and 

677 L{pygeodesy.equirectangular_}, 

678 

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

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

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

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

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

684 function L{pygeodesy.compassAngle}. 

685 

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

687 see function L{pygeodesy.equirectangular_}. 

688 

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

690 

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

692 

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

694 

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

696 and L{pygeodesy.nearestOn5}. 

697 ''' 

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

699 LatLon=self.classof, **wrap_adjust_limit) 

700 

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

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

703 coordinates. 

704 

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

706 and other keyword arguments, ignored 

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

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

709 this L{Cartesian} class or specify 

710 C{B{Cartesian}=None}. 

711 

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

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

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

715 

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

717 ''' 

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

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

720 

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

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

723 

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

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

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

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

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

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

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

731 

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

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

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

735 

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

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

738 ''' 

739 B = self.others(otherB=otherB) 

740 C = self.others(otherC=otherC) 

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

742 

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

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

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

746 

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

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

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

750 intersection} of three corresponding circles. 

751 

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

753 as B{C{radius}}). 

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

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

756 B{C{radius}}). 

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

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

759 B{C{radius}}). 

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

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

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

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

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

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

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

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

768 

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

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

771 the corresponding trilaterated points C{minPoint} and 

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

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

774 

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

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

777 

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

779 respectively largest I{radial} overlap found. 

780 

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

782 nearest respectively farthest intersection margin. 

783 

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

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

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

787 largest B{C{distance#}}. 

788 

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

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

791 no intersection or all (near-)concentric for 

792 C{B{area}=False}. 

793 

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

795 

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

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

798 or B{C{radius}}. 

799 ''' 

800 return _trilaterate5(self, distance1, 

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

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

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

804 

805 

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

807 

808 

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

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

811 (with the pointsjoined by great circle arcs). 

812 

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

814 or L{BooleanGH}). 

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

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

817 or C{None}. 

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

819 (C{bool}). 

820 

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

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

823 

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

825 

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

827 

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

829 

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

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

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

833 

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

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

836 L{ellipsoidalKarney.areaOf}. 

837 

838 @example: 

839 

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

841 >>> areaOf(b) # 8666058750.718977 

842 

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

844 >>> areaOf(c) # 6.18e9 

845 ''' 

846 if _MODS.booleans.isBoolean(points): 

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

848 

849 _at2, _t_2 = atan2, tan_2 

850 _un, _w180 = unrollPI, wrap180 

851 

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

853 p1 = p2 = Ps[0] 

854 a1, b1 = p1.philam 

855 ta1, z1 = _t_2(a1), None 

856 

857 A = Fsum() # mean phi 

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

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

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

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

862 D = Fsum() 

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

864 a2, b2 = p2.philam 

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

866 A += a2 

867 ta2 = _t_2(a2) 

868 tdb = _t_2(db, points=i) 

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

870 _1_0 + ta1 * ta2) 

871 ta1, b1 = ta2, b2 

872 

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

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

875 if z1 is not None: 

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

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

878 p1, z1 = p2, z2 

879 

880 R = abs(R * _2_0) 

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

882 R = abs(R - PI2) 

883 if radius: 

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

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

886 return float(R) 

887 

888 

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

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

891 

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

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

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

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

896 

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

898 ''' 

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

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

901 ca *= sr 

902 

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

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

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

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

907 

908 

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

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

911 # and similar logic in .ellipsoidalBaseDI._intersect3 

912 a1, b1 = s.philam 

913 

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

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

916 else: # must be a point 

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

918 hs.append(end.height) 

919 a2, b2 = end.philam 

920 if wrap: 

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

922 

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

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

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

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

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

928 -(b1 + b2) * _0_5) 

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

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

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

932 cb12 * cb21 + sb12 * sb21, 

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

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

935 

936 

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

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

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

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

941 

942 

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

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

945 '''Compute the intersections of a circle and a line given as a point and 

946 bearing or as two points. 

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 on the line (L{LatLon}). 

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

953 an other 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.