Coverage for panelaero/DLM.py: 82%
227 statements
« prev ^ index » next coverage.py v7.6.1, created at 2024-09-04 09:17 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
4import copy
5import logging
6import numpy as np
8from panelaero import VLM
10# turn off warnings (divide by zero, multiply NaN, ...) as singularities are expected to occur
11np.seterr(all='ignore')
14def calc_Qjj(aerogrid, Ma, k, method='parabolic'):
15 # calc steady contributions using VLM
16 Ajj_VLM, _ = VLM.calc_Ajj(aerogrid=copy.deepcopy(aerogrid), Ma=Ma)
17 if k == 0.0:
18 # no oscillatory / unsteady contributions at k=0.0
19 Ajj_DLM = np.zeros((aerogrid['n'], aerogrid['n']))
20 else:
21 # calc oscillatory / unsteady contributions using DLM
22 Ajj_DLM = calc_Ajj(aerogrid=copy.deepcopy(aerogrid), Ma=Ma, k=k, method=method)
23 Ajj = Ajj_VLM + Ajj_DLM
24 Qjj = -np.linalg.inv(Ajj)
25 return Qjj
28def calc_Qjjs(aerogrid, Ma, k, xz_symmetry=False):
29 # allocate memory
30 Qjj = np.zeros((len(Ma), len(k), aerogrid['n'], aerogrid['n']), dtype='complex') # dim: Ma,k,n,n
31 # Consideration of XZ symmetry like in VLM.
32 if xz_symmetry:
33 n = aerogrid['n']
34 aerogrid = VLM.mirror_aerogrid_xz(aerogrid)
36 # loop over mach number and freq.
37 for im, Ma_i in enumerate(Ma):
38 # calc steady contributions using VLM
39 Ajj_VLM, _ = VLM.calc_Ajj(aerogrid=copy.deepcopy(aerogrid), Ma=Ma_i)
40 for ik, k_i in enumerate(k):
41 if k_i == 0.0:
42 # no oscillatory / unsteady contributions at k=0.0
43 Ajj_DLM = np.zeros((aerogrid['n'], aerogrid['n']))
44 else:
45 # calc oscillatory / unsteady contributions using DLM
46 Ajj_DLM = calc_Ajj(aerogrid=copy.deepcopy(aerogrid), Ma=Ma_i, k=k_i, method='parabolic')
47 Ajj = Ajj_VLM + Ajj_DLM
48 Ajj_inv = -np.linalg.inv(Ajj)
49 if xz_symmetry:
50 Qjj[im, ik] = Ajj_inv[0:n, 0:n] - Ajj_inv[n:2 * n, 0:n]
51 else:
52 Qjj[im, ik] = Ajj_inv
53 return Qjj
56def calc_Ajj(aerogrid, Ma, k, method='parabolic'):
57 # Calculates one unsteady AIC matrix (Qjj = -Ajj^-1) at given Mach number and frequency
58 #
59 # l_2
60 # 4 o---------o 3
61 # | |
62 # u --> b_1 | l k j | b_2
63 # | |
64 # 1 o---------o 2
65 # y l_1
66 # |
67 # z.--- x
68 #
69 # M = Mach number
70 # k = omega/U, the "classical" definition, not Nastran definition!
71 # Nomencalture with receiving (r), minus (-e), plus (e), sending (s/0) point and semiwidth e following Rodden 1968
72 Pr = aerogrid['offset_j'] # receiving (r)
73 Pm = aerogrid['offset_P1'] # minus (-e)
74 Pp = aerogrid['offset_P3'] # plus (e)
75 Ps = aerogrid['offset_l'] # sending (s/0)
76 e = np.absolute(np.repeat(np.array(0.5 * ((Pp[:, 2] - Pm[:, 2]) ** 2.0 + (Pp[:, 1] - Pm[:, 1]) ** 2.0) ** 0.5, ndmin=2),
77 aerogrid['n'], axis=0)) # semiwidth
78 e2 = e ** 2.0
79 e3 = e ** 3.0
80 e4 = e ** 4.0
81 chord = np.repeat(np.array(aerogrid['l'], ndmin=2), aerogrid['n'], axis=0)
83 # Catch panels which are not defined from left to right and issue a warning.
84 # Not sure with purely vertical panels though (bottom to top vs. top to bottom)...
85 if np.any(aerogrid['N'][:, 2] < 0.0):
86 logging.warning('Detected upside down / flipped aerodynamic panels! \n'
87 'User action: Always define panels from left to right.')
89 # cartesian coordinates of receiving points relative to sending points
90 xsr = np.array(Pr[:, 0], ndmin=2).T - np.array(Ps[:, 0], ndmin=2)
91 ysr = np.array(Pr[:, 1], ndmin=2).T - np.array(Ps[:, 1], ndmin=2)
92 zsr = np.array(Pr[:, 2], ndmin=2).T - np.array(Ps[:, 2], ndmin=2)
94 # dihedral angle gamma = arctan(dz/dy) and sweep angle lambda = arctan(dx/dy)
95 sinGamma = (Pp[:, 2] - Pm[:, 2]) / (2.0 * e)
96 cosGamma = (Pp[:, 1] - Pm[:, 1]) / (2.0 * e)
97 tanLambda = (Pp[:, 0] - Pm[:, 0]) / (2.0 * e)
98 gamma = np.arcsin(sinGamma)
99 # relative dihedral angle between receiving point and sending boxes
100 gamma_sr = np.array(gamma, ndmin=2) - np.array(gamma, ndmin=2).T
102 # local coordinates of receiving point relative to sending point
103 ybar = ysr * cosGamma + zsr * sinGamma
104 zbar = zsr * cosGamma - ysr * sinGamma
106 # pre-calculate some values which will be used a couple of times
107 ybar2 = ybar ** 2.0
108 ybar4 = ybar ** 4.0
109 zbar2 = zbar ** 2.0
110 zbar4 = zbar ** 4.0
111 ratio = 2.0 * e * np.abs(zbar) / (ybar2 + zbar2 - e2)
112 L = np.log(((ybar - e) ** 2.0 + zbar2) / ((ybar + e) ** 2.0 + zbar2)) # Checked with Nastran idf1.f
114 # Condition 1, planar, the value of 0.001 is taken from Nastran
115 i0 = np.abs(zbar) / e <= 0.001
116 # Condition 2, co-planar / close-by
117 ia = (np.abs(ratio) <= 0.3) & (np.abs(zbar) / e > 0.001)
118 # Condition 3, the rest / further away
119 ir = (np.abs(ratio) > 0.3) & (np.abs(zbar) / e > 0.001)
120 # check that all conditions are captured: np.all(i0 + ia + ir) == True
122 alpha = np.zeros(e.shape)
123 funny_series = 0.0
124 for n in range(2, 8):
125 funny_series += (-1.0) ** n / (2.0 * n - 1.0) * ratio[ia] ** (2.0 * n - 4.0)
126 # Rodden 1971, eq 33, Rodden 1972, eq 31b and Rodden 1998, eq 25
127 alpha[ia] = 4.0 * e[ia] ** 4.0 / (ybar2[ia] + zbar2[ia] - e2[ia]) ** 2.0 * funny_series
129 if method == 'parabolic':
130 # Rodden et at. 1971 and 1972
132 # Initial values
133 Fparabolic = np.zeros(e.shape)
134 # Condition 1, planar
135 Fparabolic[i0] = 2.0 * e[i0] / (ybar2[i0] - e2[i0]) # OK, idf1.f
136 # Condition 2, co-planar / close-by
137 # Rodden 1971, eq 32, Checked with Nastran idf1.f
138 Fparabolic[ia] = 2.0 * e[ia] / (ybar2[ia] + zbar2[ia] - e2[ia]) * (1.0 - alpha[ia] * zbar2[ia] / e2[ia])
139 # Condition 3, the rest / further away
140 # Rodden 1971, eq 31b, Checked with Nastran idf1.f
141 Fparabolic[ir] = 1.0 / np.abs(zbar[ir]) * np.arctan2(2.0 * e[ir] * np.abs(zbar[ir]), (ybar2[ir] + zbar2[ir] - e2[ir]))
142 # Extra Condition found in Nastran idf2.f, line 26
143 # Fparabolic[np.abs(1.0/ratio) <= 0.0001] = 0.0
145 # call the kernel function with Laschka approximation
146 P1m, P2m = kernelfunction(xsr, ybar, zbar, gamma_sr, tanLambda, -e, k, Ma, method='Laschka')
147 P1p, P2p = kernelfunction(xsr, ybar, zbar, gamma_sr, tanLambda, +e, k, Ma, method='Laschka')
148 P1s, P2s = kernelfunction(xsr, ybar, zbar, gamma_sr, tanLambda, 0, k, Ma, method='Laschka')
150 # define terms used in the parabolic approximation, Nastran incro.f
151 A1 = (P1m - 2.0 * P1s + P1p) / (2.0 * e2) # Rodden 1971, eq 28
152 B1 = (P1p - P1m) / (2.0 * e) # Rodden 1971, eq 29
153 C1 = P1s # Rodden 1971, eq 30
155 A2 = (P2m - 2.0 * P2s + P2p) / (2.0 * e2) # Rodden 1971, eq 37
156 B2 = (P2p - P2m) / (2.0 * e) # Rodden 1971, eq 38
157 C2 = P2s # Rodden 1971, eq 39
159 # The "planar" part
160 # -----------------
161 # normalwash matrix, Rodden 1971, eq 34
162 D1rs = chord / (np.pi * 8.0) \
163 * (((ybar2 - zbar2) * A1 + ybar * B1 + C1) * Fparabolic
164 + (0.5 * B1 + ybar * A1) * np.log(((ybar - e) ** 2.0 + zbar2) / ((ybar + e) ** 2.0 + zbar2))
165 + 2.0 * e * A1) # Checked with Nastran idf1.f & incro.f
167 # The "nonplanar" part
168 # --------------------
169 D2rs = np.zeros(e.shape, dtype='complex')
171 # Condition 1, similar to above but with different boundary, Rodden 1971 eq 40
172 ib = (np.abs(1.0 / ratio) <= 0.1) & (np.abs(zbar) / e > 0.001)
173 D2rs[ib] = chord[ib] / (16.0 * np.pi * zbar2[ib]) \
174 * (((ybar2[ib] + zbar2[ib]) * A2[ib] + ybar[ib] * B2[ib] + C2[ib]) * Fparabolic[ib]
175 + 1.0 / ((ybar[ib] + e[ib]) ** 2.0 + zbar2[ib])
176 * (((ybar2[ib] + zbar2[ib]) * ybar[ib] + (ybar2[ib] - zbar2[ib]) * e[ib])
177 * A2[ib] + (ybar2[ib] + zbar2[ib] + ybar[ib] * e[ib]) * B2[ib] + (ybar[ib] + e[ib]) * C2[ib])
178 - 1.0 / ((ybar[ib] - e[ib]) ** 2.0 + zbar2[ib])
179 * (((ybar2[ib] + zbar2[ib]) * ybar[ib] - (ybar2[ib] - zbar2[ib]) * e[ib])
180 * A2[ib] + (ybar2[ib] + zbar2[ib] - ybar[ib] * e[ib]) * B2[ib] + (ybar[ib] - e[ib]) * C2[ib])
181 ) # Checked with Nastran idf2.f
183 # Condition 2, Rodden 1971 eq 41
184 ic = (np.abs(1.0 / ratio) > 0.1) & (np.abs(zbar) / e > 0.001)
185 # reconstruct alpha from eq 32, NOT eq 33!
186 alpha[i0] = ((2.0 * e2[i0]) / (ybar2[i0] - e2[i0])) ** 2.0 # Nastran idf2.f, line 75
187 alpha[ir] = (1.0 - Fparabolic[ir] * (ybar2[ir] + zbar2[ir] - e2[ir]) / (2.0 * e[ir])) / zbar2[ir] * e2[ir]
189 D2rs[ic] = chord[ic] * e[ic] / (8.0 * np.pi * (ybar2[ic] + zbar2[ic] - e2[ic])) \
190 * ((2.0 * (ybar2[ic] + zbar2[ic] + e2[ic]) * (e2[ic] * A2[ic] + C2[ic]) + 4.0 * ybar[ic] * e2[ic] * B2[ic])
191 / (((ybar[ic] + e[ic]) ** 2.0 + zbar2[ic]) * ((ybar[ic] - e[ic]) ** 2.0 + zbar2[ic]))
192 - alpha[ic] / e2[ic] * ((ybar2[ic] + zbar2[ic]) * A2[ic] + ybar[ic] * B2[ic] + C2[ic])
193 ) # Checked with Nastran idf2.f
195 elif method == 'quartic':
196 # Rodden et al. 1998
197 # Why chooses Rodden in 1972 and 1998 a more complicated formulation with d1,2?
198 # Only to place the tangens into the right quadrant? --> Is there no arctan2 in Fortran?!?
199 # Still, we have to use that formulation as d1,2 and epsilon will be used later in eq 34.
201 # Note that there is a difference and/or mistake (?) in eq. 23 in Roddel et al. 1998
202 # compared to eq. 30b in Roddel et al. 1972. The following values appear to be correct:
203 d1 = np.zeros(e.shape)
204 d2 = np.zeros(e.shape)
205 i1 = (ybar2 + zbar2 - e2) > 0.0
206 d1[i1] = 1.0
207 d2[i1] = 0.0
208 i2 = (ybar2 + zbar2 - e2) == 0.0
209 d1[i2] = 0.0
210 d2[i2] = 0.5
211 i3 = (ybar2 + zbar2 - e2) < 0.0
212 d1[i3] = 1.0
213 d2[i3] = 1.0
215 # Rodden 1998, eq 24 and 25
216 epsilon = np.zeros(e.shape)
217 epsilon[i0] = 2.0 * e[i0] / (ybar2[i0] - e2[i0])
218 epsilon[ia] = alpha[ia] # Rodden 1998, eq 25
219 epsilon[ir] = e2[ir] / zbar2[ir] * (1.0 - 1.0 / ratio[ir] * np.arctan(ratio[ir])) # Rodden 1998, eq 24
220 iar = ia + ir
221 # Initial values
222 Fquartic = np.zeros(e.shape)
223 # Rodden 1998, eq. 22 without terms including z because z==0
224 Fquartic[i0] = d1[i0] * 2.0 * e[i0] / (ybar2[i0] - e2[i0])
225 # Rodden 1998, eq. 22
226 Fquartic[iar] = d1[iar] * 2.0 * e[iar] / (ybar2[iar] + zbar2[iar] - e2[iar]) \
227 * (1.0 - epsilon[iar] * zbar2[iar] / e2[iar]) + d2[iar] * np.pi / np.abs(zbar[iar])
228 # check: np.allclose(Fparabolic, Fquartic)
230 # Rodden 1998
231 # call the kernel function with Desmarais approximation
232 P1m, P2m = kernelfunction(xsr, ybar, zbar, gamma_sr, tanLambda, -e, k, Ma, method='Desmarais')
233 P1mh, P2mh = kernelfunction(xsr, ybar, zbar, gamma_sr, tanLambda, -e / 2.0, k, Ma, method='Desmarais')
234 P1p, P2p = kernelfunction(xsr, ybar, zbar, gamma_sr, tanLambda, +e, k, Ma, method='Desmarais')
235 P1ph, P2ph = kernelfunction(xsr, ybar, zbar, gamma_sr, tanLambda, +e / 2.0, k, Ma, method='Desmarais')
236 P1s, P2s = kernelfunction(xsr, ybar, zbar, gamma_sr, tanLambda, 0.0, k, Ma, method='Desmarais')
238 # define terms used in the quartic approximation
239 A1 = -1.0 / (6.0 * e2) * (P1m - 16.0 * P1mh + 30.0 * P1s - 16.0 * P1ph + P1p) # Rodden 1998, eq 15
240 B1 = +1.0 / (6.0 * e) * (P1m - 8.0 * P1mh + 8.0 * P1ph - P1p) # Rodden 1998, eq 16
241 C1 = P1s # Rodden 1998, eq 17
242 D1 = -2.0 / (3.0 * e3) * (P1m - 2.0 * P1mh + 2.0 * P1ph - P1p) # Rodden 1998, eq 18
243 E1 = +2.0 / (3.0 * e4) * (P1m - 4.0 * P1mh + 6.0 * P1s - 4.0 * P1ph + P1p) # Rodden 1998, eq 19
245 A2 = -1.0 / (6.0 * e2) * (P2m - 16.0 * P2mh + 30.0 * P2s - 16.0 * P2ph + P2p) # Rodden 1998, eq 28
246 B2 = +1.0 / (6.0 * e) * (P2m - 8.0 * P2mh + 8.0 * P2ph - P2p) # Rodden 1998, eq 29
247 C2 = P2s # Rodden 1998, eq 30
248 D2 = -2.0 / (3.0 * e3) * (P2m - 2.0 * P2mh + 2.0 * P2ph - P2p) # Rodden 1998, eq 31
249 E2 = +2.0 / (3.0 * e4) * (P2m - 4.0 * P2mh + 6.0 * P2s - 4.0 * P2ph + P2p) # Rodden 1998, eq 32
251 # The "planar" part
252 # -----------------
253 # normalwash matrix, Rodden 1998, eq 20
254 D1rs = chord / (np.pi * 8.0) \
255 * (((ybar2 - zbar2) * A1 + ybar * B1 + C1 + ybar * (ybar2 - 3.0 * zbar2) * D1
256 + (ybar4 - 6.0 * ybar2 * zbar2 + zbar4) * E1) * Fquartic
257 + (0.5 * B1 + ybar * A1 + 0.5 * (3.0 * ybar2 - zbar2) * D1 + 2.0 * ybar * (ybar2 - zbar2) * E1) * L
258 + 2.0 * e * (A1 + 2.0 * ybar * D1 + (3.0 * ybar2 - zbar2 + 1.0 / 3.0 * e2) * E1)
259 )
260 # The "nonplanar" part
261 # --------------------
262 D2rs = np.zeros(e.shape, dtype='complex')
264 # Condition 1, similar to above but with different boundary, Rodden 1998 eq 33
265 ib = (np.abs(1.0 / ratio) <= 0.1) & (np.abs(zbar) / e > 0.001)
266 D2rs[ib] = chord[ib] / (16.0 * np.pi * zbar2[ib]) \
267 * (Fquartic[ib]
268 * ((ybar2[ib] + zbar2[ib]) * A2[ib]
269 + ybar[ib] * B2[ib]
270 + C2[ib]
271 + ybar[ib] * (ybar2[ib] + 3.0 * zbar2[ib]) * D2[ib]
272 + (ybar4[ib] + 6.0 * ybar2[ib] * zbar2[ib] - 3.0 * zbar4[ib]) * E2[ib]
273 )
274 + 1.0 / ((ybar[ib] + e[ib]) ** 2.0 + zbar2[ib])
275 * (((ybar2[ib] + zbar2[ib]) * ybar[ib] + (ybar2[ib] - zbar2[ib]) * e[ib]) * A2[ib]
276 + (ybar2[ib] + zbar2[ib] + ybar[ib] * e[ib]) * B2[ib]
277 + (ybar[ib] + e[ib]) * C2[ib]
278 + (ybar4[ib] - zbar4[ib] + (ybar2[ib] - 3.0 * zbar2[ib]) * ybar[ib] * e[ib]) * D2[ib]
279 + ((ybar4[ib] - 2.0 * ybar2[ib] * zbar2[ib] - 3.0 * zbar4[ib]) * ybar[ib]
280 + (ybar4[ib] - 6.0 * ybar2[ib] * zbar2[ib] + zbar4[ib]) * e[ib]) * E2[ib]
281 )
282 - 1.0 / ((ybar[ib] - e[ib]) ** 2.0 + zbar2[ib])
283 * (((ybar2[ib] + zbar2[ib]) * ybar[ib] - (ybar2[ib] - zbar2[ib]) * e[ib]) * A2[ib]
284 + (ybar2[ib] + zbar2[ib] - ybar[ib] * e[ib]) * B2[ib]
285 + (ybar[ib] - e[ib]) * C2[ib]
286 + (ybar4[ib] - zbar4[ib] - (ybar2[ib] - 3.0 * zbar2[ib]) * ybar[ib] * e[ib]) * D2[ib]
287 + ((ybar4[ib] - 2.0 * ybar2[ib] * zbar2[ib] - 3.0 * zbar4[ib]) * ybar[ib]
288 - (ybar4[ib] - 6.0 * ybar2[ib] * zbar2[ib] + zbar4[ib]) * e[ib]) * E2[ib]
289 )
290 + (zbar2[ib] * L[ib]) * D2[ib]
291 + 4.0 * zbar2[ib] * (e[ib] + ybar[ib] * L[ib]) * E2[ib]
292 )
294 # Condition 2, Rodden 1998 eq 34
295 ic = (np.abs(1.0 / ratio) > 0.1) & (np.abs(zbar) / e > 0.001)
296 D2rs[ic] = chord[ic] * e[ic] / (8.0 * np.pi * (ybar2[ic] + zbar2[ic] - e2[ic])) \
297 * (1.0 / (((ybar[ic] + e[ic]) ** 2.0 + zbar2[ic]) * ((ybar[ic] - e[ic]) ** 2.0 + zbar2[ic]))
298 * (2.0 * (ybar2[ic] + zbar2[ic] + e2[ic]) * (e2[ic] * A2[ic] + C2[ic])
299 + 4.0 * ybar[ic] * e2[ic] * B2[ic]
300 + 2.0 * ybar[ic] * (ybar4[ic] - 2.0 * e2[ic] * ybar2[ic] + 2.0 * ybar2[ic] * zbar2[ic] + 3.0 * e4[ic]
301 + 2.0 * e2[ic] * zbar2[ic] + zbar4[ic]) * D2[ic]
302 + 2.0 * (3.0 * ybar[ic] ** 6.0 - 7.0 * e2[ic] * ybar4[ic] + 5.0 * ybar4[ic] * zbar2[ic]
303 + 6.0 * e4[ic] * ybar2[ic] + 6.0 * e2[ic] * ybar2[ic] * zbar2[ic]
304 - 3.0 * e2[ic] * zbar4[ic] - zbar[ic] ** 6.0 + ybar2[ic] * zbar4[ic]
305 - 2.0 * e4[ic] * zbar2[ic]) * E2[ic]
306 )
307 - (d1[ic] * epsilon[ic] + e2[ic] / zbar2[ic] * (1.0 - d1[ic] - d2[ic] * np.pi / ratio[ic])) / e2[ic]
308 * ((ybar2[ic] + zbar2[ic]) * A2[ic]
309 + ybar[ic] * B2[ic]
310 + C2[ic]
311 + ybar[ic] * (ybar2[ic] + 3.0 * zbar2[ic]) * D2[ic]
312 + (ybar4[ic] + 6.0 * ybar2[ic] * zbar2[ic] - 3.0 * zbar4[ic]) * E2[ic]
313 )
314 ) \
315 + chord[ic] / (8.0 * np.pi) * (D2[ic] / 2.0 * L[ic] + 2.0 * (e[ic] + ybar[ic] * L[ic]) * E2[ic])
316 else:
317 logging.error('Method {} not implemented!'.format(method))
319 # add planar and non-planar parts, # Rodden eq 22
320 # the steady part D0 has already been subtracted inside the kernel function
321 Drs = D1rs + D2rs
322 return Drs
325def kernelfunction(xbar, ybar, zbar, gamma_sr, tanLambda, ebar, k, M, method='Laschka'):
326 # This is the function that calculates "the" kernel function(s) of the DLM.
327 # K1,2 are reformulated in Rodden 1971 compared to Rodden 1968 and include new
328 # conditions, e.g. for co-planar panels.
329 # Note that the signs of K1,2 are switched in Rodden 1971, in this implementation
330 # we stay with the 1968 convention and we don't want to mess with the steady part
331 # from the VLM. Also, we directly subtract the steady parts K10 and K20, as the
332 # steady contribution will be added later from the VLM.
333 # Note: Rodden has the habit of leaving out some brackets in his formulas. This
334 # applies to eq 11, 7 and 8 where it is not clear which parts belong to the denominator.
336 r1 = ((ybar - ebar) ** 2.0 + zbar ** 2.0) ** 0.5 # Rodden 1971, eq 4
337 beta2 = 1.0 - (M ** 2.0) # Rodden 1971, eq 9
338 R = ((xbar - ebar * tanLambda) ** 2.0 + beta2 * r1 ** 2.0) ** 0.5 # Rodden 1971, eq 10
339 u1 = (M * R - xbar + ebar * tanLambda) / (beta2 * r1) # Rodden 1971, eq 11
340 k1 = k * r1 # Rodden 1971, eq 12 with k = w/U
341 j = 1j # imaginary number
342 ejku = np.exp(-j * k1 * u1) # pre-multiplication
344 # direction cosine matrices
345 T1 = np.cos(gamma_sr) # Rodden 1971, eq 5
346 T2 = zbar * (zbar * np.cos(gamma_sr) + (ybar - ebar) * np.sin(gamma_sr)) # Rodden 1971, eq 21a: T2_new = T2_old*r1^2
348 # Approximation of intergrals I1,2, Rodden 1971, eq 13+14
349 I1, I2 = get_integrals12(u1, k1, method)
351 # Formulation of K1,2 by Landahl, Rodden 1971, eq 7+8
352 K1 = -I1 - ejku * M * r1 / R / (1 + u1 ** 2.0) ** 0.5
353 K2 = 3.0 * I2 + j * k1 * ejku * (M ** 2.0) * (r1 ** 2.0) / (R ** 2.0) / (1.0 + u1 ** 2.0) ** 0.5 \
354 + ejku * M * r1 * ((1.0 + u1 ** 2.0) * beta2 * r1 ** 2.0 / R ** 2.0 + 2.0 + M * r1 * u1 / R) \
355 / R / (1.0 + u1 ** 2.0) ** 1.5
357 # This is the analytical solution for K1,2 at k=0.0, Rodden 1971, eq 15+16
358 K10 = -1.0 - (xbar - ebar * tanLambda) / R
359 K20 = 2.0 + (xbar - ebar * tanLambda) * (2.0 + beta2 * r1 ** 2.0 / R ** 2.0) / R
361 # Resolve the singularity arising when r1 = 0
362 ir0xpos = (r1 == 0) & (xbar >= 0.0)
363 ir0xneg = (r1 == 0) & (xbar < 0.0)
364 K1[ir0xpos] = -2.0
365 K2[ir0xpos] = +4.0
366 K1[ir0xneg] = 0.0
367 K2[ir0xneg] = 0.0
369 # Rodden 1971, eq 27b, check: -K1*np.exp(-j*k*xbar)*T1
370 P1 = -(K1 * np.exp(-j * k * (xbar - ebar * tanLambda)) - K10) * T1
371 # Rodden 1971, eq 36b, check: -K2*np.exp(-j*k*xbar)*T2/r1**2.0
372 P2 = -(K2 * np.exp(-j * k * (xbar - ebar * tanLambda)) - K20) * T2
374 return P1, P2
377def get_integrals12(u1, k1, method='Laschka'):
379 I1 = np.zeros(u1.shape, dtype='complex')
380 I2 = np.zeros(u1.shape, dtype='complex')
382 ipos = u1 >= 0.0
383 I1[ipos], I2[ipos] = integral_approximations(u1[ipos], k1[ipos], method)
385 ineg = u1 < 0.0
386 I10, I20 = integral_approximations(0.0 * u1[ineg], k1[ineg], method)
387 I1n, I2n = integral_approximations(-u1[ineg], k1[ineg], method)
388 I1[ineg] = 2.0 * I10.real - I1n.real + 1j * I1n.imag # Rodden 1971, eq A.5
389 I2[ineg] = 2.0 * I20.real - I2n.real + 1j * I2n.imag # Rodden 1971, eq A.9
390 return I1, I2
393def integral_approximations(u1, k1, method='Laschka'):
394 if method == 'Laschka':
395 logging.debug('Using Laschka approximation in DLM')
396 I1, I2 = laschka_approximation(u1, k1)
397 elif method == 'Desmarais':
398 logging.debug('Using Desmarais approximation in DLM')
399 I1, I2 = desmarais_approximation(u1, k1)
400 elif method == 'Watkins':
401 logging.warning('Using Watkins (not preferred!) approximation in DLM.')
402 I1, I2 = watkins_approximation(u1, k1)
403 else:
404 logging.error('Method {} not implemented!'.format(method))
405 return I1, I2
408def desmarais_approximation(u1, k1):
409 # Adapted formulas from laschka_approximation
410 a12 = [0.000319759140, -0.000055461471, 0.002726074362, 0.005749551566,
411 0.031455895072, 0.106031126212, 0.406838011567, 0.798112357155,
412 -0.417749229098, 0.077480713894, -0.012677284771, 0.001787032960]
413 m = 1.0
414 b = 0.009054814793
415 j = 1j
416 ejku = np.exp(-j * k1 * u1) # pre-multiplication
417 I0 = 0.0
418 J0 = 0.0
419 for n, a in zip(range(1, 13), a12):
420 nm = n / m
421 nmbk = (2.0 ** nm) ** 2.0 * b ** 2.0 + k1 ** 2.0
422 I0 += a * np.exp(-(2.0 ** nm) * b * u1) / nmbk * ((2.0 ** nm) * b - j * k1)
423 J0 += a * np.exp(-(2.0 ** nm) * b * u1) / nmbk ** 2.0 \
424 * ((2.0 ** nm) ** 2.0 * b ** 2.0 - k1 ** 2.0 + (2.0 ** nm) * b * u1 * nmbk
425 - j * k1 * (2.0 * (2.0 ** nm) * b + u1 * nmbk))
426 # I1 as in Rodden 1971, eq A.1
427 I1 = (1.0 - u1 / (1.0 + u1 ** 2.0) ** 0.5 - 1j * k1 * I0) * ejku
428 # I2 as in Rodden 1971, eq A.6,
429 I2 = ((2.0 + j * k1 * u1) * (1.0 - u1 / (1.0 + u1 ** 2.0) ** 0.5)
430 - u1 / (1.0 + u1 ** 2.0) ** 1.5 - j * k1 * I0 + k1 ** 2.0 * J0) * ejku / 3.0
431 return I1, I2
434def laschka_approximation(u1, k1):
435 # Approximate integral I0, Rodden 1971, eq A.4
436 # Approximate integral J0, Rodden 1971, eq A.8
437 # These are the coefficients in exponential approximation of u/(1+u**2.0)**0.5
438 # The values are difficult to read in Rodden 1971 but are also given in Blair 1992, page 89.
439 a11 = [+0.24186198, -2.7918027, +24.991079, -111.59196, +271.43549, -305.75288,
440 -41.183630, +545.98537, -644.78155, +328.72755, -64.279511]
441 c = 0.372
442 j = 1j
443 ejku = np.exp(-j * k1 * u1) # pre-multiplication
444 I0 = 0.0
445 J0 = 0.0
446 for n, a in zip(range(1, 12), a11):
447 nck = n ** 2.0 * c ** 2.0 + k1 ** 2.0
448 I0 += a * np.exp(-n * c * u1) / nck * (n * c - j * k1)
449 J0 += a * np.exp(-n * c * u1) / nck ** 2.0 * (n ** 2.0 * c ** 2.0 - k1 ** 2.0 + n * c * u1 * nck
450 - j * k1 * (2.0 * n * c + u1 * nck))
451 # I1 as in Rodden 1971, eq A.1
452 I1 = (1.0 - u1 / (1.0 + u1 ** 2.0) ** 0.5 - 1j * k1 * I0) * ejku
453 # I2 as in Rodden 1971, eq A.6,
454 # but divided by 3.0 for compatibility,
455 # and with the square bracket at the correct location ;)
456 I2 = ((2.0 + j * k1 * u1) * (1.0 - u1 / (1.0 + u1 ** 2.0) ** 0.5) - u1 / (1.0 + u1 ** 2.0) ** 1.5
457 - j * k1 * I0 + k1 ** 2.0 * J0) * ejku / 3.0
458 return I1, I2
461def watkins_approximation(u1, k1):
462 # This is the old/original approximation of integrals I1,2 as in Rodden 1968, page 3.
463 # The following code is take from the previous Matlab implementation of the DLM,
464 # has been rearranged and is used for comparison.
465 a1 = 0.101
466 a2 = 0.899
467 a3 = 0.09480933
468 b1 = 0.329
469 b2 = 1.4067
470 b3 = 2.90
471 j = 1j
472 ejku = np.exp(-j * k1 * u1) # pre-multiplication
473 # evaluate with numpy, slower
474 i1 = a1 * np.exp((-b1 - (j * k1)) * u1) / (b1 + (j * k1)) \
475 + a2 * np.exp((-b2 - (j * k1)) * u1) / (b2 + (j * k1))
476 i2 = (a3 / (((b3 + (j * k1)) ** 2.0) + (np.pi ** 2.0))) * (((b3 + (j * k1)) * np.sin(np.pi * u1))
477 + (np.pi * np.cos(np.pi * u1))) \
478 * np.exp((-b3 - (j * k1)) * u1)
479 I1_temp = i1 + i2
480 I1 = ((1 - (u1 / (1 + u1 ** 2.0) ** 0.5)) * ejku) - (j * k1 * I1_temp)
482 I2_1 = I1
483 # evaluate with numpy, slower
484 I2_2_1 = a1 * np.exp(-(b1 + j * k1) * u1) / ((b1 + j * k1) ** 2.0) \
485 + a2 * np.exp(-(b2 + j * k1) * u1) / ((b2 + j * k1) ** 2.0) \
486 + ((a3 * np.exp(-(b3 + j * k1) * u1) / (((b3 + j * k1) ** 2.0 + np.pi ** 2.0) ** 2.0))
487 * (np.pi * ((np.pi * np.sin(np.pi * u1)) - ((b3 + j * k1) * np.cos(np.pi * u1)))
488 - ((b3 + j * k1) * (np.pi * np.cos(np.pi * u1) + ((b3 + j * k1) * np.sin(np.pi * u1))))))
489 I2_2 = (ejku * (u1 ** 3.0) / ((1 + u1 ** 2.0) ** (3.0 / 2.0)) - I2_1
490 - ejku * u1 / (1 + u1 ** 2.0) ** 0.5) / 3.0 - (k1 * k1 * I2_2_1 / 3.0)
491 I2 = I2_1 + I2_2
492 return I1, I2