1

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

3

4 u'''(INTERNAL) Private spherical base classes C{CartesianSphericalBase} and

5 C{LatLonSphericalBase} for L{sphericalNvector} and L{sphericalTrigonometry}.

6

7 A pure Python implementation of geodetic (lat-/longitude) functions,

8 transcoded in part from JavaScript originals by I{(C) Chris Veness 2011-2016}

9 and published under the same MIT Licence**, see

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

11'''

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

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

14

15 from pygeodesy.basics import _copysign, isbool, isinstanceof, map1

16 from pygeodesy.cartesianBase import CartesianBase, Bearing2Tuple

17 from pygeodesy.constants import EPS, EPS0, PI, PI2, PI_2, R_M, \

18 _0_0, _0_5, _1_0, _180_0, _360_0, \

19 _over, isnear0, isnon0

20 from pygeodesy.datums import Datums, _earth_ellipsoid, _spherical_datum

21 from pygeodesy.errors import IntersectionError, _ValueError, \

22 _xattr, _xError

23 from pygeodesy.fmath import favg, fdot, hypot, sqrt_a

24 from pygeodesy.interns import _COMMA_, _concentric_, _datum_, _distant_, \

25 _exceed_PI_radians_, _name_, _near_, \

26 _radius_, _too_

27 from pygeodesy.latlonBase import LatLonBase, _trilaterate5 # PYCHOK passed

28 from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS

29# from pygeodesy.namedTuples import Bearing2Tuple # from .cartesianBase

30 from pygeodesy.nvectorBase import NvectorBase, Fmt, _xattrs

31 from pygeodesy.props import deprecated_method, property_doc_, property_RO, \

32 _update_all

33# from pygeodesy.streprs import Fmt, _xattrs # from .nvectorBase

34 from pygeodesy.units import Bearing, Bearing_, _isRadius, Radians_, Radius, \

35 Radius_, Scalar_, _100km

36 from pygeodesy.utily import acos1, asin1, atan2b, atan2d, degrees90, \

37 degrees180, sincos2, sincos2d, _unrollon, \

38 tanPI_2_2, wrapPI

39

40 from math import cos, fabs, log, sin, sqrt

41

42__all__ = _ALL_LAZY.sphericalBase

43__version__ = '24.05.31'

44

45

46 class CartesianSphericalBase(CartesianBase):

47 '''(INTERNAL) Base class for spherical C{Cartesian}s.

48 '''

49 _datum = Datums.Sphere # L{Datum}

50

51 def intersections2(self, rad1, other, rad2, radius=R_M):

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

53 by a center point and a radius.

54

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

56 see B{C{radius}}).

57 @arg other: Center of the other circle (C{Cartesian}).

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

59 see B{C{radius}}).

60 @kwarg radius: Mean earth radius (C{meter} or C{None} if both

61 B{C{rad1}} and B{C{rad2}} are given in C{radians}).

62

63 @return: 2-Tuple of the intersection points, each C{Cartesian}.

64 For abutting circles, the intersection points are the

65 same C{Cartesian} instance, aka the I{radical center}.

66

67 @raise IntersectionError: Concentric, antipodal, invalid or

68 non-intersecting circles.

69

70 @raise TypeError: If B{C{other}} is not C{Cartesian}.

71

72 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}} or B{C{radius}}.

73

74 @see: U{Calculating intersection of two Circles

75 <https://GIS.StackExchange.com/questions/48937/

76 calculating-intersection-of-two-circles>} and method

77 or function C{trilaterate3d2}.

78 '''

79 x1, x2 = self, self.others(other)

80 r1, r2, x = _rads3(rad1, rad2, radius)

81 if x:

82 x1, x2 = x2, x1

83 try:

84 n, q = x1.cross(x2), x1.dot(x2)

85 n2, q1 = n.length2, (_1_0 - q**2)

86 if n2 < EPS or isnear0(q1):

87 raise ValueError(_near_(_concentric_))

88 c1, c2 = cos(r1), cos(r2)

89 x0 = x1.times((c1 - q * c2) / q1).plus(

90 x2.times((c2 - q * c1) / q1))

91 n1 = _1_0 - x0.length2

92 if n1 < EPS:

93 raise ValueError(_too_(_distant_))

94 except ValueError as x:

95 raise IntersectionError(center=self, rad1=rad1,

96 other=other, rad2=rad2, cause=x)

97 n = n.times(sqrt(n1 / n2))

98 if n.length > EPS:

99 x1 = x0.plus(n)

100 x2 = x0.minus(n)

101 else: # abutting circles

102 x1 = x2 = x0

103

104 return (_xattrs(x1, self, _datum_, _name_),

105 _xattrs(x2, self, _datum_, _name_))

106

107 @property_RO

108 def sphericalCartesian(self):

109 '''Get this C{Cartesian}'s spherical class.

110 '''

111 return type(self)

112

113

114 class LatLonSphericalBase(LatLonBase):

115 '''(INTERNAL) Base class for spherical C{LatLon}s.

116 '''

117 _datum = Datums.Sphere # spherical L{Datum}

118 _napieradius = _100km

119

120 def __init__(self, latlonh, lon=None, height=0, datum=None, wrap=False, **name):

121 '''Create a spherical C{LatLon} point frome the given lat-, longitude and

122 height on the given datum.

123

124 @arg latlonh: Latitude (C{degrees} or DMS C{str} with N or S suffix) or

125 a previous C{LatLon} instance provided C{B{lon}=None}.

126 @kwarg lon: Longitude (C{degrees} or DMS C{str} with E or W suffix) or

127 C(None), indicating B{C{latlonh}} is a C{LatLon}.

128 @kwarg height: Optional height above (or below) the earth surface (C{meter},

129 same units as the datum's ellipsoid axes or radius).

130 @kwarg datum: Optional, spherical datum to use (L{Datum}, L{Ellipsoid},

131 L{Ellipsoid2}, L{a_f2Tuple}) or earth radius in C{meter},

132 conventionally).

133 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{lat}} and B{C{lon}}

134 (C{bool}).

135 @kwarg name: Optional C{B{name}=NN} (C{str}).

136

137 @raise TypeError: If B{C{latlonh}} is not a C{LatLon} or B{C{datum}} not

138 spherical.

139 '''

140 LatLonBase.__init__(self, latlonh, lon=lon, height=height, wrap=wrap, **name)

141 if datum not in (None, self.datum):

142 self.datum = datum

143

144 def bearingTo2(self, other, wrap=False, raiser=False):

145 '''Return the initial and final bearing (forward and reverse

146 azimuth) from this to an other point.

147

148 @arg other: The other point (C{LatLon}).

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

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

151

152 @return: A L{Bearing2Tuple}C{(initial, final)}.

153

154 @raise TypeError: The B{C{other}} point is not spherical.

155

156 @see: Methods C{initialBearingTo} and C{finalBearingTo}.

157 '''

158 # .initialBearingTo is inside .-Nvector and .-Trigonometry

159 i = self.initialBearingTo(other, wrap=wrap, raiser=raiser) # PYCHOK .initialBearingTo

160 f = self.finalBearingTo( other, wrap=wrap, raiser=raiser)

161 return Bearing2Tuple(i, f, name=self.name)

162

163 @property_doc_(''' this point's datum (L{Datum}).''')

164 def datum(self):

165 '''Get this point's datum (L{Datum}).

166 '''

167 return self._datum

168

169 @datum.setter # PYCHOK setter!

170 def datum(self, datum):

171 '''Set this point's datum I{without conversion} (L{Datum}, L{Ellipsoid},

172 L{Ellipsoid2}, L{a_f2Tuple}) or C{scalar} spherical earth radius).

173

174 @raise TypeError: If B{C{datum}} invalid or not not spherical.

175 '''

176 d = _spherical_datum(datum, name=self.name, raiser=_datum_)

177 if self._datum != d:

178 _update_all(self)

179 self._datum = d

180

181 def finalBearingTo(self, other, wrap=False, raiser=False):

182 '''Return the final bearing (reverse azimuth) from this to

183 an other point.

184

185 @arg other: The other point (spherical C{LatLon}).

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

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

188

189 @return: Final bearing (compass C{degrees360}).

190

191 @raise TypeError: The B{C{other}} point is not spherical.

192 '''

193 p = self.others(other)

194 if wrap:

195 p = _unrollon(self, p, wrap=wrap)

196 # final bearing is the reverse of the other, initial one

197 b = p.initialBearingTo(self, wrap=False, raiser=raiser) + _180_0

198 return b if b < 360 else (b - _360_0)

199

200 def intersecant2(self, circle, point, other, radius=R_M, exact=False, # PYCHOK signature

201 height=None, wrap=False):

202 '''Compute the intersections of a circle and a (great circle) line

203 given as two points or as a point and bearing.

204

205 @arg circle: Radius of the circle centered at this location (C{meter},

206 same units as B{C{radius}}) or a point on the circle

207 (this C{LatLon}).

208 @arg point: A point on the (great circle) line (this C{LatLon}).

209 @arg other: An other point I{on} (this {LatLon}) or the bearing at

210 B{C{point}} I{of} the (great circle) line (compass

211 C{degrees}).

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

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

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

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

216 circle} methods.

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

218 conventionally) or C{None} for interpolated heights.

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

220 B{C{circle}}, B{C{point}} and/or B{C{other}} (C{bool}).

221

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

223 an instance of the B{C{point}} class. Both points are the same

224 instance if the (great circle) line is tangent to the circle.

225

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

227

228 @raise TypeError: If B{C{point}} is not this C{LatLon} or B{C{circle}}

229 or B{C{other}} invalid.

230

231 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{radius}},

232 B{C{exact}}, B{C{height}} or B{C{napieradius}}.

233 '''

234 p = self.others(point=point)

235 try:

236 return _intersecant2(self, circle, p, other, radius=radius, exact=exact,

237 height=height, wrap=wrap)

238 except (TypeError, ValueError) as x:

239 raise _xError(x, center=self, circle=circle, point=point, other=other,

240 radius=radius, exact=exact, height=height, wrap=wrap)

241

242 def maxLat(self, bearing):

243 '''Return the maximum latitude reached when travelling on a great circle

244 on given bearing from this point based on Clairaut's formula.

245

246 The maximum latitude is independent of longitude and the same for all

247 points on a given latitude.

248

249 Negate the result for the minimum latitude (on the Southern hemisphere).

250

251 @arg bearing: Initial bearing (compass C{degrees360}).

252

253 @return: Maximum latitude (C{degrees90}).

254

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

256 '''

257 r = acos1(fabs(sin(Bearing_(bearing)) * cos(self.phi)))

258 return degrees90(r)

259

260 def minLat(self, bearing):

261 '''Return the minimum latitude reached when travelling on a great circle

262 on given bearing from this point.

263

264 @arg bearing: Initial bearing (compass C{degrees360}).

265

266 @return: Minimum latitude (C{degrees90}).

267

268 @see: Method L{maxLat} for more details.

269

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

271 '''

272 return -self.maxLat(bearing)

273

274 def _mpr(self, radius=R_M, exact=None): # meter per radian

275 if exact and not _isRadius(radius): # see .rhumb.ekx.Rhumb._mpr

276 radius = _earth_ellipsoid(radius)._Lpr

277 return radius

278

279 @property_doc_(''' the I{Napier} radius to apply spherical trigonometry.''')

280 def napieradius(self):

281 '''Get the I{Napier} radius (C{meter}, conventionally).

282 '''

283 return self._napieradius

284

285 @napieradius.setter # PYCHOK setter!

286 def napieradius(self, radius):

287 '''Set this I{Napier} radius (C{meter}, conventionally) or C{0}.

288

289 In methods L{intersecant2} and L{rhumbIntersecant2}, I{Napier}'s

290 spherical trigonometry is applied if the circle radius exceeds

291 the I{Napier} radius, otherwise planar trigonometry is used.

292

293 @raise UnitError: Invalid B{C{radius}}.

294 '''

295 self._napieradius = Radius(napieradius=radius or 0)

296

297# def nearestTo(self, point, other, **radius_exact_height_wrap): # PYCHOK signature

298# p = self.others(point=point)

299# try:

300# p, q = _intersecant2(self, p, p, other, **radius_exact_height_wrap)

301# except (TypeError, ValueError) as x:

302# raise _xError(x, this=self, point=point, other=other, **radius_exact_height_wrap)

303# return p.midpointTo(q)

304

305 def parse(self, strllh, height=0, sep=_COMMA_, **name):

306 '''Parse a string representing a similar, spherical C{LatLon}

307 point, consisting of C{"lat, lon[, height]"}.

308

309 @arg strllh: Lat, lon and optional height (C{str}), see function

310 L{pygeodesy.parse3llh}.

311 @kwarg height: Optional, default height (C{meter}).

312 @kwarg sep: Optional separator (C{str}).

313 @kwarg name: Optional C{B{name}=NN} (C{str}), overriding this name.

314

315 @return: The similar point (spherical C{LatLon}).

316

317 @raise ParseError: Invalid B{C{strllh}}.

318 '''

319 llh = _MODS.dms.parse3llh(strllh, height=height, sep=sep)

320 return self.classof(*llh, **name)

321

322 @property_RO

323 def _radius(self):

324 '''(INTERNAL) Get this sphere's radius.

325 '''

326 return self.datum.ellipsoid.equatoradius

327

328 def _rhumbs3(self, other, wrap, r=False): # != .latlonBase._rhumbx3

329 '''(INTERNAL) Rhumb_ helper function.

330

331 @arg other: The other point (spherical C{LatLon}).

332 '''

333 p = self.others(other, up=2)

334 if wrap:

335 p = _unrollon(self, p, wrap=wrap)

336 a2, b2 = p.philam

337 a1, b1 = self.philam

338 # if |db| > 180 take shorter rhumb

339 # line across the anti-meridian

340 db = wrapPI(b2 - b1)

341 dp = _logPI_2_2(a2, a1)

342 da = a2 - a1

343 if r:

344 # on Mercator projection, longitude distances shrink

345 # by latitude; the 'stretch factor' q becomes ill-

346 # conditioned along E-W line (0/0); use an empirical

347 # tolerance to avoid it

348 q = (da / dp) if fabs(dp) > EPS else cos(a1)

349 da = hypot(da, q * db) # angular distance radians

350 return da, db, dp

351

352 def rhumbAzimuthTo(self, other, radius=R_M, exact=False, wrap=False, b360=False):

353 '''Return the azimuth (bearing) of a rhumb line (loxodrome) between

354 this and an other (spherical) point.

355

356 @arg other: The other point (spherical C{LatLon}).

357 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},

358 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).

359 @kwarg exact: If C{True}, use I{Elliptic, Krüger} L{Rhumb} (C{bool}),

360 default C{False} for backward compatibility.

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

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

363 @kwarg b360: If C{True}, return the azimuth in the bearing range.

364

365 @return: Rhumb azimuth (compass C{degrees180} or C{degrees360}).

366

367 @raise TypeError: The B{C{other}} point is incompatible or

368 B{C{radius}} is invalid.

369 '''

370 if exact: # use series, always

371 z = LatLonBase.rhumbAzimuthTo(self, other, exact=False, # Krüger

372 radius=radius, wrap=wrap, b360=b360)

373 else:

374 _, db, dp = self._rhumbs3(other, wrap)

375 z = (atan2b if b360 else atan2d)(db, dp) # see .rhumbBase.RhumbBase.Inverse

376 return z

377

378 @deprecated_method

379 def rhumbBearingTo(self, other): # unwrapped

380 '''DEPRECATED, use method C{.rhumbAzimuthTo}.'''

381 return self.rhumbAzimuthTo(other, b360=True) # [0..360)

382

383 def rhumbDestination(self, distance, azimuth, radius=R_M, height=None, exact=False):

384 '''Return the destination point having travelled the given distance from

385 this point along a rhumb line (loxodrome) of the given azimuth.

386

387 @arg distance: Distance travelled (C{meter}, same units as B{C{radius}}),

388 may be negative if C{B{exact}=True}.

389 @arg azimuth: Azimuth (bearing) of the rhumb line (compass C{degrees}).

390 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},

391 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}) if

392 C{B{exact}=True}.

393 @kwarg height: Optional height, overriding the default height (C{meter}.

394 @kwarg exact: If C{True}, use I{Elliptic, Krüger} L{Rhumb} (C{bool}),

395 default C{False} for backward compatibility.

396

397 @return: The destination point (spherical C{LatLon}).

398

399 @raise ValueError: Invalid B{C{distance}}, B{C{azimuth}}, B{C{radius}}

400 or B{C{height}}.

401 '''

402 if exact: # use series, always

403 r = LatLonBase.rhumbDestination(self, distance, azimuth, exact=False, # Krüger

404 radius=radius, height=height)

405 else: # radius=None from .rhumbMidpointTo

406 if radius in (None, self._radius):

407 d, r = self.datum, radius

408 else:

409 d = _spherical_datum(radius, raiser=_radius_) # spherical only

410 r = d.ellipsoid.equatoradius

411 r = _m2radians(distance, r, low=-EPS) # distance=0 from .rhumbMidpointTo

412

413 a1, b1 = self.philam

414 sb, cb = sincos2(Bearing_(azimuth)) # radians

415

416 da = r * cb

417 a2 = a1 + da

418 # normalize latitude if past pole

419 if fabs(a2) > PI_2:

420 a2 = _copysign(PI, a2) - a2

421

422 dp = _logPI_2_2(a2, a1)

423 # q becomes ill-conditioned on E-W course 0/0

424 q = cos(a1) if isnear0(dp) else (da / dp)

425 b2 = b1 if isnear0(q) else (b1 + r * sb / q)

426

427 h = self._heigHt(height)

428 r = self.classof(degrees90(a2), degrees180(b2), datum=d, height=h)

429 return r

430

431 def rhumbDistanceTo(self, other, radius=R_M, exact=False, wrap=False):

432 '''Return the distance from this to an other point along

433 a rhumb line (loxodrome).

434

435 @arg other: The other point (spherical C{LatLon}).

436 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},

437 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}) if

438 C{B{exact}=True}.

439 @kwarg exact: If C{True}, use I{Elliptic, Krüger} L{Rhumb} (C{bool}),

440 default C{False} for backward compatibility.

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

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

443

444 @return: Distance (C{meter}, the same units as B{C{radius}}

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

446

447 @raise TypeError: The B{C{other}} point is incompatible.

448

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

450 '''

451 if exact: # use series, always

452 r = LatLonBase.rhumbDistanceTo(self, other, exact=False, # Krüger

453 radius=radius, wrap=wrap)

454 if radius is None: # angular distance in radians

455 r = r / self._radius # /= chokes PyChecker

456 else:

457 # see <https://www.EdWilliams.org/avform.htm#Rhumb>

458 r, _, _ = self._rhumbs3(other, wrap, r=True)

459 if radius is not None:

460 r *= Radius(radius)

461 return r

462

463 def rhumbIntersecant2(self, circle, point, other, radius=R_M, exact=True, # PYCHOK signature

464 height=None, wrap=False):

465 '''Compute the intersections of a circle and a rhumb line given as two

466 points and as a point and azimuth.

467

468 @arg circle: Radius of the circle centered at this location (C{meter},

469 same units as B{C{radius}}) or a point on the circle

470 (this C{LatLon}).

471 @arg point: The rhumb line's start point (this C{LatLon}).

472 @arg other: An other point (this I{on} C{LatLon}) or the azimuth I{of}

473 (compass C{degrees}) the rhumb line.

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

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

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

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

478 circle} methods.

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

480 conventionally) or C{None}.

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

482 B{C{circle}}, B{C{point}} and/or B{C{other}} (C{bool}).

483

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

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

486 points are the same instance, wrapped or I{normalized}.

487

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

489

490 @raise TypeError: If B{C{point}} is not this C{LatLon} or B{C{circle}}

491 or B{C{other}} invalid.

492

493 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{radius}},

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

495 '''

496 m = LatLonBase.rhumbIntersecant2 if exact else \

497 LatLonSphericalBase.intersecant2

498 return m(self, circle, point, other, radius=radius, exact=exact,

499 height=height, wrap=wrap)

500

501 def rhumbMidpointTo(self, other, height=None, radius=R_M, exact=False,

502 fraction=_0_5, wrap=False):

503 '''Return the (loxodromic) midpoint on the rhumb line between

504 this and an other point.

505

506 @arg other: The other point (spherical LatLon).

507 @kwarg height: Optional height, overriding the mean height (C{meter}).

508 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},

509 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).

510 @kwarg exact: If C{True}, use I{Elliptic, Krüger} L{Rhumb} (C{bool}),

511 default C{False} for backward compatibility.

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

513 be negative if C{B{exact}=True}.

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

515 point (C{bool}).

516

517 @return: The (mid)point at the given B{C{fraction}} along the rhumb

518 line (spherical C{LatLon}).

519

520 @raise TypeError: The B{C{other}} point is incompatible.

521

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

523 '''

524 if exact: # use series, always

525 r = LatLonBase.rhumbMidpointTo(self, other, exact=False, # Krüger

526 radius=radius, height=height,

527 fraction=fraction, wrap=wrap)

528 elif fraction is not _0_5:

529 f = Scalar_(fraction=fraction) # low=_0_0

530 r, db, dp = self._rhumbs3(other, wrap, r=True) # radians

531 z = atan2b(db, dp)

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

533 r = self.rhumbDestination(r * f, z, radius=None, height=h)

534

535 else: # for backward compatibility, unwrapped

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

537 a1, b1 = self.philam

538 a2, b2 = self.others(other).philam

539

540 if fabs(b2 - b1) > PI:

541 b1 += PI2 # crossing anti-meridian

542

543 a3 = favg(a1, a2)

544 b3 = favg(b1, b2)

545

546 f1 = tanPI_2_2(a1)

547 if isnon0(f1):

548 f2 = tanPI_2_2(a2)

549 f = f2 / f1

550 if isnon0(f):

551 f = log(f)

552 if isnon0(f):

553 f3 = tanPI_2_2(a3)

554 b3 = fdot(map1(log, f1, f2, f3),

555 -b2, b1, b2 - b1) / f

556

557 d = self.datum if radius in (None, self._radius) else \

558 _spherical_datum(radius, name=self.name, raiser=_radius_)

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

560 r = self.classof(degrees90(a3), degrees180(b3), datum=d, height=h)

561 return r

562

563 @property_RO

564 def sphericalLatLon(self):

565 '''Get this C{LatLon}'s spherical class.

566 '''

567 return type(self)

568

569 def toNvector(self, Nvector=NvectorBase, **Nvector_kwds): # PYCHOK signature

570 '''Convert this point to C{Nvector} components, I{including

571 height}.

572

573 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}}

574 keyword arguments, ignored if

575 C{B{Nvector} is None}.

576

577 @return: An B{C{Nvector}} or a L{Vector4Tuple}C{(x, y, z, h)}

578 if B{C{Nvector}} is C{None}.

579

580 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}.

581 '''

582 return LatLonBase.toNvector(self, Nvector=Nvector, **Nvector_kwds)

583

584

585 def _intersecant2(c, r, p, b, radius=R_M, exact=False, height=None, wrap=False):

586 # (INTERNAL) Intersect a circle and line, see L{intersecant2}

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

588

589 if wrap:

590 p = _unrollon(c, p, wrap=wrap)

591 nonexact = exact is None

592

593 if not isinstanceof(r, c.__class__, p.__class__):

594 r = Radius_(circle=r)

595 elif nonexact:

596 r = c.distanceTo(r, radius=radius, wrap=wrap)

597 elif isbool(exact):

598 r = c.rhumbDistanceTo(r, radius=radius, exact=exact, wrap=wrap)

599 else:

600 raise _ValueError(exact=exact)

601

602 if not isinstanceof(b, c.__class__, p.__class__):

603 b = Bearing(b)

604 elif nonexact:

605 b = p.initialBearingTo(b, wrap=wrap)

606 else:

607 b = p.rhumbAzimuthTo(b, radius=radius, exact=exact, wrap=wrap,

608 b360=True)

609

610 d = p.distanceTo(c, radius=radius) if nonexact else \

611 p.rhumbDistanceTo(c, radius=radius, exact=exact)

612 if d > EPS0:

613 n = _xattr(c, napieradius=0)

614 a = p.initialBearingTo(c) if nonexact else \

615 p.rhumbAzimuthTo(c, radius=radius, exact=exact, b360=True)

616 s, c = sincos2d(b - a) # Napier's sin(A), cos(A)

617 if r > n:

618 # Napier's right spherical triangle rules (R2) and (R1)

619 # <https://WikiPedia.org/wiki/Spherical_trigonometry>

620 m = p._mpr(radius=radius, exact=exact) # meter per radian

621 if fabs(c) > EPS0:

622 d = d / m # /= chokes PyChecker

623 a = asin1(sin(d) * fabs(s)) # Napier's a

624 c = _copysign(cos(a), c)

625 d = acos1(cos(d) / c) * m

626 a *= m # meter

627 else: # point and chord center coincident

628 a, d = d, 0

629 c = cos(a / m)

630 h = (acos1(cos(r / m) / c) * m) if a < r else 0

631 else: # distance from the chord center to ...

632 a = fabs(d * s) # ... the cicle center ...

633 d *= c # ... and to the point

634 h = sqrt_a(r, a) if a < r else 0 # half chord length

635 if a > r:

636 raise IntersectionError(_too_(Fmt.distant(a)))

637 else:

638 d, h = 0, r # point and circle center coincident

639

640 _intersecant1, kwds = (p.destination, {}) if nonexact else \

641 (p.rhumbDestination, dict(exact=exact))

642 kwds.update(radius=radius, height=height)

643 t = (_intersecant1(d + h, b, **kwds),)

644 if h:

645 t += (_intersecant1(d - h, b, **kwds),)

646 else: # same instance twice

647 t *= 2

648 return t

649

650

651 def _logPI_2_2(a2, a1):

652 '''(INTERNAL) C{log} of C{tanPI_2_2}'s quotient.

653 '''

654 return log(_over(tanPI_2_2(a2), tanPI_2_2(a1)))

655

656

657 def _m2radians(distance, radius, low=EPS): # PYCHOK in .spherical*

658 '''(INTERNAL) Distance in C{meter} to angular distance in C{radians}.

659

660 @raise UnitError: Invalid B{C{distance}} or B{C{radius}}.

661 '''

662 r = float(distance)

663 if radius:

664 r = r / Radius_(radius=radius) # /= chokes PyChecker

665 if low is not None:

666 # small near0 values from .rhumbDestination not exact OK

667 r = _0_0 if low < 0 and r < 0 else Radians_(r, low=low)

668 # _0_0 if low < 0 and low < r < 0 else Radians_(r, low=low)

669 return r

670

671

672 def _radians2m(rad, radius):

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

674 '''

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

676 rad *= R_M if radius is R_M else Radius(radius)

677 return rad

678

679

680 def _rads3(rad1, rad2, radius): # in .sphericalTrigonometry

681 '''(INTERNAL) Convert radii to radians.

682 '''

683 r1 = Radius_(rad1=rad1)

684 r2 = Radius_(rad2=rad2)

685 if radius is not None: # convert radii to radians

686 r1 = _m2radians(r1, radius)

687 r2 = _m2radians(r2, radius)

688

689 x = r1 < r2

690 if x:

691 r1, r2 = r2, r1

692 if r1 > PI:

693 raise IntersectionError(rad1=rad1, rad2=rad2,

694 txt=_exceed_PI_radians_)

695 return r1, r2, x

696

697

698__all__ += _ALL_DOCS(CartesianSphericalBase, LatLonSphericalBase)

699

700# **) MIT License

701#

702# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

703#

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

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

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

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

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

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

710#

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

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

713#

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

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

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

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

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

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

720# OTHER DEALINGS IN THE SOFTWARE.