1

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

3

4 u'''(INTERNAL) Private L{geodesicx} base class, functions and constants.

5

6 Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2008-2023)

7 and licensed under the MIT/X11 License. For more information, see the

8 U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.

9'''

10

11 from pygeodesy.basics import isodd, _MODS

12 from pygeodesy.constants import _EPSmin as _TINY, _0_0, _100_0

13 from pygeodesy.errors import _or, _xkwds_item2

14 from pygeodesy.fmath import hypot as _hypot

15 from pygeodesy.karney import _CapsBase, GeodesicError, _2cos2x, \

16 _norm2, _sincos2d, _sum2_

17# from pygeodesy.lazily import _MODS, printf # .basics, _MODS

18

19 from math import fabs, ldexp as _ldexp

20

21__all__ = ()

22__version__ = '24.06.16'

23

24# valid C{nC4}s and C{C4order}s, see _xnC4 below

25_nC4s = {24: 2900, 27: 4032, 30: 5425}

26# assert (_TINY * EPS) > 0 and (_TINY + EPS) == EPS # underflow guard

27

28

29 class _GeodesicBase(_CapsBase): # in .geodsolve

30 '''(INTERNAL) Base class for C{[_]Geodesic*Exact}.

31 '''

32# def toRepr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature

33# '''Return this C{GeodesicExact*} items string.

34#

35# @kwarg prec: The C{float} precision, number of decimal digits (0..9).

36# Trailing zero decimals are stripped for B{C{prec}} values

37# of 1 and above, but kept for negative B{C{prec}} values.

38# @kwarg sep: Separator to join (C{str}).

39#

40# @return: C{GeodesicExact*} (C{str}).

41# '''

42# return Fmt.PAREN(self.named, self.toStr(prec=prec, sep=sep))

43 pass

44

45

46 class _Gfloats(dict):

47 '''(INTERNAL) Uniquify floats.

48 '''

49 n = 0 # total number of floats

50 nC4 = 0

51

52 def __init__(self, nC4): # PYCHOK signature

53 self.nC4 = nC4

54

55 def __call__(self, fs):

56 '''Return a tuple of "uniquified" floats.

57 '''

58 self.n += len(fs)

59 _f = self.setdefault

60 return tuple(_f(f, f) for f in map(float, fs)) # PYCHOK as attr

61

62 def prints(self):

63 n, u = self.n, len(self.keys())

64 d = (n - u) * _100_0 / n

65 _MODS.lazily.printf('_CX_%d: n=%d, u=%d, d=%.1f%%', self.nC4, n, u, d) # XXX

66

67

68 def _cosSeries(c4s, sx, cx): # PYCHOK shared .geodesicx.gx and -.gxline

69 '''(INTERNAL) I{Karney}'s cosine series expansion using U{Clenshaw

70 summation<https://WikiPedia.org/wiki/Clenshaw_algorithm>}.

71 '''

72 ar = _2cos2x(cx, sx)

73 y0 = t0 = y1 = t1 = _0_0

74 c4 = list(c4s)

75 _c4 = c4.pop

76 if isodd(len(c4)):

77 y0 = _c4()

78 _s2 = _sum2_

79 while c4:

80 # y1 = ar * y0 - y1 + c4.pop()

81 # y0 = ar * y1 - y0 + c4.pop()

82 y1, t1 = _s2(ar * y0, ar * t0, -y1, -t1, _c4())

83 y0, t0 = _s2(ar * y1, ar * t1, -y0, -t0, _c4())

84 # s = (y0 - y1) * cx

85 s, _ = _s2(cx * y0, _0_0, cx * t0, -cx * y1, -cx * t1)

86 return s

87

88

89_f = float # in _f2 and .geodesicx._C4_24, _27 and _30

90

91

92 def _f2(hi, lo): # in .geodesicx._C4_24, _27 and _30

93 '''(INTERNAL) For C{_coeffs}.

94 '''

95 return _ldexp(_f(hi), 52) + _f(lo)

96

97

98 def _sincos12(sin1, cos1, sin2, cos2, sineg0=False):

99 '''(INTERNAL) Compute the sine and cosine of angle

100 M{ang12 = atan2(sin2, cos2) - atan2(sin1, cos1)}.

101

102 Negate C{sin1} to get C{sin12} and C{cos12} of the sum

103 M{ang12 = atan2(sin2, cos2) + atan2(sin1, cos1)}.

104

105 @kwarg sineg0: If C{True}, make negative C{sin12} zero (C{bool}).

106

107 @return: 2-Tuple C{(sin12, cos12)}.

108 '''

109 s = sin2 * cos1 - sin1 * cos2

110 c = cos2 * cos1 + sin1 * sin2

111 if sineg0 and s < 0:

112 s = _0_0 # max(s, _0_0) or NEG0?

113 return s, c

114

115

116 def _sin1cos2(sin1, cos1, sin2, cos2):

117 '''(INTERNAL) Compute the C{sin1 * cos2} sine and its cosine.

118

119 @return: 2-Tuple C{(sin1 * cos2, hypot(sin1 * sin2, cos1)}.

120 '''

121 s = sin1 * cos2

122 c = _hypot(sin1 * sin2, cos1)

123 return s, c # _norm2(s, c)

124

125

126 def _sinf1cos2d(lat, f1):

127 '''(INTERNAL) See C{GeodesicExact} and C{_GeodesicLineExact}.

128 '''

129 sbet, cbet = _sincos2d(lat)

130 # ensure cbet1 = +epsilon at poles; doing the fix on beta means

131 # that sig12 will be <= 2*tiny for two points at the same pole

132 sbet, cbet = _norm2(sbet * f1, cbet)

133 return sbet, (cbet if fabs(cbet) > _TINY else _TINY)

134

135

136 def _xnC4(**name_nC4):

137 '''(INTERNAL) Validate C{C4order}.

138 '''

139 n, nC4 = _xkwds_item2(name_nC4)

140 if nC4 not in _nC4s or not isinstance(nC4, int):

141 raise GeodesicError(n, nC4, txt_not_=_or(*map(str, _nC4s)))

142 return _nC4s[nC4]

143

144

145# **) MIT License

146#

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

148#

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

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

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

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

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

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

155#

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

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

158#

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

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

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

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

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

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

165# OTHER DEALINGS IN THE SOFTWARE.