1

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

3

4 u'''Various utility functions.

5

6 After I{(C) Chris Veness 2011-2015} published under the same MIT Licence**, see

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

8 U{Vector-based geodesy<https://www.Movable-Type.co.UK/scripts/latlong-vectors.html>}.

9'''

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

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

12

13 from pygeodesy.basics import copysign0, isint, isstr

14 from pygeodesy.constants import EPS, EPS0, INF, NAN, NEG0, NINF, PI, PI2, PI_2, R_M, \

15 _float as _F, _isfinite, isnan, isnear0, isneg0, _M_KM, \

16 _M_NM, _M_SM, _0_0, _1__90, _0_5, _1_0, _N_1_0, _2__PI, \

17 _10_0, _90_0, _N_90_0, _180_0, _N_180_0, _360_0, _400_0

18 from pygeodesy.errors import _ValueError, _xkwds, _xkwds_get

19 from pygeodesy.interns import _edge_, _radians_, _semi_circular_, _SPACE_

20 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

21 from pygeodesy.units import Degrees, Feet, Float, Lam, Lam_, Meter, Meter2, Radians

22

23 from math import acos, asin, atan2, cos, degrees, fabs, radians, sin, tan # pow

24

25__all__ = _ALL_LAZY.utily

26__version__ = '23.09.19'

27

28# read constant name "_M_Unit" as "meter per Unit"

29_M_CHAIN = _F( 20.1168) # yard2m(1) * 22

30_M_FATHOM = _F( 1.8288) # yard2m(1) * 2 or _M_NM * 1e-3

31_M_FOOT = _F( 0.3048) # Int'l (1 / 3.2808398950131 = 10_000 / (254 * 12))

32_M_FOOT_GE = _F( 0.31608) # German Fuss (1 / 3.1637560111364)

33_M_FOOT_FR = _F( 0.3248406) # French Pied-du-Roi or pied (1 / 3.0784329298739)

34_M_FOOT_USVY = _F( 0.3048006096012192) # US Survey (1200 / 3937)

35_M_FURLONG = _F( 201.168) # 220 * yard2m(1) == 10 * m2chain(1)

36# _M_KM = _F(1000.0) # kilo meter

37# _M_NM = _F(1852.0) # nautical mile

38# _M_SM = _F(1609.344) # statute mile

39_M_TOISE = _F( 1.9490436) # French toise, 6 pieds (6 / 3.0784329298739)

40_M_YARD_UK = _F( 0.9144) # 254 * 12 * 3 / 10_000 == 3 * ft2m(1) Int'l

41

42

43 def acos1(x):

44 '''Return C{math.acos(max(-1, min(1, B{x})))}.

45 '''

46 return acos(x) if fabs(x) < _1_0 else (PI if x < 0 else _0_0)

47

48

49 def acre2ha(acres):

50 '''Convert acres to hectare.

51

52 @arg acres: Value in acres (C{scalar}).

53

54 @return: Value in C{hectare} (C{float}).

55

56 @raise ValueError: Invalid B{C{acres}}.

57 '''

58 # 0.40468564224 == acre2m2(1) / 10_000

59 return Float(ha=Float(acres) * 0.40468564224)

60

61

62 def acre2m2(acres):

63 '''Convert acres to I{square} meter.

64

65 @arg acres: Value in acres (C{scalar}).

66

67 @return: Value in C{meter^2} (C{float}).

68

69 @raise ValueError: Invalid B{C{acres}}.

70 '''

71 # 4046.8564224 == chain2m(1) * furlong2m(1)

72 return Meter2(Float(acres) * 4046.8564224)

73

74

75 def asin1(x):

76 '''Return C{math.asin(max(-1, min(1, B{x})))}.

77 '''

78 return asin(x) if fabs(x) < _1_0 else (PI_2 if x > 0 else -PI_2) # not PI3_2!

79

80

81 def atand(y_x):

82 '''Return C{atan(B{y_x})} angle in C{degrees}.

83

84 @see: Function L{pygeodesy.atan2d}.

85 '''

86 return atan2d(y_x, _1_0)

87

88

89 def atan2b(y, x):

90 '''Return C{atan2(B{y}, B{x})} in degrees M{[0..+360]}.

91

92 @see: Function L{pygeodesy.atan2d}.

93 '''

94 d = atan2d(y, x)

95 if d < 0:

96 d += _360_0

97 return d

98

99

100 def atan2d(y, x, reverse=False):

101 '''Return C{atan2(B{y}, B{x})} in degrees M{[-180..+180]},

102 optionally I{reversed} (by 180 degrees for C{azi2}).

103

104 @see: I{Karney}'s C++ function U{Math.atan2d

105 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Math.html>}.

106 '''

107 if fabs(y) > fabs(x) > 0:

108 if y < 0: # q = 3

109 d = degrees(atan2(x, -y)) - _90_0

110 else: # q = 2

111 d = _90_0 - degrees(atan2(x, y))

112 elif x < 0: # q = 1

113 d = copysign0(_180_0, y) - degrees(atan2(y, -x))

114 elif x > 0: # q = 0

115 d = degrees(atan2(y, x)) if y else _0_0

116 elif isnan(x) or isnan(y):

117 return NAN

118 else: # x == 0

119 d = _N_90_0 if y < 0 else (_90_0 if y > 0 else _0_0)

120 return _azireversed(d) if reverse else d

121

122

123 def _azireversed(azimuth): # in .rhumbBase

124 '''(INTERNAL) Return the I{reverse} B{C{azimuth}}.

125 '''

126 return azimuth + (_N_180_0 if azimuth > 0 else _180_0)

127

128

129 def chain2m(chains):

130 '''Convert I{UK} chains to meter.

131

132 @arg chains: Value in chains (C{scalar}).

133

134 @return: Value in C{meter} (C{float}).

135

136 @raise ValueError: Invalid B{C{chains}}.

137 '''

138 return Meter(Float(chains=chains) * _M_CHAIN)

139

140

141 def circle4(earth, lat):

142 '''Get the equatorial or a parallel I{circle of latitude}.

143

144 @arg earth: The earth radius, ellipsoid or datum

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

146 L{Datum} or L{a_f2Tuple}).

147 @arg lat: Geodetic latitude (C{degrees90}, C{str}).

148

149 @return: A L{Circle4Tuple}C{(radius, height, lat, beta)}

150 instance.

151

152 @raise RangeError: Latitude B{C{lat}} outside valid range and

153 L{pygeodesy.rangerrors} set to C{True}.

154

155 @raise TypeError: Invalid B{C{earth}}.

156

157 @raise ValueError: B{C{earth}} or B{C{lat}}.

158 '''

159 E = _MODS.datums._spherical_datum(earth).ellipsoid

160 return E.circle4(lat)

161

162

163 def cot(rad, **error_kwds):

164 '''Return the C{cotangent} of an angle in C{radians}.

165

166 @arg rad: Angle (C{radians}).

167 @kwarg error_kwds: Error to raise (C{ValueError}).

168

169 @return: C{cot(B{rad})}.

170

171 @raise ValueError: L{pygeodesy.isnear0}C{(sin(B{rad})}.

172 '''

173 s, c = sincos2(rad)

174 if isnear0(s):

175 raise _valueError(cot, rad, **error_kwds)

176 return c / s

177

178

179 def cot_(*rads, **error_kwds):

180 '''Return the C{cotangent} of angle(s) in C{radiansresection}.

181

182 @arg rads: One or more angles (C{radians}).

183 @kwarg error_kwds: Error to raise (C{ValueError}).

184

185 @return: Yield the C{cot(B{rad})} for each angle.

186

187 @raise ValueError: See L{pygeodesy.cot}.

188 '''

189 try:

190 for r in rads:

191 yield cot(r)

192 except ValueError:

193 raise _valueError(cot_, r, **error_kwds)

194

195

196 def cotd(deg, **error_kwds):

197 '''Return the C{cotangent} of an angle in C{degrees}.

198

199 @arg deg: Angle (C{degrees}).

200 @kwarg error_kwds: Error to raise (C{ValueError}).

201

202 @return: C{cot(B{deg})}.

203

204 @raise ValueError: L{pygeodesy.isnear0}C{(sin(B{deg})}.

205 '''

206 s, c = sincos2d(deg)

207 if isnear0(s):

208 raise _valueError(cotd, deg, **error_kwds)

209 return c / s

210

211

212 def cotd_(*degs, **error_kwds):

213 '''Return the C{cotangent} of angle(s) in C{degrees}.

214

215 @arg degs: One or more angles (C{degrees}).

216 @kwarg error_kwds: Error to raise (C{ValueError}).

217

218 @return: Yield the C{cot(B{deg})} for each angle.

219

220 @raise ValueError: See L{pygeodesy.cotd}.

221 '''

222 try:

223 for d in degs:

224 yield cotd(d)

225 except ValueError:

226 raise _valueError(cotd_, d, **error_kwds)

227

228

229 def degrees90(rad):

230 '''Convert radians to degrees and wrap M{[-270..+90]}.

231

232 @arg rad: Angle (C{radians}).

233

234 @return: Angle, wrapped (C{degrees90}).

235 '''

236 return _wrap(degrees(rad), _90_0, _360_0)

237

238

239 def degrees180(rad):

240 '''Convert radians to degrees and wrap M{[-180..+180]}.

241

242 @arg rad: Angle (C{radians}).

243

244 @return: Angle, wrapped (C{degrees180}).

245 '''

246 return _wrap(degrees(rad), _180_0, _360_0)

247

248

249 def degrees360(rad):

250 '''Convert radians to degrees and wrap M{[0..+360)}.

251

252 @arg rad: Angle (C{radians}).

253

254 @return: Angle, wrapped (C{degrees360}).

255 '''

256 return _wrap(degrees(rad), _360_0, _360_0)

257

258

259 def degrees2grades(deg):

260 '''Convert degrees to I{grades} (aka I{gons} or I{gradians}).

261

262 @arg deg: Angle (C{degrees}).

263

264 @return: Angle (C{grades}).

265 '''

266 return Degrees(deg) * _400_0 / _360_0

267

268

269 def degrees2m(deg, radius=R_M, lat=0):

270 '''Convert an angle to a distance along the equator or

271 along the parallel at an other (geodetic) latitude.

272

273 @arg deg: The angle (C{degrees}).

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

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

276 L{Datum} or L{a_f2Tuple}).

277 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

278

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

280 or equatorial and polar radii) or C{0.0} for

281 near-polar B{C{lat}}.

282

283 @raise RangeError: Latitude B{C{lat}} outside valid range and

284 L{pygeodesy.rangerrors} set to C{True}.

285

286 @raise TypeError: Invalid B{C{radius}}.

287

288 @raise ValueError: Invalid B{C{deg}}, B{C{radius}} or

289 B{C{lat}}.

290

291 @see: Function L{radians2m} and L{m2degrees}.

292 '''

293 return _radians2m(Lam_(deg=deg, clip=0), radius, lat)

294

295

296 def fathom2m(fathoms):

297 '''Convert I{Imperial} fathom to meter.

298

299 @arg fathoms: Value in fathoms (C{scalar}).

300

301 @return: Value in C{meter} (C{float}).

302

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

304

305 @see: Function L{toise2m}, U{Fathom<https://WikiPedia.org/wiki/Fathom>}

306 and U{Klafter<https://WikiPedia.org/wiki/Klafter>}.

307 '''

308 return Meter(Float(fathoms=fathoms) * _M_FATHOM)

309

310

311 def ft2m(feet, usurvey=False, pied=False, fuss=False):

312 '''Convert I{International}, I{US Survey}, I{French} or I{German}

313 B{C{feet}} to C{meter}.

314

315 @arg feet: Value in feet (C{scalar}).

316 @kwarg usurvey: If C{True}, convert I{US Survey} foot else ...

317 @kwarg pied: If C{True}, convert French I{pied-du-Roi} else ...

318 @kwarg fuss: If C{True}, convert German I{Fuss}, otherwise

319 I{International} foot to C{meter}.

320

321 @return: Value in C{meter} (C{float}).

322

323 @raise ValueError: Invalid B{C{feet}}.

324 '''

325 return Meter(Feet(feet) * (_M_FOOT_USVY if usurvey else

326 (_M_FOOT_FR if pied else

327 (_M_FOOT_GE if fuss else _M_FOOT))))

328

329

330 def furlong2m(furlongs):

331 '''Convert a furlong to meter.

332

333 @arg furlongs: Value in furlongs (C{scalar}).

334

335 @return: Value in C{meter} (C{float}).

336

337 @raise ValueError: Invalid B{C{furlongs}}.

338 '''

339 return Meter(Float(furlongs=furlongs) * _M_FURLONG)

340

341

342 def grades(rad):

343 '''Convert radians to I{grades} (aka I{gons} or I{gradians}).

344

345 @arg rad: Angle (C{radians}).

346

347 @return: Angle (C{grades}).

348 '''

349 return Float(grades=Float(rad=rad) * _400_0 / PI2)

350

351

352 def grades400(rad):

353 '''Convert radians to I{grades} (aka I{gons} or I{gradians}) and wrap M{[0..+400)}.

354

355 @arg rad: Angle (C{radians}).

356

357 @return: Angle, wrapped (C{grades}).

358 '''

359 return Float(grades400=_wrap(grades(rad), _400_0, _400_0))

360

361

362 def grades2degrees(gon):

363 '''Convert I{grades} (aka I{gons} or I{gradians}) to C{degrees}.

364

365 @arg gon: Angle (C{grades}).

366

367 @return: Angle (C{degrees}).

368 '''

369 return Degrees(Float(gon=gon) * _360_0 / _400_0)

370

371

372 def grades2radians(gon):

373 '''Convert I{grades} (aka I{gons} or I{gradians}) to C{radians}.

374

375 @arg gon: Angle (C{grades}).

376

377 @return: Angle (C{radians}).

378 '''

379 return Radians(Float(gon=gon) * PI2 / _400_0)

380

381

382 def km2m(km):

383 '''Convert kilo meter to meter (m).

384

385 @arg km: Value in kilo meter (C{scalar}).

386

387 @return: Value in meter (C{float}).

388

389 @raise ValueError: Invalid B{C{km}}.

390 '''

391 return Meter(Float(km=km) * _M_KM)

392

393

394 def _loneg(lon):

395 '''(INTERNAL) "Complement" of C{lon}.

396 '''

397 return _180_0 - lon

398

399

400 def m2chain(meter):

401 '''Convert meter to I{UK} chains.

402

403 @arg meter: Value in meter (C{scalar}).

404

405 @return: Value in C{chains} (C{float}).

406

407 @raise ValueError: Invalid B{C{meter}}.

408 '''

409 return Float(chain=Meter(meter) / _M_CHAIN) # * 0.049709695378986715

410

411

412 def m2degrees(distance, radius=R_M, lat=0):

413 '''Convert a distance to an angle along the equator or

414 along the parallel at an other (geodetic) latitude.

415

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

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

418 an L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or

419 L{a_f2Tuple}).

420 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

421

422 @return: Angle (C{degrees}) or C{INF} for near-polar B{C{lat}}.

423

424 @raise RangeError: Latitude B{C{lat}} outside valid range and

425 L{pygeodesy.rangerrors} set to C{True}.

426

427 @raise TypeError: Invalid B{C{radius}}.

428

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

430 or B{C{lat}}.

431

432 @see: Function L{m2radians} and L{degrees2m}.

433 '''

434 return degrees(m2radians(distance, radius=radius, lat=lat))

435

436

437 def m2fathom(meter):

438 '''Convert meter to I{Imperial} fathoms.

439

440 @arg meter: Value in meter (C{scalar}).

441

442 @return: Value in C{fathoms} (C{float}).

443

444 @raise ValueError: Invalid B{C{meter}}.

445

446 @see: Function L{m2toise}, U{Fathom<https://WikiPedia.org/wiki/Fathom>}

447 and U{Klafter<https://WikiPedia.org/wiki/Klafter>}.

448 '''

449 return Float(fathom=Meter(meter) / _M_FATHOM) # * 0.546806649

450

451

452 def m2ft(meter, usurvey=False, pied=False, fuss=False):

453 '''Convert meter to I{International}, I{US Survey}, I{French} or

454 or I{German} feet (C{ft}).

455

456 @arg meter: Value in meter (C{scalar}).

457 @kwarg usurvey: If C{True}, convert to I{US Survey} foot else ...

458 @kwarg pied: If C{True}, convert to French I{pied-du-Roi} else ...

459 @kwarg fuss: If C{True}, convert to German I{Fuss}, otherwise to

460 I{International} foot.

461

462 @return: Value in C{feet} (C{float}).

463

464 @raise ValueError: Invalid B{C{meter}}.

465 '''

466 # * 3.2808333333333333, US Survey 3937 / 1200

467 # * 3.2808398950131235, Int'l 10_000 / (254 * 12)

468 return Float(feet=Meter(meter) / (_M_FOOT_USVY if usurvey else

469 (_M_FOOT_FR if pied else

470 (_M_FOOT_GE if fuss else _M_FOOT))))

471

472

473 def m2furlong(meter):

474 '''Convert meter to furlongs.

475

476 @arg meter: Value in meter (C{scalar}).

477

478 @return: Value in C{furlongs} (C{float}).

479

480 @raise ValueError: Invalid B{C{meter}}.

481 '''

482 return Float(furlong=Meter(meter) / _M_FURLONG) # * 0.00497096954

483

484

485 def m2km(meter):

486 '''Convert meter to kilo meter (Km).

487

488 @arg meter: Value in meter (C{scalar}).

489

490 @return: Value in Km (C{float}).

491

492 @raise ValueError: Invalid B{C{meter}}.

493 '''

494 return Float(km=Meter(meter) / _M_KM)

495

496

497 def m2NM(meter):

498 '''Convert meter to nautical miles (NM).

499

500 @arg meter: Value in meter (C{scalar}).

501

502 @return: Value in C{NM} (C{float}).

503

504 @raise ValueError: Invalid B{C{meter}}.

505 '''

506 return Float(NM=Meter(meter) / _M_NM) # * 5.39956804e-4

507

508

509 def m2radians(distance, radius=R_M, lat=0):

510 '''Convert a distance to an angle along the equator or

511 along the parallel at an other (geodetic) latitude.

512

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

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

515 an L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or

516 L{a_f2Tuple}).

517 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

518

519 @return: Angle (C{radians}) or C{INF} for near-polar B{C{lat}}.

520

521 @raise RangeError: Latitude B{C{lat}} outside valid range and

522 L{pygeodesy.rangerrors} set to C{True}.

523

524 @raise TypeError: Invalid B{C{radius}}.

525

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

527 or B{C{lat}}.

528

529 @see: Function L{m2degrees} and L{radians2m}.

530 '''

531 m = circle4(radius, lat).radius

532 return INF if m < EPS0 else Radians(Float(distance=distance) / m)

533

534

535 def m2SM(meter):

536 '''Convert meter to statute miles (SM).

537

538 @arg meter: Value in meter (C{scalar}).

539

540 @return: Value in C{SM} (C{float}).

541

542 @raise ValueError: Invalid B{C{meter}}.

543 '''

544 return Float(SM=Meter(meter) / _M_SM) # * 6.21369949e-4 == 1 / 1609.344

545

546

547 def m2toise(meter):

548 '''Convert meter to French U{toises<https://WikiPedia.org/wiki/Toise>}.

549

550 @arg meter: Value in meter (C{scalar}).

551

552 @return: Value in C{toises} (C{float}).

553

554 @raise ValueError: Invalid B{C{meter}}.

555

556 @see: Function L{m2fathom}.

557 '''

558 return Float(toise=Meter(meter) / _M_TOISE) # * 0.513083632632119

559

560

561 def m2yard(meter):

562 '''Convert meter to I{UK} yards.

563

564 @arg meter: Value in meter (C{scalar}).

565

566 @return: Value in C{yards} (C{float}).

567

568 @raise ValueError: Invalid B{C{meter}}.

569 '''

570 return Float(yard=Meter(meter) / _M_YARD_UK) # * 1.0936132983377078

571

572

573 def NM2m(nm):

574 '''Convert nautical miles to meter (m).

575

576 @arg nm: Value in nautical miles (C{scalar}).

577

578 @return: Value in meter (C{float}).

579

580 @raise ValueError: Invalid B{C{nm}}.

581 '''

582 return Meter(Float(nm=nm) * _M_NM)

583

584

585 def _passarg(arg): # in .auxilats.auxLat, .formy

586 '''(INTERNAL) Helper, no-op.

587 '''

588 return arg

589

590

591 def _passargs(*args): # in .formy

592 '''(INTERNAL) Helper, no-op.

593 '''

594 return args

595

596

597 def radians2m(rad, radius=R_M, lat=0):

598 '''Convert an angle to a distance along the equator or

599 along the parallel at an other (geodetic) latitude.

600

601 @arg rad: The angle (C{radians}).

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

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

604 L{Datum} or L{a_f2Tuple}).

605 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

606

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

608 or equatorial and polar radii) or C{0.0} for

609 near-polar B{C{lat}}.

610

611 @raise RangeError: Latitude B{C{lat}} outside valid range and

612 L{pygeodesy.rangerrors} set to C{True}.

613

614 @raise TypeError: Invalid B{C{radius}}.

615

616 @raise ValueError: Invalid B{C{rad}}, B{C{radius}} or

617 B{C{lat}}.

618

619 @see: Function L{degrees2m} and L{m2radians}.

620 '''

621 return _radians2m(Lam(rad=rad, clip=0), radius, lat)

622

623

624 def _radians2m(rad, radius, lat):

625 '''(INTERNAL) Helper for C{degrees2m} and C{radians2m}.

626 '''

627 m = circle4(radius, lat).radius

628 return _0_0 if m < EPS0 else (rad * m)

629

630

631 def radiansPI(deg):

632 '''Convert and wrap degrees to radians M{[-PI..+PI]}.

633

634 @arg deg: Angle (C{degrees}).

635

636 @return: Radians, wrapped (C{radiansPI})

637 '''

638 return _wrap(radians(deg), PI, PI2)

639

640

641 def radiansPI2(deg):

642 '''Convert and wrap degrees to radians M{[0..+2PI)}.

643

644 @arg deg: Angle (C{degrees}).

645

646 @return: Radians, wrapped (C{radiansPI2})

647 '''

648 return _wrap(radians(deg), PI2, PI2)

649

650

651 def radiansPI_2(deg):

652 '''Convert and wrap degrees to radians M{[-3PI/2..+PI/2]}.

653

654 @arg deg: Angle (C{degrees}).

655

656 @return: Radians, wrapped (C{radiansPI_2})

657 '''

658 return _wrap(radians(deg), PI_2, PI2)

659

660

661 def _sin0cos2(q, r, sign):

662 '''(INTERNAL) 2-tuple (C{sin(r), cos(r)}) in quadrant C{0 <= B{q} <= 3}

663 and C{sin} zero I{signed} with B{C{sign}}.

664 '''

665 if r < PI_2:

666 s, c = sin(r), cos(r)

667 t = s, c, -s, -c, s

668 else: # r == PI_2

669 t = _1_0, _0_0, _N_1_0, _0_0, _1_0

670# else: # r == 0, testUtility failures

671# t = _0_0, _1_0, _0_0, _N_1_0, _0_0

672# q &= 3

673 s = t[q] or (NEG0 if sign < 0 else _0_0)

674 c = t[q + 1] or _0_0

675 return s, c

676

677

678 def sincos2(rad):

679 '''Return the C{sine} and C{cosine} of an angle in C{radians}.

680

681 @arg rad: Angle (C{radians}).

682

683 @return: 2-Tuple (C{sin(B{rad})}, C{cos(B{rad})}).

684

685 @see: U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/

686 classGeographicLib_1_1Math.html#sincosd>} function U{sincosd

687 <https://SourceForge.net/p/geographiclib/code/ci/release/tree/

688 python/geographiclib/geomath.py#l155>} and C++ U{sincosd

689 <https://SourceForge.net/p/geographiclib/code/ci/release/tree/

690 include/GeographicLib/Math.hpp#l558>}.

691 '''

692 if _isfinite(rad):

693 q = int(rad * _2__PI) # int(math.floor)

694 if q < 0:

695 q -= 1

696 t = _sin0cos2(q & 3, rad - q * PI_2, rad)

697 else:

698 t = NAN, NAN

699 return t

700

701

702 def sincos2_(*rads):

703 '''Return the C{sine} and C{cosine} of angle(s) in C{radians}.

704

705 @arg rads: One or more angles (C{radians}).

706

707 @return: Yield the C{sin(B{rad})} and C{cos(B{rad})} for each angle.

708

709 @see: function L{sincos2}.

710 '''

711 for r in rads:

712 s, c = sincos2(r)

713 yield s

714 yield c

715

716

717 def sincos2d(deg, **adeg):

718 '''Return the C{sine} and C{cosine} of an angle in C{degrees}.

719

720 @arg deg: Angle (C{degrees}).

721 @kwarg adeg: Optional correction (C{degrees}).

722

723 @return: 2-Tuple (C{sin(B{deg_})}, C{cos(B{deg_})}, C{B{deg_} =

724 B{deg} + B{adeg}}).

725

726 @see: U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/

727 classGeographicLib_1_1Math.html#sincosd>} function U{sincosd

728 <https://SourceForge.net/p/geographiclib/code/ci/release/tree/

729 python/geographiclib/geomath.py#l155>} and C++ U{sincosd

730 <https://SourceForge.net/p/geographiclib/code/ci/release/tree/

731 include/GeographicLib/Math.hpp#l558>}.

732 '''

733 if _isfinite(deg):

734 q = int(deg * _1__90) # int(math.floor)

735 if q < 0:

736 q -= 1

737 d = deg - q * _90_0

738 if adeg:

739 t = _xkwds_get(adeg, adeg=_0_0)

740 d = _MODS.karney._around(d + t)

741 t = _sin0cos2(q & 3, radians(d), deg)

742 else:

743 t = NAN, NAN

744 return t

745

746

747 def SinCos2(x):

748 '''Get C{sin} and C{cos} of I{typed} angle.

749

750 @arg x: Angle (L{Degrees}, L{Radians} or C{radians}).

751

752 @return: 2-Tuple (C{sin(B{x})}, C{cos(B{x})}).

753 '''

754 return sincos2d(x) if isinstance(x, Degrees) else (

755 sincos2(x) if isinstance(x, Radians) else

756 sincos2(float(x))) # assume C{radians}

757

758

759 def sincos2d_(*degs):

760 '''Return the C{sine} and C{cosine} of angle(s) in C{degrees}.

761

762 @arg degs: One or more angles (C{degrees}).

763

764 @return: Yield the C{sin(B{deg})} and C{cos(B{deg})} for each angle.

765

766 @see: Function L{sincos2d}.

767 '''

768 for d in degs:

769 s, c = sincos2d(d)

770 yield s

771 yield c

772

773

774 def sincostan3(rad):

775 '''Return the C{sine}, C{cosine} and C{tangent} of an angle in C{radians}.

776

777 @arg rad: Angle (C{radians}).

778

779 @return: 3-Tuple (C{sin(B{rad})}, C{cos(B{rad})}, C{tan(B{rad})}).

780

781 @see: Function L{sincos2}.

782 '''

783 rad %= PI2 # see ._wrap comments

784 if rad:

785 s, c = sincos2(rad)

786 t = (s / c) if c else (NINF if s < 0

787 else (INF if s > 0 else s))

788 else: # sin(-0.0) == tan(-0.0) = -0.0

789 c = _1_0

790 s = t = NEG0 if isneg0(rad) else _0_0

791 return s, c, t

792

793

794 def SM2m(sm):

795 '''Convert statute miles to meter (m).

796

797 @arg sm: Value in statute miles (C{scalar}).

798

799 @return: Value in meter (C{float}).

800

801 @raise ValueError: Invalid B{C{sm}}.

802 '''

803 return Meter(Float(sm=sm) * _M_SM)

804

805

806 def tan_2(rad, **semi): # edge=1

807 '''Compute the tangent of half angle.

808

809 @arg rad: Angle (C{radians}).

810 @kwarg semi: Angle or edge name and index

811 for semi-circular error.

812

813 @return: M{tan(rad / 2)} (C{float}).

814

815 @raise ValueError: If B{C{rad}} is semi-circular

816 and B{C{semi}} is given.

817 '''

818 # .formy.excessKarney_, .sphericalTrigonometry.areaOf

819 if semi and isnear0(fabs(rad) - PI):

820 for n, v in semi.items():

821 break

822 n = _SPACE_(n, _radians_) if not isint(v) else \

823 _SPACE_(_MODS.streprs.Fmt.SQUARE(**semi), _edge_)

824 raise _ValueError(n, rad, txt=_semi_circular_)

825

826 return tan(rad * _0_5)

827

828

829 def tand(deg, **error_kwds):

830 '''Return the C{tangent} of an angle in C{degrees}.

831

832 @arg deg: Angle (C{degrees}).

833 @kwarg error_kwds: Error to raise (C{ValueError}).

834

835 @return: C{tan(B{deg})}.

836

837 @raise ValueError: If L{pygeodesy.isnear0}C{(cos(B{deg})}.

838 '''

839 s, c = sincos2d(deg)

840 if isnear0(c):

841 raise _valueError(tand, deg, **error_kwds)

842 return s / c

843

844

845 def tand_(*degs, **error_kwds):

846 '''Return the C{tangent} of angle(s) in C{degrees}.

847

848 @arg degs: One or more angles (C{degrees}).

849 @kwarg error_kwds: Error to raise (C{ValueError}).

850

851 @return: Yield the C{tan(B{deg})} for each angle.

852

853 @raise ValueError: See L{pygeodesy.tand}.

854 '''

855 for d in degs:

856 yield tand(d, **error_kwds)

857

858

859 def tanPI_2_2(rad):

860 '''Compute the tangent of half angle, 90 degrees rotated.

861

862 @arg rad: Angle (C{radians}).

863

864 @return: M{tan((rad + PI/2) / 2)} (C{float}).

865 '''

866 return tan((rad + PI_2) * _0_5)

867

868

869 def toise2m(toises):

870 '''Convert French U{toises<https://WikiPedia.org/wiki/Toise>} to meter.

871

872 @arg toises: Value in toises (C{scalar}).

873

874 @return: Value in C{meter} (C{float}).

875

876 @raise ValueError: Invalid B{C{toises}}.

877

878 @see: Function L{fathom2m}.

879 '''

880 return Meter(Float(toises=toises) * _M_TOISE)

881

882

883 def truncate(x, ndigits=None):

884 '''Truncate to the given number of digits.

885

886 @arg x: Value to truncate (C{scalar}).

887 @kwarg ndigits: Number of digits (C{int}),

888 aka I{precision}.

889

890 @return: Truncated B{C{x}} (C{float}).

891

892 @see: Python function C{round}.

893 '''

894 if isint(ndigits):

895 p = _10_0**ndigits

896 x = int(x * p) / p

897 return x

898

899

900 def unroll180(lon1, lon2, wrap=True):

901 '''Unroll longitudinal delta and wrap longitude in degrees.

902

903 @arg lon1: Start longitude (C{degrees}).

904 @arg lon2: End longitude (C{degrees}).

905 @kwarg wrap: If C{True}, wrap and unroll to the M{(-180..+180]}

906 range (C{bool}).

907

908 @return: 2-Tuple C{(B{lon2}-B{lon1}, B{lon2})} unrolled (C{degrees},

909 C{degrees}).

910

911 @see: Capability C{LONG_UNROLL} in U{GeographicLib

912 <https://GeographicLib.SourceForge.io/html/python/interface.html#outmask>}.

913 '''

914 d = lon2 - lon1

915 if wrap and fabs(d) > _180_0:

916 u = _wrap(d, _180_0, _360_0)

917 if u != d:

918 return u, (lon1 + u)

919 return d, lon2

920

921

922 def _unrollon(p1, p2, wrap=False): # unroll180 == .karney._unroll2

923 '''(INTERNAL) Wrap/normalize, unroll and replace longitude.

924 '''

925 lat, lon = p2.lat, p2.lon

926 if wrap and _Wrap.normal:

927 lat, lon = _Wrap.latlon(lat, lon)

928 _, lon = unroll180(p1.lon, lon, wrap=True)

929 if lat != p2.lat or fabs(lon - p2.lon) > EPS:

930 p2 = p2.dup(lat=lat, lon=wrap180(lon))

931 # p2 = p2.copy(); p2.latlon = lat, wrap180(lon)

932 return p2

933

934

935 def _unrollon3(p1, p2, p3, wrap=False):

936 '''(INTERNAL) Wrap/normalize, unroll 2 points.

937 '''

938 w = wrap

939 if w:

940 w = _Wrap.normal

941 p2 = _unrollon(p1, p2, wrap=w)

942 p3 = _unrollon(p1, p3, wrap=w)

943 p2 = _unrollon(p2, p3)

944 return p2, p3, w # was wrapped?

945

946

947 def unrollPI(rad1, rad2, wrap=True):

948 '''Unroll longitudinal delta and wrap longitude in radians.

949

950 @arg rad1: Start longitude (C{radians}).

951 @arg rad2: End longitude (C{radians}).

952 @kwarg wrap: If C{True}, wrap and unroll to the M{(-PI..+PI]}

953 range (C{bool}).

954

955 @return: 2-Tuple C{(B{rad2}-B{rad1}, B{rad2})} unrolled

956 (C{radians}, C{radians}).

957

958 @see: Capability C{LONG_UNROLL} in U{GeographicLib

959 <https://GeographicLib.SourceForge.io/html/python/interface.html#outmask>}.

960 '''

961 r = rad2 - rad1

962 if wrap and fabs(r) > PI:

963 u = _wrap(r, PI, PI2)

964 if u != r:

965 return u, (rad1 + u)

966 return r, rad2

967

968

969 def _valueError(where, x, **kwds):

970 '''(INTERNAL) Return a C{_ValueError}.

971 '''

972 x = _MODS.streprs.Fmt.PAREN(where.__name__, x)

973 return _ValueError(x, **kwds)

974

975

976 class _Wrap(object):

977

978 _normal = False # default

979

980 @property

981 def normal(self):

982 '''Get the current L{normal} setting (C{True},

983 C{False} or C{None}).

984 '''

985 return self._normal

986

987 @normal.setter # PYCHOK setter!

988 def normal(self, setting):

989 '''Set L{normal} to C{True}, C{False} or C{None}.

990 '''

991 t = {True: (_MODS.formy.normal, _MODS.formy.normal_),

992 False: (self.wraplatlon, self.wraphilam),

993 None: (_passargs, _passargs)}.get(setting, ())

994 if t:

995 self.latlon, self.philam = t

996 self._normal = setting

997

998 def latlonDMS2(self, lat, lon, **DMS2_kwds):

999 if isstr(lat) or isstr(lon):

1000 kwds = _xkwds(DMS2_kwds, clipLon=0, clipLat=0)

1001 lat, lon = _MODS.dms.parseDMS2(lat, lon, **kwds)

1002 return self.latlon(lat, lon)

1003

1004# def normalatlon(self, *latlon):

1005# return _MODS.formy.normal(*latlon)

1006

1007# def normalamphi(self, *philam):

1008# return _MODS.formy.normal_(*philam)

1009

1010 def wraplatlon(self, lat, lon):

1011 return wrap90(lat), wrap180(lon)

1012

1013 latlon = wraplatlon # default

1014

1015 def latlon3(self, lon1, lat2, lon2, wrap):

1016 if wrap:

1017 lat2, lon2 = self.latlon(lat2, lon2)

1018 lon21, lon2 = unroll180(lon1, lon2)

1019 else:

1020 lon21 = lon2 - lon1

1021 return lon21, lat2, lon2

1022

1023 def _latlonop(self, wrap):

1024 if wrap and self._normal is not None:

1025 return self.latlon

1026 else:

1027 return _passargs

1028

1029 def wraphilam(self, phi, lam):

1030 return wrapPI_2(phi), wrapPI(lam)

1031

1032 philam = wraphilam # default

1033

1034 def philam3(self, lam1, phi2, lam2, wrap):

1035 if wrap:

1036 phi2, lam2 = self.philam(phi2, lam2)

1037 lam21, lam2 = unrollPI(lam1, lam2)

1038 else:

1039 lam21 = lam2 - lam1

1040 return lam21, phi2, lam2

1041

1042 def _philamop(self, wrap):

1043 if wrap and self._normal is not None:

1044 return self.philam

1045 else:

1046 return _passargs

1047

1048 def wraphilam(self, phi, lam):

1049 return wrapPI_2(phi), wrapPI(lam)

1050

1051 def point(self, ll, wrap=True): # in .points._fractional, -.PointsIter.iterate, ...

1052 '''Return C{ll} or a copy, I{normalized} or I{wrap}'d.

1053 '''

1054 if wrap and self._normal is not None:

1055 lat, lon = ll.latlon

1056 if fabs(lon) > 180 or fabs(lat) > 90:

1057 _n = self.latlon

1058 ll = ll.copy(name=_n.__name__)

1059 ll.latlon = _n(lat, lon)

1060 return ll

1061

1062_Wrap = _Wrap() # PYCHOK singleton

1063

1064

1065 def _wrap(angle, wrap, modulo):

1066 '''(INTERNAL) Angle wrapper M{((wrap-modulo)..+wrap]}.

1067

1068 @arg angle: Angle (C{degrees}, C{radians} or C{grades}).

1069 @arg wrap: Range (C{degrees}, C{radians} or C{grades}).

1070 @arg modulo: Upper limit (360 C{degrees}, PI2 C{radians} or 400 C{grades}).

1071

1072 @return: The B{C{angle}}, wrapped (C{degrees}, C{radians} or C{grades}).

1073 '''

1074 a = float(angle)

1075 if not (wrap - modulo) <= a < wrap:

1076 # math.fmod(-1.5, 3.14) == -1.5, but -1.5 % 3.14 == 1.64

1077 # math.fmod(-1.5, 360) == -1.5, but -1.5 % 360 == 358.5

1078 a %= modulo

1079 if a > wrap:

1080 a -= modulo

1081 return a

1082

1083

1084 def wrap90(deg):

1085 '''Wrap degrees to M{[-270..+90]}.

1086

1087 @arg deg: Angle (C{degrees}).

1088

1089 @return: Degrees, wrapped (C{degrees90}).

1090 '''

1091 return _wrap(deg, _90_0, _360_0)

1092

1093

1094 def wrap180(deg):

1095 '''Wrap degrees to M{[-180..+180]}.

1096

1097 @arg deg: Angle (C{degrees}).

1098

1099 @return: Degrees, wrapped (C{degrees180}).

1100 '''

1101 return _wrap(deg, _180_0, _360_0)

1102

1103

1104 def wrap360(deg): # see .streprs._umod_360

1105 '''Wrap degrees to M{[0..+360)}.

1106

1107 @arg deg: Angle (C{degrees}).

1108

1109 @return: Degrees, wrapped (C{degrees360}).

1110 '''

1111 return _wrap(deg, _360_0, _360_0)

1112

1113

1114 def wrapPI(rad):

1115 '''Wrap radians to M{[-PI..+PI]}.

1116

1117 @arg rad: Angle (C{radians}).

1118

1119 @return: Radians, wrapped (C{radiansPI}).

1120 '''

1121 return _wrap(rad, PI, PI2)

1122

1123

1124 def wrapPI2(rad):

1125 '''Wrap radians to M{[0..+2PI)}.

1126

1127 @arg rad: Angle (C{radians}).

1128

1129 @return: Radians, wrapped (C{radiansPI2}).

1130 '''

1131 return _wrap(rad, PI2, PI2)

1132

1133

1134 def wrapPI_2(rad):

1135 '''Wrap radians to M{[-3PI/2..+PI/2]}.

1136

1137 @arg rad: Angle (C{radians}).

1138

1139 @return: Radians, wrapped (C{radiansPI_2}).

1140 '''

1141 return _wrap(rad, PI_2, PI2)

1142

1143

1144 def wrap_normal(*normal):

1145 '''Define the operation for the keyword argument C{B{wrap}=True},

1146 across L{pygeodesy}: I{wrap}, I{normalize} or I{no-op}. For

1147 backward compatibility, the default is I{wrap}.

1148

1149 @arg normal: If C{True}, I{normalize} lat- and longitude using

1150 L{normal} or L{normal_}, if C{False}, I{wrap} the

1151 lat- and longitude individually by L{wrap90} or

1152 L{wrapPI_2} respectively L{wrap180}, L{wrapPI} or

1153 if C{None}, leave lat- and longitude I{unchanged}.

1154 Do not supply any value to get the current setting.

1155

1156 @return: The previous L{wrap_normal} setting (C{bool} or C{None}).

1157 '''

1158 t = _Wrap.normal

1159 if normal:

1160 _Wrap.normal = normal[0]

1161 return t

1162

1163

1164 def yard2m(yards):

1165 '''Convert I{UK} yards to meter.

1166

1167 @arg yards: Value in yards (C{scalar}).

1168

1169 @return: Value in C{meter} (C{float}).

1170

1171 @raise ValueError: Invalid B{C{yards}}.

1172 '''

1173 return Float(yards=yards) * _M_YARD_UK

1174

1175# **) MIT License

1176#

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

1178#

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

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

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

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

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

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

1185#

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

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

1188#

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

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

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

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

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

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

1195# OTHER DEALINGS IN THE SOFTWARE.