1

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

3

4 u'''Discrete Sine Transforms (AuxDST) in Python, transcoded from I{Karney}'s C++ class

5 U{DST<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1DST.html>}

6 in I{GeographicLib version 2.2+}.

7

8 Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2022-2023) and licensed

9 under the MIT/X11 License. For more information, see the U{GeographicLib

10<https://GeographicLib.SourceForge.io>} documentation.

11

12@note: Class L{AuxDST} requires package U{numpy<https://PyPI.org/project/numpy>} to be

13 installed, version 1.16 or newer and needed for I{exact} area calculations.

14'''

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

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

17

18 from pygeodesy.auxilats.auxily import _2cos2x

19 from pygeodesy.basics import isodd, map2, neg, _reverange, _xnumpy

20 from pygeodesy.constants import PI_2, PI_4, isfinite, \

21 _0_0, _0_5, _1_0, _inf_nan

22 from pygeodesy.fsums import Fsum, property_RO

23 from pygeodesy.lazily import _ALL_DOCS

24# from pygeodesy.props import property_RO # from .fsums

25

26__all__ = ()

27__version__ = '23.08.12'

28

29

30 class AuxDST(object):

31 '''Discrete Sine Transforms (DST) for I{Auxiliary Latitudes}.

32

33 @see: I{Karney}'s C++ class U{DST

34 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1DST.html>}.

35 '''

36 _N = 0

37

38 def __init__(self, N):

39 '''New L{AuxDST} instance.

40

41 @arg N: Size, number of points (C{int}).

42 '''

43 if N > 0:

44 self._N = int(N)

45 # kissfft(N, False) # size, inverse

46

47 @staticmethod

48 def evaluate(sinx, cosx, F, *N):

49 '''Compute the Fourier sum given the sine and cosine of the angle,

50 using I{Clenshaw} summation C{sum(B{F}[i] * sin((2*i+1) * x))}

51 for C{i in range(min(len(B{F}), *B{N}))}.

52

53 @arg sinx: The sin(I{sigma}) (C{float}).

54 @arg cosx: The cos(I{sigma}) (C{float}).

55 @arg F: The Fourier coefficients (C{float}[]).

56 @arg N: Optional, (smaller) number of terms to evaluate (C{int}).

57

58 @return: Precison I{Clenshaw} sum (C{float}).

59

60 @see: Methods C{AuxDST.integral} and C{AuxDST.integral2}.

61 '''

62 a = -_2cos2x(cosx, sinx)

63 if isfinite(a):

64 Y0, Y1 = Fsum(), Fsum()

65 n = _len_N(F, *N)

66 if isodd(n):

67 n -= 1

68 Y0 -= F[n]

69 while n > 0: # Y0, Y1 negated

70 n -= 1; Y1 -= Y0 * a + F[n] # PYCHOK semicolon

71 n -= 1; Y0 -= Y1 * a + F[n] # PYCHOK semicolon

72 r = float(_Ys(Y0, -Y1, -sinx))

73 else:

74 r = _inf_nan(a)

75 return r

76

77 @property_RO

78 def _fft_numpy(self):

79 '''(INTERNAL) Get the C{numpy.fft} module, I{once}.

80 '''

81 AuxDST._fft_numpy = fft = _xnumpy(AuxDST, 1, 16).fft # overwrite property_RO

82 return fft

83

84 def _fft_real(self, data):

85 '''(INTERNAL) NumPy's I{kissfft}-like C{transform_real} function,

86 taking C{float}[:N] B{C{data}} and returning C{complex}[:N*2].

87 '''

88 # <https://GitHub.com/mborgerding/kissfft/blob/master/test/testkiss.py>

89 return self._fft_numpy.rfftn(data)

90

91 def _ffts(self, data, cIV):

92 '''(INTERNAL) Compute the DST-III or DST-IV FFTransforms.

93

94 @arg data: Elements DST-III[0:N+1] or DST-IV[0:N] (C{float}[])

95 with DST_III[0] = 0.

96 @arg cIV: If C{True} DST-IV, otherwise DST-III.

97

98 @return: FFTransforms (C{float}[0:N]).

99 '''

100 t, N = (), self.N

101 if N > 0:

102 N2 = N * 2

103 d = list(data)

104 # assert len(d) == N + (0 if cIV else 1)

105

106 if cIV: # DST-IV

107 from cmath import exp as _cexp

108

109 def _cF(c, j, d=-PI_4 / N):

110 return c * _cexp(complex(0, d * j))

111

112 i = 0

113 else: # DST-III

114 i = 1

115 # assert d[0] == _0_0

116

117 def _cF(c, unused): # PYCHOK redef

118 return c

119

120 d += list(reversed(d[i:N])) # i == len(d) - N

121 d += list(map(neg, d[:N2]))

122 c = self._fft_real(d) # complex[0:N*2]

123 n2 = float(-N2)

124 t = tuple(_cF(c[j], j).imag / n2 for j in range(1, N2, 2))

125 return t

126

127 def _ffts2(self, data, F):

128 '''(INTERNAL) Doubled FFTransforms.

129

130 @arg data: Grid centers (C{float}[N]).

131 @arg F: The transforms (C{float}[N])

132

133 @return: Doubled FFTransforms (C{float}[N*2]).

134 '''

135 def _dmF_2(d, F):

136 return (d - F) * _0_5

137

138 def _dpF_2(d, F):

139 return (d + F) * _0_5

140

141 # Copy DST-IV order N transform to d[0:N]

142 d = self._ffts(data, True)

143 N = self._N

144 # assert len(d) >= N and len(F) >= N

145 # (DST-IV order N - DST-III order N) / 2

146 m = map2(_dmF_2, d[:N], F[:N])

147 # (DST-IV order N + DST-III order N) / 2

148 p = map2(_dpF_2, d[:N], F[:N])

149 return p + tuple(reversed(m))

150

151 @staticmethod

152 def integral(sinx, cosx, F, *N):

153 '''Compute the integral of Fourier sum given the sine and cosine

154 of the angle using I{Clenshaw} summation C{-sum(B{F}[i] / (2*i+1) *

155 cos((2*i+1) * x))} for C{i in range(min(len(B{F}), *B{N}))}.

156

157 @arg sinx: The sin(I{sigma}) (C{float}).

158 @arg cosx: The cos(I{sigma}) (C{float}).

159 @arg F: The Fourier coefficients (C{float}[]).

160 @arg N: Optional, (smaller) number of terms to evaluate (C{int}).

161

162 @return: Precison I{Clenshaw} intergral (C{float}).

163

164 @see: Methods C{AuxDST.evaluate} and C{AuxDST.integral2}.

165 '''

166 a = _2cos2x(cosx - sinx, cosx + sinx)

167 if isfinite(a):

168 Y0, Y1 = Fsum(), Fsum()

169 for r in _reverscaled(F, *N):

170 Y1 -= Y0 * a + r

171 Y1, Y0 = Y0, -Y1

172 r = float(_Ys(Y1, Y0, cosx))

173 else:

174 r = _inf_nan(a)

175 return r

176

177 @staticmethod

178 def integral2(sin1, cos1, sin2, cos2, F, *N): # PYCHOK no cover

179 '''Compute the integral of Fourier sum given the sine and cosine

180 of the angles at the end points using I{Clenshaw} summation

181 C{integral(siny, cosy, F) - integral(sinx, cosx, F)}.

182

183 @arg sin1: The sin(I{sigma1}) (C{float}).

184 @arg cos1: The cos(I{sigma1}) (C{float}).

185 @arg sin2: The sin(I{sigma2}) (C{float}).

186 @arg cos2: The cos(I{sigma2}) (C{float}).

187 @arg F: The Fourier coefficients (C{float}[]).

188 @arg N: Optional, (smaller) number of terms to evaluate (C{int}).

189

190 @return: Precison I{Clenshaw} intergral (C{float}).

191

192 @see: Methods C{AuxDST.evaluate} and C{AuxDST.integral}.

193 '''

194 # 2 * cos(y - x)*cos(y + x) -> 2 * cos(2 * x)

195 a = _2cos2x(cos2 * cos1, sin2 * sin1)

196 # -2 * sin(y - x)*sin(y + x) -> 0

197 b = -_2cos2x(sin2 * cos1, cos2 * sin1)

198 if isfinite(a) and isfinite(b):

199 Y0, Y1 = Fsum(), Fsum()

200 Z0, Z1 = Fsum(), Fsum()

201 for r in _reverscaled(F, *N):

202 Y1 -= Y0 * a + Z0 * b + r

203 Z1 -= Y0 * b + Z0 * a

204 Y1, Y0 = Y0, -Y1

205 Z1, Z0 = Z0, -Z1

206 r = float(_Ys(Y1, Y0, cos2 - cos1) +

207 _Ys(Z1, Z0, cos2 + cos1))

208 else:

209 r = _inf_nan(a, b)

210 return r

211

212 @property_RO

213 def N(self):

214 '''Get this DST's size, number of points (C{int}).

215 '''

216 return self._N

217

218 def refine(self, f, F):

219 '''Double the number of sampled points on a Fourier series.

220

221 @arg f: Single-argument function (C{callable(sigma)} with

222 C{sigma = PI_4 * j / N for j in range(1, N*2, 2)}.

223 @arg F: The initial Fourier series coefficients (C{float}[:N]).

224

225 @return: Fourier series coefficients (C{float}[:N*2]).

226 '''

227 def _data(_f, N): # [:N]

228 if N > 0:

229 d = PI_4 / N

230 for j in range(1, N*2, 2):

231 yield _f(d * j)

232

233 return self._ffts2(_data(f, self.N), F)

234

235 def reset(self, N):

236 '''Reset this DST.

237

238 @arg N: Size, number of points (C{int}).

239

240 @return: The new size (C{int}, non-negative).

241 '''

242 self._N = N = max(0, N)

243 # kissfft.assign(N*2, False) # "reset" size, inverse

244 return N

245

246 def transform(self, f):

247 '''Compute C{N + 1} terms in the Fourier series.

248

249 @arg f: Single-argument function (C{callable(sigma)} with

250 C{sigma = PI_2 * i / N for i in range(1, N + 1)}.

251

252 @return: Fourier series coefficients (C{float}[:N + 1]).

253 '''

254 def _data(_f, N): # [:N + 1]

255 yield _0_0 # data[0] = 0

256 if N > 0:

257 d = PI_2 / N

258 for i in range(1, N + 1):

259 yield _f(d * i)

260

261 return self._ffts(_data(f, self.N), False)

262

263

264 def _len_N(F, *N):

265 # Adjusted C{len(B{F})}.

266 return min(len(F), *N) if N else len(F)

267

268

269 def _reverscaled(F, *N):

270 # Yield F[:N], reversed and scaled

271 for n in _reverange(_len_N(F, *N)):

272 yield F[n] / (n * 2 + _1_0)

273

274

275 def _Ys(X, Y, s):

276 # Return M{(X - Y) * s}, overwriting X

277 X -= Y

278 X *= s

279 return X

280

281

282__all__ += _ALL_DOCS(AuxDST)

283

284# **) MIT License

285#

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

287#

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

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

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

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

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

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

294#

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

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

297#

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

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

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

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

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

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

304# OTHER DEALINGS IN THE SOFTWARE.