1# -*- coding: utf-8 -*-

2

3 u'''Single-instance C{float} and C{int} constants across C{pygeodesy}

4 modules and related functions L{pygeodesy.float_}, L{pygeodesy.isclose},

5 L{pygeodesy.isfinite}, L{pygeodesy.isinf}, L{pygeodesy.isint0},

6 L{pygeodesy.isnan}, L{pygeodesy.isnear0}, L{pygeodesy.isnear1},

7 L{pygeodesy.isnear90}, L{pygeodesy.isneg0}, L{pygeodesy.isninf},

8 L{pygeodesy.isnon0} and L{pygeodesy.remainder}.

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 _0_0, _copysign, isbool, iscomplex, isint

14 from pygeodesy.errors import _ALL_LAZY, _xError, _xError2

15 from pygeodesy.interns import _INF_, _NAN_, _sets_, _UNDER_

16# from pygeodesy.lazily import _ALL_LAZY # from .errors

17# from pygeodesy.streprs import Fmt # from .unitsBase

18 from pygeodesy.unitsBase import Float, Fmt, Int, Radius

19

20 from math import fabs, isinf, isnan, pi as _pi, sqrt

21 try:

22 from math import inf as _inf, nan as _nan # PYCHOK Python 3+

23 except ImportError: # Python 2-

24 _inf, _nan = float(_INF_), float(_NAN_)

25

26__all__ = _ALL_LAZY.constants

27__version__ = '23.03.19'

28

29

30 def _Float(**name_arg):

31 '''(INTERNAL) New named, cached C{Float}.

32 '''

33 n, arg = name_arg.popitem()

34 return Float(_float(arg), name=n)

35

36

37 def _Radius(**name_arg):

38 '''(INTERNAL) New named, cached C{Radius}.

39 '''

40 n, arg = name_arg.popitem()

41 return Radius(_float(arg), name=n)

42

43

44 def float_(*fs, **sets):

45 '''Get scalars as C{float} or I{intern}'ed C{float}.

46

47 @arg fs: One more values (C{scalar}), all positional.

48 @kwarg sets: Use keyword argument C{B{sets}=True} to

49 C{intern} each B{C{fs}}, otherwise don't.

50

51 @return: Single C{float} if only one B{C{fs}} is given,

52 otherwise a tuple of C{float}s.

53

54 @raise TypeError: Some B{C{fs}} not C{scalar}.

55 '''

56 _f = _floats.setdefault if sets.get(_sets_, False) else \

57 _floats.get

58 fl = []

59 try:

60 for i, f in enumerate(fs):

61 f = float(f)

62 fl.append(_f(f, f))

63 except Exception as x:

64 _E, t = _xError2(x)

65 _fs_i = Fmt.SQUARE(fs=i)

66 raise _E(_fs_i, f, txt=t)

67 return fl[0] if len(fl) == 1 else tuple(fl)

68

69

70 def _float(f): # in .datums, .ellipsoids, ...

71 '''(INTERNAL) Cache initial C{float}s.

72 '''

73 f = float(f)

74 return _floats.setdefault(f, f) # PYCHOK del _floats

75

76

77 def float0(*xs):

78 '''Yield C{B{x}s} as a non-NEG0 C{float}.

79 '''

80 for x in xs:

81 yield float(x) or _0_0

82

83

84 def _float0(f): # in .resections, .vector3dBase, ...

85 '''(INTERNAL) Return C{float(B{f})} or C{INT0}.

86 '''

87 if f:

88 f = float(f)

89 f = _floats.get(f, f)

90 elif f is not INT0:

91 f = _0_0

92 return f

93

94

95 def _floatuple(*fs):

96 '''(INTERNAL) Cache a tuple of C{float}s.

97 '''

98 return tuple(map(_float, fs))

99

100

101 def _1_over(x):

102 '''(INTERNAL) Return reciprocal C{1 / B{x}}.

103 '''

104 return _float(_1_0 / float(x))

105

106

107_floats = {} # PYCHOK floats cache, in .__main__

108# _float = float # PYCHOK expected

109# del _floats # XXX zap floats cache never

110

111_0_0 = _float( _0_0) # PYCHOK expected

112_0_0_1T = _0_0, # PYCHOK 1-tuple

113_0_0_9T = _0_0_1T * 9 # PYCHOK 9-tuple

114_0_0001 = _float( 0.0001) # PYCHOK expected

115_0_001 = _float( 0.001) # PYCHOK expected

116_0_01 = _float( 0.01) # PYCHOK expected

117_0_1 = _float( 0.1) # PYCHOK expected

118_0_125 = _float( 0.125) # PYCHOK expected

119_0_25 = _float( 0.25) # PYCHOK expected

120_0_26 = _float( 0.26) # PYCHOK expected

121_0_5 = _float( 0.5) # PYCHOK expected

122_1_0 = _float( 1) # PYCHOK expected

123_1_0_1T = _1_0, # PYCHOK 1-tuple

124_1_5 = _float( 1.5) # PYCHOK expected

125_1_75 = _float( 1.75) # PYCHOK expected

126_2_0 = _float( 2) # PYCHOK expected

127_3_0 = _float( 3) # PYCHOK expected

128_4_0 = _float( 4) # PYCHOK expected

129_5_0 = _float( 5) # PYCHOK expected

130_6_0 = _float( 6) # PYCHOK expected

131_8_0 = _float( 8) # PYCHOK expected

132_9_0 = _float( 9) # PYCHOK expected

133_10_0 = _float( 10) # PYCHOK expected

134_16_0 = _float( 16) # PYCHOK expected

135_32_0 = _float( 32) # PYCHOK expected

136_60_0 = _float( 60) # PYCHOK expected

137_90_0 = _float( 90) # PYCHOK expected

138_100_0 = _float( 100) # PYCHOK expected

139_180_0 = _float( 180) # PYCHOK expected

140_270_0 = _float( 270) # PYCHOK expected

141_360_0 = _float( 360) # PYCHOK expected

142_400_0 = _float( 400) # PYCHOK expected

143_720_0 = _float( 720) # PYCHOK expected

144_1000_0 = _float(1000) # PYCHOK expected

145_3600_0 = _float(3600) # PYCHOK expected

146

147_N_0_0 = float( '-0.0') # PYCHOK NOT _float!

148_N_0_5 = _float( -_0_5) # PYCHOK expected

149_N_1_0 = _float( -_1_0) # PYCHOK expected

150_N_2_0 = _float( -_2_0) # PYCHOK expected

151_N_90_0 = _float( -_90_0) # PYCHOK expected

152_N_180_0 = _float(-_180_0) # PYCHOK expected

153

154_M_KM = _1000_0 # meter per Kilo meter, see .utily

155_M_NM = _float(1852.0) # meter per Nautical Mile

156_M_SM = _float(1609.344) # meter per Statute Mile

157

158 try:

159 from sys import float_info as _f_i

160 # @see: <https://NumPy.org/doc/stable/reference/generated/numpy.finfo.html>

161 DIG = Int( DIG =_f_i.dig) # PYCHOK system's float decimal digits

162 EPS = _Float(EPS =_f_i.epsilon) # PYCHOK system's EPSilon

163 MANT_DIG = Int( MANT_DIG=_f_i.mant_dig) # PYCHOK system's float mantissa bits

164 MAX = _Float(MAX =_f_i.max) # PYCHOK system's MAX float 1.7976931348623157e+308

165# MAX_EXP = Int( MAX_EXP =_f_i.max_exp) # PYTHON system's max base 2 exponent

166 MIN = _Float(MIN =_f_i.min) # PYCHOK system's MIN float 2.2250738585072014e-308

167# MIN_EXP = Int( MIN_EXP =_f_i.min_exp) # PYTHON system's min base 2 exponent

168# RADIX = Int( RADIX =_f_i.radix) # PYTHON system's float base

169 del _f_i

170 except ImportError: # PYCHOK no cover

171 DIG = Int( DIG =15) # PYCHOK system's 64-bit float decimal digits

172 EPS = _Float(EPS =2.220446049250313e-16) # PYCHOK EPSilon 2**-52, M{EPS +/- 1 != 1}

173 MANT_DIG = Int( MANT_DIG=53) # PYCHOK float mantissa bits ≈ 53 (C{int})

174 MAX = _Float(MAX =pow(_2_0, 1023) * (_2_0 - EPS)) # PYCHOK ≈ 10**308

175# MAX_EXP = Int( MAX_ESP =1024) # 308 base 10

176 MIN = _Float(MIN =pow(_2_0, -1021)) # PYCHOK ≈ 10**-308

177# MIN_EXP = Int(MIN_EXP =-1021) # -307 base 10

178# RADIX = Int(Radix =2) # base

179

180 EPS0 = _Float( EPS0 = EPS**2) # PYCHOK near-/non-zero comparison 4.930381e-32, or EPS or EPS_2

181 EPS02 = _Float( EPS02 = EPS**4) # PYCHOK near-zero-squared comparison 2.430865e-63

182 EPS_2 = _Float( EPS_2 = EPS / _2_0) # PYCHOK ≈ 1.110223024625e-16

183 EPS1 = _Float( EPS1 =_1_0 - EPS) # PYCHOK ≈ 0.9999999999999998

184 EPS2 = _Float( EPS2 = EPS * _2_0) # PYCHOK ≈ 4.440892098501e-16

185 EPS4 = _Float( EPS4 = EPS * _4_0) # PYCHOK ≈ 8.881784197001e-16

186# _1EPS = _Float(_1EPS =_1_0 + EPS) # PYCHOK ≈ 1.0000000000000002

187_1_EPS = _Float(_1_EPS =_1_0 / EPS) # PYCHOK = 4503599627370496.0

188# _2_EPS = _Float(_2_EPS =_2_0 / EPS) # PYCHOK = 9007199254740992.0

189_EPS2e4 = _Float(_EPS2e4 = EPS2 * 1.e4) # PYCHOK ≈ 4.440892098501e-12

190_EPS4e8 = _Float(_EPS4e8 = EPS4 * 1.e8) # PYCHOK ≈ 8.881784197001e-08

191_EPSmin = _Float(_EPSmin = sqrt(MIN)) # PYCHOK = 1.49166814624e-154

192_EPSqrt = _Float(_EPSqrt = sqrt(EPS)) # PYCHOK = 1.49011611938e5-08

193_EPStol = _Float(_EPStol =_EPSqrt * _0_1) # PYCHOK = 1.49011611938e5-09 == sqrt(EPS * _0_01)

194

195_89_999_ = _Float(_89_999_= EPS1 * _90_0) # just below 90.0

196# <https://Numbers.Computation.Free.FR/Constants/Miscellaneous/digits.html>

197_1__90 = _Float(_1__90 =_1_0 / _90_0) # PYCHOK = 0.011_111_111_111_111_111_111_111_111_111_111_111_111_111_111_11111

198_2__PI = _Float(_2__PI =_2_0 / _pi) # PYCHOK = 0.636_619_772_367_581_343_075_535_053_490_057_448_137_838_582_96182

199

200_1_16th = _Float(_1_16th =_1_0 / _16_0) # PYCHOK in .ellipsoids, .karney

201_1_64th = _Float(_1_64th =_1_0 / 64) # PYCHOK in .elliptic, pow(2.0, -6)

202_1_3rd = _Float(_1_3rd =_1_0 / _3_0) # PYCHOK in .fmath

203_2_3rd = _Float(_2_3rd =_2_0 / _3_0) # PYCHOK in .fmath

204

205_K0_UTM = _Float(_K0_UTM = 0.9996) # PYCHOK in .etm, .ktm, .utm, UTM scale at central meridian

206# sqrt(2) <https://WikiPedia.org/wiki/Square_root_of_2>

207# 1.414213562373095_048_801_688_724_209_698_078_569_671_875_376_948_073_176_679_737_99

208# _1SQRT2= _Float(_1SQRT2 =sqrt(_2_0) + 1)

209_SQRT2_2 = _Float(_SQRT2_2=sqrt(_0_5)) # PYCHOK = 0.707106781186547_6 == sqrt(2) / 2

210

211 INF = Float(INF =_inf) # PYCHOK INFinity, see function L{isinf}, L{isfinite}, NOT _Float!

212 INT0 = Int( INT0= 0) # PYCHOK unique int(0) instance, see .fsums, useZ=False

213 NAN = Float(NAN =_nan) # PYCHOK Not-A-Number, see function L{isnan}, NOT _Float!

214 NEG0 = Float(NEG0=_N_0_0) # PYCHOK NEGative 0.0, see function L{isneg0}, NOT _Float!

215 NINF = Float(NINF=-INF) # PYCHOK Negative INFinity, NOT _Float!

216

217 PI = _Float(PI =_pi)

218 PI2 = _Float(PI2 =_pi * _2_0) # PYCHOK Two PI, M{PI * 2} aka I{Tau}

219 PI_2 = _Float(PI_2 =_pi / _2_0) # PYCHOK Half PI, M{PI / 2}

220 PI3 = _Float(PI3 =_pi * _3_0) # PYCHOK Three PI, M{PI * 3}

221 PI3_2 = _Float(PI3_2=_pi * _1_5) # PYCHOK PI and a half, M{PI * 3 / 2}

222 PI_3 = _Float(PI_3 =_pi / _3_0) # PYCHOK One third PI, M{PI / 3}

223 PI4 = _Float(PI4 =_pi * _4_0) # PYCHOK Four PI, M{PI * 4}

224 PI_4 = _Float(PI_4 =_pi / _4_0) # PYCHOK Quarter PI, M{PI / 4}

225

226 R_MA = _Radius(R_MA=6378137.0) # PYCHOK equatorial earth radius (C{meter}), WGS84, EPSG:3785

227 R_MB = _Radius(R_MB=6356752.3) # PYCHOK polar earth radius (C{meter}), WGS84, EPSG:3785

228 R_M = _Radius(R_M =6371008.771415) # PYCHOK mean, spherical earth radius (C{meter})

229 R_KM = _Radius(R_KM=R_M / _M_KM) # PYCHOK mean, spherical earth radius (C{KM}, kilo meter)

230 R_NM = _Radius(R_NM=R_M / _M_NM) # PYCHOK mean, spherical earth radius (C{NM}, nautical miles)

231 R_SM = _Radius(R_SM=R_M / _M_SM) # PYCHOK mean, spherical earth radius (C{SM}, statute miles)

232# See <https://www.EdWilliams.org/avform.htm>, <https://www.DTIC.mil/dtic/tr/fulltext/u2/a216843.pdf>

233# and <https://GitHub.com/NASA/MultiDop/blob/master/src/share/man/man3/geog_lib.3> based on the

234# International Standard Nautical Mile of 1,852 meter (1' latitude)

235 R_FM = _Radius(R_FM=6371000.0) # PYCHOK former FAI Sphere earth radius (C{meter})

236 R_GM = _Radius(R_GM=6371230.0) # PYCHOK avg. radius, distance to geoid surface (C{meter})

237# <http://Wiki.GIS.com/wiki/index.php/Ellipsoidal_quadratic_mean_radius>

238 R_QM = _Radius(R_QM=6372797.560856) # PYCHOK earth' quadratic mean radius (C{meter})

239# Rtri= _Radius(Rtri=6372797.5559594) # PYCHOK Rtriaxial quadratic mean radius (C{meter}), WGS84

240# Rbi = _Radius(Rbi =6367453.6345163) # PYCHOK Rbiaxial quadratic mean radius (C{meter}), WGS84

241 R_VM = _Radius(R_VM=6366707.0194937) # PYCHOK aViation/naVigation earth radius (C{meter})

242# R_AU= Meter( R_AU=149597870700.0) # PYCHOK <https://WikiPedia.org/wiki/Astronomical_unit>

243

244_INF_NAN_NINF = INF, NAN, NINF

245

246

247 def _0_0s(n):

248 '''(INTERNAL) Return an C{B{n}-tuple} of C{_0_0} zeros.

249 '''

250 return _0_0_9T[:n] if 0 <= n <= len(_0_0_9T) else (_0_0_1T * n)

251

252

253 try:

254 from math import isclose as _isclose

255 except ImportError: # Python 3.4-

256

257 def _isclose(a, b, rel_tol=1e-9, abs_tol=0):

258 '''Mimick Python 3.5+ C{math.isclose}.

259 '''

260 t, d = abs_tol, fabs(a - b)

261 if d > t:

262 r = max(fabs(a), fabs(b)) * rel_tol

263 t = max(r, t)

264 return d <= t

265

266

267 def isclose(a, b, rel_tol=1e-12, abs_tol=EPS0):

268 '''Like C{math.isclose}, but with defaults such

269 that C{isclose(0, EPS0)} is C{True} by default.

270 '''

271 return _isclose(a, b, rel_tol=rel_tol, abs_tol=abs_tol)

272

273

274 try:

275 from math import isfinite as _isfinite # in .ellipsoids, .fsums, .karney

276 except ImportError: # Python 3.1-

277

278 def _isfinite(x):

279 '''Mimick Python 3.2+ C{math.isfinite}.

280 '''

281 return not (isinf(x) or isnan(x))

282

283

284 def isfinite(obj):

285 '''Check a finite C{scalar} or C{complex} value.

286

287 @arg obj: Value (C{scalar} or C{complex}).

288

289 @return: C{False} if B{C{obj}} is C{INF}, C{NINF}

290 or C{NAN}, C{True} otherwise.

291

292 @raise TypeError: Non-scalar and non-complex B{C{obj}}.

293 '''

294 try:

295 return (obj not in _INF_NAN_NINF) and _isfinite(obj)

296 except Exception as x:

297 if iscomplex(obj): # _isfinite(complex) thows TypeError

298 return isfinite(obj.real) and isfinite(obj.imag)

299 raise _xError(x, Fmt.PAREN(isfinite.__name__, obj))

300

301

302 def isint0(obj, both=False):

303 '''Check for L{INT0} or C{int(0)} value.

304

305 @arg obj: The object (any C{type}).

306 @kwarg both: If C{true}, also check C{float(0)} (C{bool}).

307

308 @return: C{True} if B{C{obj}} is L{INT0}, C{int(0)} or

309 C{float(0)}, C{False} otherwise.

310 '''

311 return (obj is INT0 or obj is int(0) or bool(both and

312 (not obj) and isint(obj, both=True))) and not isbool(obj)

313

314

315 def isnear0(x, eps0=EPS0):

316 '''Is B{C{x}} near zero within a tolerance?

317

318 @arg x: Value (C{scalar}).

319 @kwarg eps0: Near-zero tolerance (C{EPS0}).

320

321 @return: C{True} if C{abs(B{x}) < B{eps0}},

322 C{False} otherwise.

323

324 @see: Function L{isnon0}.

325 '''

326 return bool(eps0 > x > -eps0)

327

328

329 def isnear1(x, eps1=EPS0):

330 '''Is B{C{x}} near one within a tolerance?

331

332 @arg x: Value (C{scalar}).

333 @kwarg eps1: Near-one tolerance (C{EPS0}).

334

335 @return: C{isnear0(B{x} - 1, eps0=B{eps1})}.

336

337 @see: Function L{isnear0}.

338 '''

339 return bool(eps1 > (x - _1_0) > -eps1)

340

341

342 def isnear90(x, eps90=EPS0):

343 '''Is B{C{x}} near one within a tolerance?

344

345 @arg x: Value (C{scalar}).

346 @kwarg eps90: Near-90 tolerance (C{EPS0}).

347

348 @return: C{isnear0(B{x} - 90)}.

349

350 @see: Function L{isnear0}.

351 '''

352 return isnear0(x - _90_0, eps0=eps90)

353

354

355 def isneg0(x):

356 '''Check for L{NEG0}, negative C{0.0}.

357

358 @arg x: Value (C{scalar}).

359

360 @return: C{True} if B{C{x}} is C{NEG0} or C{-0.0},

361 C{False} otherwise.

362 '''

363 return x in (_0_0, NEG0) and _copysign(1, x) < 0

364# and str(x).startswith(_MINUS_)

365

366

367 def isninf(x):

368 '''Check for L{NINF}, negative C{INF}.

369

370 @arg x: Value (C{scalar}).

371

372 @return: C{True} if B{C{x}} is C{NINF} or C{-inf},

373 C{False} otherwise.

374 '''

375 return x is NINF or (x < 0 and not isfinite(x))

376

377

378 def isnon0(x, eps0=EPS0):

379 '''Is B{C{x}} non-zero with a tolerance?

380

381 @arg x: Value (C{scalar}).

382 @kwarg eps0: Non-zero tolerance (C{EPS0}).

383

384 @return: C{True} if C{abs(B{x}) > B{eps0}},

385 C{False} otherwise.

386

387 @see: Function L{isnear0}.

388 '''

389 return not bool(eps0 > x > -eps0) # not isnear0

390

391

392 def _off90(lat):

393 '''(INTERNAL) Off 90.0 for .gars and .wgrs.

394 '''

395 return max(min(lat, _89_999_), -_89_999_)

396

397

398 try:

399 from math import remainder

400 except ImportError: # Python 3.6-

401 from math import fmod as _fmod

402

403 def remainder(x, y):

404 '''Mimick Python 3.7+ C{math.remainder}.

405 '''

406 if isnan(y):

407 x = NAN

408 elif x and not isnan(x):

409 y = fabs(y)

410 x = _fmod(x, y)

411 h = _0_5 * y

412 if x >= h:

413 x -= y

414 elif x < -h:

415 x += y

416 return x # keep signed 0.0

417

418

419 def _umod_360(deg):

420 '''(INTERNAL) Non-negative C{deg} modulo 360, basic C{.utily.wrap360}.

421 '''

422 return (deg % _360_0) or _0_0

423

424

425 def _umod_PI2(rad):

426 '''(INTERNAL) Non-negative C{rad} modulo PI2, basic C{.utily.wrapPI2}.

427 '''

428 return (rad % PI2) or _0_0

429

430

431 if __name__ == '__main__':

432

433 from pygeodesy.errors import itemsorted

434 from pygeodesy.lazily import printf

435

436 t = n = v = []

437 for n, v in itemsorted(locals()):

438 if isinstance(v, (Float, Int, Radius)):

439 printf('%9s: %r', n, v.toRepr(std=False))

440 if v.name != n:

441 raise AssertionError('%r != %r' % (n, v))

442 if v.name is not n:

443 raise AssertionError('%r is not %r' % (n, v))

444 if not n.startswith(_UNDER_):

445 t.append(n)

446 t.append(float_.__name__)

447 printf('__all__ = %r', tuple(t))

448

449# **) MIT License

450#

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

452#

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

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

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

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

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

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

459#

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

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

462#

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

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

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

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

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

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

469# OTHER DEALINGS IN THE SOFTWARE.